A&A 468, 1-18 (2007)
DOI: 10.1051/0004-6361:20067044

Properties and stability of freely propagating nonlinear density waves in accretion disks

S. Fromang - J. Papaloizou

Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge, CB3 0WA, UK

Received 29 December 2006 / Accepted 2 March 2007

Abstract
Aims. In this paper, we study the propagation and stability of nonlinear sound waves in accretion disks.
Methods. Using the shearing box approximation, we derive the form of these waves using a semi-analytic approach and go on to study their stability. The results are compared to those of numerical simulations performed using finite difference approaches such as employed by ZEUS as well as Godunov methods.
Results. When the wave frequency is between $\Omega$ and $2 \Omega$ (where $\Omega$ is the disk orbital angular velocity), it can couple resonantly with a pair of linear inertial waves and thus undergo a parametric instability. Neglecting the disk vertical stratification, we derive an expression for the growth rate when the amplitude of the background wave is small. Good agreement is found with the results of numerical simulations performed both with finite difference and Godunov codes. During the nonlinear phase of the instability, the flow remains well organised if the amplitude of the background wave is small. However, strongly nonlinear waves break down into turbulence. In both cases, the background wave is damped and the disk eventually returns to a stationary state. Finally, we demonstrate that the instability also develops when density stratification is taken into account and so is robust.
Conclusions. This destabilisation of freely propagating nonlinear sound waves may be important for understanding the complicated behaviour of density waves in disks that are unstable through the effects of self-gravity or magnetic fields and is likely to affect the propagation of waves that are tidally excited by objects such as a protoplanet or companion perturbing a protoplanetary disk. The nonlinear wave solutions described here as well as their stability properties were also found to be useful for testing and comparing the performance of different numerical codes.

Key words: accretion, accretion disks - methods: numerical - methods: analytical - hydrodynamics

   
1 Introduction

Waves and oscillations are expected to be found in accretion disks as they can be excited by a very large variety of phenomena. For example, tidal excitation can be caused by the influence of a binary companion (Artymowicz & Lubow 1994) or by a planet orbiting inside the disk (Goldreich & Tremaine 1980,1979; Nelson et al. 2000). Spiral density waves can be excited as a result of gravitational instability when the disk mass is large enough (Lin & Shu 1964; Shu 1992; Toomre 1981). It has also been suggested that waves are responsible for the puzzling phenomenon of QPOs (Kato 2001; Arras et al. 2006). Recent numerical simulations of MHD turbulence resulting from the magnetorotational instability (MRI; Balbus & Hawley 1998) have also indicated that density waves could be excited by the turbulent fluctuations of the flow (Gardiner & Stone 2005; Papaloizou et al. 2004).

There have been many studies of the properties of the large variety of waves that can propagate in disks. Most of the early work was done using simplifying assumptions and restricting the analysis to the linear case (Lubow & Pringle 1993; Lubow & Ogilvie 1998; Ogilvie & Lubow 1999; Lin et al. 1990a; Korycansky & Pringle 1995; Lin et al. 1990b). Some other work extended these results to the nonlinear regime, either semi-analytically, but retaining terms up to second order in an expansion in powers of the amplitude (Larson 1990), or numerically in order to study wave dissipation (Bate et al. 2002; Torkelsson et al. 2000).

Despite all these efforts, a complete semi-analytic solution to the wave equation, valid in the nonlinear regime for large amplitude, is still lacking. This would not only be useful for studying the properties of waves, but would also provide a valuable test for numerical codes that are more and more extensively used in the accretion disk community. The goal of this paper is to derive such a solution and to study its properties and stability. Of course, the need to derive a solution that can be found using semi-analytic methods imposes strong constraints: we use an isothermal equation of state, adopt the shearing box approximation (Goldreich & Lynden-Bell 1965) and restrict the analysis to the case where there is no variation in the direction of rotation or axisymmetry. But all these approximations simplify the analysis which can then be done including nonlinearities.

In this paper we are concerned with parametric instability occurring for non linear density waves treated in a local shearing box approximation. However, parametric instability appears also to be generic for global non axisymmetric, non magnetic disk flows such as may occur when the disk is either eccentric (Ogilvie 2001; Goodman 1993; Papaloizou 2005a,b) or warped (Torkelsson et al. 2000). It is also related to the elliptical instability observed in rotating flows that is apparently connected to turbulent breakdown (Kerswell 2002). Such turbulent breakdown does seem to occur for some of the large amplitude non circular flows studied here. All of this makes the phenomenon an important one for study in the general context of the dynamics of accretion disks. Being a local instability that occurs through the interaction of a dense spectrum of inertial waves with a primary flow, modes with arbitrarily small length scales, depending ultimately on the dissipative properties of the system may be excited, providing a challenge for numerical simulations.

The plan of the paper is as follows. In Sect. 2, we derive the equations describing the propagation of nonlinear density waves in the shearing box approximation. The solutions to these equations are compared with the results of 1D numerical solutions. In Sect. 3, we study their stability. When the frequency of the wave lies in the range $[\Omega,2\Omega]$, where $\Omega$ is the orbital angular velocity, it can interact resonantly with two linear inertial waves (whose individual frequencies, being bounded by $\Omega,$ combine to equal the wave frequency) and be unstable to a parametric instability (Torkelsson et al. 2000; Goodman 1993; Papaloizou 2005a). Neglecting vertical density stratification in this section, we derive an expression for the growth rate of the instability and describe the properties of the eigenmodes. These results are compared in Sect. 4 with 2D numerical simulations of vertically unstratified shearing boxes. Both finite difference schemes as used in ZEUS and NIRVANA (see below) as well as Godunov methods are shown to give results that agree very well with each other and with the results of the linear analysis. Finally, in Sect. 5, we investigate the effect of vertical stratification, showing the occurrence of parametric instability in that case and relating the results to the unstratified case. Finally we discuss our results and summarize conclusions in Sect. 6.

   
2 Free waves

2.1 Basic equations

In common with many studies of the local dynamics of accretion disks orbiting about a central mass, we adopt the shearing box approximation (Goldreich & Lynden-Bell 1965). In this approximation, a reference frame rotating with a fixed angular velocity $\Omega,$ which corresponds to the Keplerian angular velocity at a specified radius R₀,is adopted. The origin, at the centre of the box, is located at some specified point on the circle of radius R₀ and a system of Cartesian coordinates (x,y,z)with associated unit vectors  $( \mbox{\boldmath {$i$ }} , \mbox{\boldmath {$j$ }} , \mbox{\boldmath {$k$ }} )$ is used. Gravity is due to the central object only and an expansion of the gravitational acceleration correct to first order in x and z is adopted. Then the equations of continuity and momentum conservation can be written in the form

 

$$\frac{\partial \rho}{\partial t} + \mbox{\boldmath {$\nabla$ }} \mbox{\boldmath {$\cdot$ }} (\rho \mbox{\boldmath {$v$ }} )$$ = 0 , (1)

$$\frac{\partial \mbox{\boldmath {$v$ }} }{\partial t} + ( \mbox{\b... ...box{\boldmath {$\Omega$ }} \mbox{\boldmath {$\times$ }} \mbox{\boldmath {$v$ }} -\frac{1}{\rho} \mbox{\boldmath {$\nabla$ }} P + 3 \Omega^2 x \mbox{\boldmath {$i$ }} - \Omega^2 z \mbox{\boldmath {$k$ }} ,$$ = (2)

where $\rho$ is the density, $\mbox{\boldmath {$v$ }}$ the velocity and P is the pressure which is related to the density by an equation of state (EOS). In this paper we restrict consideration to an isothermal EOS, for which

$$P=\rho c^2 ,$$ (3)

where $c = \Omega H$ is the isothermal sound speed with H being the scale height. In the absence of waves, the equilibrium solution takes the form

$$\rho = \rho_0 \exp\left(-\frac{z^2}{2H^2}\right) ,$$ (4)

$$\mbox{\boldmath {$v$ }} = -\frac{3}{2} \Omega x \mbox{\boldmath {$j$ }} ,$$ (5)

where $\rho_0,$ being the midplane density, is a constant taken to be unity in the following. In some cases we neglect the vertical component of the gravitational acceleration, $-\Omega^2 z~ {\vec {\hat k}},$ in Eq. (2). Then vertical stratification is also neglected and $\rho = \rho_0$ throughout.

   
2.2 The wave equation

Initially, for simplicity, we neglect vertical stratification and assume that the state variables do not vary along the direction of rotation or the vertical direction and thus depend only on x and t.To derive the equations that describe wave propagation in this system, we adopt a Lagrangian description. We define X(X₀,t) to be the position, or x coordinate, of a fluid element that would be at the fixed location x=X₀ in the absence of any disturbance. The y component of the momentum conservation equation can be written

$$\frac{{\rm D} v_y}{{\rm D} t} + 2 \Omega \frac{{\rm D} X}{{\rm D} t} = 0 ,$$ (6)

where the convective derivative is

$$\frac{\rm D}{{\rm D}t}=\frac{\partial }{\partial t} + v_x \frac{\partial }{\partial x} ,$$ (7)

and we have used the fact that

$$v_x=\frac{{\rm D} X}{{\rm D}t}\cdot$$ (8)

Equation (6) can be integrated to give

$$v_y+2\Omega X = f(X_0) .$$ (9)

where f(X₀) is a function of X₀ to be determined.

Applying the same procedure to the equation of motion in the x direction, we find

$$\frac{{\rm D} v_x}{{\rm D} t} + \Omega^2 X = -\frac{1}{\rho} \frac{\partial P}{\partial x} + 2 \Omega f(X_0) ,$$ (10)

where we have used Eq. (9) to eliminate v_y. In the absence of disturbance, v_x=0 and the pressure is uniform. This means we must have

$$f(X_0)=\frac{1}{2} \Omega X_0 .$$ (11)

Thus Eq. (10) becomes

$$\frac{{\rm D}^2 X}{{\rm D} t^2} + \Omega^2 (X-X_0) = -\frac{c^2}{\rho} \frac{\partial \rho}{\partial x} \cdot$$ (12)

In order to express the density in terms of X, we use the conservation of mass. This gives

$$\rho \frac{\partial X}{\partial X_0} \equiv \rho\left( 1+{\partial \xi \over \partial X_0}\right)= \rho_0 ,$$ (13)

where we have introduced the displacement $\xi(X_0,t)$ through the relation $X(X_0,t)=X_0+\xi(X_0,t)$.

2.2.1 Travelling waves

Since the natural boundary conditions for the shearing box are periodic in shearing coordinates, we consider travelling wave solutions that depend only on the phase $\Phi=X_0-Ut$, where U is a constant the phase velocity. Then we have $\xi=\xi(\Phi)$. Using Eq. (13) to eliminate the density, we obtain

$$U^2 \frac{{\rm d}^2 \xi}{{\rm d}\Phi^2} + \Omega^2 \xi = \fra... ...}\xi}{{\rm d}\Phi} )^2} \frac{{\rm d}^2 \xi}{{\rm d}\Phi^2}\ ,$$ (14)

