A&A 399, 395-407 (2003)
DOI: 10.1051/0004-6361:20021783
R. Speith^{1} - W. Kley^{2}
1 - University of Leicester, Department of Physics and
Astronomy, Leicester LE1 7RH, UK
2 -
Universität Tübingen,
Institut für Astronomie und Astrophysik,
Abt. Computational Physics,
Auf der Morgenstelle 10, 72076 Tübingen, Germany
Received 26 July 2002 / Accepted 2 December 2002
Abstract
We study analytically and numerically the stability
of the pressure-less, viscously spreading accretion ring.
We show that the ring is unstable to small non-axisymmetric
perturbations. To perform the perturbation analysis of the ring
we use a stretching transformation of the time coordinate.
We find that to 1st order, one-armed spiral
structures, and to 2nd order additionally two-armed spiral features may
appear. Furthermore, we identify a dispersion relation determining the
instability of the ring. The theoretical results are confirmed in
several simulations, using two different numerical methods.
These computations prove independently the existence of a secular
spiral instability driven by viscosity, which evolves into
persisting leading and trailing spiral waves.
Our results settle the question whether the
spiral structures found in earlier simulations of the spreading ring
are numerical artifacts or genuine instabilities.
Key words: accretion, accretion discs - hydrodynamics - methods: numerical
In the theory of accretion discs, the idealised problem of a viscously spreading, pressure-less ring orbiting a central point mass on Keplerian orbits is often used to exemplify the main features of an evolving thin accretion disc, i.e. the inward transport of mass and outward transport of angular momentum (Pringle 1981; Frank et al. 2002). With the approximation of a small kinematic viscosity that is independent of the surface density, an analytic solution exists for this problem, stated originally by Lüst (1952) and later by Lynden-Bell & Pringle (1974).
This axisymmetric analytic solution of the viscous dust ring is frequently used as test problem for numerical methods developed to simulate accretion discs. The tested algorithms cover particle methods like smoothed particle hydrodynamics (SPH) (e.g. Flebbe et al. 1994; Speith & Riffert 1999) as well as grid-based codes like RH2D (e.g. Kley 1999).
However, numerical simulations of the evolving ring often show the appearance of additional structures in the disc. In the majority of cases, these structures consist of non-axisymmetric one-armed spirals, but sometimes also narrow concentric rings appear. Since their first discovery, it is controversially debated whether these structures are numerical artifacts or mathematical instabilities. While Maddison et al. (1996) show that for the axisymmetric case the concentric rings found in SPH simulations are numerical artifacts of the method, we demonstrate in this paper that the non-axisymmetric spiral structures result from a genuine physical instability of the problem.
Recently, Ogilvie (2001) analysed radially extended accretion discs in a work presenting a general model of eccentric accretion discs, of which the viscous ring is a simple special case. However, because his main objective was much wider, he did not elaborate on the stability properties of the ring.
The remainder of the paper is organised as follows. In the next section, we first summarise the basic hydrodynamic equations used for our analysis of the evolving disc, and present for reference the analytic solution of the viscously spreading dust ring.
In Sect. 3 we use a special perturbation technique, the time-stretching approach, which allows to obtain solutions of the spreading ring at different orders of a small expansion parameter. We derive a linear evolution equation for an eccentricity function E. It is shown that to 2nd order, the general solution allows the development of one-armed as well as two-armed spiral features.
In Sect. 4 we perform a local stability analysis of the final equation for E and deduce a dispersion relation which indicates an instability for some viscosity laws. These results are compared with numerical simulations we present in Sect. 5 using two different numerical methods, one grid-based and the other particle based. We find in general very good agreement of both numerical methods with the predictions of the local analysis. Thus, the simulations confirm the existence of a physical instability in the viscously evolving ring problem. Finally, in Sect. 6 we give our conclusions.
As we only consider cold (pressure-less) discs, we neglect the pressure throughout the rest of this paper.
For the time evolution of the viscously spreading ring, we assume
Then, with initial condition
The azimuthal velocity of the ring is equal to the Keplerian velocity,
To study the stability properties of the ring solution (8), we perform a perturbation analysis of the hydrodynamic Eqs. (1), (2) and (3).
According to the idealised problem, we assume that the averaged
kinematic viscosity
is very small and depends only on R.
To take this into account, we replace
(13) |
(14) |
To generalise the approach, we do not restrict the stability analysis
to a perturbation of the ring solution (8) but we start with
less determined unperturbed functions. The initial condition of the
problem shall consist of an arbitrary axisymmetric surface density
distribution
rotating around a central mass
.
Because self gravity of the disc can be neglected, the initial
azimuthal velocity shall be equal to the Keplerian velocity (9),
.
Then,
the following approach is made for the expansion
The approaches (16), (17) and (18) are inserted in (1), (2) and (3) considering the stretching transformation (15) and (12), and the terms of the resulting equations are ordered according to powers of .
There are three main differences of the present approach compared to previous analyses like those by Ogilvie (2001). The kinematic viscosity is assumed to be a function of radius only, while Ogilvie (2001) considered to be a function of surface density. Both dynamic and viscous timescale are taken into account, while Ogilvie (2001) assumed at the outset that non of the variables varies on the dynamic timescale. And the expansion is carried out up to second order, so that non-linear terms appear in the perturbation. Usually, a linear perturbation is performed.
In
,
the continuity equation reads
Inserting (22) into (19) and (20) results in
In
,
the continuity equation has the form
(28) |
(31) |
(34) |
Then, (35) can be solved for
,
(37) |
= | (38) | ||
= | 0 . | (39) |
The condition for the disappearance of the secular term in brackets in (35) is
That is, the surface density of the analytic solution is equal to the 0th order term of the perturbation. The azimuthal velocity of the analytic solution is also equal to the according 0th order term of the perturbation, i.e., equal to the Keplerian velocity (24). And the radial velocity of the analytic solution is equal to the axisymmetric part of the according 1st order term of the perturbation, i.e., relations (11) and (33) are identical.
Because of (40), all first order quantities are independent
of the dynamic time t. Therefore, the first order velocities take
the form
Defining a function
similar to
Ogilvie (2001)
with
Consider the 0th order results (22) and (24) and the
1st order results (42), (43) and (48). Then, in
one yields for the continuity equation
Solving (51) for
and inserting in (52) yields
(59) |
(60) |
(62) |
(63) |
There exists a second secular term in (50), which is
appearing already in (56) and (58) and which is
proportional to
.
Therefore, U_{1} has to vanish,
(65) |
To summarise, according to (64), all quantities up to second order are independent of the dynamic time t. The 2nd order quantities depend on the azimuthal angle in a linear combination of harmonic oscillations with frequencies , and , where and . Therefore, if a non-axisymmetric perturbation occurs, it has to be a superposition of one-armed and two-armed spiral structures in second order. This is an effect of the non-linearity of the analysis.
The solution of (66) rules the stability or instability of the
viscously spreading ring. Assuming that short wavelengths become
unstable first (which is supported by the numerical results), we make
the approach
(67) |
(71) |
Assume
Then, the dispersion relation (68) takes the form
(72) |
(73) |
(74) |
(76) |
(77) |
Figure 1: Evolution of the surface density of the viscous dust ring in an SPH simulation. Parameters of the disc are , disc mass , central mass , and kinematic viscosity . The initial density distribution at viscous time is shown in the top left plot. The symmetric structures seen in the distribution at may be of numerical nature due to relaxation effects of the initial particle distribution. Later, the development of the spiral instability can be seen very clearly. | |
Open with DEXTER |
The main limitation of the presented analysis is the omission of pressure. In particular, the more general analysis of Ogilvie (2001) shows that, in a disc with significant pressure, the eccentricity vector precesses rapidly and differentially due to pressure, in addition to its viscous growth or decay. When pressure dominates over viscosity, the afore mentioned characteristic radial scale where the analysis may break down is replaced by the semi-thickness of the disc. Furthermore, the eccentric instability can be enhanced or suppressed by allowing for the kinematic viscosity to depend on the surface density, by introducing a bulk viscosity, by allowing for stress relaxation, or by including three-dimensional effects.
In our simulations we did not use the artificial viscosity approach (see e.g. Monaghan & Gingold 1983) usually applied to model viscid fluids with SPH, but we chose a version of SPH that includes the entire viscous stress tensor (Flebbe et al. 1994; Riffert et al. 1995). This implementation allows for the proper treatment of viscous shear flows. Additionally, in compliance with the assumption of a cold disc, pressure was completely neglected in all the simulations. A detailed description of the application of this SPH approach to the viscously evolving ring can be found in Speith & Riffert (1999).
The simulation presented below was performed with the following parameters. The initial SPH particle distribution represents the surface density (8) of the analytic solution at the viscous time , with , disc mass , and central mass . The kinematic viscosity was set to throughout the whole simulation. The initial distribution consists of 40 000 particles which are arranged symmetrically to make sure that at the beginning the centre of mass lay exactly at the origin. Due to the stochastic nature of the SPH method, no further seed perturbations are required to bring potential instabilities to grow.
Figure 1 shows the evolution of the viscous ring during the SPH simulation in a sequence of plots of the surface density distribution. The top left plot represents the initial density distribution at viscous time . Immediately after the start of the simulation, symmetric structures develop, which can still be seen in the distribution at in the top middle panel. These narrow ring-like structures may be numerical artifacts due to the relaxation of the initial particle distribution, or they may be caused by the radial viscous overstability discussed e.g. in Kley et al. (1993), but they vanish in the further course of the simulation. Instead, the development of a spiral pattern begins, in this case at the inner edge of the disc, until finally a distinct spiral structure is covering the whole face of the disc, as can clearly be seen in the distribution at in the bottom right panel of Fig. 1.
To determine the nature of the spiral structure, a Fourier analysis of
the evolving ring was performed during the simulation.
Figure 2 gives the time evolution of the first four
azimuthal Fourier modes (m = 1 to m = 4) at radius R = R_{0}. For
small viscous times, the even modes m = 4 and especially m = 2 are
high, while the odd modes m = 1 and m = 3 are low. This effect is
due to the symmetric initial particle distribution and disappears
eventually when the symmetry of the particle distribution is lost.
Figure 2: Time evolution of the first four Fourier modes at R = R_{0} of the simulation shown in Fig. 1. Because of the symmetric initial particle distribution, even modes are high and odd modes are low at the beginning. Later, all modes decrease except the (m=1)-mode which increases exponentially as expected. The (m=3)- and (m=4)-mode merge with the noise background of the simulation, while the (m=2)-mode stays slightly above the noise level. The dashed line in the top left panel gives an approximated fit of the growth of the first Fourier mode to . | |
Open with DEXTER |
The results of the Fourier analysis, i.e., the development of a dominant one-armed spiral structure with an additional weaker component of a two-armed spiral structure, are in complete agreement with the theoretical results of the first and second order perturbation analysis as derived in Sects. 3.4 and 3.5.
Further predictions of the perturbation analysis in Sect. 3, that can be tested easily, are the relations between the first order velocities v_{R1} and . According to (45), their amplitudes have a ratio of 2, and they obey a phase shift of . To verify this, we more thoroughly studied the most evolved particle distribution of the SPH simulation presented in Fig. 1, i.e. the ring at viscous time , whose density distribution is shown in the bottom right panel of Fig. 1. In Fig. 3 the relation is plotted over radius. Here, and v_{R} denote the simulation results, which are averaged azimuthally using the SPH formalism, and and are the velocities according to the analytic solutions (9) and (11) of the viscous ring model. Taking and , the curve in diagram 3 matches the value sufficiently to satisfy relation (45). Accordingly, Fig. 4 shows the phase shift between and , and again the plot matches the expected value of satisfactorily.
Finally, we want to compare the local stability analysis and the resulting dispersion relation established in Sect. 4 with the SPH simulation results. From the top left panel of Fig. 2, the growth rate at radius R = R_{0} of the first Fourier mode may be estimated to , the latter measured in units of the angular velocity (10) at R = R_{0}. With dispersion relation (75), the corresponding wave number can be estimated to k = 42 R_{0}^{-1}, which results in a wavelength of . Comparing that result with the surface density distributions of Fig. 1 shows good agreement. That can be seen in particular in Fig. 5, where the surface density distribution of the ring at viscous time is plotted in polar coordinates, and where additionally lines of constant with k = 42 R_{0}^{-1} are drawn in for the region around . The lines match well the slopes and the wavelengths of the spiral structures. Note by the way that in this presentation of the surface density clearly can be seen that in this case the spiral structure consists of a leading (outer part of the disc) and a trailing (inner part of the disc) spiral.
Figure 3: Ratio of the amplitudes of the errors of azimuthal and radial velocities at for the disc shown in the bottom right panel of Fig. 1. Plotted is the azimuthally averaged relation over radius, where and v_{R}are the simulation results and and are the velocities (9) and (11) of the analytic solution. The curve roughly matches the expected value of (dotted line). | |
Open with DEXTER |
Figure 4: Absolute value of the phase shift between the errors of azimuthal and radial velocities, and , for the disc at shown in the bottom right panel of Fig. 1, azimuthally averaged and plotted over radius. The curve roughly matches the expected value of (dotted line). | |
Open with DEXTER |
Figure 5: Surface density distribution of the disc at as shown in the bottom right panel of Fig. 1 but plotted in polar coordinates (radius R over azimuth ) instead of Cartesian coordinates. Additionally, lines of constant with k = 42 R_{0}^{-1} are drawn in around . Note that the spiral structure consists of a leading and a trailing spiral. | |
Open with DEXTER |
(78) |
All the grid-based models presented in this paper use the same computational
domain
represented by
with
,
,
and
.
The ring
is covered in the basic model by
grid cells. However, to analyse resolution and numerical
effects we use different griddings ranging from
to
.
For our basic reference model we use for the constant dimensionless
viscosity
,
which is in fact identical
to the value used in the particle simulations.
Figure 6: Grey-scale plot of the density in polar coordinates at different evolutionary times for a model with a resolution. Times are indicated in orbital periods P_{0} at R=1 R_{0}. To compare with Fig. 1: one viscous time refers to 278 P_{0}. The time evolution of the system proceeds similar to the particle simulations. In the initial phase axisymmetric perturbations are seen first (t=17.45). They vanish later, and the system returns to axisymmetry (t=30.54). Then, from the inner parts the instability develops as a trailing spiral. At the end of the simulation a leading spiral appears in the outer region of the ring. | |
Open with DEXTER |
The initial setup consists of an axisymmetric density profile following the analytic density as given in Eq. (8). Because the original -function distribution is hard to represent numerically, we start, as above, with an initial density profile given here by the analytic solution at the viscous time , i.e. . The slight difference to the value 0.018, used in the SPH calculations, is of no importance for the subsequent evolution. For the standard viscosity this refers to the initial (dimensionless) time . Stated differently, one viscous evolution time corresponds to 278 orbital periods P_{0} for the given viscosity. The total dimensionless mass M in the disc is normalised to unity. The initial radial velocity is set to zero and the angular viscosity to .
The code is written such that axial symmetry is preserved exactly for an explicit viscosity-solver, i.e. purely axisymmetric initial conditions remain exactly axisymmetric. Hence, the initial density is disturbed randomly by typically 1 or to supply seed perturbations, and the subsequent longterm integration has to follow the evolution on a viscous timescale which is typically about 100 orbital periods for the standard model.
In Fig. 7 a radial cut through the density at the
fixed angle
(crosses) is shown for the standard model
(
grid points), for t= 52.4.
Figure 7: Radial density distribution at (crosses) and azimuthally averaged (triangles). The averaged distribution follows exactly the analytical zero order result of the diffusion equation , Eq. (8). Time is given in orbits at R=1 R_{0}. | |
Open with DEXTER |
Figure 8: Time evolution of the first 5 azimuthal modes for the standard model measured at the radius R=1. | |
Open with DEXTER |
In Fig. 8 the growth rates for the first 5
modes m=1 upto m=5 are displayed. The vertical axis refers to
the decimal logarithm of the amplitude. Clearly seen is the random initial
perturbation at the level of
which is reflected in a
start of all modes at the level of 10^{-3}. During the initial
evolution the rings spreads slowly, the amplitudes decline until
at later times (
t = 15-20) the m=1 mode begins to grow
exponentially with time. From the plot and additional runs
we may approximately determine the growth rate
for this m=1 mode
of the standard model to be
Of course, since the growth rate
depends on the
azimuthal wave number k and thus on the grid size as well,
we do expect a dependence on the numerical resolution.
Figure 9: Time evolution of the logarithm of m=1 mode (measured at R=1 R_{0}) for models with varying resolution in the angular direction. The radial number of grid points is fixed to 128. | |
Open with DEXTER |
In the next set of models the number of radial grid points was allowed to vary as well. The time evolution of the total radial kinetic energy is displayed in Fig. 10 for models with different radial (and azimuthal) grid resolutions. For very small radial resolutions N_{R}=64 there is very little or no growth. Then, for an increasing number of radial grid points the growth rates increase as well. This effect of faster growth for higher resolution is not an artifact of the numerical analysis but is to be expected from our analytical estimates. It was shown that the growth rates depend on the radial wave number k: the smaller the wavelength of the perturbation the faster the growth rates. But the smallest wavelength that can be resolved by numerical computations depends naturally on the resolution. Notice, that the wave crests in such a simulation are always resolved by only a few grid cells (see Fig. 7).
Figure 10: Time evolution of the total radial kinetic energy (dimensionless units) in the system for different numerical resolution, as indicated by the labels. The physical parameter of the models are those of the standard model. | |
Open with DEXTER |
To study the influence of the kinematic viscosity we varied from the standard value to higher and smaller values.
In Fig. 11 results are shown for several models.
A measure of the global deviation
of the density
from the axisymmetric structure has been plotted, with
defined as
It can be seen that upon increasing the viscosity the growth rates are also increasing, which is again in agreement with our analytical estimates. The initial increase of for the two models with the lowest viscosities ( and ) is related to the nearly ring like disturbances described before. These decline first, and then later the growth of the spiral disturbance sets in. The exact behaviour of the curves depends on the initial random perturbation as well.
Finally we display the obtained growth rates for different viscosity coefficients in Fig. 12, where the solid line denotes the results computed by the finite difference method. For larger viscosities the growth rates are also larger which is in agreement with Eq. (75). To compare exactly with the analytical results, the wavelength of the most unstable mode needs to be known as well. However, the wave number is a function of radius, making a detailed comparison not as easy. It is found in general that the obtained growth rates (as given for example in Fig. 12) are in the right order with the first term of Eq. (75) but are consistently lower. The radial variation of the quantities may easily account for this slight discrepancy.
Figure 11: Time evolution of the asymmetry , Eq. (79), for different kinematic viscosities. Starting from the reference value (standard, ) we lowered and increased the value of by the factors given in the key. | |
Open with DEXTER |
Figure 12: The growth rates (in units of at R=1) versus kinematic viscosity in units of the standard viscosity ( ) for different kinematic viscosities and for both numerical methods. | |
Open with DEXTER |
As a result of the perturbation analysis we found that the ring may
develop one-armed spiral structures to 1st order and one- and two-armed
spiral features to 2nd order. A local stability analysis of an
eccentricity function E provided a dispersion relation that led
for a constant kinematic viscosity to a
growth rate for secular spiral instabilities of
As a consequence for numerical simulations of accretion discs, the spiral instability has to be taken into account when using the spreading ring for instance as test problem for code development. Otherwise, occasionally emerging spiral features may be mistaken for numerical instabilities of the algorithms used for the calculation of the viscous forces.
The physical implications of our results are not as obvious. Because the kinematic viscosity is held axisymmetric throughout the whole stability analysis, the discovered non-axisymmetric instabilities may be of importance mainly in accretion discs where the viscosity is determined uniformly by external effects like irradiation.
Acknowledgements
We would like to thank the late Harald Riffert who contributed in many fruitful discussions substantially to this work and who suggested the use of the stretching transformation for the perturbation analysis. We also want to thank the referee, G. Ogilvie, for many helpful suggestions to improve and clarify this paper.Research in theoretical astrophysics at the University of Leicester is supported by a PPARC rolling grant. RS gratefully acknowledges funding by this PPARC rolling grant.
RS wishes to acknowledge the support of Advanced Micro Devices (AMD) for the Leicester Theoretical Astrophysics Group's Linux cluster on which some of the calculations were performed. Parts of his simulations were also performed on the GRAND supercomputer funded by a PPARC HPC grant.