A Planet-Driven Spiral Wave Instability

The Spiral Wave Instability Induced by a Giant Planet: I. Particle Stirring in the Inner Regions of Protoplanetary Disks


We have recently shown that spiral density waves propagating in accretion disks can undergo a parametric instability by resonantly coupling with and transferring energy into pairs of inertial waves (or inertial-gravity waves when buoyancy is important). In this paper, we perform inviscid three-dimensional global hydrodynamic simulations to examine the growth and consequence of this instability operating on the spiral waves driven by a Jupiter-mass planet in a protoplanetary disk. We find that the spiral waves are destabilized via the spiral wave instability (SWI), generating hydrodynamic turbulence and sustained radially-alternating vertical flows that appear to be associated with long wavelength inertial modes. In the interval , where denotes the semi-major axis of the planetary orbit (assumed to be 5 au), the estimated vertical diffusion rate associated with the turbulence is characterized by . For the disk model considered here, the diffusion rate is such that particles with sizes up to several centimeters are vertically mixed within the first pressure scale height. This suggests that the instability of spiral waves launched by a giant planet can significantly disperse solid particles and trace chemical species from the midplane. In planet formation models where the continuous local production of chondrules/pebbles occurs over Myr time scales to provide a feedstock for pebble accretion onto these bodies, this stirring of solid particles may add a time constraint: planetary embryos and large asteroids have to form before a gas giant forms in the outer disk, otherwise the SWI will significantly decrease the chondrule/pebble accretion efficiency.

hydrodynamics, instabilities, planets and satellites: formation, planet-disk interaction, waves

1 Introduction

When they form, planets leave traces of their presence in the disks they reside in. The most suggestive signatures include spiral arms, of which we might already have observational snapshots (e.g., Muto et al., 2012; Garufi et al., 2013; Grady et al., 2013; Currie et al., 2014; Benisty et al., 2015; Garufi et al., 2016; Stolker et al., 2016), though the physical origin of the observed spiral arms as well as the presence of planet(s) in the systems are yet to be confirmed.

In our recent work, we have shown that propagating spiral density waves, such as the ones that could be excited by a planet embedded in a gaseous protoplanetary disk, or by gravitational instability (GI), can become unstable to a spiral wave instability (SWI; Bae et al. 2016). The instability arises because the periodic forcing due to the spiral waves resonantly couples with, and transfers energy into, pairs of inertial waves (or inertial-gravity waves) at the expense of the spiral wave itself, such that the spiral waves partially dissipate. When the instability tends towards nonlinear saturation, the flow breaks down into hydrodynamic turbulence, which in turn can act as an efficient source of vertical mixing of trace chemical species and solid particles (Bae et al., 2016).

The level of turbulence in protoplanetary disks plays an important role in determining the ability of solid particles to grow from ISM sizes to eventually become planets. Solid particles grow to millimeters to centimeters in size through direct sticking collisions, after which point the so-called “bouncing barrier” limits further growth by hit-and-stick coagulation (see review by Testi et al. 2014 and references therein). If decimeter-sized particles are able to form in abundance then streaming instabilities can concentrate these “pebbles” and “boulders” via aerodynamic drag, and the resulting particle clumps can collapse gravitationally to form planetesimals of between 25 to 200 km in size (Youdin & Goodman, 2005; Johansen et al., 2007, 2015). Once large planetesimals are present, gas drag acting on the remaining pebbles can produce rapid accretion – pebble accretion (Johansen & Lacerda, 2010; Ormel & Klahr, 2010; Lambrechts & Johansen, 2012; Morbidelli & Nesvorny, 2012). Pebble accretion involving millimeter-sized chondrules can help to build up planetary embryos to the size of Mars, and can also help to grow 100 km sized asteroids that form as a result of streaming instabilities into significantly larger bodies that resemble the largest asteroids in the asteroid belt (e.g., Johansen et al., 2015). While the former process – the streaming instability – has been shown to operate in the presence of moderately strong turbulence with (Johansen et al., 2007), where denotes the canonical Shakura & Sunyaev (1973) stress parameter, for the latter process – pebble accretion – to work efficiently, it is necessary to concentrate pebbles into a thin layer at the disk midplane such that the scale height of the pebbles is smaller than the Hill radius of pebble-accreting planetesimals (Lambrechts & Johansen, 2012). Under typical disk conditions, this requires the scale height of pebbles to be less than of that of gas, which can only be satisfied with very low turbulence to avoid stirring up the pebbles. A turbulence level of may result in a too long timescale for the formation of planetary embryos or large asteroids (e.g., Johansen et al., 2015).

In the present paper, we perform three-dimensional global hydrodynamic simulations to demonstrate that the SWI operates in the presence of spiral waves excited by a Jovian mass giant planet, and to estimate the level of SWI-induced turbulence and the consequent vertical mixing of solid particles in the terrestrial body-forming and asteroid belt regions. Our simulations show that the SWI develops for the spiral waves on timescales on the order of the planetary orbital time. When the instability is fully saturated, the vertical diffusion rate estimated from gas motions in the interval , where is the semi-major axis of the planetary orbit (assumed to be 5 au in this work), is such that particles with sizes up to a few centimeters are vertically mixed within the first pressure scale height of the gas disk. This result suggests that the instability acting on the spiral waves from a gas giant can have significant influence on the dust dynamics and chemical mixing in protoplanetary disks. In particular, if accretion of chondrules/pebbles dominates the growth of terrestrial embryos in the terrestrial planet region, or large asteroids in the asteroid belt, then we suggest that the SWI can have a strong influence by limiting chondrule/pebble accretion efficiency when a gas giant planet such as Jupiter forms in the outer disk.

This paper is organized as follows. In Section 2, we discuss the expected wavelengths of the unstable inertial waves and the regions of a disk in which the SWI operates, based on the WKBJ dispersion relations. We describe our computational setup in Section 3, with a special emphasis on the choice of numerical resolution that suits for capturing unstable inertial modes. In Section 4, we begin by describing numerical results obtained with perturbations having monochromatic azimuthal modes with and , with which the spiral arms are evenly spaced in azimuth. We then describe results obtained with a full planetary potential. We present the analysis on the vertical mixing of particles arising from the SWI and discuss implications to the formation and growth of terrestrial bodies and large asteroids in Section 5, and draw conclusions in Section 6.

2 Theoretical Expectations

2.1 Monochromatic spiral waves

Before we consider spiral density waves launched by planets, which consist of a superposition of various azimuthal components (Goldreich & Tremaine, 1979), we briefly discuss theoretical expectation for the SWI driven by monochromatic and spiral waves to gain insight into the instability. These modes are chosen because they appear prominently in our simulations with Jupiter-mass planets, as discussed below. We focus mainly on spiral waves that propagate in the disk regions that lie interior to the planet in this paper, and leave discussion of outward propagating spiral waves for future work. Throughout the paper, we use the term monochromatic modes/waves to denote the spiral waves that are evenly spaced in azimuth. For a more detailed discussion of the theoretical background, we refer readers to Section 2 of Bae et al. (2016) and references therein.

In a differentially rotating disk, the WKBJ dispersion relation for local disturbances can be written as (Goodman, 1993)


In the dispersion relation, is the mode frequency, is the epicyclic frequency, is the vertical Brunt-Väisälä frequency, and are the radial and vertical wave numbers associated with the wave vector , and is the sound speed.

The WKBJ dispersion relation for a spiral wave with azimuthal mode number is given by


when self-gravity is neglected and is assumed. Interior to the inner Lindblad resonance (ILR), this gives a radial wave number of


