Spiral Wave Instability in Disks

Self-Destructing Spiral Waves: Global Simulations of a Spiral Wave Instability in Accretion Disks

Jaehan Bae11affiliation: Department of Astronomy, University of Michigan, 1085 S. University Ave., Ann Arbor, MI 48109, USA , Richard P. Nelson22affiliation: Astronomy Unit, Queen Mary University of London, Mile End Road, London E1 4NS, UK , Lee Hartmann11affiliation: Department of Astronomy, University of Michigan, 1085 S. University Ave., Ann Arbor, MI 48109, USA , Samuel Richard22affiliation: Astronomy Unit, Queen Mary University of London, Mile End Road, London E1 4NS, UK jaehbae@umich.edu, r.p.nelson@qmul.ac.uk, lhartm@umich.edu, samuel.richard@qmul.ac.uk

We present results from a suite of three-dimensional global hydrodynamic simulations which show that spiral density waves propagating in circumstellar disks are unstable to the growth of a parametric instability that leads to break-down of the flow into turbulence. This spiral wave instability (SWI) arises from a resonant interaction between pairs of inertial waves, or inertial-gravity waves, and the background spiral wave. The development of the instability in the linear regime involves the growth of a broad spectrum of inertial modes, with growth rates on the order of the orbital time, and results in a nonlinear saturated state in which turbulent velocity perturbations are of a similar magnitude to those induced by the spiral wave. The turbulence induces angular momentum transport, and vertical mixing, at a rate that depends locally on the amplitude of the spiral wave (we obtain a stress parameter in our reference model). The instability is found to operate in a wide-range of disk models, including those with isothermal or adiabatic equations of state, and in viscous disks where the dimensionless kinematic viscosity . This robustness suggests that the instability will have applications to a broad range of astrophysical disk-related phenomena, including those in close binary systems, planets embedded in protoplanetary disks (including Jupiter in our own Solar System) and FU Orionis outburst models. Further work is required to determine the nature of the instability, and to evaluate its observational consequences, in physically more complete disk models than we have considered in this paper.

Subject headings:
accretion disks, hydrodynamics, instabilities, waves

1. Introduction

Spiral density waves are excited in circumstellar disks by a variety of processes, and play an important role in the dynamical evolution and observational appearance of these systems. They arise as normal modes in self-gravitating disks (Lin & Shu, 1964), where they provide an efficient mechanism for the transport of mass and angular momentum through gravitational torques and an advective wave flux (Lynden-Bell & Kalnajs, 1972; Papaloizou & Savonije, 1991; Laughlin & Bodenheimer, 1994). Spiral shocks can be excited by stellar companions in close binary systems, including cataclysmic variables, X-ray binaries and T Tauri binaries (Lin & Papaloizou, 1979; Sawada et al., 1986; Artymowicz & Lubow, 1994). The excitation of nonlinear spiral waves by giant planets also plays an important role in structuring protoplanetary disks, and in driving the migration of embedded planets (Bryden et al., 1999; Kley, 1999; Nelson et al., 2000). Models for the origin of FU Orionis outbursts have been presented where spiral waves, originating in the gravitationally unstable outer regions of protostellar disks during the early infall phase, propagate into the inner disk regions and trigger the magnetorotational instability by heating and ionizing the disk gas there (Gammie, 1999; Armitage et al., 2001; Zhu et al., 2010; Bae et al., 2014). Recent work has also shown that the infall of low angular momentum material onto a protostellar disk during the early infall phase can also generate global spiral waves (Lesur et al., 2015).

In this paper, we show that spiral waves in circumstellar disks are unstable to the growth of a parametric instability that leads to the disk flow becoming turbulent. We came across this phenomenon while extending the two-dimensional hydrodynamic simulations of triggered FU Orionis outbursts in Bae et al. (2014) to three dimensions. The instability arises because pairs of inertial waves, or inertial-gravity waves, couple resonantly to the spiral wave, leading to the extraction of energy and angular momentum from it and the growth of a broad spectrum of inertial waves. Similar parametric instabilities, leading to the growth of inertial waves and the generation of small scale turbulence, have been reported to arise in numerous circumstances where disk fluid elements are subjected to periodic forcing. Goodman (1993) and Ryu & Goodman (1994) showed that a parametric instability arises due to the elliptical distortion of a disk caused by an external orbiting companion. Warped disks have been shown to be subject to parametric instability (Gammie et al., 2000; Ogilvie & Latter, 2013), as have globally eccentric disks (Papaloizou, 2005a, b; Barker & Ogilvie, 2014). In a study similar to the one that we present here, Fromang & Papaloizou (2007) used shearing box simulations and analytical calculations to show that axisymmetric, nonlinear sound waves traveling in disks are also subject to parametric instability.

In this paper, we examine the stability of propagating spiral waves in circumstellar disks, by carrying out three-dimensional global hydrodynamic simulations at high resolution. We consider a broad range of simplified disk models and input physics, including both isothermal and adiabatic equations of state, viscous and inviscid disks, and vertically stratified and non-stratified density profiles. Simulations are computed with four independent numerical codes, and the outcomes are found to be in good agreement with each other. Our main result is that the spiral density waves are unstable to the growth of the aforementioned parametric instability, leading to the development of turbulence in the disk models. We refer to the instability as the spiral wave instability (SWI). The SWI is found to be robust across the full range of models that we considered, suggesting that it will be influential in the dynamical evolution of disks that contain nonlinear spiral density waves of whatever origin.

This paper is organized as follows. In Section 2, we discuss the theoretical background to the SWI to set the scene for the numerical simulations. We describe our numerical methods in Section 3, including the basic equations solved and the disk models examined. In Section 4, we present results from cylindrical disk models as a demonstration of the instability in the simplest global model with a non-stratified, globally isothermal initial setup. Results from vertically stratified disk models with isothermal and adiabatic equations of state are discussed in Sections 5 and 6, respectively. We discuss the application of the SWI to various astrophysical disk phenomena, with a particular emphasis on protoplanetary disks, in Section 7, and we provide concluding remarks in Section 8. A reader who is mainly interested in the physical manifestation of the instability, and its potential applications, can safely skip Sections 24 and read from Section 5 onwards.

2. Theoretical Background

Working in the context of a cylindrical coordinate system (, , ), we consider the situation where a spiral density wave, with azimuthal mode number , propagates inwards in a quasi-Keplerian disk. The disk is assumed to be isothermal in the vertical direction (i.e. temperature independent of ). We consider models in which the thermodynamic response of the gas is adiabatic, such that , and models where the disk gas responds isothermally such that . Under these conditions, the sound speed of the gas is independent of , and according to linear theory a spiral density wave that is launched from an inner Lindblad resonance (ILR) and has no vertical structure should propagate inwards maintaining its two-dimensional nature (Lin et al., 1990). Here, the absence of vertical structure simply means that perturbed quantities associated with the wave are independent of height.

A parametric instability may arise when an oscillator is subjected to a particular form of forcing at twice the natural frequency of the oscillator (e.g. a swinging pendulum with a moment of inertia that varies sinusoidally at twice the pendulum’s mean natural frequency). In the context of astrophysical disks, parametric instabilities have been examined that arise from the excitation of pairs of inertial modes in disks that are elliptical due to tidal distortion by an external binary companion (Goodman, 1993) or due to the propagation of axisymmetric, nonlinear sound waves (Fromang & Papaloizou, 2007), where the two excited inertial waves have frequencies that satisfy the resonance condition , where the frequency of the forcing disturbance is . Viewed in a frame corotating with a fluid element in the disk, the frequency associated with a spiral wave that propagates interior to its inner Lindblad resonance is given by , where is the local orbital angular frequency and is the pattern speed of the spiral wave measured in the inertial frame. In the specific case of an external companion on a circular orbit, would be the angular frequency of the companion’s orbit.

The WKBJ dispersion relation for local disturbances in a differentially rotating disk is given by (Goodman, 1993)


where is the mode frequency, is the epicyclic frequency defined as


is the Brunt-Väisälä frequency defined as


and are the wave numbers associated with the wave vector and is the sound speed. In the high frequency limit, , , Equation (1) supports acoustic waves where . In the low frequency limit, where , the dispersion relation becomes


where the angle between the wave vector and the -axis . Equation (4) governs the behavior of inertial-gravity modes. These correspond to fluid elements undergoing epicyclic motions, with the fluid motions being confined to planes that are tilted by an angle with respect to the disk midplane, and which are perpendicular to the wave vector, . Coriolis and buoyancy forces provide the restoring forces in the horizontal and vertical directions, respectively. In the absence of vertical buoyancy, as in a disk with an isothermal response to perturbations ( when ), the dispersion relation simplifies to one that describes inertial waves


