SELF-DESTRUCTING SPIRAL WAVES: GLOBAL SIMULATIONS OF A SPIRAL-WAVE INSTABILITY IN ACCRETION DISKS

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 ${\omega }_{{\rm{i,1}}}+{\omega }_{{\rm{i,2}}}={\omega }_{{\rm{s}}}$ 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 m (in the absence of self-gravity) is given by (assuming k_Z = 0)

$$\begin{eqnarray}&&{m}^{2}{({\rm{\Omega }}-{{\rm{\Omega }}}_{{\rm{p}}})}^{2}={\kappa }^{2}+{c}_{s}^{2}{k}_{R,{\rm{s}}}^{2},\end{eqnarray} \tag{ 9 }$$

giving a radial wave number interior to the ILR

$$\begin{eqnarray}&&{k}_{R,{\rm{s}}}^{2}=\displaystyle \frac{{m}^{2}{({\rm{\Omega }}-{{\rm{\Omega }}}_{{\rm{p}}})}^{2}-{\kappa }^{2}}{{c}_{s}^{2}}.\end{eqnarray} \tag{ 10 }$$

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) ${k}_{R,{\rm{i}}}={{nk}}_{R,{\rm{s}}}$, where n 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

$$\begin{eqnarray}&&\displaystyle \frac{{k}_{Z,{\rm{i}}}^{2}}{{\omega }_{{\rm{i}}}^{2}-{N}^{2}}+\displaystyle \frac{{k}_{R,{\rm{i}}}^{2}}{{\omega }_{{\rm{i}}}^{2}-{\kappa }^{2}}-\displaystyle \frac{{\omega }_{{\rm{i}}}^{2}/{c}_{s}^{2}}{{\omega }_{{\rm{i}}}^{2}-{N}^{2}}=0,\end{eqnarray} \tag{ 11 }$$

giving

$$\begin{eqnarray}&&{k}_{Z,{\rm{i}}}^{2}=-\displaystyle \frac{{n}^{2}{k}_{R,{\rm{s}}}^{2}}{{\omega }_{{\rm{i}}}^{2}-{\kappa }^{2}}({\omega }_{{\rm{i}}}^{2}-{N}^{2})+\displaystyle \frac{{\omega }_{{\rm{i}}}^{2}}{{c}_{s}^{2}}.\end{eqnarray} \tag{ 12 }$$

Substituting the resonance condition ω_i = ω_s/2 and ${k}_{R,{\rm{s}}}^{2}=({\omega }_{{\rm{s}}}^{2}-{\kappa }^{2})/{c}_{s}^{2}$ into Equation (12) gives

$$\begin{eqnarray}&&{k}_{Z,{\rm{i}}}^{2}=\displaystyle \frac{{({\omega }_{{\rm{s}}}/2)}^{2}}{{c}_{s}^{2}}-{n}^{2}\displaystyle \frac{{\omega }_{{\rm{s}}}^{2}-{\kappa }^{2}}{{c}_{s}^{2}}\displaystyle \frac{{({\omega }_{{\rm{s}}}/2)}^{2}-{N}^{2}}{{({\omega }_{{\rm{s}}}/2)}^{2}-{\kappa }^{2}}.\end{eqnarray} \tag{ 13 }$$

First, we consider the behavior of ${k}_{Z,{\rm{i}}}$ for disks in which N² = 0:

$$\begin{eqnarray}&&{k}_{Z,{\rm{i}}}^{2}=\displaystyle \frac{{({\omega }_{{\rm{s}}}/2)}^{2}}{{c}_{s}^{2}}\left(1-{n}^{2}\displaystyle \frac{{\omega }_{{\rm{s}}}^{2}-{\kappa }^{2}}{{({\omega }_{{\rm{s}}}/2)}^{2}-{\kappa }^{2}}\right).\end{eqnarray} \tag{ 14 }$$