where is the Doppler-shifted frequency of the incoming waves, and is the pattern speed of the spiral wave measured in the inertial frame.

Inertial modes excited by a spiral wave must have radial length scales that are similar to or smaller than the wavelength of the incoming wave. Thus, assuming that the excited inertial modes have very similar spatial structure and frequencies, which is appropriate to the high wave number limit, we can relate the wave number of the excited inertial waves to as


where is an integer. As noted in Bae et al. (2016), this integer relationship holds for waves in periodic shearing boxes, and is introduced for the global models that we consider here for the purpose of choosing a simple relation between the wavelengths of the spiral and inertial waves. Given the radial wave number and the frequency of the spiral wave (in the local fluid frame), we can estimate the vertical wave number of the excited inertial waves using the dispersion relation in Equation (1) and the resonance condition in the high wave number limit (; Fromang & Papaloizou, 2007). This results in


when the Brunt-Väisälä frequency, , is important or


when is negligible. In order for to be real for arbitrary values of from Equation (6), such that inertial waves can be excited, one requires


to be satisfied.

In rotating disks, inertial modes have frequencies in the interval . This, together with Equation (7), implies that for an arbitrary perturbation with azimuthal mode number the SWI operates in the radial region of


in the inner disk and


in the outer disk. Focusing on the inner disk, there will always be inertial mode pairs available for a resonant interaction with spiral waves having a Doppler-shifted frequency inside the ILR. On the other hand, for waves, the SWI cannot operate in disk regions where , because the spiral wave frequency at these radii is too large, so no inertial modes can participate in a resonant interaction.

Figure 1: (a) Contour plots of the vertical wavelength of the unstable inertial modes , in units of local scale height, calculated with Equation (5) and the disk model that will be introduced in Section 3.2. (b) Same as (a) but for the radial wavelength . (c) The ratio of vertical to radial wavelength of the inertial modes. In all panels, the black and white dashed lines indicate where and . The black regions are where the dispersion relation does not have a physical solution, and thus the SWI is believed to be forbidden. The red dashed curves indicate where the local buoyancy frequency equals to a half of the doppler-shifted frequency of the spiral waves: . The upper half of each panel presents the wavelengths of unstable inertial modes in response to monochromatic waves, whereas the lower half of each panel shows those in response to monochromatic waves. We use as a representative example purely for the purposes of illustration, but the numbers on the colorbars can be linearly scaled to other modes as noted in the labels.

In Figure 1, we present two-dimensional maps showing the vertical and radial wavelengths of the inertial modes that satisfy the resonance conditions for monochromatic and spiral waves, assuming an adiabatic response of gas with and the disk model that will be introduced in Section 3.2. The wavelengths are calculated using Equations (3), (4), and (5).

We point out two important features from the figure. First, as we have shown above, there is a limited region where the SWI can be triggered for a given monochromatic perturbation. At the midplane, where the Brunt-Väisälä frequency becomes zero, the SWI can operate in the entire region inside of the ILR for waves (), whereas it can only operate in the region for waves. Second, when buoyancy is considered, the SWI-permitted regions are confined towards the midplane at the radii where the SWI can operate because the SWI can operate only in disk regions where is satisfied (Bae et al., 2016). The vertical extent of this SWI-permitted region is above and below the midplane, and this is where we will focus on in this paper. While modes can trigger the SWI near the surface region () at because of buoyancy, and this may have important implications for the morphology of spiral arms traced by observations of infrared scattered light, we do not consider this further here as we are mainly concerned with the mixing of solid particles located near the midplane.

2.2 Planet-induced spiral waves

A Jovian planet on a circular orbit will excite a spiral wave pattern that is a superposition of components with different azimuthal mode numbers, but with each component displaying a pattern speed, , that is equal to the Keplerian angular velocity of the planet. Nonlinear effects may cause the relation between the different components to vary as a function of radius. As shown later in the paper via a Fourier analysis, the spiral wave pattern at any radius interior to the Lindblad resonances due to a Jovian planet can be decomposed into a sum of modes where the dominant contributors are the and modes. At any point in time, one can assume that a protoplanetary disk will support the whole spectrum of possible inertial modes, which, in the absence of a strong excitation mechanism, will be present with very low amplitudes. The presence of a spiral wave pattern consisting of different azimuthal components, each with different Doppler-shifted frequencies as observed by fluid elements orbiting in the disk, will lead to excitation of those inertial modes that can resonantly interact with any of the azimuthal mode components that make up the incoming spiral wave. As such, we expect that a giant planet-induced spiral wave will lead to the excitation of inertial modes in a manner that is similar to a superposition of the response to both the and spiral modes described above. For lower mass planets, where the strengths of the spiral wave components will not be so strongly concentrated towards the and modes, the SWI may also be excited by the higher- modes.

Although our focus in this paper is on planets with circular orbits, it is worthwhile briefly discussing how the above picture changes if orbital eccentricity is considered. If the planet is on an eccentric Keplerian orbit, then -fold spiral waves are excited at Lindblad resonances with pattern speeds , where is a positive integer and is the mean motion of the planet (Goldreich & Tremaine, 1980). When considering planet-disk interactions at Lindblad resonances, the strength of the torque exerted on the disk by the planet, , scales as , where is the orbital eccentricity. Therefore, increasing the planet eccentricity gives rise to new spiral waves, associated with higher values of , being excited, with their amplitudes increasing as the eccentricity increases. Considering Lindblad resonances located interior to the planet that launch inward propagating spiral waves, these additional spiral waves can have higher pattern speeds than the equivalent -fold spirals associated with planets on circular orbits. For example, taking and , we have a wave that is excited with pattern speed instead of as is the case for a planet on a circular orbit. Clearly, introducing the possibility of eccentric orbits increases the range of inertial mode frequencies that can be resonantly excited.

Figure 2: Two-dimensional distributions of (left) optical depth and (right) the dimensionless cooling time for the initial disk. The red dashed curves in the right panel indicate where . The main body of the disk is expected to behave adiabatically, since the cooling timescale is much longer than the dynamical timescale in the region.

3 Numerical Methods

3.1 Basic Equations

We solve the hydrodynamic equations for mass, momentum, and internal energy conservation in the three-dimensional spherical coordinates :


In the above equations, is the mass density, is the velocity, is the pressure, is the gravitational potential from the central star, is the external potential (see below), is the internal energy per unit volume, and is the cooling rate (see below).

In the case of monochromatic spiral waves, we implement the following potential form:


Here, determines the spiral wave amplitude which is assumed to be constant over time , is the azimuthal mode number, is the pattern speed, is the cylindrical radius, is the radius about which the potential is centered, and is the radial width of the potential. We assume that the pattern speed is the local Keplerian frequency at its central position and that .

When a planetary companion is considered, its potential is included as


where is the planetary mass, and are the radius vectors of the center of grid cells in question and of the planet, and is the smoothing length. We increase the planetary mass from zero to its full mass (i.e., ) over ten orbital times. In three-dimensional calculations, the smoothing length is used only to avoid singularities in the potential on the grid scale. We thus adopt the cell diagonal size at the position of the planet as the smoothing length. With the disk model introduced in Section 3.2, the smoothing length corresponds to about of a scale height, or about of the Hill radius when the planet is fully grown to . In the planet run, we also include the indirect potential that arises because the origin of the coordinate system is based on the central star and not the center of the mass of the system.

We make use of an adiabatic equation of state. The gas pressure and the internal energy are thus related through , with an adiabatic index adopted. To realize the radiative cooling of the disk we implement a simple, but physically motivated, cooling scheme. This assumes relaxation of the internal energy towards the background disk temperature at each location on the cooling timescale . Then, the cooling rate can be written as