To simplify the discussion, we use the term inertial modes/waves to denote inertial-gravity waves (supported by Coriolis and buoyancy forces) and inertial waves (supported only by the Coriolis force) from now on.

We might expect parametric instability to arise locally, through the coupling of a pair of inertial waves to an incoming spiral wave with doppler-shifted frequency , when the two inertial waves have frequencies that obey the relation . Working in the high wave number limit where (Fromang & Papaloizou, 2007), we have


For the specific case of an spiral wave, as considered in this paper, this becomes


We see that the frequencies of inertial waves are independent of wave number, and simply depend on the orientation of the wave vector. For a disk that is finite in both radius and height, the wave numbers, and , and associated wave frequencies, will be discrete. Considering a situation where the radial and vertical wavelengths are given by and , respectively (where is the scale height), such that , and (, ) are integers, then we have


It is obvious from Equation (8) that if we consider arbitrarily large and , then can be matched with arbitrary precision to the value of , where . In other words, inertial waves have a dense spectrum, and the resonance condition can always be matched to high precision for large enough . Fromang & Papaloizou (2007), considering a disk with an isothermal equation of state, comment on the fact that the density of inertial modes is greater in a disk with vertical density stratification than in a non-stratified disk because of the larger range of length scales. Here, we also note that an adiabatic disk with also increases the density of inertial modes, because scans a range of frequencies as a function of height, and so increases the probability that the resonant condition for the parametric instability can be met at any given radius (i.e. one combination of and may satisfy the resonance condition near the midplane where , and another combination may satisfy the resonance condition away from the midplane where ).

2.1. Wavelengths of parametrically excited inertial modes

Obtaining a full and detailed understanding of which specific pairs of inertial modes will couple to an incoming spiral wave by satisfying the resonance condition would require determination of the radial and vertical eigenmode structure, and associated eigenfrequencies, for the global disk models that we consider. This goes beyond the scope of this paper (but see Barker & Latter (2015) for a study of a similar problem in the context of the vertical shear instability). Instead, we use simple arguments, based on a local WKBJ picture, to consider the approximate wave numbers associated with unstable inertial modes to help interpret and understand the simulations presented in later sections.

The WKBJ dispersion relation for a spiral wave with azimuthal mode number (in the absence of self-gravity) is given by (assuming )


giving a radial wave number interior to the ILR


Inertial modes that are excited by the spiral wave must occur on radial length scales similar to or smaller than the incoming wavelength, so for simplicity we write the wave number of the excited inertial waves (assuming they have very similar spatial structure and frequencies, appropriate to the high wave number limit) , where is an integer. 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 from Equation (1), written using the modified notation




Substituting the resonance condition and into Equation (12) gives


First, we consider the behavior of for disks in which


For to be real, such that inertial waves can be excited, we require that and , and these conditions are both satisfied at all radii interior to the inner Lindblad resonance of an spiral wave in a Keplerian disk. (These conditions, however, are not satisfied everywhere for a spiral wave that propagates outwards from an outer Lindblad resonance, as we discuss later in Section 7.) Considering the influence of buoyancy in Equation (13), we see that can become imaginary when , indicating that inertial modes will not be excited at certain locations above the midplane where the buoyancy frequency is large. Taking the large limit, the transition occurs at . Given that the resonance condition for the excited inertial waves is , this indicates that inertial modes cannot be excited where . Later in the paper, we use Equation (13) to construct maps of stable and unstable regions of our disks for comparison with the simulation results, and to consider the influence of numerical resolution on our ability to resolve unstable inertial waves that arise on small scales.

As we have emphasized, there are a number of caveats that should be considered before applying the above arguments to the simulation results. The most obvious is the fact that a local WKBJ analysis, that applies to modes where and , is not strictly applicable to modes with , as may be excited by spiral waves in a global disk model. Furthermore, for simplicity we have assumed that there is an integer relation between the wavelengths of the excited inertial waves and the incoming spiral modes, and this is not precisely the situation that should arise because the linear inertial modes supported in the global disk model should be quantized to fit within the radial boundaries of the disk, and hence are not expected to have an integer relation with the spiral waves. A full analysis of the modes that can be excited should take account of the full structure of the eigenfunctions in both radial and vertical directions, and their frequencies, which represents a significant computational problem (e.g. Barker & Latter, 2015), and goes beyond the scope of this paper. It is noteworthy, however, that Lubow & Pringle (1993) compute the vertical structure of inertial-gravity waves in vertically isothermal disks that support buoyancy forces, and demonstrate that inward traveling waves have their wave energy increasingly confined towards the midplane as they enter disk radii where large fractions of the disk vertical domain have , since this defines the region of space where the waves can (or cannot) propagate. Similar forbidden zones arise in the above WKBJ analysis for essentially the same reason, as we have described above.

2.2. Nonlinear wave propagation in vertically stratified disks

In linear theory, a two-dimensional spiral wave without any vertical structure, launched at an inner Lindblad resonance in a vertically isothermal disk with sound speed independent of , should propagate towards the star with wave fronts remaining perpendicular to the midplane. As described in the next section, the underlying axisymmetric disk model will normally have midplane density varying according to , with , and variation with height scaling as , where is the density scale height. For such a model, the wave front associated with an inward propagating wave will see an increasing midplane density as it moves in, but at high altitudes the wave front will see a sharply decreasing density. We might then expect that nonlinear effects experienced by the wave to be greater at high altitudes than near the midplane. These effects include distortion of the wave form and an enhanced propagation speed due to advection of the wave by the perturbed radial velocity. We expect this will result in curvature of the spiral wave fronts as these waves propagate towards the star. An important consequence of this is that vertical hydrostatic equilibrium will not be maintained at the wave fronts, resulting in the generation of vertical motions as the wave moves inwards, and an oscillating non-axisymmetric corrugation of the disk surface.

3. Numerical Methods

3.1. Basic Equations

We solve the hydrodynamic continuity, momentum, and internal energy equations:


where is the mass density, is the velocity, is the pressure, is the gravitational potential due to the central object, is the spiral potential (see Section 3.3), is the viscous stress tensor, and is the internal energy per unit volume.

We make use of both cylindrical and spherical coordinates. We consider cylindrical disk models where we neglect the vertical component of gravity from the central object, so the gravitational potential is , where is the gravitational constant, is the mass of the central object and is the cylindrical radius. These models are computed using cylindrical coordinates. The purpose of using non-stratified, cylindrical models is to illustrate the nature of the instability as clearly as possible, by removing/minimizing complications that may arise from spherical geometry, concave wave structure (see Section 5.1), etc. Furthermore, adopting a cylindrical disk model allows us to run models with periodic boundary conditions in the vertical direction, and hence allows us to demonstrate that the existence of the instability is independent of the chosen boundary conditions.

We also consider models where the density is vertically stratified. The gravitational potential is then , where is the spherical radius. These models are computed using spherical coordinates.

We make use of two forms of equation of state. In the isothermal calculations, we use an isothermal equation of state , where is the isothermal sound speed, and do not solve the internal energy equation. In the adiabatic calculations, we relate the gas pressure and the internal energy through , where is the adiabatic index, and solve the internal energy equation. We do not include the effects of cooling in the internal energy equation. Instead, we explore the effects of cooling by changing the adiabatic index value. The idea is that a disk with very rapid cooling acts as if it is isothermal, and hence has an effective , whereas a disk with very slow cooling acts as if it is adiabatic, and thus its effective is equal to the actual ratio of specific heats in the disk, here taken to be . Intermediate values of the effective correspond to intermediate cooling rates. We take this approach because determining the Brunt-Väisälä frequency is not trivial in the presence of cooling, making comparison between numerical results and theoretical predictions difficult.

We include physical and artificial viscosity (Stone & Norman, 1992). Lyra et al. (2016) has recently considered the local generation of entropy through strong spiral shocks induced by massive planets, that could in turn drive local convective motions. We have confirmed, however, that the shock dissipation heating is negligible in our adiabatic calculations because the spiral waves considered in this work are much weaker.