which is the governing equation that describes travelling wave propagation. It has been derived assuming no magnetic field is present. We comment here that a vertical or an azimuthal field (or both) that depend only on $\Phi$ can easily be incorporated into the above analysis. However, these details are beyond the scope of this paper. We comment also that there are some similarities of the above analysis to that of Larson (1990). He considered travelling waves in a disc though not explicitly in a shearing box. However, his analysis focused on spiral waves and only terms of second order or lower in the wave amplitude were retained. This is inexact and inapplicable to the large amplitude waves considered here.

2.3 First integral of motion

By considering a first integral of Eq. (14), it is possible to interpret the problem as corresponding to the motion of a particle in a potential well. To see this, let p be defined through

$$p=\frac{{\rm d}\xi}{{\rm d}\Phi} \cdot$$ (15)

Using this, Eq. (14) can be written

$$\left(U^2 -\frac{c^2}{(1+p)^2}\right) p\frac{{\rm d}p}{{\rm d}\xi}+\Omega^2 \xi =0 .$$ (16)

which can be integrated with respect to $\xi$ to give

$$\frac{1}{2}\Omega^2\xi^2 = \frac{1}{2\Omega^2}\left(U^2 -\fra... ...frac{{\rm d}p}{{\rm d}\Phi}\right)^2 =U_0^2- U^2\tilde{V}(p) ,$$ (17)

where the dimensionless function $\tilde{V}(p)$ is given by

$$\tilde{V}(p)= \frac{p^2}{2}+(c/U)^2\left(\frac{p}{1+p}-\ln(1+p)\right)$$ (18)

and U₀ is a constant. Equation (17) shows by analogy with classical mechanics that the evolution of p as a function of $\Phi$ can be periodic, consisting of oscillations in the wells of the function $\tilde{V}(p),$ should they exist. The bounds of such oscillations are determined by the condition $\tilde{V}(p) = (U_0/U)^2$. By definition, $\xi$ and p should oscillate around p=0 which should be a minimum of V and be such that p > -1. Such a situation is only possible if $c / U \le 1.$ The form of $\tilde{V}(p)$ is shown in Fig. 1. The different curves correspond to (c/U)²=0.005(solid line), 0.02 (dashed line), 0.1 (dotted line), 0.25 (dotted-dashed line) and 0.5 (dotted-dotted-dotted-dashed line). The form of $\tilde{V}(p)$ is similar for all these values of (c/U)².It increases for -1 < p < -1+c/U, decreases for -1+c/U < p < 0 and increases again for p>0. The well allows p to oscillates around 0. Figure 1 also shows that the amplitude of the oscillations of p is bounded such that |p|<1-c/U.

2.4 The linear limit

When |p| is small, the wave becomes linear. This happens when $\vert p\vert \ll 1-c/U$. In this limit the potential is parabolic, being given by

$$\tilde{V}(p) \sim\frac{p^2}{2} \left( 1-(c/U)^2 \right) .$$ (19)

The system then behaves as an harmonic oscillator. In this linear limit, Eq. (14) becomes

$$(U^2-c^2)\frac{{\rm d}^2\xi}{{\rm d}\Phi^2}+\Omega^2\xi=0 .$$ (20)

The solution is

$$\xi=\xi_0 \cos (kX_0-\omega t) + \xi_1 \sin (kX_0-\omega t) ,$$ (21)

with $k = \Omega / \sqrt{U^2-c^2},$and $\omega = kU$ which together imply the standard dispersion relation for linear sound waves modified by rotation in the form $\omega^2=k^2 c^2 + \Omega^2.$

	Figure 1: Profile of the function $\tilde{V}(p)$ for different values of (c/U)2. The solid, dotted, dashed, dotted-dashed and dotted-dotted-dotted-dashed curves respectively correspond to (c/U)2=0.005,0.02,0.1,0.25,0.5.
Open with DEXTER

2.5 The effect of vertical stratification

The same equations as above apply when the disk is vertically stratified. The nonlinear wave has exactly the same properties in that case. Both the perturbed and unperturbed density are simply multiplied by a factor $\exp \left(-\frac{z^2}{2H^2}\right)$which contains the entire vertical dependence.

   
2.6 Solutions of the wave equation

Figure 2: Left panel: radial profile of the density for a wave of velocity (c/U)2=0.3, for different wave amplitudes: $\rho _{\rm max}/\rho _0 =1.1$ ( dashed line), 1.4 ( dotted line) and 1.8 ( solid line). The solution shows a cusp around its maximum as the amplitude of the wave is increased. Right panel: same as the left panel, but for (c/U)2=0.9 and $\rho _{\rm max}/\rho _0=1.01$ ( dashed line), 1.03 ( dotted line) and 1.052. Note the smaller wave amplitude and the smaller periods in that case compared to the case (c/U)2=0.3.

Open with DEXTER

Using Eq. (15), the wave Eq. (14) can be written as a system of two coupled first order differential equations:

  

$$\frac{{\rm d}\xi}{{\rm d}\Phi} = p .$$ (22)

$$\frac{{\rm d}p}{{\rm d}\Phi} = - \left( \frac{\Omega}{c} \right)^2 \frac{(c/U)^2 \xi}{1-\frac{(c/U)^2}{ ( 1+p )^2}}\cdot$$ (23)

Once the phase $\Phi$ and displacement $\xi$ are normalized by the natural shearing box length scale $\Omega/c,$ a solution of the wave equation is specified by the parameter c/U and the maximum value of the density in units of $\rho_0,$ $\rho_{\rm max}/\rho_0.$ The latter fixes the initial conditions for the integration, which is performed until the solution reaches the next maximum. This procedure gives the spatial period T_w of the wave, or, equivalently, its wavenumber $k_{x,w}=2\pi/T_w.$

The solutions we obtained are illustrated in Fig. 2. The left panel shows the structure of the wave when (c/U)²=0.3 for different values of the relative wave amplitude $\rho _{\rm max}/\rho _0 =1.1$, 1.4 and 1.8 (we note that the maximum possible free wave amplitude is $\rho_{\rm max}/\rho_0 = U/c = 1.825$for these parameters). For these three values, we respectively found T_w=9.58H, 9.43H and 9.23H. As the amplitude of the wave is increased, nonlinearity becomes more and more important and a cusp eventually forms around the location of the density maximum. The right panel shows the structure of the wave when (c/U)²=0.9. Again, three different wave amplitudes illustrate the transition from linear to nonlinear wave profiles. These are given by $\rho _{\rm max}/\rho _0=1.01$,1.03 and 1.052. We note that similar waveforms occur for both values of c/U. However, for corresponding waveforms, when (c/U)²=0.9, the amplitudes and periods of the wave are much smaller than when (c/U)²=0.3.

   
2.7 Numerical tests

The results presented in Sect. 2 can be used to provide 1D tests for hydrodynamic codes. Here we present such test calculations performed with two finite difference codes, ZEUS (Hawley & Stone 1995) and NIRVANA (Ziegler & Yorke 1997), and a Godunov code, RAMSES (Teyssier 2002). The implementation of shearing box simulations using a Godunov code is described in detail in Gardiner & Stone (2005) and was subsequently applied using RAMSES (Fromang et al. 2006).