where denotes the heat capacity at constant volume.

The cooling timescale is calculated every time step for each grid cell, using the optical depth and the temperature. We follow the approach used in Lyra et al. (2016) and previously in Lyra et al. (2010) and Horn et al. (2012) in a vertically-integrated two-dimensional approximation. To briefly summarize, the cooling timescale is estimated to be the radiative timescale that is defined as


where and with and being the Stefan-Boltzmann constant and the effective optical depth. The effective optical depth is given by


to take account of both optically thick and thin limits (Hubeny, 1990; D’Angelo et al., 2003).

The integration in Equation (16) is done over a sphere that has a radius equal to the local pressure scale height. This results in a cooling time of


To obtain the optical depth , we first calculate the optical depth due to the material above and below each grid cell as




In practice, the integration is done in , instead of , for simplicity. Then, we calculate the optical depth as to give a correct midplane cooling rate (Lyra et al., 2016). For the opacity in the above equations, we use the Rosseland mean opacity of Zhu et al. (2009).

The optical depth and the cooling time of the initial disk are presented in Figure 2. As shown, the main body of the disk is optically thick, with the cooling time much longer than the dynamical time: . Therefore, we expect that the main body of the disk behaves essentially adiabatically, and hence will be stable against the growth of the vertical shear instability (Nelson et al., 2013).

Our calculations are inviscid, but artificial viscosity and the associated heating are included in the momentum and internal energy equations (Stone & Norman, 1992). As discussed in Bae et al. (2016), the effective dimensionless kinematic viscosity associated with the numerical diffusion for FARGO3D (Benítez-Llambay & Masset, 2016) operates with a value that is below .

3.2 Disk Models

We begin with an initial radial power-law temperature distribution in the disk that is independent of height:


where is the temperature at the location of the perturber . The isothermal sound speed is related to the temperature by , where is the gas constant and is the mean molecular weight. Thus, Equation (21) corresponds to the radial sound speed distribution given by


In all models, we adopt and such that the entire disk has the aspect ratio of 0.05.

The initial density and azimuthal velocity profiles are constructed to satisfy hydrostatic equilibrium (e.g., Nelson et al., 2013):




We choose the initial density distribution in such a way that the vertically-integrated surface density becomes , where is the minimum mass solar nebular (MMSN) model of Hayashi (1981). A value of is adopted accordingly, to have the surface density power-law slope that matches to the slope of the MMSN model, .

The initial radial and meridional velocities are set to zero, but uniformly distributed random perturbations are added as white noise, at the level of , to the initial meridional velocity field in order to seed the instability.

3.3 Computational Setup

We expect the main body of the disk to behave adiabatically, since the disk cooling timescale is much longer than the dynamical timescale in the region (see Figure 2). Thus, the regions in which the SWI is allowed to operate can be inferred using the linear analysis described in Section 2 and depicted in Figure 1, assuming an adiabatic equation of state with .

As seen in Figure 1, the SWI-permitted regions are confined around the midplane, with a thickness that is as small as above and below the midplane. Therefore, we aim to have a numerical resolution with which unstable inertial modes with can be well captured. At the midplane, the ratio of the vertical to radial wavelengths of unstable inertial modes increases over radius, except near the ILR where a singularity appears. The largest ratio of for perturbations with is (see Figure 1), and thus the radial wavelength that we aim to resolve will be .

Since inertial modes can have arbitrarily small length scales, smaller scale unstable modes can only be captured with higher resolutions than the one used in the present work. On the other hand, it seems that the modes with the largest length scale contain the largest amount of kinetic energy (Bae et al., 2016). In the simulations presented in Section 4 as well as the ones in Bae et al. (2016), we observe that smaller scale modes grow at earlier times, but only with small velocity amplitudes. These modes are then masked by larger scale modes that grow at later times with larger velocity amplitudes, supporting the conjecture that the modes with the largest length scales contain most of the kinetic energy as the system approaches the nonlinear saturated state. One effect of under-resolving the small scale modes, however, is that the development of a realistic breakdown into hydrodynamic turbulence may be prevented, as the cascade of energy from large to small scales in the fully saturated nonlinear regime will be prevented, affecting the nature of the flow and the statistics of the quasi-turbulent flow that develops due to the SWI.

Our simulation domain extends from to in radius, from to in the meridional direction (covering 4 scale heights above and below the midplane), and from 0 to in azimuth. We have tested with various numerical resolutions and find that 12 grid cells per wavelength is required to properly resolve unstable inertial modes. To achieve , or equivalently at the midplane, we adopt 726 logarithmically-spaced radial grid cells. We adopt 144 and 754 uniformly-spaced grid cells in the meridional and the azimuthal directions, respectively, with which choice .

Figure 3: (Left Top) Two-dimensional distribution of density perturbation in a plane for model. The two white curves near the midplane indicate where the local Brunt-Väisälä frequency equals to one half of the doppler-shifted wave frequency of spiral waves (), between which region is expected to be unstable to the SWI. (Left Bottom) Contour plots of the meridional velocity normalized by the local sound speed, for the dotted rectangle region in the left top panel. The and axes are drawn isotropically so that the radial and vertical wavelengths of unstable inertial modes can be compared with each other from the figure. The black dashed lines indicate where . Note that the SWI has grown in the entire disk region inside the ILR. (Right) Time evolution of the meridional kinetic energy at various radius bins, normalized by in each bin. Note that the exponential growth of the instability appears at earlier times for smaller radii.

The radial boundary condition is chosen to have a zero gradient for all variables. We further implement a wave damping zone (de Val-Borro et al., 2006) in the intervals and , since reflected waves at the radial boundaries may affect triggering and growth of the SWI, as they have the same pattern speed as the incoming waves. Periodic boundary conditions are used in azimuth since the simulation domain covers . At the meridional boundaries, we use an outflow boundary condition such that all velocity components in the ghost zones have the same values as the last active zones, but the meridional velocity is set to 0 if directed toward the disk midplane. The temperature in the ghost zones is set to have the same value as in the last active zones. The density in the ghost zone is then obtained by solving the hydrostatic equilibrium in the meridional direction:


We make use of FARGO3D (Benítez-Llambay & Masset, 2016), which runs on both central processing units (CPUs) and graphics processing units (GPUs). In order to deal with the high numerical resolution, our calculations use a cluster of GPUs, which allows us to perform the calculations within a reasonable time frame thanks to a large speed up with respect to CPUs. With four NVIDIA Tesla K20x GPUs on the University of Michigan high-performance computing cluster4, it takes about one month for the planet run in Section 4.2 to evolve for 200 orbits. We enable the FARGO (Fast Advection in Rotating Gaseous Objects) orbital advection module (Masset, 2000).

In the following sections we will use the orbital distance of the planet (or perturbing potential) as the length unit (assumed to be 5 au), and the orbital time at as the time unit.

Since we adopt physical units to consider cooling, the models are not scalable in principle. However, we believe that the results can be qualitatively applicable as far as the main body of a disk is sufficiently optically thick, because the disk response to spiral waves is not very sensitive to the cooling timescale when (Zhu et al., 2015, and also from our test runs during early phases of the present work).

Figure 4: Same as Figure 3, but for model. Note that the instability has grown beyond , but not inward of the radius as no inertial modes there have an adequate frequency to satisfy the resonant condition with the incoming spiral waves. Also, the exponential growth in appears for and but not at the inner two radius bins, supporting the fact that the SWI does not operate at . The increase in in the inner two radius bins at later times () is presumably because the spiral waves propagate through the SWI-turbulent region.