Run Label Code Numerical Resolution Kinematic
( or or ) Viscosity
R512 FARGO3D 0 1.0
R128 FARGO3D 0 1.0
R256 FARGO3D 0 1.0
R768 FARGO3D 0 1.0
AMP0.625 FARGO3D 0 1.0
AMP1.25 FARGO3D 0 1.0
AMP2.5 FARGO3D 0 1.0
AMP10 FARGO3D 0 1.0
V6 FARGO3D 1.0
V5 FARGO3D 1.0
V4 FARGO3D 1.0
GAM01 FARGO3D 0 1.01
GAM05 FARGO3D 0 1.05
GAM1 FARGO3D 0 1.1
GAM2 FARGO3D 0 1.2
GAM3 FARGO3D 0 1.3
GAM4 FARGO3D 0 1.4
Table 1Model Parameters

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 . In all models presented, is chosen such that the ratio of disk scale height to radius at is . The isothermal sound speed is related to the temperature by , where is the gas constant and is the mean molecular weight. Thus, Equation (18) corresponds to the radial sound speed distribution given by


For the main simulation set we assume , so that the disk has the same temperature everywhere. This temperature structure might not be very realistic, especially when the simulation domain extends more than an order of magnitude in radius. On the other hand, one can safely ignore the vertical shear instability in these models, because the instability requires a non-zero vertical gradient in the disk angular velocity which arises when there is a radial temperature gradient (Nelson et al., 2013). In the Appendix, we present models with a radial temperature gradient (). In those models, we find that the spiral wave instability and the vertical shear instability coexist in the absence of kinematic viscosity, and thus identifying intrinsic features arising from the SWI is challenging. The vertical shear instability can be suppressed with a non-zero kinematic viscosity, but adding kinematic viscosity not only suppresses the vertical shear instability but also damps the growth of small scale unstable modes excited by the SWI. To summarize, by assuming the globally isothermal temperature structure, we aim to illustrate the nature of the SWI as clearly as possible, in the absence of any other hydrodynamic instabilities.

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




The midplane density at , , is chosen such that the total disk mass is of the central object’s mass. When vertical stratification of the disk is ignored in the cylindrical disk models, the above equations simplify to:




We fix the radial power-law index of the midplane gas density to for all models. The initial radial and vertical/meridional velocity are set to zero, but uniformly distributed random perturbations are added as white noise, at the level of , to the initial vertical/meridional velocity in order to seed the instability. We have tested the effects of changing the amplitudes of the random perturbations, and found that the growth rate and the saturation level of the instability are not sensitive to the initial noise level.

Our simulation domain extends from to 10 in radius and from 0 to in azimuth. In cylindrical models, the vertical domain extends from to 0.1. In spherical models, the meridional domain covers above and below the midplane. The model parameters are summarized in Table 1.

3.3. Spiral Potential

Figure 1.— (Left) The radial velocity in the disk midplane induced by the spiral potential with amplitudes of (solid) from model AMP0.625 and (dashed) from model R512. The arrows indicate the location of the inner and outer Lindblad resonances (ILR and OLR). (Right) The azimuthal distribution of the midplane density at , normalized by the azimuthal average. As in the left panel, the solid and dashed curves are when and , respectively. The data are taken after the spiral waves are well established over the entire disk, but before the spiral wave instability perturbs the wave structures.

In order to excite spiral waves, we adopt the following potential :


where is the amplitude, assumed to be constant over time, is the azimuthal mode number, is the pattern speed, is the radius about which the potential is centered, and is the radial width of the potential. We adopt the values and , and assume that the pattern speed is the local Keplerian frequency at its central position: . In this work, we only focus on spiral potentials with .

Our fiducial model assumes potential amplitude of . With this amplitude we intend to produce velocity perturbations at that are about of the local sound speed111As we mentioned in the introduction, the spiral wave instability was found while studying spiral wave propagation as part of a study of accretion outbursts in FU Orionis systems. This spiral amplitude produces velocity and density fluctuations in the inner disk similar to the ones produced by the gravitationally unstable outer disk, with which an accretion outburst was triggered (see Figure 5b of Bae et al. 2014).. The effect of varying the spiral potential strength on the spiral wave instability is investigated in Section 5.3.

Prior to discussing the results of simulations in any detail, we briefly illustrate the effects of changing the spiral wave amplitude on the wave forms and propagation speeds of the waves in the midplane of a stratified disk. In Figure 1, we present the radial velocity profiles in the disk midplane of a stratified disk model, driven by two spiral potential amplitudes of and . As shown, the overall structure of the spiral waves driven in the numerical simulations agrees well with theoretical expectations: from the WKBJ dispersion relation for a non self-gravitating disk one expects to see waves propagating interior and exterior to the inner and outer Lindblad resonances, with the waves being evanescent in the region between the resonances (Lin & Shu, 1964; Goldreich & Tremaine, 1978, 1979). The wave with the small potential amplitude is nearly sinusoidal. In the inner disk, the velocity amplitude is about a few per cent of the sound speed and is nearly constant over radius. When the spiral potential amplitude is increased to , on the other hand, the waves are rather strong, generating sharp wave fronts that have radial velocities up to about of the sound speed. At , the two spiral potential amplitudes generate relative density enhancements of () and () above the azimuthal average. In addition to these differences in perturbed velocities and densities, we also note that the propagation speed of the higher amplitude wave is slightly faster than that of the lower amplitude wave. This is a nonlinear effect arising from the fact that the magnitude of the perturbed velocity contributes to the wave speed through the advection term in the momentum equation. The magnitude of this effect as a function of height in the disk plays an important role in determining the shape of the spiral wave fronts in the stratified disk models, as discussed in Section 2 and demonstrated in Section 5.

3.4. Boundary Conditions

One needs to take particular care with the boundary conditions when simulating an instability, in order to ensure that what is being observed is not an artifact generated by the boundaries. The requirement that models with vertical stratification are able to maintain hydrostatic equilibrium has forced us to pay particular attention to the meridional boundary conditions. The radial boundary condition is chosen to have a zero gradient for all variables in all models, and a wave damping zone is implemented (de Val-Borro et al., 2006) in the intervals and . We choose one local orbital time for the wave damping timescale. As the simulation domain covers in azimuth, periodic boundary conditions are used in azimuth.

Vertical periodic boundary conditions are used in models where we neglect vertical stratification. In the stratified models, we use an outflow meridional 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 directing toward the disk midplane. In adiabatic calculations, 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:


In order to investigate the influence of the meridional boundaries in the stratified models, we tested various boundary conditions including the zero-gradient boundary condition, the standard outflow boundary condition, the reflecting boundary condition, and the wave-damping boundary conditions. We find that the choice of the meridional boundary condition does not affect the triggering of the parametric instability. For the Godunov codes, however, we find that enforcing hydrostatic equilibrium in the meridional ghost zones using Equation (25) maintains the low-density surface region more stably than with other boundary conditions, allowing us to better observe the early development of the instability.

3.5. Codes

We use four independent grid-based codes: FARGO3D, NIRVANA, PLUTO, and INABA3D. These include two finite difference codes and two Godunov codes, to ensure the robustness of the results.

FARGO3D is a finite difference code developed with special emphasis on disk simulations (Benítez-Llambay & Masset, 2016). An orbital advection algorithm is implemented as in its two-dimensional predecessor FARGO (Fast Advection in Rotating Gaseous Objects; Masset 2000) code. We have tested with and without the FARGO algorithm and found that the result is not dependent on the use of the FARGO algorithm.

Similar to FARGO3D, NIRVANA uses an algorithm very similar to that used in the ZEUS code to solve the equations of ideal MHD (Ziegler & Yorke, 1997; Stone & Norman, 1992). This scheme uses operator splitting, dividing the governing equations into source and transport terms. Advection is performed using the second-order monotonic transport scheme (van Leer, 1977).

PLUTO is a general-purpose Godunov code (Mignone et al., 2007). We employ the piece-wise parabolic reconstruction and third-order Runge-Kutta time integration. We note that lower-order schemes, e.g. the piece-wise linear reconstruction and second-order Runge-Kutta time integration, do not produce a converged Riemann solution for the problem considered in this work. The FARGO orbital advection module was enabled for the calculations with PLUTO.

INABA3D is a finite volume code using an unsplit MUSCL Hancock scheme, which is a second order method in time and space. The physical values on the faces of each cells are calculated through an exact Riemann solver. The code is based on the 2D code described in Inaba et al. (2005), and the third dimension has been implemented as in Richard et al. (2013).