We consider only the case (c/U)²=0.3 here since this is the case we will study in more detail in Sect. 4. We perform different 1D numerical simulations of wave propagation for different wave amplitudes. In each case, the size of the box is set equal to the wavelength of the wave. We used periodic boundary conditions in shearing coordinates and followed the wave over many cycles.

The point of these tests is to demonstrate that we can successfully follow free propagating non linear waves. In particular we need to verify that this can be done under the conditions adopted in Sect. 4 where we study the hydrodynamic stability of the waves. For that reason, we present results obtained with a single resolution having 320 grid cells in the radial direction. In these simulations, the timestep $\Delta t$ is defined by

$$\Delta t=C_0\frac{\Delta x}{\textrm{max}\vert v_x\vert+c}$$ (24)

where $\Delta x$ stands for a grid cell size and C₀ is the Courant number, which should be positive and smaller than one for the CFL condition to be satisfied (Press et al. 1986).

2.7.1 Results

We considered waves with two amplitudes. These are given by $\rho _{\rm max}/\rho _0=1.05$ and $\rho _{\rm max}/\rho _0=1.8$. The results obtained with ZEUS and RAMSES are illustrated in Fig. 3. Results obtained using NIRVANA are very similar to the former. Each of the panels in Fig. 3 shows the radial profile of the density at different times measured in wave periods. The upper panels illustrate results obtained for the low amplitude case. The lower panels illustrate results obtained for the high amplitude case. In the former case, these are shown for t=0, t=250 and t=500. In the latter case, results are shown for t=0, t=30 and t=60. The left panels show results obtained with ZEUS with C₀=0.5, the central panels show results obtained with ZEUS using C₀=0.1 and the right panels show results obtained with RAMSES with C₀=0.7.

Figure 3: Results of the numerical tests illustrating the propagation of 1D nonlinear waves in the shearing box. The parameters are such that (c/U)2=0.3 and $\rho _{\rm max}/\rho _0=1.05$ ( upper row) and $\rho _{\rm max}/\rho _0=1.8$ ( lower row). The resolution is with 320 grid points in the radial direction. Upper and lower left panels show results obtained with ZEUS using a Courant number C0=0.5, the upper and lower central panels are obtained with ZEUS using C0=0.1 The upper and lower right panels are obtained with RAMSES using C0=0.7. For the upper panels, the plots are for times (measured in cycles or wave periods, with one cycle or period being the time for the wave to cross the box) t=0 ( solid line), t=250 ( dotted line) and t=500 ( dashed line). For the lower panels, the plots are for times t=0 ( solid line), t=30 ( dotted line) and t=60 ( dashed line). Note that the profiles at t=0 are computed by solving the wave equation as described in Sect. 14.

Open with DEXTER

2.7.2 Effect of Courant number

Firstly, it is apparent that the ZEUS results obtained with C₀=0.5 are never acceptable for long simulations. When the wave amplitude is small, we observed that its profile is dispersed on timescales of hundreds of periods. When the wave amplitude is large, small scale oscillations indicative of numerical instability, appear on timescales of tens of periods.

As indicated in the central panels, both problems are cured by decreasing the Courant number. The wave drift relative to the centre of the box seen at large times in both cases is due to the modification of the wave period as its amplitude is slowly damped by numerical dissipation.

Finally, as shown in the right hand side panels, the results obtained with RAMSES using C₀=0.7 are very stable and less dissipative than the ZEUS results. In conclusion, despite differences in the codes, we found that both Godunov and finite difference codes are able to follow the evolution of these waves for a long enough timescale to perform the simulations we will present in Sect. 4 provided the Courant number is chosen appropriately. For ZEUS or NIRVANA this means that $C_0 \leq 0.1.$

   
3 Stability of free waves

Having found exact nonlinear solutions corresponding to freely propagating density waves, an important issue is their stability. A possibility is that such waves undergo parametric instability through interaction with lower frequency inertial or gravity waves (Goodman 1993; Papaloizou 2005a). In addition to the production of these additional waves, such instabilities may result in the production of turbulence and the decay of the initial wave amplitude and so affect its ability to propagate over large distances.

Accordingly we now study the stability of non linear density waves to parametric decay into inertial waves (gravity waves being absent in our model). This can be done while maintaining axisymmetry, ie. we assume dependence of the state variables on x and z only.

3.1 Stability analysis

To study the stability of the nonlinear waves described above we start from the equations of motion in Lagrangian form. Following the discussion of Sect. 2, when the vertical component of the gravitational acceleration, and thus vertical stratification is neglected, these are seen to take the form

$$\frac{{\rm D}^2 X}{{\rm D} t^2} + \Omega^2 (X-X_0) = -\frac{1}{\rho} \frac{\partial P}{\partial x}$$ (25)

$$\frac{{\rm D}^2 Z}{{\rm D} t^2} = -\frac{1}{\rho} \frac{\partial P}{\partial z} \cdot$$ (26)

We now want to study the evolution of a small perturbation of the state defined by the unperturbed wave. To do so, we write

$$X = X_0+\xi_x + \tilde{\xi}_x$$ (27)

$$Z = Z_0 + \tilde{\xi}_z ,$$ (28)

where $\mbox{\boldmath {$\tilde \xi$ }} = (\tilde{\xi}_x , \tilde{\xi}_z)$ has a small magnitude in comparison to that of $\xi_x \equiv \xi$ the displacement associated with the wave used above. We therefore have $\Delta X=\tilde{\xi}_x$ and $\Delta Z=\tilde{\xi}_z$, where $\Delta$ stands for the Lagrangian perturbation. Upon noting that $\Delta$ and ${\rm D}/{\rm D}t$ commute (Lynden-Bell & Ostriker 1967), we obtain

  

$$\frac{{\rm D}^2 \tilde{\xi}_x}{{\rm D} t^2} + \Omega^2 \tilde{\xi... ...rtial P}{\partial x}\right) = -\Delta\left(\frac{\partial h}{\partial x}\right)$$ (29)

$$\frac{{\rm D}^2 \tilde{\xi}_z}{{\rm D} t^2} = -\Delta\left(\frac{... ...ial P}{\partial z}\right) = -\Delta\left( \frac{\partial h}{\partial z}\right),$$ (30)

where $h=c^2 \ln \rho$ is the enthalpy. Writing the Eulerian variation of a quantity f by f', we have the relation $\Delta f= f' + \mbox{\boldmath {$\tilde{\xi}$ }} . \mbox{\boldmath {$\nabla$ }} f$ (Lynden-Bell & Ostriker 1967), from which we obtain

$$\Delta\left(\frac{\partial h}{\partial x}\right) = \frac{\partial h'}{\partial x}+\tilde{\xi}_x \frac{\partial^2 h}{\partial x^2}$$ (31)

$$\Delta\left(\frac{\partial h}{\partial z}\right) = \frac{\partial h'}{\partial z}\cdot$$ (32)

Using the continuity equation $\rho' + \mbox{\boldmath {$\nabla$ }} (\rho \mbox{\boldmath {$\tilde{\xi}$ }} )=0$, we write

$$h'=c^2 \frac{\rho'}{\rho}=-\frac{c^2}{\rho}\frac{\partial \rh... ...^2}{\rho}\frac{\partial \rho \tilde{\xi}_z}{\partial z} \cdot$$ (33)

Combining Eqs. (29), (30) and (33), we get

  

$$\left(\frac{\partial}{\partial t}+v_x\frac{\partial}{\partial x}\... ...rho} \frac{\partial \rho}{\partial x} \frac{\partial \tilde{\xi}_x}{\partial x}$$ (34)

$$\left(\frac{\partial}{\partial t}+v_x\frac{\partial}{\partial x}\... ...rho} \frac{\partial \rho}{\partial x} \frac{\partial \tilde{\xi}_x}{\partial z}$$ (35)

where we used Eq. (7) that defines the convective derivative. Adopting a reference frame that moves with the speed U of the background wave (discussed in Sects. 14-2.6), we also have $\rho=\rho_0/(1+p)$ and v_x=-U(1+p) where $p=\partial \xi_x/\partial x$. These, along with the relation $v_y=-\frac{3}{2}\Omega X - \frac{1}{2}\Omega \xi$, fully determine the system. Finally, note that dissipation is neglected in our analysis. In Appendix A, we describe briefly how to incorporate it and the changes it produces. However, in the following, we assume the evolution is inviscid as it simplifies the analysis.

   
3.2 The case p = 0

In the absence of a background wave, p=0 and the previous system reduces to the simpler form

  

$$\left(\frac{\partial}{\partial t}-U\frac{\partial}{\partial x}\ri... ...}_x}{\partial x^2} + c^2 \frac{\partial^2 \tilde{\xi}_z}{\partial x \partial z}$$ (36)

$$\left(\frac{\partial}{\partial t}-U\frac{\partial}{\partial x}\ri... ...\partial z^2} + c^2 \frac{\partial^2 \tilde{\xi}_x}{\partial x \partial z}\cdot$$ (37)

In this case, the radial displacement of a fluid element can be taken to have the form