3.4 Diagnostics

In order to examine the growth and the saturation of the SWI, we compute the volume-integrated meridional kinetic energy :


While the velocity perturbation driven by spiral waves is larger at high altitude than near the midplane in general, traces perturbations near the midplane region because of the vertically stratified density structure.

In the figures presenting the meridional kinetic energy, we normalize by which we define as the volume-integrated kinetic energy of gas parcels assuming that they move at the local sonic velocity:


We also compute the Shakura-Sunyaev stress parameter as


where is a density-weighted mean pressure at radius , in order to examine the sustained level of angular momentum transport driven by the turbulent flow at saturation of the instability (although we note that a contribution to also arises from the propagating spiral waves that is difficult to disentangle from that arising from the SWI-induced turbulence).

4 Results

4.1 Results with Monochromatic Waves

Before we describe our main results with a planet, we introduce models with monochromatic perturbations. The perturbing waves are imposed using Equation (13) with azimuthal mode numbers (Section 4.1.1; hereafter model) or (Section 4.1.2; hereafter model). The potential amplitudes are chosen so that the linear phase of the instability spans several orbital times, and therefore the growth of individual unstable modes can be captured: for the model and for the model.


In Figure 3, we present the density distribution for model, where the brackets denote an azimuthal average, along with the meridional velocity distribution near the midplane where the SWI is expected to operate. The density perturbation is about in the midplane and about at the surface when the spiral waves are fully established. The larger perturbation towards the surface is presumably because of a nonlinear effect due to faster advection of the wave at higher altitudes in the disk (Bae et al., 2016, see also Zhu et al. 2015).

The checkerboard pattern shown in the meridional velocity distribution is a generic feature of the linear growth phase of the SWI (Bae et al., 2016), and was also observed by Fromang & Papaloizou (2007) in their shearing box simulations of nonlinear, axisymmetric sound waves propagating in astrophysical disks. As expected from the dispersion relations, the regions in which the SWI develops are confined toward the midplane, because the large Brunt-Väisälä frequency near the surface does not allow any inertial modes there to resonantly interacting with the incoming waves. During the linear phase, the magnitude of the perturbed meridional velocity increases over time, but remains a few to about ten percent of the local sound speed. The wavelengths of the dominant unstable inertial modes varies over radius: at the midplane, at , whereas and at . The ratio of vertical to radial wavelengths of the unstable inertial modes increases toward larger radii and height. Note that these unstable mode properties are well described by the linear analysis introduced in Section 2.

Also shown in Figure 3 is the time evolution of the volume-integrated meridional kinetic energy for various radius bins. The exponential growth of , for example during at , indicates the development of the SWI. At smaller radii the instability becomes apparent at earlier times because of the shorter interaction period between inertial modes and the spiral waves.


Figure 4 displays the density distribution for the model. The density perturbation is about in the midplane and about at the surface when the spiral waves are fully established. As seen in the figure, the instability is developing beyond at . The SWI does not operate inward of this radius because the Doppler-shifted perturbation frequency from the incoming waves is too large to excite any inertial modes, which is in very good agreement with the expectation from the dispersion relations. As in the model, the ratio of vertical to radial wavelength of unstable inertial modes increases towards larger radii and height.

The time evolution of is presented in Figure 4. The exponential growth of is only observed in and , supporting the conjecture that the instability is forbidden at smaller radii. The increase in at later times (e.g., ) in the inner disk regions of is probably because the spiral waves propagate through the turbulent region of , and hence generate some vertical motion there due to the curvature of the spiral wave fronts in the meridional plane discussed in Bae et al. (2016).

4.2 Results with a Perturbing Planet

We now describe the results for the planet run. We first discuss the overall evolution in Section 4.2.1 and then focus on the development of the SWI in Section 4.2.2.

Figure 5: Radial distributions of the azimuthally averaged surface density at , 40, 80, 120, 160, and .
Figure 6: Two-dimensional distribution of in the midplane, taken when the planetary mass is fully grown to at . The planet is located at . The -axis is plotted using a logarithmic scale to stretch the inner disk.

Overall Evolution

We evolve the planet run for . In Figure 5, we present the azimuthally averaged surface density distributions at every . The planet opens a gap around its orbit and creates a mild density bump at the inner gap edge. In the inner disk (), the disk looses material through the meridional boundary because of vertical flows induced by the spiral waves.

As inferred from the figure, the depth of the gap does not reach full saturation by the end of the run at . Recently, Fung & Chiang (2016) carried out three-dimensional calculations to examine gap opening by planets in a locally isothermal disk. Their results indicate that the saturation of the gap depth requires about 1000 planetary orbits (or even longer depending on the planetary mass) with a disk viscosity of . On the other hand, in their case, the density perturbation driven by spiral waves remains roughly constant after 100 planetary orbits (Fung, D., private communication). It is therefore reasonable to assume that the amplitudes of the spiral waves in our simulations have also saturated.

In Figure 6, we present the density distribution in the midplane at , at which time the planetary mass is fully grown to . The planet excites three distinguishable spiral arms. This multi-armed feature can be understood as a consequence of non-linear mode-coupling (e.g., Artymowicz & Lubow, 1992; Lee, 2016) and has been shown to arise in recent numerical simulations (e.g., Fung & Dong, 2015; Juhász et al., 2015; Zhu et al., 2015). Looking in more details at the spiral arms, the first arm is connected to the planet, and the second and third arms start to appear at and at , respectively. The azimuthal separation between the spiral arms, as well as the relative strength of the arms, vary over radius. At , for instance, the second arm is ahead of the first arm and the third arm is behind the first arm. The density perturbation by the second arm at this radius is about twice as strong as the perturbation driven by the other two arm. On the other hand, at , the second arm is ahead of the first arm and the third arm is behind the first arm. The density perturbation by the second arm is still the strongest, but only stronger than the perturbation driven by the first arm.

Figure 7: Time evolution of the Fourier amplitudes for various azimuthal modes () in different radius bins. Note that the amplitudes remain nearly constant after .

In order to more quantitatively examine the gas response to the spiral waves, in Figure 7 we plot the time evolution of the Fourier amplitudes of the density perturbation at different radii. As seen in the figure, spiral waves driven by a planet are not monochromatic, but consist of various azimuthal components that are superimposed on each other. We emphasize that, most importantly for our purpose, the Fourier amplitudes of modes remain nearly constant after at all radii, despite the fact that the gap depth has increased by about a factor of two in between and . Another important feature is that the strongest mode shifts toward smaller azimuthal wave numbers at smaller radii. The mode has the largest Fourier amplitude at , the and modes have a comparable amplitude at , and the dominates at . For completeness, we note that there are no spiral modes that propagate interior to the planet because formally the ILR is not present in the disk. The appearance of a finite amplitude Fourier component arises presumably because the disk generates an perturbation due to the appearance of a low amplitude vortex, or because the disk becomes mildly eccentric.

Figure 8: Two-dimensional distributions of the meridional velocity in a plane at (left) and (right) . The two sets of white curves indicate where the local Brunt-Väisälä frequency equals to one half of the Doppler-shifted wave frequency, for monochromatic waves (the ones close to the midplane) and waves (the ones close to the surface). The vertical dashed line indicates the radial location inward of which the SWI is forbidden for monochromatic waves. The black dashed lines denotes where . The two axes are drawn isotropically.
Figure 9: (Left) Time evolution of the meridional kinetic energy at various radius bins, normalized by in each bin. (Right) Time evolution of the Reynolds Stress at various radius bins. In both panels, we add vertical offsets of 0.001, 0.002, 0.003, and 0.004 to the blue, green, yellow, and red curves for illustration purposes.