As shown in Table 1, most simulations in this paper were computed using the FARGO3D code. We note that FARGO3D and PLUTO use a logarithmically spaced radial grid, whereas NIRVANA and INABA3D use a uniform mesh. The innermost grid cell size in the FARGO3D and PLUTO runs was , and this value was used for the uniform meshes in the NIRVANA and INABA3D runs.

3.6. Code Units

We choose code units such that . Since we do not include cooling, all the calculations are scalable to any system of interest. In the following sections we will use the orbital time at as the time unit. We will denote this quantity as . One orbital time at the spiral potential center , , corresponds to .

Figure 2.— Time evolution of the integrated vertical kinetic energy for non-stratified models with different codes (CYL-F, CYL-N, CYL-P, CYL-I). The black dashed line corresponds to the growth rate of .

3.7. Diagnostics

In order to examine the growth and evolution of the spiral wave instability we calculate the volume-integrated vertical kinetic energy, , in non-stratified disk models (as these are computed using cylindrical coordinates) and the volume-integrated meridional kinetic energy, , in the stratified disk models (as these use spherical polar coordinates):


The volume integration is performed in between and , in order to capture growth in a local region of the disk to better measure exponential growth rates.

Figure 3.— Distributions of the vertical velocity in a plane () interior to at various times for model CYL-F. The upper panels show the linear growth phase of the SWI. The instability saturates when the perturbed vertical velocity of the unstable modes is comparable to the radial velocity of the imposed spiral waves.
Figure 4.— Same as Figure 3, but for the radial velocity.

4. Non-Stratified Disk Models

We start discussion of the simulations by describing results of the non-stratified disk models. In Figure 2, we present the time evolution of the integrated vertical kinetic energy , where we see that this quantity exponentially increases at early times () before saturating to a quasi-steady value at late times, indicating the presence of an instability.

Figure 5.— The predicted (left) vertical and (right) radial wavelengths of the unstable modes and as a function of radius, calculated from the relation given in Equation (13) and Equation (10) with , where is an integer, and and are the wave numbers of the spiral wave and the excited inertial modes, respectively. The over-plotted symbols indicate the measured vertical and radial wavelengths of the unstable inertial modes at the innermost spiral arm location at (; diamonds) and 30 (; triangles). In the right panel, the solid curves present the radial wavelengths and the dashed curves show the number of radial grid cells in one inertial mode wavelength ().

In Figures 3 and 4, we display the two-dimensional distributions of the vertical and radial velocities at some selected times. Since vertical density stratification of the disk is neglected, and the sound speed is independent of height, the spiral waves propagate purely radially before the SWI sets in, maintaining wave fronts that are perpendicular to the disk midplane, as seen in the first two panels of Figure 4. Hence, the spiral waves in themselves are not expected to generate any vertical motion in the disk. We see that the inertial modes grow fastest in the inner disk regions, as their frequencies and growth rates are larger there, and the instability spreads to larger radii as time progresses. At , the inertial modes that are excited by the spiral wave form a checkerboard pattern in the vertical velocities, as seen in Figure 3. Close visual inspection of the plot, combined with a Fourier analysis, shows that in the disk inner regions around , the vertical wavelength associated with the excited inertial modes and the radial wavelength . These wavelengths correspond to approximately 3 and 12 grid cell spacings in the vertical and radial directions, respectively, indicating that the fastest growing modes have the shortest wavelengths that can be represented on the computational grid, presumably limited by the vertical resolution in this case. This is in agreement with Fromang & Papaloizou (2007), who showed that short wavelength inertial waves excited by axisymmetric, nonlinear sound waves have higher growth rates than longer wavelength modes.

Figure 3 shows that at slightly later times, , the length scales associated with the most prominent perturbations have increased substantially, presumably indicating that these longer wavelength modes grow on longer time scales, but may also contain more energy than the shorter wavelength modes seen to grow earlier. This latter point may, however, be affected by numerical diffusion damping the smallest scale modes, so some caution is required when interpreting this result. Inspection of the figure, and Fourier analysis, indicate that these modes have wavelengths and . We have plotted these and the earlier wavelength values in Figure 5, which shows the vertical and radial wavelengths of the inertial modes that are predicted to be excited by the spiral wave using the simple local analysis presented in Section 2. Noting the caveats associated with the analysis discussed earlier in the paper, we see that the early growing modes have radial wavelengths that are that of the incoming spiral wave, and the later growing modes have radial wavelengths very similar to the spiral wave. The numbers of grid cells covered by the different wavelengths are also indicated in this figure, demonstrating that the earlier growing modes occur on length scales close to the grid scale, whereas the later growing modes are more comfortably represented on the grid.

Figure 6.— Distribution of the normalized perturbed density in the plane at from the CYL-F model.
Figure 7.— Distributions of the vertical velocity in a vertical plane obtained with (left) NIRVANA, (middle) PLUTO, and (right) INABA3D. The upper panels present the distribution at to compare the early evolution of the instability with CYL-F model. For the CYL-I model, we present the data taken at since the growth of the instability starts at a later time (see Figure 2). The lower panels present the distribution at the end of calculations ().

The growth rate measured during the exponential growth phase shown in Figure 2 is measured to be . While this represents the growth associated with the superposition of numerous growing inertial modes between the radii , it is likely to be dominated by the longer wavelength modes. The instability saturates at about . As shown in Figure 3 and 4, the maximum perturbed vertical and radial velocities of the unstable inertial modes at saturation are on the order of the radial velocity of the background wave. This is presumably an indication of the fact that the linear instability saturates when the perturbed velocity grows large enough to disrupt the structure of the background wave, as also argued by Fromang & Papaloizou (2007). During the non-linear phase of the instability, the checkerboard shaped cells merge together, creating alternating vertical flows. At , the vertical flows are locally as fast as of the sound speed.

In Figure 6, we present the perturbed density distribution in a horizontal (, ) plane at , where the angled brackets denote an azimuthal average and . As expected from the snapshots in the vertical plane, the spiral arms are highly perturbed, showing fragmentary structures. Moving along the azimuth at any given radius shows a clear asymmetry. We will further discuss the appearance of the disrupted spiral structures in Section 7.2.3, along with their observational implications.

4.1. Code Comparison

In order to demonstrate the robustness of the instability we have run the same calculation with the four independent codes introduced in Section 3.5. As shown in Figure 2, clear exponential growth of the instability is observed with all codes. The growth time measured during the linear phase ranges between . The saturated perturbed energy levels also show good agreement with each other, being within a factor at .

We note, however, that the early evolution of the perturbed kinetic energy () seems to be more code dependent than the later evolution. In particular, the exponential growth seen at in the CYL-F model is not observed with the other codes. As discussed earlier, this is when the small scale modes, with wavelengths comparable to grid cell size, grow in the models. The growth rates of such small wavelength modes are inevitably dependent on the dissipative properties of each code.

In Figure 7, we present contour plots of the vertical velocity obtained with NIRVANA, PLUTO and INABA3D. Although the detailed evolution of the instability may slightly differ, it is clear that the same instability is being captured by all codes – the instability exhibits a checkerboard pattern during the linear phase, and alternating vertical flows when saturated.

As mentioned in Section 3.2, we initially add small, random, vertical velocities to seed the instability. Adding noise puts energy into the inertial modes in the disk. When the initial random perturbation is not added, we note that the cylindrical model does not trigger the SWI with the finite volume codes PLUTO and INABA3D. These codes do an excellent job of maintaining the symmetry associated with the initial conditions because the conservative properties of these codes are enforced by the pairwise exchange of fluid properties between neighboring grid cells. In the finite difference codes, however, we find that the SWI still develops when starting without the initial perturbations, presumably because machine precision level errors are introduced by the use of finite differences.

Figure 8.— Time evolution of the integrated meridional kinetic energy . Results with four different resolutions are plotted. The plot indicates that the numerical results more or less converge at and beyond. We run the two low resolution models for a longer time to see if the saturation level eventually converges toward the value with higher resolutions, but find that the perturbed energy level does not further increase after in both and run.
Figure 9.— Contour plots presenting two-dimensional distributions of the meridional kinetic energy density in a vertical plane (). The snapshots are taken (upper-left) before the spiral wave instability is triggered, (upper-right) at the beginning of the linear growth phase, (lower-left) at the middle of the linear growth phase, and (lower-right) at the end of the calculation when the instability is fully saturated. The white dashed contours in the upper-left panel show the surface. Note that the wave fronts are curved towards the midplane due to nonlinear advection of the wave. The arrows in the upper-right panel indicate the four locations where the spiral arms intersect the plane at that time.

5. Vertically Stratified, Isothermal Disk Models