For ${k}_{Z,{\rm{i}}}$ to be real, such that inertial waves can be excited, we require that ${\omega }_{{\rm{s}}}^{2}\gt {\kappa }^{2}$ and ${({\omega }_{{\rm{s}}}/2)}^{2}\lt {\kappa }^{2}$, and these conditions are both satisfied at all radii interior to the ILR of an m = 2 spiral wave in a Keplerian disk. (These conditions, however, are not satisfied everywhere for a spiral wave that propagates outward from an outer Lindblad resonance, as we discuss later in Section 7.) Considering the influence of buoyancy in Equation (13), we see that ${k}_{Z,{\rm{i}}}$ can become imaginary when ${({\omega }_{{\rm{s}}}/2)}^{2}\lt {N}^{2}$, indicating that inertial modes will not be excited at certain locations above the midplane where the buoyancy frequency is large. Taking the large n limit, the transition occurs at ${({\omega }_{{\rm{s}}}/2)}^{2}={N}^{2}$. Given that the resonance condition for the excited inertial waves is ω_i = (ω_s/2), this indicates that inertial modes cannot be excited where ω_i < N. 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, which applies to modes where ${k}_{Z}H\gg 1$ and ${k}_{R}H\gg 1$, is not strictly applicable to modes with k_R ∼ 1/H, 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 toward the midplane as they enter disk radii where large fractions of the disk vertical domain have ${N}^{2}\gt {\omega }_{{\rm{i}}}^{2}$, 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.

In linear theory, a two-dimensional spiral wave without any vertical structure, launched at an ILR in a vertically isothermal disk with sound speed independent of Z, should propagate toward the star with wave fronts remaining perpendicular to the midplane. As described in the next section, the underlying axisymmetric disk model will normally have a midplane density varying according to $\rho {(R,0)=\rho ({R}_{0},0)(R/{R}_{0})}^{p}$, with p < 0, and a variation with height scaling as $\rho (R,Z)\sim \rho (R,0)\exp (-{Z}^{2}/(2{H}^{2}))$, where H 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 toward 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.

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho v)=0,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\rho \left(\displaystyle \frac{\partial v}{\partial t}+v\cdot {\rm{\nabla }}v\right)=-{\rm{\nabla }}P-\rho {\rm{\nabla }}({{\rm{\Phi }}}_{* }+{{\rm{\Phi }}}_{p})+{\rm{\nabla }}\cdot {\rm{\Pi }},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{\rm{\nabla }}\cdot ({ev})=-P{\rm{\nabla }}\cdot v,\end{eqnarray} \tag{ 17 }$$

where ρ is the mass density, v is the velocity, P is the pressure, Φ_* is the gravitational potential due to the central object, Φ_p is the spiral potential (see Section 3.3), Π is the viscous stress tensor, and e is the internal energy per unit volume.

We make use of both cylindrical (R, ϕ, Z) and spherical (r, θ, ϕ) coordinates. We consider cylindrical disk models where we neglect the vertical component of gravity from the central object, so the gravitational potential is Φ_* = GM_*/R, where G is the gravitational constant, M_* is the mass of the central object, and R 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 Φ_* = GM_*/r, where $r=\sqrt{{R}^{2}+{Z}^{2}}$ 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 $P=\rho {c}_{s}^{2}$, where c_s 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 $P=(\gamma -1)e$, 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 γ ∼ 1, 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 γ = 1.4. 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.

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