Spiral Wave Instability

The SWI starts to grow from very early times even before the planetary mass is fully grown () at , because the perturbation is strong enough there even with a sub-Jovian mass planet. The weak, wave-like inter-arm structures at and the break-up of the second arm at seen in Figure 6, are some of the indications of the instability.

In Figure 8, we display the meridional velocity field in a vertical plane at and . In the snapshot taken at , one can see the checkerboard pattern suggesting that the SWI is growing there, although the individual unstable modes are not as clearly identifiable as in the monochromatic case due to the complex nature of planet-driven waves (i.e., superposition of various azimuthal modes). However, we note that there is no signature of the SWI in at this time epoch. This is probably because (1) the perturbation in the region is small compared with the region (Figure 6) and/or (2) the component, with which the SWI is forbidden in this region, has a comparable power to the mode during the early evolution (e.g., region in Figure 7) so the growth rate of the instability could be reduced. At , checkerboard patterns appear in the inner disk region of . The steep increase in the meridional kinetic energy around , particularly for the region, presented in Figure 9 supports the suggestion that the SWI is developing in the region.

Figure 10: (Top) Three-dimensional global view of the meridional velocity field in the bottom half of the simulations domain: the upper surface shows the disk midplane. The quantity is plotted in units of the local sound speed: . on the X-axis, on the Y-axis in units of radians, and on the Z-axis in units of radians. Note that the instability creates radially-alternating vertical flows when saturated. The vertical flows have a magnitude of order of a few tens of percent at the midplane. (Bottom) A volume-rendering view of the density in a logarithmic scale. The snapshots are taken at the end of the simulation (). The planet is located at .

As shown in Figure 9, the overall perturbed vertical kinetic energy and the stress associated with the combined action of the propagating spiral waves and the SWI (we have not sought to disentangle the contributions from these in this paper) reach quasi-steady state after . This adds further evidence that the strength of the SWI is probably not very sensitive to the gap depth after this point as the spiral wave amplitudes have reached steady values. We note that increases after for . However, we do not find any significant changes in the strength of spiral arms or background disk structures to explain this rise. One possible explanation is that some large-scale inertial modes, which have to have smaller growth rates than small-scale modes (Fromang & Papaloizou, 2007), become unstable at this time, although we were not able to identify any corresponding individual mode growing as the disk has already developed strong vertical flows by this time, which we will discuss now.

In Figure 10, we present a three-dimensional global view of the meridional velocity and the density at the end of the simulation (). As shown, the checkerboard meridional velocity patterns merge to create radially alternating vertical flows when the SWI saturates. Such vertical flows have been pointed out in the non-stratified, isothermal cylindrical models of Bae et al. (2016), in which other sources capable of inducing vertical motion, such as the vertical shear instability that is known to cause the excitation of corrugation modes (Nelson et al., 2013), are absent. We thus believe that the alternating vertical flow is a generic outcome of the SWI. Visual inspection of the upper panel of figure 10 suggests the azimuthal structure of these modes is of low degree, corresponding to perhaps a mixture of corrugation-type modes and/or warp or tilt modes, which are inertial modes for which the vertical velocity perturbations have no nodes in the vertical direction. In the case of corrugation modes, the whole disk column at a given radius, including the midplane, moves up and down in a coherent manner. For tilt modes the disk acts as a series of narrow rings each of which tilts rigidly. As these disturbances arise in the form of waves, they propagate through the disk as corrugation or warping disturbances. The vertical length scale of the vertical flows is similar to the thickness of the disk region in which the SWI is permitted, suggesting that it is set by the wavelength of the largest unstable inertial modes that are allowed to grow. So these might be modes which are excited directly by the SWI, or they may be the result of nonlinear mode coupling which allows energy to seep into these modes from other modes that are excited by the SWI. Being coherent and quite large scale, once excited these corrugation/tilt modes become a prominent feature of the flow because they do not damp efficiently. We find that the associated vertical flows can have perturbed vertical velocities that are in the range of the local sound speed. One caveat, however, which we have already alluded to in Section 3.3 and which causes us to be cautious in interpreting our results, is that a resolution of cells per scale height in the vertical and radial directions does not allow for a turbulent cascade to develop efficiently. While some transfer of energy to small scales undoubtedly arises in our simulations, leading to dissipation on the grid scale, it may also be the case that the low resolution also favors the development of the SWI in such a manner that the large scale modes are more prominent than they would be in a more highly resolved simulation. Testing the outcome of the SWI as a function of resolution is a task that will be undertaken in the future when available computational resources allow such a study to be conducted.

5 Discussion

5.1 Particle Stirring Induced by the SWI

We measure the vertical diffusion coefficient to estimate the rate of vertical mixing of dust particles induced by the SWI. In order to do this, we restart the planet run described in Section 4.2 from , with outputs of the meridional velocity in three dimensions every . We calculate the vertical diffusion coefficient using the approximation , where is the correlation time of the vertical velocity fluctuations, (Fromang & Papaloizou, 2006). In practice, we use the meridional velocity rather than . The quantity represents the ensemble and time average of the mean velocities calculated at all grid cells at the cylindrical radii listed below. Using this definition for the diffusion coefficient implicitly assumes that the turbulence properties are uniform at all heights for each value of .

In obtaining an estimate for , we generate a time series for at various radii in the disk ( and 0.7), within . We then compute the autocorrelations of the time series and fit the computed autocorrelations with a form of


following Nelson & Gressel (2010). Here, indicates the relative strength of the sinusoidal feature in the autocorrelation function, is the frequency associated with the sinusoidal component, and is the correlation time. Figure 11 shows an example of the fit at . For this example, the correlation time is or . While it varies over radius, note that the correlation time is of order of , which is longer than the ones obtained for MHD turbulence in protoplanetary disk simulations (; Fromang & Papaloizou, 2006; Yang et al., 2009; Nelson & Gressel, 2010). We suspect that this longer correlation time, combined with the prominent oscillatory component displayed by the autocorrelation function, may arise because of the contribution of the coherent radially-alternating flows that arise in the simulations as described above, such that the resulting flow is a mixture of coherent vertical motion and smaller scale turbulence. The correlation time obtained for various radii is listed in Table 1.

With the vertical velocity fluctuation that is of the local sound speed, the vertical diffusion coefficient is in natural units between and . Then, the vertical mixing time can be estimated as


Using the prescription


the vertical diffusion rate is characterized by values of in the range , which is comparable to or greater than the Reynolds stress measured at these radii (see Figure 9).

() () () () (cm) (cm)
0.3 0.050 0.43 2.62 0.0065 0.62 2.49 4.63
0.4 0.059 0.55 2.17 0.0075 0.90 4.64 3.20
0.5 0.055 0.51 1.44 0.0043 2.06 5.96 1.75
0.6 0.112 0.46 0.99 0.0124 0.95 9.28 3.99
0.7 0.054 0.41 0.70 0.0020 7.34 8.71 0.82
Table 1: Planet Run Results

Since we are interested in the maximum size of particle that can be vertically mixed by the SWI-induced turbulence, we calculate the settling time of particles and compare it with the vertical diffusion time. The settling time can be obtained from the terminal velocity by equating the drag force and the vertical component of the stellar gravity. The drag force depends on the size of the particle , relative to the mean free path of the gas , where is the mean molecular weight of gas, is the mass of a Hydrogen atom, is the gas mass density, and (Chapman & Cowling, 1970) is the molecular collisional cross section. When , the Epstein law is applicable and the drag force can be written as