$$\mbox{\boldmath {$\tilde{\xi}$ }} \equiv \mbox{\boldmath {$\t... ...\propto {\rm e}^{{\rm i}(\omega_{I,n} t +n k_{x,w} x -k_z z)}$$ (38)

where n, as required by the radial periodicity condition, is an integer. Then Eqs. (36) and (37) lead to the conditions

  

$$\left[(\omega_{I,n} - n\omega)^2-\Omega^2 -n^2k_{x,w}^2c^2\right]\tilde{\xi}_x^n + n k_{x,w}k_z c^2\tilde{\xi}_z^n = 0$$ (39)

$$n k_{x,w} k_z c^2 \tilde{\xi}_x^n +\left[(\omega_{I,n}-n \omega)^2-k_z^2 c^2\right]\tilde{\xi}_z^n = 0.$$ (40)

From these a dispersion relation can be derived in the form

$$[\omega_{I,n}-n \omega]^4-\left((n^2 k_{x,w}^2+k_z^2)c^2+\Omega^2\right)[\omega_{I,n}-n \omega]^2+\Omega^2 k_z^2 c^2=0,$$ (41)

where, here, $\omega= k_{x,\omega}$. The low frequency branch of this dispersion relation corresponds to inertial waves for which the frequency is given by

$$(\omega_{I,n}-n \omega)^2=\omega_n^2=\frac{k^2c^2+\Omega^2}{2... ...ac{4k_z^2\Omega^2c^2}{(k^2c^2+\Omega^2)^2} \right]} \right) ,$$ (42)

where k²=n² k_x,w²+k_z². In the large wavenumber limit, this reduces to $\omega_{I,n}-n \omega= \pm \omega_n = \pm \frac{k_z}{k} \Omega$. The frequency of these inertial waves is thus smaller than $\Omega$ in the original non-comoving frame.

   
3.3 The case p $\ne$ 0

When $p \ne 0$, but is in the small wave amplitude limit (i.e. $\vert p\vert \ll 1$), terms of higher than first order in p may be neglected, so that the system of equations becomes

  

$$\left[\frac{\partial}{\partial t}-U(1+p)\frac{\partial}{\partial x}\right]^2 \tilde{\xi}_x + \Omega^2 \tilde{\xi}_x c^2 \frac{\partial^2 \tilde{\xi}_x}{\partial x^2} + c^2 \frac{\pa... ...} - c^2 \frac{\partial p}{\partial x} \frac{\partial \tilde{\xi}_x}{\partial x}$$ = (43)

