A&A 468, 1-18 (2007)
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
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 and (where 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
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 , where is the orbital angular velocity, it can interact resonantly with two linear inertial waves (whose individual frequencies, being bounded by 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.
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(X0,t) to be the position, or x coordinate,
of a fluid element that would be at the fixed location x=X0 in the absence of any
disturbance. The y component of the momentum
conservation equation can be written
Applying the same procedure to the equation of motion in the x direction, we find
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
When |p| is small, the wave becomes linear. This happens when
In this limit the potential is parabolic, being given by
|Figure 1: Profile of the function 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|
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 which contains the entire vertical dependence.
|Figure 2: Left panel: radial profile of the density for a wave of velocity (c/U)2=0.3, for different wave amplitudes: ( 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 ( 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:
The solutions we obtained are illustrated in Fig. 2. The left panel shows the structure of the wave when (c/U)2=0.3 for different values of the relative wave amplitude , 1.4 and 1.8 (we note that the maximum possible free wave amplitude is for these parameters). For these three values, we respectively found Tw=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)2=0.9. Again, three different wave amplitudes illustrate the transition from linear to nonlinear wave profiles. These are given by ,1.03 and 1.052. We note that similar waveforms occur for both values of c/U. However, for corresponding waveforms, when (c/U)2=0.9, the amplitudes and periods of the wave are much smaller than when (c/U)2=0.3.
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)2=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
is defined by
|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 ( upper row) and ( 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|
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 C0=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
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.
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
In the absence of a background wave, p=0 and the previous system
reduces to the simpler form
but is in the small wave amplitude limit (i.e. ),
terms of higher than first order in p may be neglected, so that
the system of equations becomes
equations relating the total of eight components
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
|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 . 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)2=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 kx,w and the wave frequency are then fixed. The instability couples two inertial waves with radial wavenumbers nkx,w and (n+1)kx,w. These should satisfy the resonant condition, given by Eq. (54), and the dispersion relations, given by Eq. (41) and its counterpart with For given values of n, and kx,w; , and kz can be determined from these algebraic conditions. At first we shall ignore any additional constraints on kz 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 ()
the solution is
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)2=0.3 and . The structure of this mode will be compared in the following section with the results of 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)2=0.3.
|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 .|
|Open with DEXTER|
We first study the case of a background wave with a small amplitude such that . 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 Tw=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 and 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 [-Tw/2,Tw/2] in x and 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 (Nx,Nz)=(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 +kz and -kz 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 +kz by adding a
velocity perturbation of magnitude
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
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 vz 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
|Figure 6: The time history of the maximum value of the vertical velocity ( left panel) and of the background wave amplitude ( 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 . The first two panels, both obtained with ZEUS, respectively show the density (normalised by the mean density ) 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|
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.
|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 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 . 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 of the instability measured during the linear phase for wave amplitudes such that , 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 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 . 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 . 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 (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 and (c/U)2=0.3. The upper figure plots 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 =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 as compared to 100 wave periods when (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 ( ) 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 ( ) 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. vz 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.
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.
The properties of linear inertial waves in stratified isothermal disks have been
studied by Lubow & Pringle (1993). We recall these below for
in their notation, as this is the case of interest
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
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 .
|Figure 12: Time history of the vertical velocity in the stratified model, for the small amplitude wave ( left panel) and for the large amplitude wave ( 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 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)2=0.3, and we consider only two wave amplitudes defined such that and . In our standard runs, the resolution is kept the same in the xdirection, Nx=320, as is the extent of the computational domain, . 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 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 kz and -kz 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, , and the right panel to the large amplitude wave . 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 and . 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, , while the right hand side panels refer to the large amplitude wave . 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 . 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
to 10. As an example, the structure of the flow
resulting from the mode n=9 is illustrated in
Fig. 14 (note that kz is determined
as the value that leads to the satisfaction of
the parametric resonance condition
and both modes with wavenumber +kz and
-kz 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
in rough agreement with
the simulation. However, it predicts a growth rate
larger than the numerically determined value
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
which we define through
Table 1: Modes ranked according to the precision (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, is given in the last column, while the second and third columns respectively give the values of n and nz. The parameters of the model are such that (c/U)2=0.3 and the background wave has an amplitude .
|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|
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 , 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 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.
S.F. acknowledges the important contribution of Romain Teyssier to the adaption of the numerical scheme of RAMSES to the shearing box equations.
Here we consider the effect of a small constant kinematic viscosity
the linear growth of the parametric instability. To do so, we study
the properties of unperturbed inertial waves (see
Sect. 3.2) when
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
Another application of this result is to estimate a numerical Reynolds number Re for our simulations. Re is given by