where is the gas velocity and is the thermal velocity of gas molecules. When , we can apply the Stokes law and the drag force can be written as


where is the drag coefficient. We use (Probstein & Fassio, 1969; Whipple, 1972), since the Reynolds number is measured to be always smaller than unity in the region and for the particle sizes of interest. Here, and is the gas molecular viscosity.

By setting , where is the mass of a solid particle and is the internal density of the dust particles, one can obtain the settling time . For the Epstein regime, this gives


whereas for the Stokes regime,


It is immediately obvious from Equations (34) and (35) that the equations result in the same settling time when . Also, note that Equations (34) can be written as , or , where the stopping time for the Epstein regime is defined by (Weidenschilling, 1977) and .

Then, the maximum particle size that is significantly mixed by the turbulence (i.e., ) can be computed by setting Equation (30) to either Equation (34) or (35). This leads to


As seen in Table 1, in the interval , we find that solid particles  cm in size can be vertically mixed within the first pressure scale height.

Figure 11: The normalized autocorrelation function (ACF) at from the planet run with diamond symbols, and the fit obtained with Equation (29) with solid curves.

We note that the diffusion coefficient obtained from our simulation is for gas molecules. For solid particles that are not perfectly coupled to gas, the diffusion coefficient is (Youdin & Lithwick, 2007), where is the Stokes number. In the disk assumed here, however, centimeter-sized particles have within the first scale height throughout the entire inner disk so . The correlation time, vertical diffusion coefficient, the mixing time, and the maximum particle size that can be vertically mixed by the turbulence obtained for the planet run are summarized in Table 1.

The above calculations relating to particle mixing assume that the SWI generates homogeneous hydrodynamic turbulence at each radius considered, but as discussed already the simulations also show strong evidence for there being coherent large scale vertical motions that might also contribute to lofting particles away from the midplane. Particles can only be advected within the coherent vertical flows if the vertical drag force that they feel due to the moving gas exceeds the vertical force of gravity. The maximum particle size for which this is true can be estimated by equating these forces – –, assuming that the vertical gas velocity is equal to its characteristic value measured in the simulations, which is of the local sound speed. Evaluating this at , corresponding to 3 au in physical units, we obtain a maximum particle size of  cm that can be lofted, consistent with the value obtained from the calculations above within a factor of few.

5.2 Implications for the Growth of Large Asteroids and Terrestrial Planet Embryos

Through a combination of the streaming instability producing planetesimals with sizes up to few  km (Youdin & Goodman, 2005; Johansen et al., 2007, 2015) and subsequent pebble accretion onto the larger planetesimals (Johansen & Lacerda, 2010; Ormel & Klahr, 2010; Lambrechts & Johansen, 2012; Morbidelli & Nesvorny, 2012), it seems plausible that, at 5 au in a proto-solar nebula-like gaseous disk, planetesimals can grow to become 10 Earth-mass () solid cores within about one million years (Bitsch et al., 2015; Levison et al., 2015; Morbidelli et al., 2015). This presumably allows sufficient time for such cores to accrete a substantial amount of gas to eventually form gas giants, assuming that another million or so years is required for completion of the gas accretion phase (e.g., Pollack et al., 1996; Movshovitz et al., 2010), and that protoplanetary disks around T Tauri stars dissipate in a few million years (Haisch et al., 2001; Hernández et al., 2007; Mamajek, 2009). While the details of how gas giant planets accrete their gaseous envelopes remain uncertain, it is clear that they must form before protoplanetary disks lose a significant fraction of their gas content (but see Tanigawa & Tanaka 2016).

In the Solar System, Jupiter is generally thought to have been the first planet that finished accretion. By the time Jupiter has formed, it is commonly assumed that terrestrial bodies have grown to become Moon or Mars mass embryos (see review by Dauphas & Chaussidon 2011 and Morbidelli et al. 2012 and references therein). Planetary embryos then complete their growth during the gas-free debris phase through collisions with planetesimals and/or other embryos, which takes several tens to hundreds million years after gas has removed.

More recently, models of the growth of terrestrial planet embryos and large asteroids through the accretion of pebbles and/or chondrules that are continuously generated over multi-Myr time scales have been presented (Johansen et al., 2015; Levison et al., 2015). In the work of Johansen et al. (2015) it was shown that the formation of planetesimals via high resolution streaming instability simulations leads to the formation of bodies with sizes up to  km, and that subsequent chondrule accretion onto these bodies over time scales of  Myr is needed to explain the presence of large asteroids such as Ceres and Vesta in the asteroid belt. Similarly, a model of planetary embryo growth that involves only the collisional accretion of planetesimals formed by the streaming instability was shown to be unable to form Mars-size embryos within the required time frame, and that the addition to the model of chondrule accretion over Myr time scales led to a dramatic increase in the efficacy of embryo formation. Levison et al. (2015) have presented a model in which embryos in the terrestrial planet region grow through pebble accretion over Myr time scales, with the pebbles being continuously generated during this time period, and show that such a model is able to explain certain features of the terrestrial planet system. The model introduces the outer solar system planets at the end of the embryo formation epoch, and subsequent evolution through giant impacts leads to final planetary systems that appear to be a good fit to the basic structure of the inner solar system, including a small mass for Mars.

The results of our present work suggest that a Jupiter-mass planet forming within the first few Myr in a protoplanetary disk can produce turbulence and vertical stirring of solid particles interior to its orbit. The influence seems to be quite significant, such that solid particles with sizes up to several centimeters can be vertically well mixed. While pebble/chondrule accretion takes advantage of the large accretion cross section produced by aerodynamic drag, which can be orders of magnitude larger than the geometric cross section of the target planetesimals, for solid particles having a scale height () that is larger than the size of the planetesimal Hill radius () the accretion rate will be reduced by a factor of (Lambrechts & Johansen, 2012). For example, assuming the planetesimal mass of  g with which , the accretion rate for particles that have same scale height as the gas () will therefore drop by a factor of . This probably results in too long a timescale to form planetary embryos and large asteroids through pebble accretion within the typical lifetime of protoplanetary disks. While this argument considers the decrease in the number density of pebbles near the midplane only, the larger relative velocity between pebbles and planetesimals as well as the non-zero orbital inclination and eccentricity of planetesimals that might arise in the SWI-driven turbulence can make the situation even worse. Thus, if pebble accretion is the dominant process by which terrestrial planet embryos and large asteroids gain their mass, the requirement for efficient pebble accretion may put a time and/or space constraint on giant planet formation: planetary embryos have to form before a gas giant forms in the outer disk, otherwise the strong stirring induced by the SWI will significantly drop the chondrule/pebble accretion efficiency. It is interesting to note that isotope measurements of Martian meteorites indicate that Mars, as a stranded planetary embryo (e.g., Chambers & Wetherill, 1998), grew rapidly within about 2 Myr after the birth of the solar system (Dauphas & Pourmand, 2011), and then halted accretion. It is interesting to speculate that the formation of Jupiter at or close to its current location may have occurred at around this time, making pebble accretion onto Mars inefficient and effectively halting its growth and that of the large asteroids.