$$\left[\frac{\partial}{\partial t}-U(1+p)\frac{\partial}{\partial x}\right]^2\tilde{\xi}_z \hspace{1.2cm} c^2 \frac{\partial^2 \tilde{\xi}_z}{\partial z^2} + c^2 \frac{\pa... ...^2 \frac{\partial p}{\partial x} \frac{\partial \tilde{\xi}_x}{\partial z}\cdot$$ = (44)

In this limit, the background wave is linear and p can be written $p=\epsilon \cos (k_{x,w}x),$ with $\epsilon$ being a small parameter. If the frequency of the background wave, $\omega$, lies between $\Omega$ and $2 \Omega$, then it is possible for two inertial waves with wavenumbers n k_x,w and (n+1)k_x,w to interact with each other by coupling through the original wave. Working in the frame comoving with the wave, noting that the time dependence is separable, we write the perturbed eigenvector as a linear combination of inertial waves in the form

$$\mbox{\boldmath {${\tilde \xi}$ }} = ( \mbox{\boldmath {$\til... ...\delta \sigma_{n+1}) t + {\rm i}(n+1) k_{x,w}x - {\rm i}k_z z}$$ (45)

where the wave frequencies are respectively shifted by $\delta \sigma_n$ and $\delta \sigma_{n+1}$ which are both taken to be of order $\epsilon.$ For the time dependence to be separable, so allowing the modes to interact, we require that

$$\omega_{I,n}+\delta \sigma_n=\omega_{I,n+1}+\delta \sigma_{n+1}.$$ (46)

Using Eq. (42), if $\omega_{I,n}=\omega_n +n \omega$and $\omega_{I,n+1}=-\omega_{n+1} +(n+1) \omega$, this is equivalent to

$$\omega_n+\omega_{n+1}=\omega+\delta \sigma_{n+1}-\delta \sigma_n=\omega+\delta \omega ~ .$$ (47)

Here $\omega_n, \omega_{n+1}, \omega$ and thus $\delta \omega$ are given, while (47) gives an interaction condition involving $\delta \sigma_n$ and $\delta \sigma_{n+1}.$A further condition on the latter two, enabling the completion of the solution, is obtained as follows:

Four linear equations relating the total of eight components of $\mbox{\boldmath {$\tilde \xi$ }} ^n , \mbox{\boldmath {$\tilde \xi$ }} ^{n+1}, \mbox{\boldmath {$\tilde \eta$ }} ^n$ and $\mbox{\boldmath {$\tilde \eta$ }} ^{n+1}$ can be obtained after substituting into Eqs. (43) and (44) and requiring that the resulting Fourier coefficients corresponding to the terms in Eq. (45) vanish. Setting $D_j = w_j^2-\Omega^2-j^2 k_{x,w}^2 c^2$these are

    

$$\epsilon \left[ D_n\tilde{\eta}_x^n + n k_{x,w} k_z c^2 \tilde{\e... ..._{n+1}-\frac{\omega^2}{2}-\frac{k_{x,w}^2c^2}{2}\right] \tilde{\xi}_x^{n+1} = 0$$ (48)

$$\epsilon \left[ D_{n+1}\tilde{\eta}_x^{n+1} + (n+1) k_{x,w} k_z c... ...a \omega_n-\frac{\omega^2}{2}-\frac{k_{x,w}^2c^2}{2}\right] \tilde{\xi}_x^n = 0$$ (49)

$$\epsilon \left[ (w_n^2-k_z^2 c^2)\tilde{\eta}_z^n + n k_{x,w} k_z... ...tilde{\xi}_z^{n+1} - \frac{\epsilon}{2} k_z k_{x,w} c^2 \tilde{\xi}_x^{n+1} = 0$$ (50)

$$\epsilon \left[ (w_{n+1}^2-k_z^2 c^2)\tilde{\eta}_z^{n+1} + (n+1)... ...right] \tilde{\xi}_z^n + \frac{\epsilon}{2} k_z k_{x,w} c^2 \tilde{\xi}_x^n = 0$$ (51)

Eqs. (39) and (40) and their counterparts for $n \rightarrow n+1$ give an additional four linear equations. Being linear, in order that the total of eight equations be soluble, their determinant must vanish. This leads to the requirement that

$$\delta \sigma_n \delta \sigma_{n+1} = - \frac{\epsilon^2}{4} ... ...{C_n C_{n+1}} \frac{1}{\omega_n \omega_{n+1}} = -\sigma^2 ~ ,$$ (52)

where

$$\left\{ \begin{array}{ccl} A_n &=& \left[(n+1)(k_z^2 c^2 -\omega_... ...2 + k_z^2 c^2 + (n+1)^2 k_{x,w}^2 c^2 - 2 \omega_{n+1}^2 ~. \end{array} \right.$$

Upon noting that $\delta \sigma_{n+1}=\delta \sigma_n+\delta \omega$, we find

$$\delta \sigma_n=\frac{\delta \omega}{2} \pm \sqrt{\frac{(\delta \omega)^2}{4}-\sigma^2}.$$ (53)

When $\delta \omega<2\sigma$, there is a non zero imaginary part to both $\delta \sigma_n$ and $\delta \sigma_{n+1}$ and the amplitude of the two coupled inertial modes grows: this is a parametric instability. The maximum growth rate $\sigma$ is obtained when $\delta \omega=0$, which occurs when the resonant condition satisfied by the inertial waves is:

$$\omega_n+\omega_{n+1}=\omega .$$ (54)

When $\delta \omega$ does not vanish, the growth rate of the instability is reduced, as also noted by Goodman (1993). We note that on account of symmetry with respect to reflection in the midplane, the growth rate does not change when $k_z \rightarrow -k_z$. This enables solutions to be classified according to their symmetry with respect to reflection in the midplane.

   
3.4 Unstable mode properties

Figure 4: Properties of the unstable modes of the parametric instability for a wave characterized by (c/U)2=0.3 and an amplitude such that the maximum density $\rho _{\rm max}/\rho _0=1.05$. The two panels respectively show the vertical wavelength, in units of the scale height, of the unstable mode and the growth rate in units of the wave period as a function of the radial wavenumber n. In both panels, the solid lines represent the large n limit given by Eqs. (55) and (56). Note that the mode vertical wavelengths become dense at large nsuch that the requirement of being equal to the length of a prescribed vertical domain such as the scale height when multiplied by some integer can be approached arbitrarily closely.

Open with DEXTER

In this section, we study the properties of modes that undergo parametric instability for the special case (c/U)²=0.3. We focus here on modes that occur when the resonant condition given by Eq. (54) is satisfied exactly.

For a given background wave amplitude, the spatial period of the background wave is calculated as described in Sect. 2.6. The radial wavenumber of the background wave k_x,w and the wave frequency $\omega$ are then fixed. The instability couples two inertial waves with radial wavenumbers nk_x,w and (n+1)k_x,w. These should satisfy the resonant condition, given by Eq. (54), and the dispersion relations, given by Eq. (41) and its counterpart with $n \rightarrow n+1.$For given values of n, $\omega$ and k_x,w; $\omega_n$, $\omega_{n+1}$ and k_z can be determined from these algebraic conditions. At first we shall ignore any additional constraints on k_z that may arise from boundary conditions such as the requirement that the vertical wavelength should be the scale height divided by an integer.

We note that in the large wavenumber limit ($n \gg 1$) the solution is given by

$$\omega_n=\omega_{n+1}=\frac{\omega}{2} ~ , ~ ~ k_z=n k_{x,w}\sqrt{\frac{\omega^2}{4\Omega^2 - \omega^2}} \cdot$$ (55)

In the same limit, Eq. (53) for the growth rate gives

$$\sigma \sim \frac{\epsilon}{8} \Omega \left(\frac{c}{U}\right... ...{x,w} H)^2}\frac{4\Omega^2 -\omega^2}{4\Omega^2}\right) \cdot$$ (56)

To illustrate the properties of the unstable modes, in Fig. 4 we plot the growth rate $\tilde{\sigma}=\sigma (2\pi/\omega)$ and vertical wavelength as a function of n for a wave whose amplitude is such that $\rho _{\rm max}/\rho _0=1.05$. It is seen that results obtained assuming the large n limit applies are always in good agreement with the actual results. We comment that in the limit of large n the modes become dense in the sense that the requirement that the product of the vertical wave length and some integer be the extent of a vertical domain or the scale height can be satisfied with arbitrary accuracy. This follows from Eq. (55) and the fact that the rationals are dense.

Finally, the structure of the eigenmode is illustrated in a particular case in Fig. 5. This shows a snapshot of the vertical velocity in the (x,z) plane for a mode with n=4, in the case of a background wave for which (c/U)²=0.3 and $\rho _{\rm max}/\rho _0=1.05$. The structure of this mode will be compared in the following section with the results of numerical simulations.

   
4 Numerical simulations

In this section, we present the results of numerical simulations performed with ZEUS, NIRVANA and RAMSES. The goals are first to compare the properties of the instability during the linear phase of the parametric instability with the analytical theory presented above. Then we investigate the stability of strongly nonlinear waves. We also study the evolution of the instability during its nonlinear phase. In the next section we focus on waves for which (c/U)²=0.3.

   
4.1 Small amplitude waves

	Figure 5: Snapshot of the nondimensional vertical velocity in the (x,z) plane showing the structure of the unstable mode of the parametric instability when n=4. The parameters of the background wave are (c/U)2=0.3 and $\rho _{\rm max}/\rho _0=1.05$.
Open with DEXTER

We first study the case of a background wave with a small amplitude such that $\rho _{\rm max}/\rho _0=1.05$. The properties of the unstable modes expected in that case were derived in Sect. 3.4. From the solution of Eq. (14) we find that T_w=9.6H. The analytical results (see Fig. 4) show that the most unstable modes have the largest wavenumbers. However, because of their small scale structure, they will be difficult to represent numerically. In order to illustrate the properties of the parametric instability, we isolate the mode for which n=4. This is achieved by adjusting the extent of the vertical domain appropriately. The results of Sect. 3.4 indicate that $\lambda_z=2.77H$ and $\tilde{\sigma}=0.06$ for that mode.

To follow its evolution, we performed 2D numerical simulations in the (x,z) plane with ZEUS, NIRVANA and RAMSES. For the n=4 mode to develop, we chose a computational box that covers the range [-T_w/2,T_w/2] in x and $[-\lambda_z/2,\lambda_z/2]$ in z (so that the radial extent matches the wavelength of the background wave while the vertical extent matches the vertical wavelength of the unstable mode). We used periodic boundary conditions in both x and z. Of course, other modes, compatible with the boundary conditions (i.e. those which have a ratio of box size to vertical wavelength that is a large integer) could in principle grow as well. However, very small scales modes will be strongly damped by numerical diffusion (see Appendix A). In practise, we found that they did not develop at the resolutions we investigated. We also found that an accurate description of the n=4mode requires roughly 64 cells per wavelength. This translates into a resolution (N_x,N_z)=(320,64). For smaller resolutions, the growth rate is reduced, although the evolution is qualitatively similar.

At time t=0, the density and velocity components are set equal to the solutions of the system of differential equations given by Eqs. (22) and (23). A perturbation is also added to each component of the velocity. If this perturbation was random, then a combination of modes with vertical wavenumbers +k_z and -k_z would grow with relative amplitudes depending on the initial conditions. This would make the comparison with the linear theory less illustrative. Rather, we chose to add a perturbation that is close to the mode with wavenumber +k_z by adding a velocity perturbation of magnitude $\delta v$ such that

$$\delta v=\delta v_0 [-\sin(n k_{x,w}x-k_z z)+\sin((n+1)k_{x,w} x-k_z z)],$$ (57)

where we took $\delta v_0=10^{-4} c$. In practise, we found this procedure prevented the mode with wavenumber -k_z from growing during the course of our simulations.

The time history of the maximum value of the vertical velocity is plotted in the left panel of Fig. 6. The solid line shows the results obtained with ZEUS, while the dot-dashed and dashed line respectively correspond to the results obtained with NIRVANA and RAMSES. All curves indicate an exponential growth at early times. The evolution expected from linear theory is also plotted with the dotted line. There is very good agreement between this and the numerical results, which themselves are in very good agreement. Results from all three codes show that the linear instability saturates after about 100 wave periods, at which point the maximum vertical velocity has reached about $10\%$ of the sound speed. This is of the order of the radial velocity associated with the background wave. This indicates that the linear instability saturates when the amplitude of the velocities of the unstable inertial modes becomes large enough to perturb the structure of the background wave. During the nonlinear phase, the maximum value of v_z at particular time is first found to oscillate before it saturates (as obtained with ZEUS) or slowly decays (as obtained with RAMSES). The oscillation is due to the linear instability recurring quasi periodically between t=100 and t=400, before the structure of the underlying wave is sufficiently modified. Finally, we comment that the different evolution of the simulations in the different codes is due to their different dissipative properties. In order to illustrate the effect of the instability on the wave structure, we plot the time history of its amplitude on the right hand side panel of Fig. 6 (with the same conventions as used on the left hand side panel). Once again, there is a good agreement between the codes during the linear phase of the evolution and during the recurring instability phase. RAMSES displays a larger amplitude than ZEUS and NIRVANA during that time. An additional run performed with NIRVANA using twice the resolution gave similar results to those of RAMSES, suggesting that this difference is due to the larger numerical dissipation in the finite difference schemes. At late times, all codes find that the wave amplitude is reduced to smaller and smaller amplitudes, although there are differences due to the different small scale dissipative properties.

	Figure 6: The time history of the maximum value of the vertical velocity ( left panel) and of the background wave amplitude $(\rho _{\rm max}/\rho _0 -1)$ ( right panel) as obtained with ZEUS ( solid line), NIRVANA ( dot-dashed line) and RAMSES ( dashed line). The dotted line shown in the left panel indicates the behaviour expected from linear theory.
Open with DEXTER

Figure 7: Snapshots showing the structure of the parametric instability in the (x,z) plane in the case of a linear background wave whose parameters are (c/U)2=0.3 and $\rho _{\rm max}/\rho _0=1.05$. The first two panels, both obtained with ZEUS, respectively show the density (normalised by the mean density $\rho _0$) and vertical velocity distribution (normalised by the sound speed). The third and forth panels show snapshots of the vertical velocity obtained with NIRVANA and RAMSES respectively. Both are very similar to the ZEUS results. All the plots are obtained at t=50.

Open with DEXTER

Figure 7 gives plots of the structure of the eigenmode obtained with the three different codes. The top panel shows the density distribution, normalised by the mean density, in the (x,z) plane after 50wave periods, as obtained with ZEUS. Being in the linear phase, it is hardly modified compared to the initial profile. The next three panels, from top to bottom, represent the vertical velocity (normalised by the sound speed c) respectively with ZEUS, NIRVANA and RAMSES. They look almost identical to each other and extremely similar to the expected theoretical result shown in Fig. 5.

Taken together, all the results presented here demonstrate a very good agreement between the numerical and analytical results. In the next section we investigate the properties of the parametric instability that destabilizes waves of larger amplitude.

   
4.2 Larger amplitude waves

Figure 8: Growth rate of the n=4 mode of the parametric instability for different background wave amplitudes that are characterized by (c/U)2=0.3. The stars are determined by fitting the results of numerical simulations obtained with ZEUS. The circles are obtained using the results of the linear analysis (using Eq. (52), the parameter $\epsilon$ is taken to be the amplitude of the Fourier mode whose wavelength equals that of the wave).

Open with DEXTER

	Figure 9: Same as Fig. 7, but for a large amplitude background wave such that $\rho _{\rm max}/\rho _0=1.8$. The snapshots are obtained at time t=10. Note the distorted structure of the eigenmode close to the density maximum. It results from the strong nonlinearities of the background wave at this location.
Open with DEXTER

Using ZEUS, we followed the evolution of waves of different amplitude. We used the same set up (resolution, initial perturbation) as for the small amplitude wave described above. In each case, the size of the computational box is tuned to match the background wave period in the radial direction and the wavelength of the n=4 mode in the vertical direction. The results are illustrated in Fig. 8. We show the growth rates $\tilde{\sigma}$ of the instability measured during the linear phase for wave amplitudes such that $\rho _{\rm max}/\rho _0=1.05$, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7 and 1.8. They can be compared with the growth rates computed according to Eq. (52). To calculate these, the amplitude $\epsilon$ that appears in Eq. (52) is taken to be equal to the amplitude of the Fourier mode with wavelength equal to the radial width of the box. Figure 8 shows that there is very good agreement between the results of the numerical simulation and the analytical results for wave amplitudes smaller than $\rho_{\rm max}/\rho_0=1.5$. For larger amplitudes, the analytical estimates tend to underestimate the growth rate slightly because of the strongly nonlinear nature of the wave. This nonlinearity modifies the structure of the eigenmode, as shown in Fig. 9 for the case $\rho _{\rm max}/\rho _0=1.8$. This figure is similar to Fig. 7: the top panel shows the density distribution in the (x,z) plane and the bottom panel represents the vertical velocity. Both were obtained with ZEUS at time t=10, well within the linear phase of the instability. A strong cusp, already noted in Fig. 2, is seen almost undisturbed on the upper snapshot. The lower one shows the structure on the eigenmode. It is similar to the plots of Fig. 7 although its pattern is distorted close to the wave cusp. This is probably affecting the growth rate of the mode as indicated in Fig. 8.

	Figure 10: Same as Fig. 6, but for a large amplitude background wave such that $\rho _{\rm max}/\rho _0=1.8$ (only the ZEUS and RAMSES results are represented, respectively with the solid and dashed lines, while NIRVANA results are very similar). Note the faster growth of the instability compared to the smaller amplitude case.
Open with DEXTER

Figure 11: Structure of the flow during the nonlinear stage of the parametric instability. Left panels correspond to a background wave amplitude $\rho _{\rm max}/\rho _0=1.05$ and (c/U)2=0.3. The upper figure plots $\rho /\rho _0$ in the (r,z) plane, while the lower figure shows vz/c. Both are plotted at t=500. The right hand side panels shows the same quantities at t=27 for a background wave amplitude $\rho _{\rm max}/\rho _0$=1.8 and (c/U)2=0.3.

Open with DEXTER

The long term evolution of the instability in this strongly nonlinear regime is further illustrated by Fig. 10, which is the equivalent of Fig. 6. ZEUS results are plotted using the solid line while RAMSES results are plotted with the dashed line. Both are almost identical and results obtained using NIRVANA displays a very similar behaviour. The evolution of the instability is qualitatively similar to the small amplitude wave case, although the timescale for the instability to saturate is much shorter, being 20 wave periods when $\rho _{\rm max}/\rho _0=1.8$ as compared to 100 wave periods when $\rho _{\rm max}/\rho _0=1.05$ (once again, the linear instability saturates when the perturbed velocities reach that of the background wave, which happens to be close to the speed of sound in the large amplitude case). In addition, no recurring instability phase is found before the amplitude of the wave starts to decrease. This is better illustrated in Fig. 11. The left panels show the density structure (upper panels) and vertical velocity (lower panels) obtained in the small amplitude case ( $\rho _{\rm max}/\rho _0=1.05$) at time t=500. Even at those late times, although smaller scale disturbances are involved, the regular pattern characteristic of the parametric instability is still visible. This is in contrast with the snapshots obtained in the large amplitude case ( $\rho _{\rm max}/\rho _0=1.8$) at time t=27 and shown in the right panels of Fig. 11. The nature of the flow is far more disorganized and looks more like turbulence in that case. Finally, we want to comment that the system eventually goes back to equilibrium in all cases, with the wave amplitude being strongly damped. This might seem to disagree with Figs. 6 and 10, where the vertical velocity saturates and shows very little decrease. This is because a vertical flow of alternating direction sets in at large times due to the periodic boundary conditions in the vertical direction and damps only very slowly. v_z being only a function of x, this flow simply superposes to the equilibrium without disturbing it. We emphasize here that this is only an artifact of the numerical setup.

   
5 Stratified case

The results of the previous section give confidence that the nature of the instability is properly captured by the analytical work presented in Sect. 3. However, this analysis neglected the disk density vertical stratification, which is the focus of the present section.

5.1 Inertial wave properties

The properties of linear inertial waves in stratified isothermal disks have been studied by Lubow & Pringle (1993). We recall these below for the case $\gamma=1$ in their notation, as this is the case of interest here. Let v'_z(x,z) be the Eulerian perturbation of the velocity in the vertical direction. Because of vertical stratification, the z dependence in v'_z(x,z) is no longer harmonic as it was in the unstratified case. Instead, it can be shown to be of the form

$$v'_z(x,y,z)=\tilde{v_z}(z) \exp [{\rm i}(nk_{x,w}x-\omega_n t)] ,$$ (58)

where v_z^'(z) is a polynomial in z of the form:

$$\tilde{v_z}(z)=\sum_{m=0}^{n_z} c_m z^m .$$ (59)

Here n_z is given by

$$n_z+1=\frac{\omega_{n,n_z}^2}{\Omega^2} \left( 1- \frac{n^2 k_{x,w}^2 c^2}{\omega_{n,n_z}^2-\Omega^2} \right) ,$$ (60)

where the frequency of inertial waves, which constitute the low frequency branch of this relation, is now denoted by $\omega_{n,n_z}$since it depends on the two integers n and n_z. In addition, the coefficients c_m satisfy a recurrence relation, given by Eq. (49) of Lubow & Pringle (1993). In the large n_z limit, $\tilde{v_z}(z)$ oscillates with a period

$$\lambda_{z,n_z} \sim \frac{2\pi H}{\sqrt{n_z}}$$ (61)

in the vertical direction in the vicinity of the equatorial plane. Because of separability, two inertial waves with the same n_z but different values of n have the same vertical structure. They can therefore interact together in the presence of the background nonlinear sound wave when they satisfy the same resonant condition as in the unstratified case:

$$\omega_{n,n_z}+\omega_{n+1,nz}=\omega .$$ (62)

As for the unstratified case, the spectral frequencies $\omega_{n,n_z}$ are dense in the interval $[0,\Omega]$(this can be demonstrated by showing that the difference $\omega_{n+1,n_z}-\omega_{n,n_z}$ tends to zero as n_z tends to infinity): there will always be an infinite number of values of n_z such that Eq. (62) is satisfied to a specified accuracy. But note that the modes are in fact denser in the stratified case than in the unstratified case. This is because n_z appears linearly in Eq. (62) whereas the corresponding expression in the unstratified case would contain n_z²(see Eq. (55)). This has the consequence that the separation between neighbouring modes is less in the former case. A physical explanation for this is that this results from there being a continuous range of density scale heights in the stratified disk model compared to only one length scale in the unstratified case provided by the size of the computational box. This makes parametric instability more readily realised in the stratified case.

In a numerical simulation, however, only those modes for which the wavelength is big enough (or n small enough) to be well resolved will have a non vanishing growth rate. This limits the number of realisable modes to those that satisfy the resonant condition to within a given accuracy and are such that $n \leq n_{\rm max}$.

5.2 Numerical simulations

Figure 12: Time history of the vertical velocity in the stratified model, for the small amplitude wave $\rho _{\rm max}/\rho _0=1.05$ ( left panel) and for the large amplitude wave $\rho _{\rm max}/\rho _0=1.8$ ( right panel). In both cases, the dotted line represents the exponential fit obtained during the linear phase of the instability. The growth rates that result are respectively $\tilde{\sigma}=0.038$ and 0.32.

Open with DEXTER

The simulations we describe in this section are the stratified analogs of the simulations presented in Sect. 4 (we only report results obtained with ZEUS here as we have shown good agreement between the different numerical methods above). The wave velocity is identical, (c/U)²=0.3, and we consider only two wave amplitudes defined such that $\rho _{\rm max}/\rho _0=1.05$ and $\rho _{\rm max}/\rho _0=1.8$. In our standard runs, the resolution is kept the same in the xdirection, N_x=320, as is the extent of the computational domain, $x \in [-T_w/2,T_w/2]$. In the vertical direction, the computational domain is extended to cover 4 scale heights on both sides of the equatorial plane and the resolution is consequently increased to 193 grid points. The initial density is set to be in hydrostatic equilibrium in the vertical direction (this gives a Gaussian vertical profile for the density) while the profile of the variables in the radial direction is unchanged compared to the unstratified case. Finally, we use periodic boundary conditions in the radial direction and outflow boundary conditions in the vertical directions. At t=0, all components of the velocity are randomly perturbed with an amplitude equal to $0.01 \%$ of the sound speed. Note that this random initial condition does not break the symmetry of the problem at t=0. Therefore, we expect modes with wavenumbers k_z and -k_z to develop at the same time.

As in the unstratified simulations, the parametric instability is seen to develop in both runs. This is illustrated in Fig. 12 which shows the time history of the maximum value of the vertical velocity in both cases (as before, time is measured in units of the wave period). The left panel corresponds to the small amplitude wave, $\rho _{\rm max}/\rho _0=1.05$, and the right panel to the large amplitude wave $\rho _{\rm max}/\rho _0=1.8$. Similarly to the unstratified case, the growth of the instability is faster in the latter case. The growth rates, measured during the linear phase and shown in Fig. 12 using a dotted line, are respectively $\tilde{\sigma}=0.038$ and $\tilde{\sigma}=0.32$. Both are comparable (although slightly smaller) than the growth rates obtained in Sect. 4.2 and plotted in Fig. 8. We also checked that the amplitude of the wave starts to decrease as the instability enters the nonlinear regime and that the system goes back to hydrostatic equilibrium. The nonzero velocities seen at large times in Fig. 12 are due to motions in the disk upper layers which are possibly long wavelengths acoustic modes excited by the parametric instability. Because of the low inertia of the gas, they take very long to damp. The structure of the flow at various stages of the simulation is illustrated in both simulations by showing snapshots of different quantities in the (x,z) plane in Fig. 13. The left hand side snapshots correspond to the small amplitude wave case and the right hand side snapshots to the large amplitude wave case (note that both plots only cover the domain [-2.5H,2.5H] in the vertical direction, ie only a fraction of the actual computational box). On both sides, the upper panel shows the distribution of the density during the linear phase, respectively at times t=150 and t=17. As in Fig. 10, the cusp in the density is apparent in the latter case. The middle panels in Fig. 13 represent the vertical velocity at the same times as the upper panels. The stripped structure of the flow characteristic of the parametric instability is obvious in both simulations. In the small amplitude wave case, the structure of the unstable eigenmode is regular while it is distorted in the neighbourhood of the cusp in the large amplitude wave case. This is again similar to the results obtained in the unstratified simulations. Also similar is the outcome of the instability during the nonlinear regime. This is illustrated by the lower panels of Fig. 13 which show the structure of the vertical velocity in both cases, respectively at times t=450 and t=45 (i.e. well into the nonlinear regime). The small amplitude wave simulation still displays a very regular structure while the entire flow broke down into turbulence in the large amplitude wave case. In conclusion, all of the results obtained in these stratified simulations are very similar to the unstratified results of Sect. 4.

Figure 13: Snapshots of the density and velocity in the (x,z) plane showing the structure of the parametric instability in a stratified disk. The panels on the left correspond to the small amplitude wave, $\rho _{\rm max}/\rho _0=1.05$, while the right hand side panels refer to the large amplitude wave $\rho _{\rm max}/\rho _0=1.8$. On both sides, the top panels show the density (normalised by the mean density) during the linear phase (respectively at t=150 and t=17), the middle panels show the vertical velocity (normalised by the sound speed) at the same times and the bottom panels show the vertical velocity during the nonlinear phase (respectively at t=450 and t=45).

Open with DEXTER

	Figure 14: Vertical velocity in the (x,z) plane resulting from the n=9 mode only, calculated according to the linear results of Sect. 3 in the case of the small amplitude wave $\rho _{\rm max}/\rho _0=1.05$. It is seen to be in good agreement with the numerical simulations (see the snapshots given by Fig. 13).
Open with DEXTER

One may wonder, however, what determines the particular mode(s) that grow in these stratified simulations. Unlike the unstratified case, no vertical scale is determined by the size of the computational box. To get some insight into this question, we now compare the results of the small amplitude wave simulation to the linear analysis (performed in the unstratified box) presented in Sect. 3.4. We expect the latter to agree with the numerical results when large wavenumber modes are excited. Indeed, the instability is truly local in that case and should not be affected by the vertical density stratification. In Fig. 13, the dominant mode that emerges in the simulations seems to be such that $n \sim 8$ to 10. As an example, the structure of the flow resulting from the mode n=9 is illustrated in Fig. 14 (note that k_z is determined as the value that leads to the satisfaction of the parametric resonance condition and both modes with wavenumber +k_z and -k_z are included). It is roughly in agreement with the left hand side snapshots of Fig. 13, although some signs of the vertical stratification are seen in the latter through an increase of the velocity at large z. Figure 13 also displays a signature of large scale variations in the radial direction, possibly indicating that a smaller n mode may be growing at the same time (see below). The linear analysis indicates that the n=9 mode should have a wavelength $\lambda_z=1.34H$, in rough agreement with the simulation. However, it predicts a growth rate $\sigma_{\rm lin}=0.062$, larger than the numerically determined value $\sigma_{\rm num}=0.038$. This discrepancy is due partly to numerical damping reducing the growth rate. It is also due to the resonant condition being actually only marginally satisfied, which reduces the growth rate (Goodman 1993). We quantify that accuracy using the parameter $\theta$ which we define through

$$\theta=\frac{\omega_{n,n_z}+\omega_{n+1,n_z}}{\omega}-1 .$$ (63)

In the notation of Sect. 3.3, $\theta=\delta \omega/\omega$. Equation (53) states that the growth rate is unchanged when $\theta \ll \tilde{\sigma}/\pi$, while the instability is suppressed when $\theta \ge \tilde{\sigma}/\pi$. According to the linear results of Sect. 3.4, $\tilde{\sigma}/\pi \sim0.01$ for the small amplitude wave parameters we are studying here. We therefore expect that only those modes that have $\theta \le 0.01$ will grow with a significant growth rate. Using Eq. (60), we searched for all the modes having $n \le 10$ and $\theta \le 0.01$. The modes satisfying these criteria are listed in Table 1. The mode with the smallest value of $\theta$ has (n,n_z)=(9,21) and probably corresponds to the mode described above and seen in the simulations. In fact, Eq. (61) gives $\lambda_{z,n_z}=1.37H$ for that mode, which is in good agreement both with the results of the linear analysis mentioned above and to the results of the numerical simulations. All the other modes, with $n \in [7,10]$ and increasing values of $\theta$ probably contribute in a complicated way to the growth of the parametric instability. Table 1 also reports the existence of a mode with (n,n_z)=(4,4) and $\theta=8.9 \times 10^{-3}$ which probably accounts for the large scale structure seen in Fig. 13. As the resolution increases, we expect more modes with smaller wavelength to become important and to contribute to the structure of the flow. This is also expected as a consequence of numerical dissipation acting on the smallest scales (see appendix A). To illustrate these points, we performed a simulation with the resolution decreased by a factor of two: (N_x,N_z)=(160,93). Everything else is kept identical. The structure of the flow in that case is illustrated on the left panel of Fig. 15. This is again a snapshot of the vertical velocity in the (x,z) plane at time t=200. In agreement with the previous discussion, the growth of the parametric instability is dominated by an eigenmode whose wavelength is larger than the high resolution run. In fact, its structure probably results from the growth of the n=4 mode listed in Table 1. To check that hypothesis, the flow is compared in the right panel with the structure of the n=4 mode as computed according to the unstratified linear analysis. The agreement between the analytical results and the numerical simulations is very good.

Table 1: Modes ranked according to the precision $\theta$ (see its definition in the text) within which they satify the resonant condition given by Eq. (62). For each mode, labelled according to the first column, $\theta$ is given in the last column, while the second and third columns respectively give the values of n and n_z. The parameters of the model are such that (c/U)²=0.3 and the background wave has an amplitude $\rho /\rho _{\rm max}=1.05$.

Figure 15: Snapshot of the vertical velocity at time t=200 in the low resolution run: (Nx,Nz)=(160,93) ( left panel). The color scale spans the interval [-0.15c,0.15c]. Note the larger scale of the unstable eigenmode that grows in that case as compared to the standard run with twice that resolution. It is compared on the right panel with the structure of the n=4 eigenmode as calculated according to the results of Sect. 3.4 (recall that both n and n+1 contribute). Good agreement is found between the numerical and analytical results

Open with DEXTER

   
6 Conclusion

In this paper, we have derived and solved a simple ODE which describes the structure of nonlinear density waves propagating in accretion disks. Such an analysis is important for understanding the properties, stability and dissipation of waves that propagate in disks in many situations: tidal excitation by a companion, waves generated by the disk self-gravity when it is massive enough or even waves generated by MHD turbulence resulting from the MRI. For the analysis to remain tractable, we adopted in the present paper the shearing box approximation, restricted the equation of state to be isothermal and the wave to have no variation in the direction of the mean shear flow i.e. to be axisymmetric. Because the waves were free and stationary in a frame moving with the phase velocity, shocks that would require a continual forcing or energy input were not considered here.

Using the solution of the ODE that describes the propagation of waves of arbitrary amplitudes to provide initial conditions, we followed their subsequent evolution in 1D using different numerical methods. A Godunov code, RAMSES, was found to give very accurate results. Finite difference codes such as ZEUS or NIRVANA were found to diverge from the correct solution or develop small scale oscillations on long timescales. These problems were cured in our case by decreasing the Courant number to C=0.1. However, we comment that reducing the timestep only delays the onset of the problems which were always found to develop provided the waves were evolved for long enough.

We also studied the stability of the waves. When their frequency is in the range $[\Omega,2\Omega]$, they can interact resonantly with two linear inertial waves and undergo a parametric instability. In this context we note that this frequency range may be increased in stratified disks if $\gamma \ne 1$ for perturbations allowing the existence of g modes (Lubow & Pringle 1993). A wider class of waves would thus become candidates for parametric instability. Neglecting vertical density stratification, we derived an expression for the growth rate when the amplitude of the background wave is small. This was found to be in good agreement with 2D numerical simulations using RAMSES and ZEUS or NIRVANA (provided that in the latter cases, the timestep was tuned as described above). Both unstratified and stratified simulations were carried out. We found that the precise structure of the mode that grows depends on the resolution we used. This is because we used no physical dissipation and rather relied on the resolution dependent numerical viscosity to damp small scale modes. Using these simulations, we could study the nonlinear evolution of the instability. For small amplitude waves the motion remained ordered but large amplitude waves broke down to a turbulent state as might be expected from consideration of the related elliptical instability in rotating flows. In all cases, we found that the amplitude of the wave was damped as the system relaxed back to a stationary state. Of course, this is because no external energy that could keep the instability going is provided. This differs from realistic situations that can occur, such as when the disk is tidally forced by a companion. Then the background wave will be constantly forced. In such cases, the external forcing will keep driving the instability which we expect to result in a complicated quasi-steady state. We will study forced waves, which also allow the incorporation of shocks, in a future paper.

Acknowledgements

S.F. acknowledges the important contribution of Romain Teyssier to the adaption of the numerical scheme of RAMSES to the shearing box equations.

   
Appendix A: The effect of a physical dissipation

Here we consider the effect of a small constant kinematic viscosity $\nu$ on the linear growth of the parametric instability. To do so, we study the properties of unperturbed inertial waves (see Sect. 3.2) when $\nu$ is nonzero. Note however that we neglect the effect of the dissipation on the background wave because its scale is much larger. In the large wavenumber limit, we may neglect the effect of possible stratification, the inertial modes are nearly incompressible and Eqs. (36) and (37) can be written as

  

$$\left(\frac{\partial}{\partial t}-U\frac{\partial}{\partial x}\right)^2 \tilde{\xi}_x + \Omega^2 \tilde{\xi}_x c^2 \frac{\partial^2 \tilde{\xi}_x}{\partial x^2} + c^2 \frac{\pa... ...(\frac{\partial}{\partial t}-U\frac{\partial}{\partial x}\right) \tilde{\xi}_x,$$ = (A.1)

$$\left(\frac{\partial}{\partial t}-U\frac{\partial}{\partial x}\right)^2\tilde{\xi}_z \hspace{1.2cm} c^2 \frac{\partial^2 \tilde{\xi}_z}{\partial z^2} + c^2 \frac{\pa... ...(\frac{\partial}{\partial t}-U\frac{\partial}{\partial x}\right) \tilde{\xi}_x,$$ = (A.2)

when viscosity is taken into account. Using Eq. (38), these two equations lead to:

  

$$\left[\omega_n^2-n^2k_{x,w}^2c^2-{\rm i}\nu k^2 \omega_n\right]\tilde{\xi}_x^n + n k_{x,w}k_z c^2\tilde{\xi}_z^n = 0,$$ (A.3)

$$n k_{x,w} k_z c^2 \tilde{\xi}_x^n +\left[\omega_n^2-k_z^2 c^2-i\nu k^2 \omega_n\right]\tilde{\xi}_z^n = 0,$$ (A.4)

where $\omega_n=\omega_{I,n} - n\omega$ (note that we assumed that $k^2c^2 \gg \Omega^2$ in writing the above equations). Neglecting terms in order $\nu^2$, these equations are in fact identical to their inviscid counterparts, Eqs. (39) and (40), provided $\omega_n$ is replaced by $\omega_{n,{\rm visc}}=\omega_n-i\nu k^2/2$. This indicates that inertial waves are damped on a timescale

$$\tau_{\rm damp}=\frac{2}{\nu k^2} \cdot$$ (A.5)

There will be a corresponding reduction in the growth rate of the parametric instability. In numerical simulations performed with Eulerian codes, we assume that the numerical dissipation has an effect similar to that describe above, with a "numerical'' viscosity  $\nu_{\rm num}$. For a given resolution, $\nu_{\rm num}$ is fixed and Eq. (A.5) shows that modes with larger wavenumbers are damped more rapidly, in agreement with our findings.

Another application of this result is to estimate a numerical Reynolds number Re for our simulations. Re is given by

$${\rm Re}=\frac{cH}{\nu_{\rm num}} \cdot$$ (A.6)

In Sect. 4.1, we found that the linear growth rate of the n=4 mode is unaffected by numerical dissipation, which indicates that

$$\tau_{\rm damp} \gg \frac{1}{\sigma} \cdot$$ (A.7)

This gives an upper limit for the numerical viscosity of the codes and a lower limit for the Reynolds number of our simulations. Introducing the parameters of the simulation in Eq. (A.7), we found $\rm Re \ge 10^3$. This is consistent with the value reported by Balbus et al. (1996) for the Reynolds number of Eulerian codes with similar resolutions.

References

• Arras, P., Blaes, O., & Turner, N. J. 2006, ApJ, 645, L65 [NASA ADS] [CrossRef]
• Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [NASA ADS] [CrossRef] (In the text)
• Balbus, S., & Hawley, J. 1998, Rev. Mod. Phys., 70, 1 [NASA ADS] [CrossRef] (In the text)
• Balbus, S. A., Hawley, J. F., & Stone, J. M. 1996, ApJ, 467, 76 [NASA ADS] [CrossRef] (In the text)
• Bate, M. R., Ogilvie, G. I., Lubow, S. H., & Pringle, J. E. 2002, MNRAS, 332, 575 [NASA ADS]
• Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Gardiner, T. A., & Stone, J. M. 2005, in Magnetic Fields in the Universe: From Laboratory and Stars to Primordial Structures, ed. E. M. de Gouveia dal Pino, G. Lugones, & A. Lazarian, AIP Conf. Proc., 784, 475
• Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125 [NASA ADS] (In the text)
• Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [NASA ADS] [CrossRef]
• Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [NASA ADS] [CrossRef]
• Goodman, J. 1993, ApJ, 406, 596 [NASA ADS] [CrossRef]
• Hawley, J., & Stone, J. 1995, Comput. Phys. Comm., 89, 127 [NASA ADS] [CrossRef] (In the text)
• Kato, S. 2001, PASPJ, 53, 1 [NASA ADS]
• Kerswell, R. R. 2002, Ann. Rev. Fluid Mech., 34, 83 [NASA ADS] (In the text)
• Korycansky, D. G., & Pringle, J. E. 1995, MNRAS, 272, 618 [NASA ADS]
• Larson, R. B. 1990, MNRAS, 243, 588 [NASA ADS] (In the text)
• Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646 [NASA ADS] [CrossRef]
• Lin, D. N. C., Papaloizou, J. C. B., & Savonije, G. J. 1990a, ApJ, 365, 748 [NASA ADS] [CrossRef]
• Lin, D. N. C., Papaloizou, J. C. B., & Savonije, G. J. 1990b, ApJ, 364, 326 [NASA ADS] [CrossRef]
• Lubow, S. H., & Ogilvie, G. I. 1998, ApJ, 504, 983 [NASA ADS] [CrossRef]
• Lubow, S. H., & Pringle, J. E. 1993, ApJ, 409, 360 [NASA ADS] [CrossRef]
• Lynden-Bell, D., & Ostriker, J. P. 1967, MNRAS, 136, 293 [NASA ADS] (In the text)
• Nelson, R., Papaloizou, J., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 [NASA ADS]
• Ogilvie, G. I. 2001, MNRAS, 325, 231 [NASA ADS]
• Ogilvie, G. I., & Lubow, S. H. 1999, ApJ, 515, 767 [NASA ADS] [CrossRef]
• Papaloizou, J. C. B. 2005a, A&A, 432, 757 [NASA ADS] [CrossRef] [EDP Sciences]
• Papaloizou, J. C. B. 2005b, A&A, 432, 743 [NASA ADS] [CrossRef] [EDP Sciences]
• Papaloizou, J. C. B., Nelson, R. P., & Snellgrove, M. D. 2004, MNRAS, 350, 829 [NASA ADS]
• Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1986, Numerical Recipes: The Art of Scientific Computing (Cambridge Univ. Press) (In the text)
• Shu, F. H. 1992, Physics of Astrophysics, Vol. II, ed. F. H. Shu (University Science Books), 476
• Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, ed. S. M. Fall, & D. Lynden-Bell, 111
• Torkelsson, U., Ogilvie, G. I., Brandenburg, A., et al. 2000, MNRAS, 318, 47 [NASA ADS]
• Ziegler, U., & Yorke, H. W. 1997, Comp. Phys. Comm., 101, 54 [NASA ADS] (In the text)