$$\begin{eqnarray}&&T(R)={T}_{0}{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{q},\end{eqnarray} \tag{ 18 }$$

where T₀ is the temperature at R₀ = 1. In all models presented, T₀ is chosen such that the ratio of disk scale height to radius at R₀ is H₀/R₀ = 0.05. The isothermal sound speed is related to the temperature by ${c}_{s}^{2}={ \mathcal R }T/\mu$, where ${ \mathcal R }$ is the gas constant and μ is the mean molecular weight. Thus, Equation (18) corresponds to the radial sound speed distribution given by

$$\begin{eqnarray}&&{c}_{s}(R)=\displaystyle \frac{{H}_{0}}{{R}_{0}}{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{q/2}.\end{eqnarray} \tag{ 19 }$$

For the main simulation set we assume q = 0, 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 (q = −1). In those models, we find that the SWI and the vertical shear instability co-exist 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):

$$\begin{eqnarray}&&\rho (R,Z)={\rho }_{0}{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{p}\exp \left(\displaystyle \frac{{{GM}}_{* }}{{c}_{s}^{2}}\left[\displaystyle \frac{1}{\sqrt{{R}^{2}+{Z}^{2}}}-\displaystyle \frac{1}{R}\right]\right)\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&{v}_{\phi }(R,Z)\\ &&\,={\left[(1+q)\displaystyle \frac{{{GM}}_{* }}{R}+(p+q){c}_{s}^{2}-q\displaystyle \frac{{{GM}}_{* }}{\sqrt{{R}^{2}+{Z}^{2}}}\right]}^{1/2}.\end{eqnarray} \tag{ 21 }$$

The midplane density at R = R₀, ρ₀, is chosen such that the total disk mass is 10% of the central object's mass. When vertical stratification of the disk is ignored in the cylindrical disk models, the above equations simplify to

$$\begin{eqnarray}&&\rho (R,Z)={\rho }_{0}{\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{p}\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&{v}_{\phi }(R,Z)={\left[\displaystyle \frac{{{GM}}_{* }}{R}+(p+q){c}_{s}^{2}\right]}^{1/2}.\end{eqnarray} \tag{ 23 }$$

We fix the radial power-law index of the midplane gas density to p = −1.5 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 ${10}^{-6}\,{c}_{s}$, 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 0.5 to 10 in radius and from 0 to 2π in azimuth. In cylindrical models, the vertical domain extends from Z = −0.1 to 0.1. In spherical models, the meridional domain covers $\pm 4{H}_{0}/{R}_{0}$ above and below the midplane. The model parameters are summarized in Table 1.

Table 1.  Model Parameters

Run Label	Code	Numerical Resolution	${ \mathcal A }$	Kinematic	γ
		(NR or Nr × Nϕ × Nz or Nθ)		Viscosity ν
CYL-F	FARGO3D	512 × 128 × 128	5.0 × 10−4	0	1.0
CYL-N	NIRVANA	3200 × 128 × 128	5.0 × 10−4	0	1.0
CYL-P	PLUTO	512 × 128 × 128	5.0 × 10−4	0	1.0
CYL-I	INABA3D	3200 × 128 × 128	5.0 × 10−4	0	1.0
R512	FARGO3D	512 × 128 × 128	5.0 × 10−4	0	1.0
R128	FARGO3D	128 × 32 × 32	5.0 × 10−4	0	1.0
R256	FARGO3D	256 × 64 × 64	5.0 × 10−4	0	1.0
R768	FARGO3D	768 × 192 × 192	5.0 × 10−4	0	1.0
AMP0.625	FARGO3D	512 × 128 × 128	6.25 × 10−5	0	1.0
AMP1.25	FARGO3D	512 × 128 × 128	1.25 × 10−4	0	1.0
AMP2.5	FARGO3D	512 × 128 × 128	2.5 × 10−4	0	1.0
AMP10	FARGO3D	512 × 128 × 128	1.0 × 10−3	0	1.0
V6	FARGO3D	512 × 128 × 128	5.0 × 10−4	10−6	1.0
V5	FARGO3D	512 × 128 × 128	5.0 × 10−4	10−5	1.0
V4	FARGO3D	512 × 128 × 128	5.0 × 10−4	10−4	1.0
GAM01	FARGO3D	512 × 128 × 128	6.25 × 10−5	0	1.01
GAM05	FARGO3D	512 × 128 × 128	6.25 × 10−5	0	1.05
GAM1	FARGO3D	512 × 128 × 128	6.25 × 10−5	0	1.1
GAM2	FARGO3D	512 × 128 × 128	6.25 × 10−5	0	1.2
GAM3	FARGO3D	512 × 128 × 128	6.25 × 10−5	0	1.3
GAM4	FARGO3D	512 × 128 × 128	6.25 × 10−5	0	1.4

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{p}(R,\phi ,t)={ \mathcal A }\cos [m(\phi -{{\rm{\Omega }}}_{p}t)]{e}^{-{(R-{R}_{p})}^{2}/{\sigma }_{p}^{2}},\end{eqnarray} \tag{ 24 }$$

where ${ \mathcal A }$ is the amplitude, assumed to be constant over time, m is the azimuthal mode number, Ω_p is the pattern speed, R_p is the radius about which the potential is centered, and σ_p is the radial width of the potential. We adopt the values R_p = 5 and ${\sigma }_{p}=1$ and assume that the pattern speed is the local Keplerian frequency at its central position: ${{\rm{\Omega }}}_{p}={({{GM}}_{* }/{R}_{p}^{3})}^{1/2}$. In this work, we only focus on spiral potentials with m = 2.

Our fiducial model assumes a potential amplitude of ${ \mathcal A }=5.0\times {10}^{-4}$. With this amplitude we intend to produce velocity perturbations at R = 1 that are about ∼20% of the local sound speed.³ The effect of varying the spiral potential strength on the SWI 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 ${ \mathcal A }=5.0\times {10}^{-4}$ and ${ \mathcal A }=6.25\times {10}^{-5}$. 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 ${ \mathcal A }=6.25\times {10}^{-5}$ is nearly sinusoidal. In the inner disk, the velocity amplitude is about a few percent of the sound speed and is nearly constant over radius. When the spiral potential amplitude is increased to ${ \mathcal A }=5\times {10}^{-4}$, on the other hand, the waves are rather strong, generating sharp wave fronts that have radial velocities up to about 40% of the sound speed. At r = 0.6, the two spiral potential amplitudes generate relative density enhancements of ∼7% (${ \mathcal A }=6.25\times {10}^{-5}$) and ∼37% (${ \mathcal A }=5\times {10}^{-4}$) 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.

Zoom In Zoom Out Reset image size

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 R = [0.5, 0.6] and R = [9, 10]. We choose one local orbital time for the wave-damping timescale. As the simulation domain covers 2π 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:

$$\begin{eqnarray}&&\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial }{\partial \theta }(\rho {c}_{s}^{2})=\displaystyle \frac{{v}_{\phi }^{2}}{\tan \theta }.\end{eqnarray} \tag{ 25 }$$

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.

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). 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 (Stone & Norman 1992; Ziegler & Yorke 1997). 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 cell 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 Δr = 0.002934, and this value was used for the uniform meshes in the NIRVANA and INABA3D runs.

We choose code units such that GM_* = R₀ = Ω(R = R₀) = 1. 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 R = 1 as the time unit. We will denote this quantity as $1\,{T}_{{\rm{orb}}}$. One orbital time at the spiral potential center R = R_p, T_p, corresponds to 11.2 T_orb.

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

$$\begin{eqnarray}&&{e}_{Z,\theta }=\displaystyle \frac{1}{2}{\int }_{V}\rho {v}_{Z,\theta }^{2}\,{dV}.\end{eqnarray} \tag{ 26 }$$

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

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 $\sim 1.2\mbox{--}1.4\,{T}_{{\rm{orb}}}^{-1}$. The saturated perturbed energy levels also show good agreement with each other, being within a factor ∼3 at t = 250.

We note, however, that the early evolution of the perturbed kinetic energy (t ≲ 30) seems to be more code dependent than the later evolution. In particular, the exponential growth seen at 5 ≲ t ≲ 15 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.

Zoom In Zoom Out Reset image size

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.

In Figure 8, we present the time evolution of the integrated meridional kinetic energy e_θ. Contour plots of the meridional kinetic energy density $\rho {v}_{\theta }^{2}$ are plotted in Figure 9. We make use of the meridional kinetic energy density $\rho {v}_{\theta }^{2}$ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As shown in Figure 8, e_θ increases rapidly very early in the simulations at 10 ≲ t ≲ 15. 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 toward the star. As discussed in Section 2, this seems to arise as a nonlinear effect due to 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 inward. 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 t ∼ 20, and e_θ saturates at t ∼ 30. 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 (t = 200), 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.

In order to test the effect of numerical resolution, we run the reference model with four different resolutions: (N_r × N_ϕ × N_θ) = (128 × 32 × 32), (256 × 64 × 64), (512 × 128 ×128), and (768 × 192 × 192). 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 $({N}_{r}\times {N}_{\phi }\times {N}_{\theta })=(512\times 128\times 128)$ 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 (512 × 128 × 128) and (768 × 192 × 192) 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 (256 × 64 × 64) and (128 × 32 × 32) 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 the R768 model at t = 22 and 200. Compared to the results obtained with (512 × 128 × 128) grid cells presented in Figure 9, one can see more fine scale structures developing at the beginning of the instability (t = 22). 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 r = 0.6, the spiral waves induce a 39% density enhancement in the midplane in R768 model, whereas the enhancement is 37% 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.

Zoom In Zoom Out Reset image size

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, ${ \mathcal A }=6.25\times {10}^{-5}$, up to the largest value ${ \mathcal A }=1.0\times {10}^{-3}$.

Figure 11 shows the time evolution of e_θ. 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 (t ≲ 20), 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 ${ \mathcal A }=6.25\times {10}^{-5}$, for example, the instability starts at about t = 50, while in the reference model the instability starts almost immediately after the initial propagation of the waves.

Zoom In Zoom Out Reset image size

From ${ \mathcal A }=6.25\times {10}^{-5}$ to ${ \mathcal A }=2.5\times {10}^{-4}$, 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 ${ \mathcal A }\gt 2.5\times {10}^{-4}$. Comparing the ${ \mathcal A }=2.5\times {10}^{-4}$ model with the ${ \mathcal A }={10}^{-3}$ model, the density enhancement increases only by about 25%, presumably as a result of nonlinear dissipation. The nonlinear dissipation is observed over a broad disk region (r ≲ 1.6), and interestingly, it is more significant at smaller radii.

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: ν = 10⁻⁶, 10⁻⁵, and 10⁻⁴. In the canonical α prescription of Shakura & Sunyaev (1973), the kinematic viscosity ν corresponds to

$$\begin{eqnarray}&&\alpha (R)=1.1\,\times \,{10}^{-3}\left(\displaystyle \frac{\nu }{{10}^{-6}}\right){\left(\displaystyle \frac{R}{0.5}\right)}^{-1.5}.\end{eqnarray} \tag{ 27 }$$

If we assume that radii are measured in units of au, then a value of ν = 10⁻⁵ corresponds to α = 3.9 × 10⁻³ at 1 au, and α = 3.5 × 10⁻⁴ 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 ν ≤ 10⁻⁵, but is completely suppressed when ν = 10⁻⁴. 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 ν = 0 and ν = 10⁻⁶ models is only marginal, suggesting that numerical diffusion operates with a value that is moderately below ν = 10⁻⁶ in the inviscid model.

Zoom In Zoom Out Reset image size

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 ν ≤ 10⁻⁵. With $\nu ={10}^{-4}$, we do not find any signatures of the instability growing throughout the duration of the calculation. Compared with the results of the reference model at t = 24 presented in Figure 9, we find that the small scale unstable modes are absent with ν = 10⁻⁶. With $\nu ={10}^{-5}$, small scale modes are further suppressed and the unstable modes that have grown in the model have noticeably larger length scales than in the ν = 0 and ν = 10⁻⁶ 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 $\sim {k}^{2}\nu$, this demonstrates that we expect the high k (small λ) modes to be preferentially damped.

Zoom In Zoom Out Reset image size

We emphasize that the instability operates with a fairly large viscosity of ν = 10⁻⁵. In terms of the canonical α description, $\nu ={10}^{-5}$ corresponds to α ∼ 0.008 at R = 0.6 (right outside of the inner damping zone) and α ∼ 0.001 at R = 2. The fact that the SWI 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.

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 r = 0.5 was Δr = 0.002934. Similarly, the run undertaken by NIRVANA, described in the Appendix, used a uniform radial mesh with the same grid spacing, corresponding to N_r = 3200 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 ν ≥ 10⁻⁵, 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 3D calculations, since the smaller the mass of the 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 180° 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 labeled as GAM4 in Table 1. Purely for the purposes of illustration, we present maps for the inertial modes predicted to arise when n = 4, where n is defined as ${k}_{R,{\rm{i}}}={{nk}}_{R,{\rm{s}}}$. However, since the second term in Equation (13) is typically smaller than the first term, the contour plots can be linearly scaled to other modes⁵ as indicated in the color bar labels.

Zoom In Zoom Out Reset image size

The maps imply that there is a resolution limit on which modes can be represented by the computational grid. For n = 4, the midplane region interior to R ∼ 0.9 is poorly resolved by ≲8 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 n 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 R = 1, the vertical and radial wavelengths of the excited modes start to exceed ∼10 grid cell spacings, suggesting that these modes can be represented on the grid.

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 ${\alpha }_{r\phi }$ for the reference model R512 according to ${\alpha }_{r\phi }(r,\theta )\equiv \langle \rho \delta {v}_{r}\delta {v}_{\phi }\rangle /\overline{P}(r)$ where $\overline{P}(r)$ is a density-weighted mean pressure at radius r. The simple arithmetic average of ${\alpha }_{r\phi }$ over θ in between r = 0.5 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 ${\alpha }_{r\phi }\sim 0.0016$. 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 inward, 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 ${\alpha }_{r\phi }\sim 5\times {10}^{-4}$. 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 break down into turbulence.

Zoom In Zoom Out Reset image size

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 magnetorotational turbulence (Nelson & Papaloizou 2004; Nelson 2005). We estimate the vertical diffusion coefficient associated with the correlated vertical velocity fluctuations using the approximation ${{ \mathcal D }}_{Z}=\langle {v}_{Z}^{2}\rangle {\tau }_{{\rm{corr}}}$, where τ_corr is the correlation time of the vertical velocity fluctuations, v_Z (Fromang & Papaloizou 2006). In obtaining an estimate for τ_corr, we generate a time series for v_Z at different locations in the disk, defined by the region $| Z| \lt 1H$ at r = 1, 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 $\exp (-t/{\tau }_{{\rm{corr}}})$, leading to an average of the measured correlation times τ_corr ≃ 0.12T_orb. Using natural units, the diffusion coefficient is estimated to be ${{ \mathcal D }}_{Z}\simeq 9\times {10}^{-6}$, such that the timescale for vertical mixing ${t}_{{\rm{mix}}}={H}^{2}/{{ \mathcal D }}_{Z}\simeq 44$ orbits.

Zoom In Zoom Out Reset image size

Particles in the Epstein drag regime have a settling time that can be expressed as ${t}_{{\rm{settle}}}={(2\pi {t}_{{\rm{s}}}{\rm{\Omega }})}^{-1}$ orbits. t_s is the stopping time defined by ${t}_{{\rm{s}}}=({\rho }_{{\rm{grain}}}a)/(\rho {v}_{{\rm{therm}}})$ (Weidenschilling 1977), where ρ_grain is the internal density of the grain particles (typically ∼2 g cm⁻³), a is the grain size, and v_therm is the thermal velocity of gas molecules. Assuming that particles are significantly mixed by the turbulence when t_mix ∼ t_settle, we expect this to occur for particles with Stokes numbers ${\rm{St}}\equiv {t}_{s}{\rm{\Omega }}\leqslant 0.004$. This corresponds to a particle size a ∼ 2 cm at 1 au and a ∼ 2 mm at 5 au in the midplane of a minimum mass solar nebula (Hayashi 1981).⁶

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 wavelengths. 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)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 the 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 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:

• 1.  
  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.
• 2.  
  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 three-dimensional (3D) simulations designed to address this very issue and in a forthcoming paper we will present a study that examines the viability of this picture.
• 3.  
  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.
• 4.  
  Formation of circumbinary planets. The situation regarding the growth of the SWI for spiral waves that propagate outward 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 0 < ω_i ≤ Ω, 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 ω_s ≃ ω_i/2. A fluid element at large radius in a disk sees a large value for ω_s, 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.
• 5.  
  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 through 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).