With our initial density profile the disk is optically thick in the midplane. The main body of the disk therefore behaves fully adiabatically, confining the SWI towards the midplane. As the disk evolves the optical depth of the disk will decrease by accretion of material and/or particle growth and settling, and therefore the disk will behave more isothermally. This will allow the SWI to operate in broader disk regions in height and possibly strengthen the spiral arms with less resistance from gas pressure, although the turbulence level and the significance of vertical mixing have to be further examined with relevant disk models. If the level of turbulence via the SWI remains relatively constant, regardless of the thermal properties of the disk, or decreases over time, this implies that will decrease as the disk loses its mass over time (see Equation 36). In this case, it is possible that planetesimals and/or planetary embryos resume accreting pebbles during the later evolution as far as the disk still has some pebbles that survived rapid migration. On the other hand, when the disk becomes optically thin the vertical shear instability and/or the MRI could operate down to the midplane and provide some turbulence that may limit pebble accretion efficiency.

If a giant planet migrates over a significant distance to reach its final position, as in the “Grand Tack” model proposed for the early Solar System (e.g., Walsh et al., 2011), the SWI could have a significant effect over a broader range of disk radii than otherwise, because the growth timescale of the instability is orders of magnitude shorter than the migration time scale.

It might be possible that gas giants form as early as during the time when their host disks still gain mass from their natal clouds, as recent observation of HL Tau may suggest (ALMA Partnership et al., 2015). Whether or not such an early giant planet formation is common is still an open question. If giant planets form at such an early time, the disk may be left with a very narrow window in time to grow terrestrial bodies in the disk if pebble accretion is the dominant mechanism.

5.3 Dependence of SWI on Planetary Mass

We have focused on the SWI arising from the spiral waves excited by a giant planet in the inner regions of a protoplanetary disk in this paper. One may wonder if there is a mass requirement for triggering the SWI, and how the strength of the SWI varies with planet mass.

First, given that the instability involves energy from the spiral waves being transferred into inertial modes, it is clear that stronger spiral waves will generate higher amplitudes for unstable inertial modes. In the nonlinear saturated state, we therefore expect that the resulting hydrodynamic turbulence will be more vigorous for larger mass planets.

Second, our general picture of how the instability operates is that a protoplanetary disk hosts the full spectrum of inertial modes that are present at low amplitude. When a source of spiral waves (e.g., a planet) is present, it excites waves with a range of azimuthal mode numbers and Doppler-shifted frequencies. Then, the inertial modes that are resonant with spiral waves with specific azimuthal mode numbers are able to grow due to the periodic forcing. In the case of a giant, gap forming planet, we have shown that the dominant spiral modes that propagate in the inner disk are the and modes, presumably because these are excited at ILRs that lie outside of the low density gap region. As the planet mass is lowered, however, the gap depth reduces and we expect that higher values of will become prominent, with the amplitudes of waves associated with lower values decreasing in proportion to the planet mass. The increasing prominence of the higher spiral modes may then allow the SWI to operate in regions closer to the planet because, at the midplane, the SWI operates in the radial intervals in the inner disk and in in the outer disk. Again, these conditions come from the fact that, as the value of increases, the synodic period between orbiting fluid elements and the planet needs to increase to match the resonance condition with inertial modes in the disk.

Based on the two points discussed above, we expect that (1) the strength of SWI-driven turbulence will be weaker for lower mass planets; and (2) the region where the instability operates will be narrower and closer to the planet, as the planet mass decreases to the point that the mode is no longer strong enough to induce the SWI in the inner disk.

In order to support these two points, we have run additional simulations with planetary masses of . For these runs, we were forced to use a numerical resolution of because of the computational cost. This is a factor of two lower than our reference model and so thus the growth of the SWI is not properly captured, although we can measure spiral wave amplitudes.

The Fourier amplitudes for different planetary masses are plotted in Figure 12 as a function of azimuthal mode number. At both and , the Fourier amplitudes are larger for more massive planets in general. It is therefore reasonable to expect that more massive planets will produce larger energy injection rates into the inertial modes, and accordingly stronger turbulence through the SWI. We see that the maximum amplitude shifts towards lower degree azimuthal modes as increases, presumably due to non-linear mode coupling (e.g., Lee, 2016) and gap formation which suppresses the amplitudes of the high- modes that are excited close to the planet. This trend is particularly apparent for . As the SWI is expected to operate for and perturbations at radii (it can also be excited by the mode between , but this is only just contained in the radial range under consideration), and for modes only at radii , we expect that the SWI operates more strongly with larger planetary masses that have most power in or in these regions (). When , the Fourier amplitudes are relatively flat for . If the SWI can operate for such a low mass planet, then we might expect that the strength of the outcome will be much weaker than for larger mass planets, and that the regions unstable to the instability will be limited to a narrow radial region around the planetary orbit because of the weakness of the spiral wave. Global disk calculations performed at high resolution will be required to confirm or refute these conjectures, and to fully characterize how the SWI operates as a function of planet mass.

Figure 12: The Fourier amplitude for various azimuthal modes at (upper) and (lower) , taken at for different planetary masses of (black) , (blue) , (green) , (yellow) , and (red) . Black open circles are for our reference model with the high numerical resolution of , whereas all the other models run with low resolution of . The arrows indicate the azimuthal mode numbers with which perturbing waves are capable of exciting the SWI at each radius bin.

5.4 Caveats and Future Work

As seen in Figure 11, the ACF of the vertical velocity does not converge but oscillates around zero with finite amplitude. This is because the vertical flows developing from the SWI not only contain repeated growth and decay of vertical motions, but also display sustained and coherent oscillations, likely due to the excitation of longer wavelength inertial waves. Due to this apparently complex superposition of turbulence and coherent oscillations, the analysis presented in Section 5.1 has to be viewed with some caution by keeping in mind the uncertainty in estimating the correlation time. Nonetheless, a factor of a few uncertainty in the correlation time would not change our conclusions about particle stirring in a qualitative sense.

At this point, it is probably worthwhile pointing out that some previous simulation studies, in which grains and pebbles are implemented as Lagrangian particles, indicate evidence of vertical mixing of solid particles through the SWI. As a potential example, we note that in Figure 13 of Zhu et al. (2014), where the authors placed a 0.65 Jupiter-mass planet around a solar-mass star, particles a and b – that are about a millimeter and a centimeter in size, respectively, when the MMSN disk model and the semi-major axis of 5 au are assumed – are vertically dispersed in the outer disk ( in the figure). In particular, the particle distributions show some wiggling morphology near the midplane, which is tempting to explain with the radially alternating vertical flows arising from the SWI. The vertical dispersal of particles is not seen in the inner disk and we speculate that this is because of the insufficient numerical resolution to properly capture the unstable inertial modes there. They used uniform radial grid cells, with which one scale height is resolved with only about 5 grid cells in both radial and vertical directions at , as opposed to with about 15 grid cells at .

In the future, more quantitative conclusions will be able to be made from high resolution gas and particle simulations, where one can obtain the vertical distribution of solid particles stirred by the SWI, both interior and exterior to the planet. Furthermore, some of the questions that we have discussed in this paper such as the outcome of the SWI as a function of planet mass, and the influence that this has on pebble accretion, can also be addressed in the context of improved disk models that contain a more sophisticated model of disk thermodynamics. And finally, hydrodynamic simulations combined with radiative transfer calculations can be used to determine the observational appearance of disks with planets, and to determine what influence the SWI has on images formed through scattered light and thermal emission.

6 Conclusion