In this section, we present results of vertically stratified models with an isothermal equation of state. The calculations are initiated with the vertically stratified density structure and the corresponding azimuthal velocity that satisfy the hydrostatic equilibrium, as introduced in Equations (20) and (21). We first introduce our reference model, for which the spiral wave amplitude and which has zero kinematic viscosity. The numerical resolution is . The effects of varying the numerical resolution, spiral potential amplitude, and viscosity will be discussed in the following sections.

Figure 10.— Contour plots of the meridional kinetic energy density in a vertical plane (), obtained with grid cells.

As discussed in Section 2, we expect that the spectrum of inertial modes that can be excited by the spiral waves will be dense, and therefore it is difficult to identify individual growing modes during the simulations. The best way of doing this would be to compute the eigenfunctions associated with these modes, such that they could be seeded in the initial conditions and their growth rates measured. We take a different approach by using white noise to perturb the initial conditions, thereby adding energy to the full spectrum of inertial modes, and ensuring the robustness of the SWI by examining the triggering of the instability under a variety of disk conditions, instead of identifying and characterizing individual unstable modes.

5.1. A Reference Run (R512)

In Figure 8, we present the time evolution of the integrated meridional kinetic energy . Contour plots of the meridional kinetic energy density are plotted in Figure 9. We make use of the meridional kinetic energy density for the stratified models, instead of the meridional velocity, in order to better illustrate the development of the instability in the main body of the disk.

As shown in Figure 8, increases rapidly very early in the simulations at . This initial increase does not represent the growth of unstable inertial modes. As shown in the first panel of Figure 9, the spiral waves are not perpendicular to the disk midplane in the stratified models, but instead have a concave shape that develops as they propagate towards the star. As discussed in Section 2, this seems to arise as a nonlinear effect due faster advection of the wave at higher altitudes in the disk because the initial planar wave fronts move through steeply decreasing density profiles there, whereas the density increases near the midplane as the wave moves inwards. The curvature results in a loss of vertical hydrostatic equilibrium in the vertical direction, giving rise to the meridional kinetic energy.

The exponential growth associated with the SWI starts at , and saturates at . The second and third panel of Figure 9 show the linear growth and saturation of unstable modes at the innermost spiral arm. The unstable modes grow faster at smaller radii as shown, because of the higher mode frequencies and growth rates closer to the central object. When the instability is saturated (), the disk ends up in a complicated quasi-steady turbulent state. We discuss the turbulence produced via the SWI and its implications for angular momentum transport and vertical mixing later in Section 7.2.1.

5.2. Effect of Numerical Resolution

In order to test the effect of numerical resolution, we run the reference model with four different resolutions: , and . Increasing the numerical resolution should allow the growth of smaller scale inertial modes.

As seen in Figure 8, the growth rate and the saturated perturbed energy are reasonably well converged at and beyond. As expected, the highest resolution run displays the highest level of perturbed energy in the saturated state on average, because of the excitation of smaller scale modes, but the difference between the and runs is modest, indicating that the modes containing the largest amount of energy are captured in the lower resolution calculation. We note that there is a significant difference between the two higher resolution cases and the and cases, with the lowest resolution run in particular showing a much lower growth rate and saturation level.

In Figure 10, we present contour plots of the meridional kinetic energy density for R768 model at and 200. Compared to the results obtained with grid cells presented in Figure 9, one can see more fine scale structures developing at the beginning of the instability (). It is also apparent that, with the higher resolution, the unstable modes have grown faster at the third and fourth innermost spiral arms. In the R768 model, we find that the spiral waves create a larger density enhancement at the wave fronts than in the reference model: at , the spiral waves induce of density enhancement in the midplane in R768 model, whereas the enhancement is in the reference model. Since the growth rate of unstable modes is an increasing function of the wave amplitude (Fromang & Papaloizou, 2007, and see below), the stronger perturbations lead to faster growth of the unstable modes, in addition to there being more modes excited at high resolution.

Figure 11.— Time evolution of the meridional kinetic energy for different potential amplitude . The percentage in the parentheses indicate the density enhancement induced by the imposed spiral waves at (refer to the right panel of Figure 1).
Figure 12.— Time evolution of the meridional kinetic energy for different values of the kinematic viscosity. Note that the spiral wave instability operates with , while the instability is suppressed when .

5.3. Effect of Spiral Potential Amplitude

In this section, we explore the effect of varying the imposed spiral potential strength. We increase the spiral potential amplitude by factors of two from the smallest value, , up to the largest value .

Figure 11 shows the time evolution of . It is immediately obvious from the figure that the instability has a higher perturbed energy level with stronger spiral waves. During the initial phase when the spiral waves propagate into the inner disk (), stronger spiral waves create more vertical motion. In addition, the linear growth phase begins at an earlier time with a stronger spiral amplitude. In the case with , for example, the instability starts at about , while in the reference model the instability starts almost immediately after the initial propagation of the waves.

From to , the density enhancement doubles when the spiral potential amplitude is doubled. However, the density enhancement induced by the spiral waves does not linearly increase with the potential strength when . Comparing model with model, the density enhancement increases only by about , presumably as a result of nonlinear dissipation. The nonlinear dissipation is observed over a broad disk region (), and interestingly, it is more significant at smaller radius.

Figure 13.— Contour plots showing two-dimensional distributions of the meridional kinetic energy density in the vertical plane with different kinematic viscosities. The upper panels show the distributions at the beginning of the linear growth phase, and the lower panels show the distributions at . Note that the instability operates with (left) and (middle) , but is suppressed with (right) .

5.4. Effect of Viscosity

We add a constant kinematic viscosity to our reference model to study the triggering and the evolution of the SWI in the presence of viscous dissipation. We test with three different kinematic viscosities: , and . In the canonical prescription of Shakura & Sunyaev (1973), the kinematic viscosity corresponds to


If we assume that radii are measured in units of AU, then a value of corresponds to at 1 AU, and at 5 AU. In Figure 12, we plot the time evolution of the integrated meridional kinetic energy. As the plot indicates, the SWI operates with , but is completely suppressed when . The linear growth of the instability begins at a later time with a larger viscosity, probably because the viscosity damps the smaller scale fastest growing modes. Also, the saturated energy level is reduced with a larger viscosity, for the same reason, although the difference between and models is only marginal, suggesting that numerical diffusion operates with a value that is moderately below in the inviscid model.

Figure 14.— Time evolution of the meridional kinetic energy for different adiabatic indices .
Figure 15.— Contour plots of the meridional kinetic energy density for adiabatic models with different values. The black curves indicate where the local buoyancy frequency equals to a half of the local Doppler-shifted frequency of the spiral waves: . The results from the run demonstrate the fact that the region where the SWI develops is confined to the regions where .

Contour plots of meridional kinetic energy density are shown in Figure 13. As inferred from the energy evolution, one can see that the SWI operates with . With , we do not find any signatures of the instability growing throughout the duration of the calculation. Comparing with the results of the reference model at presented in Figure 9, we find that the small scale unstable modes are absent with . With , small scale modes are further suppressed, and the unstable modes that have grown in the model have noticeably larger length scales than in the and models. Fromang & Papaloizou (2007) show that, except for the largest wavelength inertial modes, the growth rates are almost constant. Given that we expect modes to be viscously damped when the damping rate exceeds the growth rate, and the viscous damping rate scales as , this demonstrates that we expect the high (small ) modes to be preferentially damped.

We emphasize that the instability operates with fairly large viscosity of . In terms of the canonical description, corresponds to at (right outside of the inner damping zone) and at . The fact that the spiral wave instability operates in quite dissipative disks may be important for disks that sustain fully developed MHD turbulence, such as the accretion disks around cataclysmic variables and/or black holes. Furthermore, the growth rate of the instability depends on the amplitude of the spiral wave, so the value of for which the SWI operates will depend on the incoming wave amplitude.

6. Vertically Stratified, Adiabatic Disk Models

In this section, we adopt an adiabatic equation of state to investigate the influence of the thermodynamics on the growth and development of the SWI. We consider six different adiabatic index values: . By using different values of , we mimic a disk that experiences cooling at different rates. Reducing acts as if the disk has more efficient cooling. Using the broad range of values, we aim to demonstrate that the SWI operates in both optically thin regions (isothermal) and optically thick regions (adiabatic). As we shall show below, the disks become more susceptible to the SWI with stronger adiabatic responses (i.e. larger value). In contrast to the reference model that we have presented earlier in this paper, we use a smaller spiral potential amplitude of to allow us to capture the growth of the instability at early times.

