A&A 418, 325-335 (2004)
C. Schäfer1 - R. Speith2,1 - M. Hipp3 - W. Kley1
1 - Computational Physics, Auf der Morgenstelle 10C, 72076 Tübingen, Germany
2 - Department of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK
3 - Wilhelm-Schickard-Institut für Informatik, Sand 13, 72076 Tübingen, Germany
Received 2 July 2003 / Accepted 5 January 2004
We have performed Smoothed Particle Hydrodynamics (SPH) simulations to study the time evolution of one and two protoplanets embedded in a protoplanetary accretion disc. We investigate accretion and migration rates of a single protoplanet depending on several parameters of the protoplanetary disc, mainly viscosity and scale height. Additionally, we consider the influence of a second protoplanet in a long time simulation and examine the migration of the two planets in the disc, especially the growth of eccentricity and chaotic behaviour. One aim of this work is to establish the feasibility of SPH for such calculations considering that usually only grid-based methods are adopted. To resolve shocks and to prevent particle penetration, we introduce a new approach for an artificial viscosity, which consists of an additional artificial bulk viscosity term in the SPH-representation of the Navier-Stokes equation. This allows accurate treatment of the physical kinematic viscosity to describe the shear, without the use of artificial shear viscosity.
Key words: accretion, accretion discs - hydrodynamics - methods: numerical - stars: formation - stars: planetary systems
It is generally believed that planet formation is part of star formation. Stars are formed by gravitational collapse of a molecular cloud. During this process, angular momentum conservation leads to the formation of a circumstellar disc around a central object, the protostar. The material in the disc slowly accretes onto the protostar in the centre until the protostar becomes a star. Circumstellar discs composed of dust and gas have been found around young T Tauri stars (Beckwith & Sargent 1996). Star formation theory suggests a typical lifetime of such circumstellar discs in the range of 106 to 107 yr, which determines the time frame for planet formation (see Kallenbach et al. 2000 and references therein). It is assumed that planets form in these circumstellar (or protoplanetary) discs. Two mechanisms are proposed, either fast formation by direct fragmentation due to gravitational instabilities (Boss 2003,2001; Mayer et al. 2002), or slow formation during a process of several steps: at first the dust in the disc grows via a so-called hit-and-stick mechanism and settles in the mid-plane of the disc, where a dense layer forms (Beckwith et al. 2000). There the solid material is able to combine into larger bodies, which eventually form planetesimals. Finally, planetary-sized objects develop out of these planetesimals by collisions (Ohtsuki et al. 2002). The latter, slower scenario leads to a longer formation time. Calculations by Pollack et al. (1996) suggest a formation time of 1 to 10 million years for Jupiter and Saturn, and of 2 to 16 million years for Uranus, respectively.
This scenario leads to essentially circular orbits of the planets as they form in a keplerian disc. Due to the lack of material in the inner regions of the disc, more massive planets can only be created at larger distances to the central object. However, while these properties fit well to our solar system where the eccentricities are rather small and the massive planets are located in the outer regions, the model cannot explain the appearance of highly massive planets on significantly eccentric orbits at orbital radii down to 0.04 AU, as found among the extrasolar planets.
A possible solution to this problem is planet migration. It is generally assumed that these planets have formed further out in the disc and have migrated inwards later. This migration has its cause in tidal interactions of the planet with the disc (Goldreich & Tremaine 1980). The planet exerts tidal torques on the disc which lead to the formation of spiral density wave perturbations. Simultaneously, since the disc loses its axisymmetry, a net torque acts on the planet itself, which starts migrating in radial direction (Lin et al. 2000). The migration rate depends on different factors, such as the mass of the planet, the density distribution in the disc and the viscosity of the disc material.
Another problem lies in the high masses of the discovered planets. The tidal interaction of the planet can eventually lead to the opening of a gap at the planet's orbital region in the disc, which has a severe impact on the accretion rate. Until the numerical calculations by Artymowicz & Lubow (1996), who modelled the accretion on circumbinaries, and the simulations of a Jovian protoplanet embedded in a disc by Bryden et al. (1999); Kley (1999), it has been assumed that accretion stops after gap formation. These simulations however yield a final non-vanishing accretion rate which depends on the viscosity in the disc. Since then, further work has been done on this by various other authors (Nelson et al. 2000; Bryden et al. 2000; Lubow et al. 1999; Lin et al. 2000; Terquem et al. 2000).
While the numerical method known as Smoothed Particle Hydrodynamics (SPH) has been used for one of the first simulations that examined the system of an embedded protoplanet in a disc (Sekiya et al. 1987), later simulations were nearly exclusively performed with grid-based methods, mainly because computer resources have not been available for other high resolution simulations. This has changed only recently. In this paper we present a feasibility study for the simulation of planet-disc interactions with SPH, applying a new approach for the artificial viscosity, and we show the results of simulations with one and two embedded protoplanets in a viscous accretion disc. Very recently, also Lufkin et al. (2004) have used SPH for investigating planet-disc interaction. However, they focus on the scenario where planets form due to selfgravitating instabilities, while we consider non-selfgravitating discs in two dimensions ().
The outline of the paper is as follows. In the next section we present the physical model of the system, including an overview on gap formation and planet migration theory. Numerical issues will be discussed in Sect. 3, where especially the approach for a new artificial bulk viscosity will be described. The results of several simulations will be given in Sect. 4, followed by a discussion in Sect. 5. Finally, we draw our conclusions in Sect. 6.
We only consider a non-selfgravitating, thin accretion disc model for the protoplanetary disc. The disc is located in the z=0 plane, rotating around the z-axis. The vertical thickness H of the disc is assumed to be small in comparison to the radial extension r, . Therefore it is possible to vertically integrate the hydrodynamic equations and use vertically averaged quantities only, neglecting the vertical motion in the disc. The protostar in the centre of the disc always is assumed to have a mass of .
The equations describing a viscid flow of gas in the disc are the
continuity equation and the Navier-Stokes equation. In the Lagrangian scheme the
For computational simplicity we follow the usual approach and use an
isothermal equation of state, where the vertically integrated pressure p is related
to the surface density
Viscous processes are of great importance in the theory of accretion
discs. In order to move radially inwards, the disc material
has to lose its angular momentum. Therefore a permanent flow of angular
momentum radially outwards and of mass radially inwards takes place in the disc.
Often, theories of protostellar discs assume
viscous accretions discs, see e.g. Papaloizou et al. (1999).
We follow the prescription by Shakura & Sunyaev (1973) and assume an
effective -viscosity, where the kinematic viscosity
in the disc
is given by
Using perturbation theory, the tidal interaction of the planet and the accretion disc can be analyzed analytically if the planetary mass is much smaller than the mass of the central object, see Goldreich & Tremaine (1979,1980) for details. The planet exerts torques on the disc material (and vice versa), which can lead to gap opening at the planet's orbital region if these tidal torques exceed the internal viscous torques in the disc (Papaloizou & Lin 1984). The faster (slower) moving disc material inside (outside) the planet's orbit exerts a positive (negative) torque on the planet. The planet will begin to migrate radially as soon as these torques do not cancel out. There are two different timescales for migration depending on whether a gap has formed or not, which are called type I and type II migration respectively (Ward 1997).
As soon as the planet has grown to a sufficient high mass, the response of the disc to the perturbation becomes non-linear, and gap formation occurs.
We apply Smoothed Particle Hydrodynamics (SPH) for our simulations. SPH is a grid-free Lagrangian particle method for solving the system of hydrodynamic equations for compressible and viscous fluids first introduced by Gingold & Monaghan (1977) and Lucy (1977). SPH is especially suited for boundary-free problems with high density contrasts. Rather than being solved on a grid, the equations are solved at the positions of the so-called particles, each of which represents a volume of the fluid with its physical quantities like mass, density, temperature, etc. and moves with the flow according to the equations of motion. For a more detailed description see, e.g., Benz (1990) and Monaghan (1992).
Our SPH-code is fully parallelized to make optimal use of today's supercomputers which is necessary to achieve the required resolution and accuracy. The parallel implementation is done with the so called ParaSPH Library, a generalized set of routines to simplify and optimize the development of parallel particle codes (Bubeck et al. 1999). This library allows a clear separation between the implementation of the physical problem to be simulated and the parallelization. A programmer using the library works with parallel arrays containing the physical quantities of the SPH-particles and iterators to step through the particles and their neighbours. Domain decomposition, load balancing, nearest neighbour search and communication is hidden in the library. ParaSPH uses MPI (Message Passing Interface) for the communication and works on every parallel architecture providing an MPI library. On a Cray T3E one can achieve a parallel speedup of about 120 using 256 nodes for a problem with 360 000 particles. On the Beowulf-Cluster where we have performed the majority of our simulations - a Pentium III based Linux-Cluster with a Myrinet communication network located at the University of Tübingen - the parallel speedup is about 60 on 128 processors.
To model the physical shear viscosity we use the implementation
introduced by Flebbe et al. (1994)
and solve the Navier-Stokes equation including the entire viscous
stress tensor (6), in contrast to the usual
approach of a scalar artificial viscosity term. For the physical shear
flow, we assume a constant kinematic viscosity
and a vanishing
bulk viscosity .
Then, the acceleration due to the shear
In order to prevent the particles from mutual penetration,
we add an additional artificial bulk viscosity to the viscous stress
tensor. The bulk part of the viscous stress tensor is given by
and leads to an additional acceleration of particle i according to
In Appendix A we demonstrate the low intrinsic numerical shear stress of this new artificial bulk viscosity approach.
To study the migration of the protoplanet in detail, the
equation of motion of the protoplanet is integrated. The planet's
acceleration due to the gravitational forces reads
The initial setup of our standard model places a Jovian planet with mass at AU and a central object with mass , which gives a mass ratio of . The surface density falls off as r-3/4, and the particles with equal masses are distributed according to this density between AU and AU. The mass of the disc is set to and the initial particle velocities are vr = 0 and as in the unperturbed case. The typical number of particles at the beginning of the simulation is . The standard model uses a Mach number of , that is a scale height H/r = 0.05, and the kinematic viscosity parameter is set to , which corresponds to an value of about at the orbital radius of the planet.
As the density falls off with radius, fewer particles are located in the outer regions of the disc than near the central object. Hence, a constant number of interaction partners is more suited to the problem than a fixed smoothing length. Unless otherwise mentioned the number of interaction partners per particle is set to 100. As mentioned above, we use the XSPH-method to avoid mutual particle penetration and we set x=0.5 in all our simulations.
In this section we present the results of the simulations with various input parameters. First we consider the case of a single planet in the disc. Afterwards, the case of two planets and their mutual interactions is investigated.
In order to study the formation of a gap, we have performed several simulations of a disc with an embedded protoplanet. The disc scale height and the kinematic viscosity have been used as input parameters and have been varied accordingly. Here, the gravitational forces of the particles on the planet have not been considered, keeping it on a circular orbit around the central star. Moreover, the planet was assumed to be able to accrete material without growing in mass, its Roche lobe has been kept constant throughout the simulations.
Just after the start of the simulation the planet disturbs the density
distribution in the disc effectively, which leads to the typical spiral
density wave pattern in the disc which can be seen
in Fig. 1, where a greyscale plot of the surface
density is presented for the system after 10 orbits of the protoplanet.
In the case of this reference model, gap formation begins
immediately after the start of the simulation run.
orbits of the planet the decrease of the density at the protoplanet's
orbital region is already clearly visible. The evolution of the azimuthally
averaged surface density
during the gap formation process is shown in Fig. 2.
After one hundred orbits of the protoplanet, the decrease of density at
the orbital region is more than two orders in magnitude.
|Figure 1: Greyscale plot of the surface density in the disc in the case of the reference model after 10 orbits of the protoplanet in the disc. The gap clearing has already started, and the spiral density wave pattern is visible throughout the disc. The circle represents the Roche lobe of the protoplanet according to Eggleton's formula.|
|Open with DEXTER|
|Figure 2: Evolution of the azimuthally averaged surface density for five different times in the case of the reference model. The planet has a fixed orbit at 5.2 AU.|
|Open with DEXTER|
The radial size of the gap depends only slightly on the viscosity, as can be seen in Fig. 3, where the surface density is plotted for three simulation runs with different kinematic viscosity: (reference model), , and . A higher viscosity leads to a faster radial motion of the gas in the disc, and thus to more loss at the inner boundary, which results in a slightly lower density profile after the same number of orbits. A similar behaviour was found by Kley (1999). Because in Fig. 3 the simulation time is rather short compared with the viscous timescales, the gap structure looks very similar in each case. The differences may become more pronounced for longer simulation runs. Unlike the viscosity, the scale height has a severe impact on the depth of the gap as can be seen in Fig. 4. There, the surface density of the disc after 100 orbits of the protoplanet is plotted for the scale heights H/r = 0.1, H/r = 0.05, and H/r = 0.025. The smaller the scale height, the lower the temperature and the deeper the gap, as can be seen clearly. Additionally, the spiral density wave pattern gets tighter as the sound speed is lower for smaller scale heights.
|Figure 3: Influence of the kinematic viscosity parameter on the gap size. The azimuthally averaged surface density is plotted for three different viscosities in units of after 100 orbits of the protoplanet. The planet has a fixed orbit with a radial distance of 5.2 AU.|
|Open with DEXTER|
|Figure 4: Influence of the scale height of the disc. The greyscale plot shows the surface density in the disc after 100 orbits of the protoplanet for three different scale heights.|
|Open with DEXTER|
In order to study the migration of the protoplanet we consider again the standard reference model of a Jovian protoplanet in a circular orbit with radius . To achieve a steady-state particle distribution we proceed as in the last section and let the system evolve until the gap is nearly fully developed before we allow the gravitational force of the particles to act on the protoplanet. Then the simulation is restarted and the trajectory of the planet is integrated as well. For the reference model we expect pure type II migration on a viscous timescale. For this simulation we used nearly 360 000 particles to model a disc with mass and , . Particles that get closer to the protoplanet than 50% of its Roche radius are taken out of the simulation and assumed to be accreted by the protoplanet. In one model, the mass of the accreted particles is added to the mass of the planet, in a second, the mass of the planet is fixed for all times. The total simulation time is about 30 000 years.
The evolution of the radial distances of the planets is shown in Fig. 5. Initially the migration rate is very high and about the same for both the planet with fixed mass and the planet with increasing mass. After 3000 years the radial distances to the central object differ distinctly, with the more massive planet migrating slightly faster. Similar behaviour was found in numerical simulations by Nelson et al. (2000). The planet with a fixed mass ends up at an orbital distance of after 20 000 years, which would correspond to a steady-state migration rate of . The planet with increased mass has a radial distance of after 30 900 years, or a steady-state migration rate of . As can be seen in Fig. 5, the migration rates are not constant but decrease considerably. This is caused primarily by the mass loss of the disc, leading to less gravitational interaction between the planet and the disc.
|Figure 5: Evolution of the radial distance of the planet to the central object. The planet is initially located at a radial distance of . In one model, the accreted mass is added to the planet's mass, in the other the mass of the planet remains fixed.|
|Open with DEXTER|
|Figure 6: Evolution of the accretion rate on the planet for a planet with fixed mass and a planet with increasing mass. The planet is originally located at a radial distance of 5.2 , and migrates inwards as shown in Fig. 5.|
|Open with DEXTER|
The accretion rates of the planets for both models are presented in Fig. 6. The initial accretion rate of about drops rapidly in the first 3000 years to a value of for both the planet with constant mass and the planet with increasing mass. At that time, the inner part of the accretion disc has vanished, as the material was either accreted by the protoplanet or by the central object. After the loss of the inner disc, the accretion rate decreases more slowly. Between 3000 and 7500 years of simulation time, the planet with increasing mass has a slightly lower accretion rate. The situation changes after 10 000 years, when the planet with the fixed mass accretes less. At the end of the simulation, both planets have about the same accretion rate of , although the planet with mass accumulation has already doubled its mass. However, during that time the constant loss of material through the open boundaries and the accretion on the planet lead to a decrease in disc mass of nearly one magnitude, which has to be considered for comparisons to grid-based simulations with closed boundaries.
As it is generally assumed that the mutual interaction of gravitating objects in an accretion disc can finally lead to eccentric orbits, we have included another protoplanet into our simulations.
For the simulation with two protoplanets the setup was
changed as follows: one protoplanet is located at a distance of
the central gravitating object, the other one at twice that. The
kinematic viscosity parameter is again set to
and the disc scale height is 0.05. Now the outer radius of
the disc is
and the inner boundary is set to
number of particles in this run was 320 000 and the disc mass was set
The evolution of the density in the disc is
shown in Fig. 7. As in the case of a single Jovian
planet, gap formation starts immediately after the beginning of the
simulation. This happens for each planet individually, until the particles between the
inner and outer planet have been accreted by the inner planet, and a common
gap has formed. The innermost region of the disc has already vanished after 200 orbits of
the inner planet.
|Figure 7: Evolution of the azimuthally averaged surface density in the case of two protoplanets in the disc. The number of orbits of the inner planet is indicated. The inner planet is located at , the outer planet at twice that distance from the central object.|
|Open with DEXTER|
After 200 orbits of the inner planet, the simulation is restarted and the
trajectories of the two planets are integrated to examine their migration.
The evolution of the semi-major axis of the two planets is
shown in Fig. 8 for the next 30 000 years.
During the first 5000 years the inner planet mainly stays on its orbital
radius around the central star because it is not
surrounded anymore by disc material (see Fig. 7)
which could exert torques. On the other hand
the outer planet migrates inwards due to the action of the outer disc,
resulting in a radial approach of the two planets. As the orbital distance
of the two planets decreases, they capture each other in a 2:1 mean motion
resonance at approximately t= 7500 years. Upon capture
their eccentricities slowly increase, as
can be seen in Fig. 9.
By the end of the simulation, the eccentricity of the outer planet has
grown to 0.1 and the eccentricity of the inner planet to 0.18. In the end,
the planets have reached masses of 1.34 (inner planet)
(outer planet), respectively.
The mass of the disc changed from
at the start
of the simulation to
the end. After all, only 20% of the mass lost by the disc
have been accreted by the protoplanets.
The resonant capture and
overall evolution of this system
is nearly identical to that described in Kley (2000).
|Figure 8: Evolution of the semi-major axis of the two protoplanets. The initial parameters of the planets are the same as in Fig. 7, but now the trajectories of the planets are integrated.|
|Open with DEXTER|
|Figure 9: Evolution of the eccentricities of the two protoplanets during the migration, same simulation as presented in Fig. 8.|
|Open with DEXTER|
In the case of a single Jovian planet we obtain the same influence of the disc scale height on the gap formation and on the structure of the spiral wave pattern as found already by several other authors (see, e.g., Bryden et al. 1999; Kley 1999), and as postulated by the theory. For the parameters of the reference model, Eq. (13) yields a migration rate of for pure type II migration, for a planet located at a radial distance of . We find migration rates of the same order of magnitude by averaging over the complete simulation time. The migration rate for this simulation agrees well with migration rates found by D'Angelo et al. (2003), who examined the migration and growth of one protoplanet in a disc for various initial planet masses, using a grid-based method and closed boundary conditions.
Interestingly, the accretion rate for a planet with a fixed mass and a planet with increasing mass do not differ considerably at the end of a simulation time of 30 000 years. Presumably due to the loss of disc material at the open boundaries, the resolution is too coarse to resolve the difference correctly. After a simulation time of 15 000 years, the accretion rate onto the planets is about yr for our reference model with a kinematic viscosity coefficient of . This value agrees well with the results from simulations performed by Lubow et al. (1999), who used the ZEUS hydrodynamics code to simulate disc accretion onto high-mass planets. They find an accretion rate of yr, which is half the value that was published by Kley (1999). However, Lubow et al. (1999) use different boundary conditions, namely reflecting closed boundaries at the inner edge of the disc, and open boundaries at the outer edge. The results of the simulations by Bryden et al. (1999) provide accretion rates in the same range, although they use an initial density profile which is constant throughout the disc. Extensive simulations with three different grid codes, nirvana, rh2d, and fargo, by Nelson et al. (2000) provide comparable accretion rates for a Jovian planet.
For a single planet in the accretion disc the orbit remains circular, although the planet migrates towards the central star. This suggests that the gravitational interaction between the planet and the protoplanetary disc cannot excite a change in the eccentricity of the orbit of the planet, at least for low mass planets (Papaloizou et al. 2001). Therefore, in systems with a single extra-solar planet mostly circular orbits should be expected. This clearly appears to be in contradiction to the observations, and is why different mechanisms, e.g. resonances, for the evolution of high eccentricities in single-planet systems have to be assumed.
We have found that a number of 360 000 SPH-particles, which has been used in the simulation of migration and planet accretion, is sufficient in 2D to resolve an expected accretion rate of several yr. However, the use of an artificial bulk viscosity is essential, as extensive test calculations have shown. Without the use of this additional viscosity, the particle orbits penetrate each other, density shock waves cannot be resolved, and even the formation of the gap may be suppressed. Note that in contrast to the normally used artificial -viscosity approach (see Monaghan & Gingold 1983) which leads to an effective physical shear viscosity (Nelson et al. 1998), this artificial bulk viscosity has very little impact on the shear kinematic viscosity in the disc, as can be seen in Appendix A.
The number of SPH-particles was not high enough though to study the flow around the protoplanet within its Hill radius, as was done by D'Angelo et al. (2002) using nested-grid techniques. This would have made it possible to examine the accretion process onto the protoplanet more precisely. D'Angelo et al. (2002) find the formation of a so-called proto-Jovian disc around the protoplanet which dominates the accretion process onto the planet. However, in our simulations the spatial resolution was too low to resolve the formation of such a disc.
Generally, SPH seems to be slightly more computationally expensive for this kind of application than comparable simulations with grid-based schemes. This is mostly due to the intrinsic noise of the particle method that requires a comparatively high particle number to achieve a similar spatial resolution. The advantage of this numerical method is its Lagrangian nature, which allows an easy treatment of open boundaries and complex geometries, like multiplanetary systems.
In order to study the evolution of a protoplanetary system with several planets, we included a second protoplanet in the simulation. After the first 200 orbits of the inner planet, the inner disc has already nearly vanished leading to a lower accretion and migration rate of the inner planet in comparison to the outer planet. Due to different migration speeds of the two planets they approach each other radially and are eventually captured into a 2:1 mean motion resonance. Upon capture, the dynamical interaction between the two planets increases and their eccentricities begin to grow considerably. This resonant behaviour has been found and investigated numerically already by several authors (e.g. Kley 2000; Lee & Peale 2002; Nelson & Papaloizou 2002; Kley et al. 2004). Depending on the subsequent evolution (distance of migration and damping), the system might become unstable eventually.
Thus, the gravitational interaction between protoplanets in the disc can explain the appearance of resonant extra-solar planets with high eccentricities as found for example in GJ 876, HD 82943, or 55 Cnc.
According to our simulations, gap formation does not inhibit accretion onto the protoplanet. Our simulations indicate that a protoplanet can grow up to several Jupiter masses before it has migrated to the central star. The simulations with two protoplanets show a possible mechanism for the formation of a two-planet system with a lower-mass planet near the central star and a more massive planet with a larger semi-major axis. The results of our SPH simulations compare favourably with those of grid codes. Although the computational costs are rather high, and therefore SPH seems to be less effective to model this problem, it is well suited to verify simulation results achieved with other, grid-based numerical schemes because of its completely different numerical approach. Moreover, it may be useful for 3D-simulations or self-gravitating flows.
Part of this work was funded by the German Science Foundation (DFG) in the frame of the Collaborative Research Centre SFB 382 Methods and Algorithms for the Simulation of Physical Processes on Supercomputers. Additionally, CS wishes to thank the UK Astrophysical Fluid Facility (UKAFF) in Leicester for the kind hospitality during a visit, where parts of this work were initiated, and he would like to acknowledge the funding of this visit by the EU FP5 programme. We are also deeply grateful for the comments and suggestions of an anonymous referee which helped to clarify and improve the paper considerably.
|Figure A.1: Evolution of the azimuthally averaged surface density of the viscous ring at three different times, using only physical viscosity in the simulation. Dotted lines denote the analytic solutions.|
However, applying only physical viscosity according to Eq. (16) is insufficient to simulate the dynamical evolution of a protoplanet in a disc, because interpenetration of SPH particles occurs such that the spiral density wave pattern in the disc cannot be resolved. This problem is reliably prevented by the new artificial bulk viscosity presented in Eq. (19) in combination with the XSPH-scheme of Eq. (22). To test for any spurious numerical shear viscosity that may be introduced by the usage of these algorithms, a simulation similar to that presented in Fig. A.1 was performed with additional artificial bulk viscosity, where we set f = 0.5 and XSPH is switched on with x = 0.5. The results are shown in Fig. A.2, where the quantities equivalent to those in Fig. A.1 are plotted. As can be seen, the global evolution of the viscous ring, which is very sensitive to the effective shear viscosity in the simulation, seems not to be affected at all by the new artificial viscosity approach.
|Figure A.2: Evolution of the azimuthally averaged surface density of the viscous ring at three different times, additionally using artificial bulk viscosity and XSPH together with the physical viscosity. Dotted lines denote the analytic solutions.|
|Figure A.3: Evolution of the azimuthally averaged surface density of the viscous ring at three different times, applying only artificial bulk viscosity. The dotted line represents the initial analytic solution.|
|Figure A.4: Greyscale plots of the surface density of the inner part of the protoplanetary disc in the case of the reference model similar to Fig. 1 but in Cartesian coordinates, with ( upper panel) and without ( lower panel) using the XSPH-algorithm. Note that the planet cannot be seen on these plots, since it is located at a radial distance of .|
|Figure A.5: Evolution of the azimuthally averaged surface density of the viscous ring at three different times, applying only the usual artificial viscosity. The dotted line represents the initial analytic solution.|
This can be tested further by a simulation where only artificial bulk viscosity is used. The results are shown in Fig. A.3. Any evolution of the viscous ring away from the initial distribution, which is indicated by the dashed line, has to be due to inherent numerical shear in the simulation. It can clearly be seen that the intrinsic effective shear of the artificial bulk viscosity approach is very small. A more quantitative analysis gives a value of less than 9% of , which can be neglected in all our simulations of the protoplanetary disc. Note that the numerical effective shear viscosity is not constant throughout the ring but seems to be higher at the edges, especially at the inner edge, than in the centre.
It has turned out that it is necessary to use the XSPH-algorithm in addition to the artificial bulk viscosity to prevent mutual particle penetration more effectively. This can be seen in Fig. A.4, where two different simulations of the reference model of the protoplanetary disc are shown, in the upper panel with XSPH switched on and in the lower panel with XSPH switched off. The surface density is plotted in greyscale, similar to the result presented in Fig. 1, but only for the inner part of the disc and in Cartesian coordinates. Both plots show basically the same features, but in the case with XSPH the structures are more distinct and better resolved, while in the case without XSPH the distribution is more diffuse.
For the sake of completeness it should be mentioned, that a frequently used very efficient different approach to prevent particle interpenetration is the usual artificial viscosity introduced by Monaghan & Gingold (1983), and there especially the quadratic -term. However, this artificial viscosity suffers from a high level of spurious intrinsic shear. This becomes obvious in the simulation of the viscous ring shown in Fig. A.5, where only this artificial viscosity is used with a rather small value of . The effective physical shear is already of the order of .