We have presented inviscid three-dimensional global hydrodynamic simulations of protoplanetary disks in which we excite spiral waves either monochromatically using an imposed external potential or using a full planetary potential. Using monochromatic waves, we first show that the region of the SWI operating and the properties of the unstable inertial modes show very good agreement with the predictions made by matching resonance criteria using simple dispersion relations for spiral and inertial waves. When a Jupiter-mass planet is added in the disk, the spiral waves it launches have various azimuthal components superimposed on the dominant or modes, depending on the radial position in the disk, and the SWI is found to grow on the order of the planet orbital time. When it is saturated, the SWI generates turbulence and radially alternating vertical flows that have characteristic radial and vertical length scales of about one local pressure scale height and vertical velocities of order of few tens percent of the local sound speed. Using the alpha prescription, the associated vertical diffusion rate of gas is estimated to be characterized by in the range of radii , where is the semi-major axis of the planet. At this rate, solid particles up to a few centimeters in size can be vertically mixed within the first pressure scale height. Since protoplanetary disks are believed to remain laminar, and thus induce no or very little particle stirring, as suggested by recent magnetized wind models (Bai & Stone, 2013; Gressel et al., 2015), the SWI can be the mechanism controlling the degree of vertical settling of solid particles in the optically-thick regions of planet-hosting disks where other hydrodynamic instabilities are not thought to operate.

While more quantitative results will be obtained with future high resolution simulations that include particles, we conjecture that significant stirring of solid particles through the instability of the spiral waves excited by giant planets can have an influence on the formation and growth of large asteroids and terrestrial planet embryos in the inner disk, if these bodies are actively growing in the presence of a giant planet (such as Jupiter in the protosolar nebula). In particular, if growth of these bodies proceeds mainly through the accretion of chondrules/pebbles, then they have to form before a gas giant forms in the outer disk as the stirring of solid particles is likely to significantly decrease the chondrule/pebble accretion efficiency. Since the properties of the turbulence driven by the SWI are dependent upon the background thermal structure of the disk as well as the planetary mass, future studies that survey parameter space at high numerical resolution are needed to further examine the significance of the instability for protoplanetary disk evolution and for planet formation.

The authors thank the anonymous referee for a prompt report and helpful comments that improved the initial manuscript. J.B. thanks Edwin Bergin, Jeffrey Fung, Wing-Kit Lee, and Zhaohuan Zhu for helpful conversations. This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. The authors acknowledge the San Diego Supercomputer Center at University of California, San Diego for providing HPC resources that have contributed to the research results reported within this paper.


  1. affiliation: Department of Astronomy, University of Michigan, 1085 S. University Ave., Ann Arbor, MI 48109, USA
  2. affiliation: Astronomy Unit, Queen Mary University of London, Mile End Road, London E1 4NS, UK
  3. affiliation: Department of Astronomy, University of Michigan, 1085 S. University Ave., Ann Arbor, MI 48109, USA
  4. http://arc-ts.umich.edu/systems-and-services/flux/
  5. footnotetext: For the quantities in the brackets , we average them over azimuth within .


  1. ALMA Partnership, Brogan, C.L., Pérez, L. M., et al.  2015, ApJ, 808, L3
  2. Artymowicz, P., & Lubow, S. H. 1992, ApJ, 389, 129
  3. Bae, J., Nelson, R. P., Hartmann, L., & Richard, S. 2016, ApJ, 829, 13
  4. Bai, X.-N & Stone, J. M. 2013, ApJ, 769, 76
  5. Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6
  6. Benítez-Llambay, P., & Masset, F. 2016, ApJS, 223, 11
  7. Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112
  8. Chambers, J. E., & Wetherill, G. W. 1998, Icarus, 136, 304
  9. Chapman, S., & Cowling, T. G.  1970, The Mathematical Theory of Non-uniform Gases. An account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases (3rd ed.; Cambridge: Cambridge Univ. Press)
  10. Currie, T., Muto, T., Kudo, T., et al. 2014, ApJ, 796, L30
  11. Dauphas, N., & Chaussidon, M. 2011, AREPS, 39, 351
  12. Dauphas, N., & Pourmand, A. 2011, Natur, 473, 489
  13. D’Angelo, G., Henning, T., & Kley, W.  2003, ApJ, 599, 548
  14. de Val-Borro, Edgar, R. G., M., Artymowicz, P., et al.  2006, MNRAS, 370, 529
  15. Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751
  16. Fromang, S., & Papaloizou, J. 2007, A&A, 468, 1
  17. Fung, J., & Dong, R. 2015, ApJ, 815, L21
  18. Fung, J., & Chiang, E. 2016, , 815, L21
  19. Garufi, A., Quanz, S. P., Avenhaus, H., et al. 2013, A&A, 560, A105
  20. Garufi, A., Quanz, S.P., & Schmid, H. M. 2016, A&A, 588, A8
  21. Goldreich, P., & Tremaine, S.  1979, ApJ, 233, 857
  22. Goldreich, P., & Tremaine, S.  1980, ApJ, 241, 425
  23. Goodman, J. 1993, ApJ, 406, 596
  24. Grady, C. A., Muto, T., Hashimoto, J., et al. 2013, ApJ, 762, 48
  25. Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84
  26. Haisch, K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153
  27. Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  28. Hernández, J., Hartmann, L., Megeath, T., et al. 2007, ApJ, 662, 1067
  29. Horn, B., Lyra, W., Mac Low, M.-M., & Sándor, Z.  2012, ApJ, 750, 34
  30. Hubeny, I. 1990, ApJ, 351, 632
  31. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Natur, 448, 1022
  32. Johansen, A., & Lacerda, P. 2010, MNRAS, 404, 475
  33. Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, SciA, 1, 1500109
  34. Juhász, A., Benisty, M., Pohl, A., et al.  2015, MNRAS, 451, 1147
  35. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32
  36. Lee, W.-K.  2016, arXiv:1604.08941
  37. Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015a, Natur, 524, 322
  38. Lyra, W., Paardekooper, S.-J., & Mac Low, M.-M.  2010, ApJ, 715, L68
  39. Lyra, W., Richert, A. J. W., Boley, A., et al. 2016, ApJ, 817, 102
  40. Mamajek, E. E. 2009, in AIP Conf. Proc. 1158, Exoplanets and Disks: Their Formation and Diversity, ed. T. Usuda, M. Tamura, & M. Ishii (Melville, NY: AIP), 3
  41. Masset, F. 2000, A&AS, 141, 165
  42. Morbidelli, A., & Nesvorny, D. 2012, A&A, 546, A18
  43. Morbidelli, A., Lunine, J. I., O’Brian, D. P., Raymond, S. N., & Walsh, K. J.  2012, AREPS, 40, 251
  44. Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B.  2015, Icarus, 258, 418
  45. Movshovitz, N., Bodenheimer, P., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616
  46. Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJ, 748, L22
  47. Nelson, R. P., & Gressel, O. 2010, MNRAS, 409, 639
  48. Nelson, R. P., Gressel, O., & Umurhan, O. M.  2013, MNRAS, 435, 2610
  49. Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43
  50. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  51. Probstein, R. F., & Fassio, F.  1969, AIAA, 8, 4
  52. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  53. Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753
  54. Stolker, T., Dominik, C., Avenhaus, H., et al.  2016, arXiv:1603.00481
  55. Tanigawa, T., & Tanaka, H. 2016, ApJ, 823, 48
  56. Testi, L., Birnstiel, T., Ricci, L., et al.  2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. of Arizona Press), 339
  57. Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. M. 2011, Nature, 475, 206
  58. Whipple, F. L. 1972, From Plasma to Planet, ed. A. Evlius (NewYork:Wiley), 211
  59. Weidenschilling, S. J. 1977, MNRAS, 180, 57
  60. Yang, C.-C., Mac Low, M.-M, & Menou, K. 2009, ApJ, 707, 1233
  61. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
  62. Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
  63. Zhu, Z., Hartmann, L., & Gammie, C. 2009, ApJ, 694, 1045
  64. Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-N. 2014, ApJ, 785, 122
  65. Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description