In Figure 14, we present the time evolution of the integrated meridional kinetic energy. When the spiral waves initially propagate into the inner disk (), a more rapid and higher amplitude growth of is observed with a larger value, and it is clear that the SWI triggers at earlier times with larger values. For example, the linear growth of the instability begins at with , but at with and at with . For , it is difficult to disentangle the initial growth of the perturbed kinetic energy due to the spiral wave propagation from the growth of the instability.

As we discussed in Section 2, the Brunt-Väisälä frequency varies over a wide range of values over height in adiabatic disks, with the range of values covered being a function of the value of the effective adopted. This suggests that at each radius in an adiabatic disk, the spiral modes will have a larger number of inertial modes with which to resonantly interact, with the number being larger for larger values of . We might therefore expect to see a more rapid growth of the SWI for the larger values of , and this expectation is borne out by the simulation results.

We also pointed out in Section 2 that there can be forbidden regions for the SWI in adiabatic disks. The existence of these can be inferred from Equation (13), which shows that the vertical wave number has no physical solution in disk regions where in the large limit. Although in global models, such as those we present here, the analysis of the properties of the inertial waves should also be global and not local, the analysis of the vertical structure of inertial mode eigenfunctions by Lubow & Pringle (1993) also shows that the energy in the modes is confined near the midplane, and the modes do not propagate in regions where the inertial mode frequencies .

In Figure 15, we display contour plots of the meridional kinetic energy density at the end of calculations () with , 1.2, and 1.4. In the figure, the black curves connect the disk regions where the local Brunt-Väisälä frequency matches to a half of the Doppler-shifted spiral wave frequency: . As discussed above, the region surrounded by the black curves is where the inertial waves are expected to be confined. The forbidden region is not clearly seen in the and 1.2 models (although there is a strong hint of it in the latter model), because the region is located near the meridional boundary, and acoustic waves excited by the SWI are likely to be present there along with vertical motions induced by the spiral waves. We note, however, the strong signature of inertial mode confinement near the midplane in the model.

Finally, we note that varying the adiabatic index affects the shape and strength of spiral waves. As is increased, spiral waves are more openly wound, due to the larger sound speed, and the density enhancement at wave fronts is reduced. These, in turn, can have influence on the growth rate of the SWI and the saturated energy level, although the existence of the SWI in adiabatic disks is clearly robust. At , for example, the density perturbation before the SWI sets in is with , with , and only with . The reduction in density enhancement at spiral wave fronts for larger values of may explain why the saturated state of the model with is somewhat smaller than in the other models (as seen in Figure 14), in spite of the obvious faster growth rate associated with this run.

Figure 16.— (a) Contour plots of the vertical wavelength of the inertial modes that are subject to the resonant interaction with the imposed spiral waves assumed in this work. (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 dashed curves indicate where . The black regions are where the dispersion relation does not have a physical solution. In (a) and (b), the upper half of the disk shows the wavelength in units of disk scale height at each radius, whereas the lower half of the disk shows the wavelength in units of vertical and radial grid cell size at each position. We use the mode as a representative, but the numbers on the colorbars can be linearly scaled to other modes as noted in the labels.

7. Discussion

Our results suggest that the SWI can operate under a broad range of disk conditions. In addition, the instability presumably operates regardless of the origin of spiral waves, as long as the waves maintain a fixed period until the instability develops222We have confirmed that the SWI develops for the spiral waves driven by a companion in a disk.. This implies a potential significance of the SWI in various astrophysical disks.

7.1. Why Has the Instability Not Been Reported Previously?

At this point, one might wonder why the SWI has not been reported previously if it operates under a broad range of disk conditions as we claim.

First of all, it is possible that the SWI has been present in previous studies, but masked by turbulence driven via other processes. For instance, the recent study of nonlinear spiral waves in a disk excited by a massive planet by Lyra et al. (2016) used high resolution in an essentially inviscid disk, but the violence of the disk response to the strong nonlinear forcing, combined with the presence of convective motions, probably masked the presence of the SWI in that study. Also, we speculate the SWI might have been present in the recent 3D MHD calculation of accretion disks in cataclysmic variable (CV) systems by Ju et al. (2016). In this study, however, the SWI must be obscured by the vigorous MRI turbulence.

Considering the reference run R512, this was computed using a logarithmic radial grid, where the grid cell size at was . Similarly, the run undertaken by NIRVANA, described in the appendix, used a uniform radial mesh with the same grid spacing, corresponding to grid cells. This is higher resolution than is usually undertaken in global simulations of giant planets embedded in protoplanetary disks, for example, although high resolution in the vicinity of the planet is often achieved with local mesh refinement (e.g Kley et al., 2001; Klahr & Kley, 2006; Gressel et al., 2013; Szulágyi et al., 2014). Viscosity, with , is also often included in these models, which has the effect of partially damping the instability. We note that simulations with low-mass planets are unlikely to show growth of the SWI within the typical durations of three-dimensional calculations, since the smaller the mass of companion, the weaker the instability and the longer it takes for it to develop (see Figure 11). In addition, the primary and secondary arm separation is known to be less than for low-mass planets (e.g. Fung & Dong, 2015; Zhu et al., 2015), in which case the resonance between the spiral arms and inertial modes is not exact. In this case, the SWI may grow at a reduced rate, requiring a longer simulation duration to observe the instability develop.

In Figure 16, we present two-dimensional maps showing the vertical and radial wavelengths of the inertial modes that satisfy the resonance conditions with the imposed spiral waves. The wavelengths are calculated using Equations (10) and (13), using the background disk structure and computational domain used in the adiabatic simulation labelled as GAM4 in Table 1. Purely for the purposes of illustration, we present maps for the inertial modes predicted to arise when , where is defined as . However, since the second term in Equation (13) is typically smaller than the first term, the contour plots can be linearly scaled to other modes333We confirmed that this is invalid only in a very narrow region near the boundary between the SWI forbidden and permitted regions where . as indicated in the color bar labels.

The maps imply that there is a resolution limit on which modes can be represented by the computational grid. For , the midplane region interior to is poorly resolved by grid cells in both the vertical and radial directions, indicating that these inertial modes cannot be well represented on the mesh there. Smaller values of are hence expected to be associated with the modes seen to grow in this region, in agreement with the results described in Section 4. Moving out to the midplane region beyond , the vertical and radial wavelengths of the excited modes start to exceed grid cell spacings, suggesting that these modes can be represented on the grid.

7.2. Implications

7.2.1 Angular Momentum Transport

In order to estimate the rate of angular momentum transport induced by the turbulent flow arising from the SWI, we calculate the Shakura-Sunyaev stress parameter for the reference model R512 according to where is a density-weighted mean pressure at radius . The simple arithmetic average of over in between and 2 is plotted in Figure 17 as a function of time. We see that the advected flux of angular momentum due to the spiral waves at the beginning of the simulation gives rise to an apparent Reynolds stress . This is not the accretion stress experienced by the disk, but instead represents the negative angular momentum flux associated with the spiral wave as it propagates inwards, and only that fraction of the wave angular momentum that is deposited in the gas through wave dissipation acts to drive accretion. Once the SWI develops we see that the correlated velocity fluctuations generate a sustained accretion stress . This value is naturally smaller than that associated with the unattenuated spiral wave because the SWI operates over a range of radii in the disk, and the energy and angular momentum associated with the wave are injected into the disk matter over this radial range via the breakdown into turbulence.

Figure 17.— Time evolution of Shakura-Sunyaev parameter , averaged over in between and 2.
Figure 18.— Velocity vectors in a plane at from the R512 model. The map shows only to due to large velocity in the upper layer. Note the turbulent eddies generated by the spiral wave instability. A velocity vector with is shown on the lower-left corner.

7.2.2 Vertical Mixing

In Figure 18, we present velocity vectors in a vertical plane at the end of the reference run. As can be seen, the SWI creates a set of turbulent eddies that will be effective at inducing vertical mixing of the gas, and any dust particles that it contains. It is also the case that forces arising from the turbulent density fluctuations will induce stochastic migration and eccentricity/inclination excitation of larger bodies, such as planetesimals and protoplanets, in a protoplanetary disk in which nonlinear spiral waves propagate, similar to what is observed in disks that sustain magneto-rotational turbulence (Nelson & Papaloizou, 2004; Nelson, 2005). We estimate the vertical diffusion coefficient associated with the correlated vertical velocity fluctuations using the approximation , where is the correlation time of the vertical velocity fluctuations, (Fromang & Papaloizou, 2006). In obtaining an estimate for , we generate a time series for at different locations in the disk, defined by the region at , and filter these time series to remove the sinusoidal component induced by the spiral wave. We then compute the autocorrelations of these time series, to which we fit the function , leading to an average of the measured correlation times T. Using natural units, the diffusion coefficient is estimated to be , such that the time scale for vertical mixing orbits.

Particles in the Epstein drag regime have a settling time that can be expressed as orbits. is the stopping time defined by (Weidenschilling, 1977), where is the internal density of the grain particles (typically  g cm), is the grain size and is the thermal velocity of gas molecules. Assuming that particles are significantly mixed by the turbulence when , we expect this to occur for particles with Stokes numbers . This corresponds to a particle size cm at 1 AU and mm at 5 AU in the midplane of a minimum mass solar nebula (Hayashi, 1981).444We note that the mean free path of molecules is  cm at 1 AU, so using the Epstein drag formula here is only marginally justified for 2 cm-sized pebbles.

7.2.3 Non-axisymmetric Spiral Features

Recent near-IR scattered light observations have revealed complex spiral arm morphologies in some protoplanetary disks (e.g. Garufi et al., 2013; Benisty et al., 2015; Garufi et al., 2016), thanks to the advent of extreme adaptive optics. The scattered light is believed to trace (sub-)m-sized dust grains in the disk atmosphere since protoplanetary disks are highly optically thick at near-IR wavelength. Interestingly, even in nearly face-on systems exhibiting well-defined spiral arms, the brightness of the spiral arms at the same distance from the central object significantly differ from each other (e.g. SAO 206462; Garufi et al. 2013, MWC 758; Benisty et al. 2015).

In Figure 19, we present two-dimensional distributions of the azimuthal density variation at various heights in the disk, and in Figure 20, we show how the perturbed density varies with azimuth at different disk heights. We see that spiral waves are disturbed by the instability at all heights, but with more prominent disruption occurring higher in the disk. This will obviously have significant implications when attempting to predict the appearance of spiral waves for comparison with observations (e.g. Dong et al., 2016)

Figure 19.— Distribution of at in the (left) plane, (middle) plane, and (right) plane for R512 model.
Figure 20.— Azimuthal distributions of the perturbed density at different heights in the disk at . The solid curves show the distributions before the SWI sets in, and the dashed curves show the distributions when the instability is saturated.

7.3. Future Work

Our results suggest that the SWI is robust, and possibly operates under a variety of disk conditions. The simplified disk models assumed in the present work help to illuminate the nature of the instability, but the significance of the instability will have to be tested with more realistic disk models.

The influence of thermodynamic effects on the instability will need to be tested using a more realistic model for the thermal evolution, including processes such as radiation transport and external irradiation. We have shown that the SWI is likely to operate in both optically thin (isothermal) and thick (adiabatic) disk regions. On the other hand, the outcome of the instability is quite sensitive to the disk thermal model, as we showed by varying the adiabatic index. Real disks are likely to behave adiabatically near the midplane and isothermally near the surface, so the outcome of t-he instability in such disks will be more complicated than shown by the models we have presented here. In particular, the forbidden zones, where inertial waves cannot be excited, will depend on the detailed thermodynamics due to the dependency on the local Brunt-Väsälä frequency.

Our simulations have shown that the SWI can operate in viscous disks, suggesting that it might also operate in disks that sustain turbulence through the magnetorotational instability (MRI) (Balbus & Hawley, 1991). It certainly appears to be the case that the instability could operate in dead zones (Gammie, 1996), and regions where non ideal MHD effects such as ambipolar diffusion act to moderate the strength of MRI turbulence (Bai, 2015). Precisely how the SWI operates under these conditions will require high resolution MHD simulations of spiral waves propagating in magnetized disks.

In addition to understanding how the inclusion of additional physical processes could modify the SWI, future work is needed to examine how the SWI changes our understanding of astrophysical phenomena where spiral waves play an important role. The list of these includes the following:
– Spiral wave propagation in disks with external binary companions, such as protoplanetary disks, CV systems and X-ray binaries. Issues of interest include how efficiently the spiral waves are damped by the SWI and how this affects the wave amplitude as a function of position in the disk, and whether or not the turbulent stresses associated with the SWI modify the outburst cycles associated with CV systems. In recent work, Ju et al. (2016) showed that spiral shocks can co-exist with MRI turbulence in CV systems, so it will be interesting to investigate applications of the SWI in the systems.
– FU Orionis outbursts driven by spiral waves excited by self-gravity in disks around young stars. The idea here is that spiral waves launched in the outer disk during the infall phase can heat the inner disk sufficiently to ionize it, and drive accretion onto the star by switching on the MRI (Zhu et al., 2010; Bae et al., 2014). We came across the SWI while conducting 3-D simulations designed to address this very issue, and in a forthcoming paper we will present a study that examines the viability of this picture.
– Planet formation in disks with external binary companions. Numerous exoplanets are known to orbit one star that is a member of a binary system, perhaps the most famous of these being Cephei (Hatzes et al., 2003). The excitation of spiral waves in a protoplanetary disk by an external binary is already known to have a strong influence on planet formation (Kley & Nelson, 2008; Paardekooper et al., 2008). The excitation of turbulence by the SWI will also have an important influence, particularly on the settling of dust grains and the growth of small particles.
– Formation of circumbinary planets. The situation regarding the growth of the SWI for spiral waves that propagate outwards from an outer Lindblad resonance is different than for an inward propagating wave. Ignoring the effects of buoyancy, and noting that inertial wave frequencies must lie in the range , we infer that the resonant excitation of inertial waves by an inner companion occurs in a relatively narrow range of radii where the doppler-shifted frequency of the spiral . A fluid element at large radius in a disk sees a large value for , such that the resonance condition cannot be satisfied, unless buoyancy forces are also included in a disk where the Brunt-Väsälä frequency is large. This suggests that there may be a range of radii in circumbinary disks where the SWI operates, making it difficult for small particles to grow there, and for planets to form in situ. The details will depend on the disk thermodynamics.
– Formation of planets in the presence of a giant planet. It is generally believed that Jupiter was the first planet to form in the Solar System, and its early presence has been invoked to explain a number of physical and dynamical features that are observed today, such as the dynamical state of the asteroid belt and the small mass of Mars (Walsh et al., 2011). It will be of interest to examine how the onset of the SWI, induced by the spiral waves excited by a growing Jupiter, influenced the orbits of asteroids though stochastic forcing, and modified the growth of the terrestrial planets or their precursor embryos. The influence could be strong if this occurred mainly through chondrule/pebble accretion, as explored recently by Johansen et al. (2015) and Levison et al. (2015).

8. Conclusion

We have presented the results of high resolution, 3-D simulations of circumstellar disks that are perturbed by two-armed spiral density waves, and we have shown that a broad range of disk models are subject to a parametric instability involving the excitation of pairs of inertial waves that interact resonantly with the spiral wave. This spiral wave instability (SWI) gives rise to turbulence that transports angular momentum and causes vertical mixing. The apparent robustness of the SWI under changes to physical conditions in the disks suggests that it may arise in a broad range of astrophysical settings. Future work is required to understand and evaluate its influence on the physical evolution and observational appearance of these systems.

Authors thank the anonymous referee for a helpful report that improved the initial manuscript. J.B. thanks Fred Adams and Steve Lubow for valuable conversations, and Charles Gammie for insightful comments on the initial manuscript. This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor, and by HPC resources provided by IT Services at Queen Mary University of London. 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 and the Texas Advanced Computing Center at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC Operations grant ST/K0003259/1. DiRAC is part of the national E-Infrastructure.

Appendix A Vertically Stratified, Locally Isothermal Models

We introduce locally isothermal disk models, in which a radial temperature gradient is imposed by assuming and , and an isothermal equation of state is adopted. With this setup, the disk aspect ratio maintains a constant value at all radii. Since the radial temperature gradient produces a non-zero vertical gradient in the disk rotation (see Equation 21), this model is susceptible to the vertical shear instability (Nelson et al., 2013). We implement the density slope , so the initial disk structure is identical to the T1R-0 and T1R-0-3D models of Nelson et al. (2013). The model parameters are summarized in Table 2.

In order to check whether or not the vertical shear instability develops in the absence of imposed spiral waves, we first run calculations with the spiral potential turned off (). We vary the kinematic viscosity using the values . As inferred from the time evolution of the integrated meridional kinetic energy presented in Figure 21, the vertical shear instability develops when but is suppressed with , which is in good agreement with Nelson et al. (2013). The small fluctuations at in the VISO-V6 model and at in the VISO-V5 model are signatures of the vertical shear instability, but note that they do not grow over time and damp out because of viscous dissipation.

Having confirmed that the disk is stable to the vertical shear instability with , we now introduce spiral waves into the disk. We test with four different kinematic viscosity values: , and . As shown in Figure 21, the result is consistent with the globally isothermal models: the spiral wave instability is triggered with and is completely suppressed with . Also, the saturated energy level is smaller with larger viscosity, as expected.

We run the model with the four codes introduced in Section 3.5. As we discussed earlier, the detailed evolution is not identical because of the difference in dissipative properties of the codes, and the fact that NIRVANA and INABA3D use uniform radial meshes, whereas FARGO3D and PLUTO use logarithmically spaced grid cells. Nevertheless, the results show reasonably good agreement to each other, and all codes reproduce the SWI.

Run label Code Numerical Resolution Kinematic
( ) Viscosity
VISO-V0 FARGO3D -1 0 1.0
VISO-V6 FARGO3D -1 1.0
VISO-V5 FARGO3D -1 1.0
VISO-V4 FARGO3D -1 1.0
Table 2Locally Isothermal Model Parameters
Figure 21.— Time evolution of the meridional kinetic energy for vertically isothermal models with different kinematic viscosities. The solid curves show the results with the spiral potential amplitude of , while the dashed curves show the results with to check whether the vertical shear instability develops in the absence of spiral waves. The vertical shear instability is suppressed with , which is in good agreement with Nelson et al. (2013, see their Figure 9). The spiral wave instability is suppressed with as in the globally isothermal models.
Figure 22.— Time evolution of the meridional kinetic energy for the vertically isothermal models with different codes (VISO-V6, VISO-N, VISO-P, VISO-I).


  • Armitage et al. (2001) Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705
  • Artymowicz & Lubow (1994) Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651
  • Bae et al. (2014) Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61
  • Bai (2015) Bai, X.-N. 2015, ApJ, 798, 84
  • Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  • Barker & Latter (2015) Barker, A. J., & Latter, H. N. 2015, MNRAS, 450, 21
  • Barker & Ogilvie (2014) Barker, A. J., & Ogilvie, G. I. 2014, MNRAS, 445, 2637
  • Bate et al. (2002) Bate, M. R., Ogilvie, G. I., Lubow, S. H., & Pringle, J. E. 2002, MNRAS, 332, 575
  • Benisty et al. (2015) Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6
  • Benítez-Llambay & Masset (2016) Benítez-Llambay, P., & Masset, F. 2016, ApJS, 223, 11
  • Bryden et al. (1999) Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., & Papaloizou, J. C. B. 1999, ApJ, 514, 344
  • de Val-Borro et al. (2006) de Val-Borro, Edgar, R. G., M., Artymowicz, P., et al.  2006, MNRAS, 370, 529
  • Dong et al. (2016) Dong, R., Fung, J., & Chiang, E. 2016, arXiv:1602.04814
  • Fromang & Papaloizou (2006) Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751
  • Fromang & Papaloizou (2007) Fromang, S., & Papaloizou, J. 2007, A&A, 468, 1
  • Fung & Dong (2015) Fung, J., & Dong, R. 2015, ApJ, 815, L21
  • Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
  • Gammie (1999) Gammie, C. F. 1999, Astrophysical Discs - an EC Summer School, 160, 122
  • Gammie et al. (2000) Gammie, C. F., Goodman, J., & Ogilvie, G. I. 2000, MNRAS, 318, 1005
  • Garufi et al. (2013) Garufi, A., Quanz, S. P., Avenhaus, H., et al. 2013, A&A, 560, A105
  • Garufi et al. (2016) Garufi, A., Quanz, S.P., & Schmid, H. M. 2016, A&A, 588, A8
  • Goldreich & Tremaine (1978) Goldreich, P., & Tremaine, S.  1978, ApJ, 222, 850
  • Goldreich & Tremaine (1979) Goldreich, P., & Tremaine, S.  1979, ApJ, 233, 857
  • Goodman (1993) Goodman, J. 1993, ApJ, 406, 596
  • Gressel et al. (2013) Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59
  • Hatzes et al. (2003) Hatzes, A. P., Cochran, W. D., Endl, M., et al. 2003, ApJ, 599, 1383
  • Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  • Inaba et al. (2005) Inaba, S., Barge, P., Daniel, E., & Guillard, H.  2005, A&A, 431, 365
  • Johansen et al. (2015) Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, SciA, 1, 1500109
  • Ju et al. (2016) Ju, Wenhua, Stone, J. M., & Zhu, Z.  2016, ApJ, in press
  • Klahr & Kley (2006) Klahr, H., & Kley, W. 2006, A&A, 445, 747
  • Kley (1999) Kley, W. 1999, MNRAS, 303, 696
  • Kley et al. (2001) Kley, W., D’Angelo, G., & Henning, T. 2001, ApJ, 547, 457
  • Kley & Nelson (2008) Kley, W., & Nelson, R. P. 2008, A&A, 486, 617
  • Laughlin & Bodenheimer (1994) Laughlin, G., & Bodenheimer, P. 1994, ApJ, 436, 335
  • Lesur et al. (2015) Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9
  • Levison et al. (2015) Levison, H. F., Kretke, K. A., Walsh, K. J., & Bottke, W. F. 2015, Proceedings of the National Academy of Science, 112, 14180
  • Lin & Shu (1964) Lin, C. C., & Shu, F. H.  1964, ApJ, 140, 646
  • Lin & Papaloizou (1979) Lin, D. N. C., & Papaloizou, J. 1979, MNRAS, 186, 799
  • Lin et al. (1990) Lin, D. N. C., Papaloizou, J. C. B., & Savonije, G. J. 1990, ApJ, 364, 326
  • Lubow & Ogilvie (1998) Lubow, S. H., & Ogilvie, G. I.  1998, ApJ, 504, 983
  • Lubow & Pringle (1993) Lubow, S. H., & Pringle, J. E.  1993, ApJ, 409, 360
  • Lynden-Bell & Kalnajs (1972) Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1
  • Lyra et al. (2016) Lyra, W., Richert, A. J. W., Boley, A., et al. 2016, ApJ, 817, 102
  • Masset (2000) Masset, F. 2000, A&AS, 141, 165
  • Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al.  2007, ApJS, 170, 228
  • Nelson et al. (2000) Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18
  • Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M.  2013, MNRAS, 435, 2610
  • Nelson & Papaloizou (2004) Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849
  • Nelson (2005) Nelson, R. P. 2005, A&A, 443, 1067
  • Ogilvie & Latter (2013) Ogilvie, G. I., & Latter, H. N. 2013, MNRAS, 433, 2420
  • Paardekooper et al. (2008) Paardekooper, S.-J., Thébault, P., & Mellema, G. 2008, MNRAS, 386, 973
  • Papaloizou & Savonije (1991) Papaloizou, J. C., & Savonije, G. J. 1991, MNRAS, 248, 353
  • Papaloizou & Terquem (1995) Papaloizou, J. C. B., & Terquem, C.  1995, MNRAS, 274, 987
  • Papaloizou (2005a) Papaloizou, J. C. B.  2005a, A&A, 432, 743
  • Papaloizou (2005b) Papaloizou, J. C. B.  2005b, A&A, 432, 757
  • Richard et al. (2013) Richard, S., Barge, P., & Le Dizès, S.  2013, A&A, 559, A30
  • Ryu & Goodman (1994) Ryu, D., & Goodman, J. 1994, ApJ, 422, 269
  • Sawada et al. (1986) Sawada, K., Matsuda, T., & Hachisu, I. 1986, MNRAS, 219, 75
  • Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  • Stone & Norman (1992) Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753
  • Szulágyi et al. (2014) Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65
  • van Leer (1977) van Leer, B. 1977, Journal of Computational Physics, 23, 276
  • Walsh et al. (2011) Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. M. 2011, Nature, 475, 206
  • Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  • Ziegler & Yorke (1997) Ziegler, U., & Yorke, H. W. 1997, Computer Physics Communications, 101, 54
  • Zhu et al. (2010) Zhu, Z., Hartmann, L., Gammie, C. F., et al. 2010, ApJ, 713, 1134
  • Zhu et al. (2015) 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