A&A 431, 365-379 (2005)
DOI: 10.1051/0004-6361:20041085

A two-phase code for protoplanetary disks

S. Inaba1,2 - P. Barge1 - E. Daniel3,4 - H. Guillard4


1 - Laboratoire d'Astrophysique de Marseille, CNRS, BP 8, 13376 Marseille Cedex 12, France
2 - Department of Earth and Planetary Sciences, Tokyo Institute of Technology, 2-12-1 Oookayama, Meguro-ku, Tokyo 152-8551, Japan
3 - EPU - UMR CNRS 6595, 5 rue Enrico Fermi, 13453 Marseille Cedex 13, France
4 - ProgetSMASH - INRIA, 2004 route des Lucioles, 06902 Sophia Antipolis, France

Received 14 April 2004 / Accepted 21 October 2004

Abstract
A high accuracy 2D hydrodynamical code has been developed to simulate the flow of gas and solid particles in protoplanetary disks. Gas is considered as a compressible fluid while solid particles, fully coupled to the gas by aerodynamical forces, are treated as a pressure-free diluted second phase. The solid particles lose energy and angular momentum which are transfered to the gas. As a result particles migrate inward toward the star and gas moves outward. High accuracy is necessary to account for the coupling. Boundary conditions must account for the inward/outward motions of the two phases. The code has been tested on one and two dimensional situations. The numerical results were compared with analytical solutions in three different cases: i) the disk is composed of a single gas component; ii) solid particles migrate in a steady flow of gas; iii) gas and solid particles evolve simultaneously. The code can easily reproduce known analytical solutions and is a powerful tool to study planetary formation at the decoupling stage. For example, the evolution of an over-density in the radial distribution of solids is found to differ significantly from the case where no back reaction of the particles onto the gas is assumed. Inside the bump, solid particles have a drift velocity approximately 16 times smaller than outside which significantly increases the residence time of the particles in the nebula. This opens some interesting perspectives to solve the timescale problem for the formation of planetesimals.

Key words: hydrodynamics - planets and satellites: formation - solar system: formation

1 Introduction

Observations at various wavelengths reveal that most young stars are surrounded by circumstellar disks of gas and dust (Beckwith & Sargent 1996) which are believed to be the place where planets form. During the early stage of their evolution, protoplanetary disks are mainly composed of molecular hydrogen and micron-size solid particles (i.e., dust). The total mass of dust deduced from IR observations is in the range from $1 \times 10^{-5}$ to $1 \times 10^{-3} \;M_{\odot}$. The total mass of gas is about 100 times larger than the dust mass, a ratio similar to that found in the interstellar medium (Thi et al. 2001). IR emission from protoplanetary disks is found also to strongly decrease as a function of time, probably due to particle growth and to formation of planetesimals, on a time scale of the order of $3 \times 10^{6}$ yr (Briceño et al. 2001).

Understanding the processes responsible for particle growth requires one to investigate the dynamics of solid particles coupled to the gas disk. When gas/particle interactions are neglected, particles are known to move around the star with Keplerian velocity while gas moves at a slightly slower velocity since it is radially supported by a pressure gradient. Due to this difference in their rotational velocities, the solid particles experience a head wind which makes them lose energy and angular momentum. As a result, they spiral inward toward the star at a rate which depends on their size. The drift velocity is very slow for the largest particles, due to strong inertia, but is faster at smaller sizes. In a standard model of the nebula the drift velocity peaks at about 50 m/s or 1 AU/century for meter-size particles (Weidenschilling 1980) and is nearly independent of the distance from the star. On the other hand, the smallest particles are strongly coupled to the gas with a very slow drift inward toward the star. Such drift motions imply that growing solid particles may be lost into the star when they reach a characteristic size of the order of 1 m. The loss time scale is 103 yr, much less than the lifetime of the gas disk ( $3 \times 10^{6}$ yr), so that it seems very difficult to keep such particles growing in the gas disk before dispersal of the circumstellar material due to stellar evolution.

What happens in the vertical direction, perpendicular to the plane of the disk, should not be forgotten. Solid particles experience gas drag and sink toward the equatorial plane where they form a sub-layer of increasing density. In the pioneering scenario proposed by Safronov (1969) and by Goldreich & Ward (1973) the dust sub-layer may become gravitationally unstable and fragments into a number of clumps that condense under their own gravity and form kilometer-size planetesimals. This model is very attractive because solid particles rapidly reach comet size, avoiding the sticking problem of a continuous growth scenario. Once formed, these primordial bodies are large enough to survive much longer against losses into the star (Goldreich & Ward 1973; Safronov 1969; Sekiya 1998). However, Weidenshilling (1980) and Cuzzi et al. (1993) pointed out that the gravitational instability could be inhibited due to turbulence induced by the shear between the gas layer and the particle sub-layer. Indeed, particles in the sub-layer have a nearly Keplerian velocity whereas in the upper regions of the disk, far from the equatorial plane, gas dominates and drives the rotational velocity. So, the gas motion is sheared in the vertical direction and is thought to produce turbulence in the thickness of the disk. Turbulence tends to inhibit particle sedimentation toward the equatorial plane and, consequently, gravitational instabilities. Only particles larger than 1 m radius are still able to sink because they cannot be stirred up by turbulence (Cuzzi et al. 1993). So, if not destroyed under radial drift, a sub-layer of particles of this size would be dense enough for gravitational instability to occur.

Another point to be taken into account is the exchange of energy and angular momentum between the gas and the solid component. It is expected to play an important role in the equatorial dust-dominated sub-layer. For example solid particles are stirred up by the turbulent motions but turbulence may be also damped through interaction with the particles. So far, the self-consistent evolution of a protoplanetary disk in which the transfer of energy and angular momentum are allowed between the two components, both gas and solid, has never been investigated in detail.

The present code is developed to address the problem and will help to explore planetary formation at the decoupling stage when the gas and solid phases are competing with one another. Some authors have also addressed this problem but either incompletely (under specific assumptions) or locally (in peculiar configurations). Nakagawa et al. (1986) used a two-fluid description and, under a steady state assumption, derived analytical expressions of the drift velocities as a function of the solid/gas mass density ratio and of the radial pressure gradient. They showed that the drift velocity of a solid particle can becomes four times smaller than that given by Adachi et al. (1986) when the solid to gas mass density ratio is of the order of one. Cuzzi et al. studied a dust/gas evolution with a rough axisymmetrical model of the shear induced turbulence. More recently, Johansen et al. (2004) used local shearing box simulations to study the evolution of gas and solid particles in a vortex configuration.

In this work our goal is to study the full time evolution, outside the context of a steady state assumption, and using 2D global simulations. The numerical code includes a self-consistent description of the gas/particle interactions: gas and solid particles are treated as two phases of fluid and interact with each other through aerodynamical drag forces. The accuracy of the code is tested against a number of the analytical solutions obtained by Nakagawa et al. (1986). A 2D code is obviously easier to develop due to shorter computing times but also permits us to explore many initial conditions and a wide range of parameters. This code will be sufficient, in a first step, to study the evolution of the drift velocities in the disk. Of course the actual problem is a 3D one since the vertical evolution leads to an equatorial sub-layer strongly enriched in solid particles. Building such a code will be the next step of our work.

In Sect. 2 we show the basic equations describing the time evolution of a flow of gas and particles around a star like the Sun. The numerical method and the boundary conditions are described in Sect. 3. Tests and comparisons are presented in Sects. 4 and 5. It is shown that the numerical results agree quite well with 1D analytical solutions and a reasonable grid resolution for the code is given, allowing further exploration of the gas/solid evolution in the protoplanetary disks. Section 5 is devoted to tests in the 2D case. We apply the code to the lifetime problem of solid particles in a disk in Sect. 6. Conclusions are presented in Sect. 7.

2 Basic equations and initial conditions

The starting point is a flat disk of gas and solid particles rotating around a solar-like star and our goal is to study the time evolution of the two components under aerodynamical coupling. The problem is limited to a two-dimensional flow where all physical quantities (e.g., density) only depend on r and $\theta $, where r is the distance from the star and $\theta $ is the angle between the x-axis and the position vector.

2.1 Equations for the gas and solid phases

The protoplanetary disk is described as a two phase flow, with gas and solid particles, in nearly Keplerian motion around the star. Governing equations for the gas and solid components are mass conservation, Euler and energy conservation equations. The viscosity of the gas is completely neglected and Euler equations are used instead of Navier-Stokes equations. A second assumption is that the solid component induces no pressure because the number density and the collisional frequency of the particles are small. Further, the solid phase is assumed to be in the form of non-growing identical spheres, a simplifying assumption but also a way to study the disk evolution as a function of the particle size. Gas and solid particles are coupled by aerodynamical drag forces and can exchange energy and angular momentum.

Using cylindrical coordinates centered on the star the basic equations for the gas component are:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle \frac{\partial \sigma_... ...x{\scriptsize g}} \frac{G ~M_{\odot}}{r^2} \end{array}\right. \end{displaymath} (1)

where $\sigma_{\mbox{\scriptsize g}}$ is the gas surface density, p and $e_{\mbox{\scriptsize g}}$ are vertically integrated pressure and specific energy of the gas, respectively; $u_{\mbox{\scriptsize g}}$ and $v_{\mbox{\scriptsize g}}$ are the radial and tangential velocities of the gas; Fr and $F_{\theta}$ are the radial and tangential components of the drag forces acting on the solid particles, respectively. Specific energy of the gas is obtained with its pressure and kinetic energy:

 \begin{displaymath}\displaystyle \sigma_{\mbox{\scriptsize g}} e_{\mbox{\scripts... ...g}} (u_{\mbox{\scriptsize g}}^2 + v_{\mbox{\scriptsize g}}^2), \end{displaymath} (2)

where $\gamma$ is the ratio of specific heats at constant pressure and volume (we set up $\gamma=1.4$). The equation of state is that for an ideal gas:

 \begin{displaymath}\displaystyle p = \sigma_{\mbox{\scriptsize g}} k_{\mbox{\scriptsize B}} T_{\mbox{\scriptsize g}}/\mu, \end{displaymath} (3)

where $k_{\mbox{\scriptsize B}}$ is the Boltzmann constant, $T_{\mbox{\scriptsize g}}$ is the temperature and $\mu$ is the mean molecular weight of the gas. For a single phase of a gas disk, the third equation of Eq. (1) can be written as

 \begin{displaymath}\displaystyle \frac{\partial l_{\mbox{\scriptsize g}} } {\par... ...e g}} v_{\mbox{\scriptsize g}} + p r)} {\partial \theta} = 0, \end{displaymath} (4)

where $l_{\mbox{\scriptsize g}}$ is angular momentum of the gas per unit surface and is given by

 \begin{displaymath}\displaystyle l_{\mbox{\scriptsize g}} = r \sigma_{\mbox{\scriptsize g}} v_{\mbox{\scriptsize g}}. \end{displaymath} (5)

Equation (4) conserves angular momentum, similar to the third equation of Eq. (1) in Li et al. (2001). For the solid component, the equations are similar to the gas equations but without the pressure terms; they read:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle \frac{\partial \sigma_... ...{\mbox{\scriptsize d}}}{r} -F_{\theta}, \\ \end{array}\right. \end{displaymath} (6)

where $\sigma_{\mbox{\scriptsize d}}$ is the surface density of the solid component, $u_{\mbox{\scriptsize d}}$ and $v_{\mbox{\scriptsize d}}$ are the radial and tangential velocities of the solid phase, respectively. The gas drag acting on the particles and the "particle drag'' acting on the gas have the same modulus but opposite signs because of Newton's third law.

The expression of the drag force differs according to the ratio of the particle radius to the mean free path of the gas molecules, $r_{\mbox{\scriptsize d}}/l$. If it is larger than 1.5, the drag force is given by the Stokes formula:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle F_r = \frac{3}{2} \fr... ...iptsize d}} - v_{\mbox{\scriptsize g}}),\\ \end{array}\right. \end{displaymath} (7)

where $\rho_{\mbox{\scriptsize g}}$ is gas density and c the velocity of sound given by:

 \begin{displaymath}\displaystyle c = \sqrt{\frac{8 k_{\mbox{\scriptsize B}} T_{\mbox{\scriptsize g}}}{\pi \mu}}, \end{displaymath} (8)

$\rho_{\mbox{\scriptsize mat}}$ is the density of the solid material (a mean value $\rho_{\mbox{\scriptsize mat}}=3\;\mbox{g/cm}^3$ is assumed). The gas density is set to a mean value equal to that in the equatorial plane of the nebula:

 \begin{displaymath}\displaystyle \rho_{\mbox{\scriptsize g}} = \frac{\sigma_{\mbox{\scriptsize g}}}{\sqrt{\pi} h_{\mbox{\scriptsize g}}}, \end{displaymath} (9)

where $h_{\mbox{\scriptsize g}}$ is the scale height of an isothermal gas disk given by (Hayashi et al. 1985)

 \begin{displaymath}\displaystyle h_{\mbox{\scriptsize g}}= \frac{\sqrt{\pi} c}{2 \Omega_{\mbox{\scriptsize K}}}\cdot \end{displaymath} (10)

Note that, with the above assumptions, the gas density is set to the maximum value in the vertical direction, so the drag force is systematically overestimated. On the other hand if $r_{\mbox{\scriptsize d}}/l \leq 1.5$, the drag law is given by the Epstein formula:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle F_r = \frac{\rho_{\mbo... ...scriptsize d}} - v_{\mbox{\scriptsize g}}). \end{array}\right. \end{displaymath} (11)

The evolution of a 2D protoplanetary disk with two phases, gas and solid, is then described by solving simultaneously Eqs. (1) and (6).

2.2 Initial conditions

Initial conditions are chosen within the framework of the disk model introduced by Hayashi et al. (1985). The gas is assumed to be a mixture of hydrogen molecules and helium with a small amount of heavy elements; its mean molecular weight is $\mu \simeq 2.34~{m_{\mbox{\scriptsize H}}}$ where $m_{\mbox{\scriptsize H}}$ is the mass of a hydrogen atom. The surface density of the gas and its temperature are initially given by simple power law distributions:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle \sigma_{\mbox{\scripts... ...{r}{\mbox{1~AU}}\right)^{-1/2} \mbox{K}.\\ \end{array}\right. \end{displaymath} (12)

Similar choices are made for the solid component assuming that power law distributions have the same indices but different magnitudes.

Initial values of the velocities, for the gas and the particles, are chosen using the analytical expression derived by Nakagawa et al. (1986) in which the surface density and the pressure of the gas have a constant power index (see Eq. (12)). So, the velocity of the gas is given by:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle u_{\mbox{\scriptsize g... ...pha)^2 g^2} \eta v_{\mbox{\scriptsize K}}, \end{array}\right. \end{displaymath} (13)

where $v_{\mbox{\scriptsize K}}$ is the Keplerian velocity, $\alpha = \sigma_{\mbox{\scriptsize d}}/\sigma_{\mbox{\scriptsize g}}$ is the ratio of the surface densities of solid and gas, and g is the non dimensional friction coefficient corresponding to drag laws which reads:

 \begin{displaymath}g = \left\{ \begin{array}{l} \displaystyle \frac{3}{2} \frac... ...iptsize d}}} \qquad \mbox{Epstein regime}, \end{array}\right. \end{displaymath} (14)

and $\eta$ is the deviation from Keplerian velocity due to the gas pressure gradient:

 \begin{displaymath}\displaystyle \eta = - \frac{r \partial p/\partial r} {2 \sigma_{\mbox{\scriptsize g}} v_{\mbox{\scriptsize K}}^2}\cdot \end{displaymath} (15)

On the other hand, the velocity of the solid component is given by:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle u_{\mbox{\scriptsize d... ...)^2 g^2} \eta v_{\mbox{\scriptsize K}}.\\ \end{array}\right. \end{displaymath} (16)

Independently of the values of the parameters, the solid component has negative radial velocity which corresponds to inward drift, while the gas component has a positive radial velocity and migrates outward. In Fig. 1, the drift velocities of the solid particles are plotted for various values of the particle radius, 1 cm, 10 cm and 100 cm, and for two different cases: $\alpha =0$ and 1. The first one is the limit of the weak coupling situations ( $\alpha \ll 1$) in which there is little back reaction of the gas from the solid particles (particles look like passive tracers), therefore the migration velocity of the gas is nearly zero and the drift velocity of the solid particles has large values. The case $\alpha = 1$ corresponds to a situation in which the gas and the particles are completely coupled (particles play an active role), therefore the drift velocity of the solid particles and the migration velocity of the gas have the same amplitude but opposite signs. The rotational velocity of the gas is slightly slower than the Kepler one because the gas component is supported by the pressure against the gravity of the star. The solid component is not pressure-supported but has also an azimuthal velocity somewhat smaller than the Kepler one since it is subject to the gas drag.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{1085f1.eps}\end{figure} Figure 1: Radial velocity of the gas (dotted lines) and the solid component (solid lines) for particle radii of 1 cm, 10 cm, and 100 cm. $\alpha \ll 1$ corresponds to the case of passive particles and $\alpha = 1$ corresponds to the case of active particles, fully coupled to the gas. The size of the "optimal'' particles, those with the highest migration velocity, depends on the distance from the star.
Open with DEXTER

The gas drag is proportional to the relative velocity between the two phases, that is the small velocity difference $\delta v$ between the two phases. Therefore, study of the coupled evolution of the two phases requires numerical errors to be much smaller than the solid/gas velocity difference. High accuracy is required to compute the evolution of the velocities. This point is illustrated in Fig. 2 which shows the velocity differences between the gas and solid components for various particle sizes. The velocity of the particle relative to gas is only 10-6 times the mean rotational velocity. Numerical errors were determined from series of simulations. Small particles, with size less than 1 cm, have very small velocity differences with respect to the gas and are strongly coupled to the motion of the gas. This is why in Sect. 4 we used the relative velocities of centimeter-size particle to test the accuracy of the code and estimate the grid resolution required for subsequent simulations.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{1085f2.eps}\end{figure} Figure 2: Relative velocities $\delta v$ between the gas and the solid particles for the cases $\alpha =0$ (solid lines) and $\alpha = 1$ (dotted lines).
Open with DEXTER

In the limit of a single phase disk, only composed of gas, the radial and tangential velocities are given by:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle u_{\mbox{\scriptsize g... ...frac{2p}{\sigma_{\mbox{\scriptsize g}}}}\\ \end{array}\right. \end{displaymath} (17)

(i.e., the gas flows steadily around the star and does not migrate in the inward or outward direction).

3 Simulation method

The fluid equations for the gas and solid components, Eqs. (1) and (6) of the last section, are solved simultaneously using the finite volume method with an operator splitting procedure. The main steps of the code are described in Thevand et al. (1999). In this section we present some parts of the code, either the most important ones or the most specific to our problem. The flow of gas and particles is studied in a ring around the star with inner and outer radii denoted by $r_{\mbox{\scriptsize in}}$ and $r_{\mbox{\scriptsize out}}$, respectively. For the numerical description, the variables are discretized on an annular grid divided into nr and $n_\theta$ pieces, respectively for the radial and azimuthal directions. The grid is regularly spaced with:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle \Delta r = \frac{r_{\m... ...ta \theta = \frac{2 \pi}{n_\theta}\cdot \\ \end{array}\right. \end{displaymath} (18)

The cells are labeled with the index i in the r-direction while they are labeled with the index j in the $\theta $-direction. The position of the cell, labeled (i,j), is given by:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle r_i = r_{\mbox{\script... ...) \Delta \theta \qquad (j=0,n_\theta+1).\\ \end{array}\right. \end{displaymath} (19)

The configuration of the cell is illustrated in Fig. 3. In the r-direction the boundaries of the cell are given by $r_i+ \Delta r/2$ and $r_i- \Delta r/2$, while in the $\theta $-direction they are given by $\theta_j+\Delta \theta/2$ and $\theta_j-\Delta \theta/2$, respectively. With the finite volume method, the average values of the variables are stored for each cell, then their time evolution is computed. For example, $\sigma_{\mbox{\scriptsize g,ij}}^n$ is the surface density of the gas in the cell (i,j) at time t=tn.


  \begin{figure} \par\includegraphics[width=7cm,clip]{1085f3.eps}\end{figure} Figure 3: Configuration of the cell labeled by i and j in our 2D polar grid around the star.
Open with DEXTER

In the operator splitting procedure the changes in the variables due to the source and advective terms are computed successively at each time step; then the procedure advances to the next time step, $t^{n+1}=t^n + \Delta t$. If we denote the solution operators for the source terms and the advective terms by Ct and St, the surface density of the gas in the cell (i,j) at t=tn+1 is calculated by (Eq. (15.18) in Toro 1999)

 \begin{displaymath}\displaystyle \sigma_{\mbox{\scriptsize g,ij}}^{n+1} = S^{\D... ...\Delta t} S^{\Delta t/2} \sigma_{\mbox{\scriptsize g,ij}}^{n}. \end{displaymath} (20)

With the above calculation, we have a second order scheme as a whole if S and C have at least a second order accuracy. The source terms are computed with a fourth order Runge-Kutta scheme while the numerical flux corresponding to the advection terms is obtained using a second order MUSCL-Hancock scheme and an exact Riemann solver. Since the computation cost with the Runge-Kutta scheme is not large, the fourth order is preferred to the second order. There are also differences following the phase of the fluid (i.e., gaseous or solid). For the gas component, a "van Leer'' flux limiter helps to remove possible spurious oscillations. For the solid component, no limiter is used except if the particle surface density in a cell becomes 1.1 times larger than that in the neighboring cells; in that case we use a van Leer limiter for the solid component, as for the gas component.

The choice of appropriate boundary conditions at the inner and outer limits of the ring is also difficult. We have chosen extrapolated boundary conditions to permit connection of the simulated ring to the rest of the disk. From the initial conditions, Eq. (12), we know that surface densities (of gas and solid) and gas temperature are given by the power law distributions whose indices are assumed constant within the disk boundaries. Boundary conditions are set up using the surface densities and the pressure of three inner cells in the computational domain (i=3, 4, 5) and extrapolating them to find values inside the ghost cells (i=0, 1, 2). The number of ghost cells in our code is larger than 2, as is commonly used. Indeed, from our experience with the present code, numerical errors decrease with increasing ghost cell numbers and are steady for more than three ghost cells; this is the reason for our choice. The surface densities (for gas and solid) and the gas temperature in the inner ghost cells are given by:

 \begin{displaymath}\left\{ \begin{array}{l} \displaystyle \sigma_{\mbox{\scripts... ...}} \left(\frac{r_2}{r_3}\right)^{-3/2}.\\ \end{array}\right. \end{displaymath} (21)

Variables in the outer ghost cells ( i=nr+1, nr, nr-1) are derived in the same way as for inner ghost cells. Velocities of the gas and the particles in the ghost cells are obtained using the analytical solution, Eqs. (13) and (16). All variables are determined in the ghost cells and the numerical fluxes can easily be computed. Since only a small part of the disk is simulated, the boundary conditions are chosen to determine the fluxes (of gas and particles) with the highest accuracy. Further, we are interested in local phenomena that do not influence the boundary itself (see the next sections). The boundary conditions close the system of equations and we are able to perform numerical simulations.

4 Testing the code accuracy in one dimension

Now, we check whether the code has the required accuracy to compute the small velocity differences between the two phases and to study the gas/particle evolution. To this end numerical outputs of our code are compared with known analytical solutions.

First, we test the capacity of the code to keep a gas disk in steady rotation around the star. In other words our questions are: with what duration and with what precision does the code satisfy the initial steady state solution for the gas? In this respect it is necessary to study the numerical error as a function of grid resolution, that is, the number of cells ( $n_r, n_\theta$) used in a given simulation. Of course, the higher the resolution the smaller the numerical errors but, also, the higher the cost in terms of computing time and memory. It is therefore valuable to optimize the resolution in simulating a gas disk with adequate accuracy. A second test is to check that the code appropriately describes the evolution of solid particles embedded in a steady flow of gas and that their drift rate is consistent with known analytical solutions. As Fig. 1 indicates, the particle drift velocity depends on the radial distance from the star, so that particles migrating from the outer regions of the disk increase the density of solids in the inner regions. The last test is to check the accuracy of the code in the case of a protoplanetary disk in which the gas and the solid components are coupled.

4.1 Steady state of the gas disk

The steady state solution of a gas disk flowing around a central star is given by Eqs. (12) and (17), with inner and outer boundaries chosen at 5 and 10 AU, respectively. Simulations are performed over 1000 yr and for three different resolutions of the grid: (nr, $n_\theta$)=(100, 100), (200, 100), and (200, 200), where nr and $n_\theta$ are the numbers of cells in the r and $\theta $ directions, respectively. Initially, a steady state is assumed; numerical errors are 0 but increase with the number of computing steps. The time scale of the mechanisms is of the order of a few dynamical time scales (i.e., much shorter than the 1000 yr in our computational domain). Figure 4 shows the deviations of the numerical results from the initial steady state after an integration of 1000 yr. For example, the deviation of the tangential velocity from the initial steady state is given by $\Delta v_{\mbox{\scriptsize g}} = v_{\mbox{\scriptsize g}} - v_{\mbox{\scriptsize g,0}}$. The largest numerical errors are not found at a fixed position in the disk and can be significantly reduced by increasing the resolution in the r-direction. On the other hand, the numerical errors are nearly independent of the resolution in the $\theta $-direction; this likely comes from the axisymmetrical nature of our initial conditions.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{1085f4.eps}\end{figure} Figure 4: Deviation of the numerical results from the steady state solution of the gas for various resolutions: (100, 100) (solid lines), (200, 100) (dotted lines) and (200, 200) (dashed lines).
Open with DEXTER

To describe the coupled evolution of the two phases requires an accurate calculation with the numerical code of the velocity difference between the gas and solid components as given by Eqs. (7) and (11). The magnitude of the relative velocities are plotted in Fig. 2 as a function of the distance from the star and for various radii. It indicates the level of accuracy necessary to compute the drag forces and to build up a two-phase code for protoplanetary disks. From our experience, the numerical error of the tangential velocity of the gas is larger than that of the radial velocity whereas the tangential component of the relative velocity is smaller than its radial component (see Fig. 2). In our computational domain, the velocity of a particle relative to the gas is found to increase with the distance from the star (see Fig. 2) and comparisons are carried out with the smallest and largest values of the relative velocity at the inner and outer boundaries, respectively 5 and 10 AU. We also studied the importance of the grid resolution on the numerical outputs of our code. In Fig. 5 the numerical errors on the tangential velocity are compared to the expected relative velocity for various grid resolutions (nr=50, 100, 200, and 400). In the $\theta $-direction the grid solution is less important than in the radial direction (see the deviation from the steady state solution in Fig. 4.) and was set as $n_\theta =100$. At the lowest values (nr=50 and 100) the accuracy is not sufficient to numerically compute the drag forces but the numerical error decreases with increasing resolution. A grid resolution of ( $n_r, n_\theta$) = (200,100) seems to be a good compromise in building up a 2-phase code. The numerical error is close to a power law of the radial resolution with an index -2 (the dotted line in Fig. 5), consistent with the second order accuracy of the code. An application with a higher resolution was performed in Sect. 6.


  \begin{figure} \par\includegraphics[width=7cm,clip]{1085f5.eps}\end{figure} Figure 5: Dependence of the numerical error on the radial resolution. The resolution in the $\theta $-direction is fixed to be $n_\theta =100$. The dotted line set up the slope -2.
Open with DEXTER

4.2 Radial drift of the solid particles

Now, we have to check whether the code can correctly handle the evolution of the solid component. To this end we first assume that the gas always stays in a steady flow around the star (given by Eqs. (12) and (17)) and that the particles have no back reaction onto the gas. In this case the particles are submitted to a continuous gas drag, which makes them drift radially towards the star, and the results of the numerical simulations can be compared with known analytical solutions. This is, of course, an oversimplified situation when the density of solid particles is close to that of the gas, $\sigma_{\mbox{\scriptsize d}} = \sigma_{\mbox{\scriptsize g}}$. In the next subsection, it will be used as a convenient test case.

The analytical expressions for the time evolution of the surface density of the solid component were derived from the paper of Nakagawa et al. (1986). At distances from the star which range from 5 to 10 AU, the drag force is given by the Epstein law for a solid particle with radius of 1 cm and Eq. (16) becomes the approximate expression of the migration velocity given by:

 \begin{displaymath}\displaystyle u_{\mbox{\scriptsize d}} = - 18.4 \left( \frac{r}{\mbox{1~AU}} \right)^{3/2} \mbox{cm/s}. \end{displaymath} (22)

Centimeter-size particles have very small relative velocity with respect to the gas (much smaller than bigger particles with radii of the order of 10 cm or 100 cm) so that computing drag forces with a good accuracy may become tricky. Computations remain consistent if the gas is assumed to flow in a laminar way around the star, as assumed in the present paper. The migration velocity of the particles is given by a simple power law and it is possible to derive an analytical expression for the surface density of the solid component at any time t by solving a Cauchy problem (Youdin & Shu 2002):

 \begin{displaymath}\displaystyle \sigma_{\mbox{\scriptsize d}} = \sigma_{\mbox{\scriptsize d,0}} \left( \frac{r_0}{r} \right)^{5/2}, \end{displaymath} (23)

where $\sigma_{\mbox{\scriptsize d,0}}$ is the initial surface density of the solid component and r0 is the distance that separates the particles from the star at time t=0, given by:

 \begin{displaymath}\displaystyle r_0 = r \left( 1 + 0.5 \frac{u_{\mbox{\scriptsize d}} t}{r} \right)^{-2}\cdot \end{displaymath} (24)

(r0 is always larger than r since the drift velocity is negative). The evolution of the surface density of the solid component as a function of time has been simulated and is sketched in Fig. 6 by snapshots of the evolution at t=0, 500 yr and 1000 yr. This figure indicates that even a low resolution grid, $(n_r,n_\theta)=(50,50)$, is sufficient to correctly fit the analytical solution. Solid particles migrating from the outer region of the disk increase the surface density of the solid component in the inner disk regions. Our numerical results are in very good agreement with the analytical solutions on the whole domain studied in this paper and even after 1000 years of evolution.


  \begin{figure} \par\includegraphics[width=7.3cm,clip]{1085f6.eps}\end{figure} Figure 6: Surface density of the solid component at t=500 (panel a) and 1000 years (panel b). The results of our numerical simulations (open circles) are compared to the analytical solution (solid lines). The initial surface density of the solid component is given by the dashed lines. In this plot the results are obtained with a grid resolution of (50, 50). The surface density increases with time.
Open with DEXTER

In Fig. 7, the numerical error on the computed surface density of the solid component has been plotted for three different values of the resolution $(n_r,n_\theta)=(15,15)$, (20, 20), and (50, 50). In the outer region the numerical errors are similar for three resolutions. In the inner regions the largest errors occur for the low resolutions. However, even at the boundaries, the errors are much smaller than the absolute values of the surface density. Any resolution used in the present version of the code seems sufficient to describe the migration of the solid component.


  \begin{figure} \par\includegraphics[width=7.3cm,clip]{1085f7.eps}\end{figure} Figure 7: Numerical errors for the surface density of the solid component. The errors correspond to simulations made with three differentresolutions: (50, 50) (solid lines), (20, 20) (dotted lines), and (15, 15) (dashed lines). Errors are large near the boundaries of the disk.
Open with DEXTER

4.3 Simultaneous evolution of the gas and solid components

In this subsection, the gas and solid components are coupled and can evolve as a function of time. The solid particles drift towards the central star, losing energy and angular momentum which is transfered to the gas moving outward. The analytical solution derived by Nakagawa et al. (1986) provides approximate expressions of the velocities of the two phases, if the distributions of gas and solids are slowly evolving in a disk rotation. Initial conditions were chosen similar to that used in the previous subsection, with $\sigma_{\mbox{\scriptsize d}} = \sigma_{\mbox{\scriptsize g}}$; this is in order to know better the effect of the exchange of energy and angular momentum between the two phases. In this case, which corresponds to choosing $\alpha = 1$ in Eq. (16), the two components are coupled and the drift velocity of the solid particles is found to be four times smaller than in the case with no coupling ($\alpha =0$) discussed above. Using this reduced drift velocity instead of Eq. (22), it is possible to analytically estimate the surface density of the solid component at 1000 years. A comparison is made in Fig. 8 of this approximate solution with the solid surface density arising from the complete numerical simulation. In spite of differences near the boundaries of the computational domain the semi-analytical calculations quite well agree with the numerical computations.


  \begin{figure} \par\includegraphics[width=7.1cm,clip]{1085f8.eps}\end{figure} Figure 8: Surface density of the solid component at t=1000 years resulting from the evolution of a two-phase disk, when gas and solid are fully coupled. a) The numerical results (dashed lines) compared to the approximate solutions(solid lines). The dotted line refers to the initial surface density. The surface density does not increase with time due to the strong coupling. b) Surface density deviations: difference between the numerical results and the approximate solution (solid line), difference between the numerical results and the initial distribution (dashed line). The resolution of the grid is (200, 100).
Open with DEXTER

5 Testing the code accuracy in two dimensions

Now, the code is tested on a full 2D situation looking at the evolution of non-axisymmetric initial conditions. For this test we again find some of the results obtained by Li et al. (2001) on the Rossby wave instability (RWI) and make comparisons with solutions of the linear analysis (Li et al. 2000). Unfortunately, to our knowledge, there is no solution for non-axisymmetric problems with two phases of fluid using linear analysis; therefore, we concentrate on tests with a single fluid.

Li et al. (2000) performed a linear analysis of the RWI and showed that a disk with an annular overpressure (a few times the scale height) can become unstable to non-axisymmetric perturbations. They consider a disk with a Gaussian shaped bump given by:

 \begin{displaymath}\displaystyle \sigma_{\mbox{\scriptsize g}} = 1 + (A-1) \exp \left[ -\frac12 \left( \frac{r-1}{\Delta r} \right)^2 \right], \end{displaymath} (25)

where A and $\Delta r$ are the amplitude and the width of the bump, respectively. Starting from the above disk, they found the growth rate of the perturbation and the corresponding eigenfunctions of the unstable mode. For the test we use the same initial conditions and eigenfunctions as those of Li et al. (2000) and perform a simulation with the resolution (nr, $n_\theta$) = (200, 100). The evolution of the radial velocities at r=1 and 0.7 obtained with the numerical simulations is compared to evolution corresponding to the growth rate of perturbation derived by the linear analysis (see Fig. 9). The numerical results are found in good agreement with the solution of the linear analysis.


  \begin{figure} \par\includegraphics[width=8cm,clip]{1085f9.eps}\end{figure} Figure 9: Evolution of the radial velocities at r=1 and 0.7, and $\theta =0$ in the linear growth stage. Numerical results (dashed lines) are compared with solutions issued from the linear analysis (solid lines). The amplitude and width of the pressure bump are 1.35 and 0.05, respectively. The azimuthal mode number of the perturbation is 3. This case corresponds to the T5 case in Li et al. (2001).
Open with DEXTER

Li et al. (2001) also studied the nonlinear evolution of the disk and showed that vortices form in the region initially occupied by the pressure bump. We succeed in simulating the same nonlinear evolution, reproducing the formation of vortices. Figure 10 shows snapshots of the pressure at 2.7 orbital revolution obtained from the numerical simulations with nr = 200 and $n_\theta =100$ and 200. A chain of three vortices forms because the initial perturbation has an azimuthal mode number of 3. The overall structure of the vortices is the same as in Li's simulations although the structure is a bit different at smaller scale, for example at the edges of the vortices. The number of meshes in the $\theta $-direction is necessary to display the small scale structure. In contrast the overall structure was unchanged when increasing the azimuthal ($\theta $) direction. These results clearly confirm that our code is relevant for two dimensional problems and the study of disk evolution.

6 Application to the lifetime problem

The lifetime of the solid particles in protoplanetary disks is generally estimated neglecting the back reaction of the particles onto the gas. Under this assumption the drift velocity in the meter-size range is so large that particles are lost into the star in only a few thousand years, a much shorter time scale than the few million years for gas survival around stars. So, how are planetesimals and planets formed if growing particles are progressively swept into the star?

  \begin{figure} \par\includegraphics[width=8cm]{1085f10-CMJN.eps}\end{figure} Figure 10: Snapshots of the pressure in units of 10-3 after 2.7 orbital revolutions at r=1. Numerical results for two different azimuthal resolutions: $n_\theta =100$ (panel a) and 200 (panel b). The number of meshes in the radial direction is set up to 200. The inner and outer boundaries of the disk are 0.4 and 2, respectively.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{1085f11.eps}\end{figure} Figure 11: Evolution of a two-phase disk when perturbed by a Gaussian bump in the distribution of the solid particles. The dotted lines correspond to the initial conditions. The solid lines correspond to the the surface densities and the radial velocities (gas and solids) after a 100 yr evolution. Note that the surface density of the solid component is plotted on a log scale. The resolution of the grid is (400, 200).
Open with DEXTER

We are exploring the possibility that the coupling of the two phases, with the exchange of energy and angular momentum, could help to solve the paradox. Indeed, an important result of Sect. 4 is that "active'' solid particles migrate inward toward the star more slowly than "passive'' particles; the angular momentum is transfered to the gas which migrates outward. The coupling of the two phases can change the gas pressure and, in relation, the spatial distribution of the gas component. Consequently, the radial pressure gradient, the parameter $\eta$ (Eq. (15)) and the particle drift velocity are also modified. For example, a vanishing pressure gradient is able to stop particle drift and gas migration. Such a situation could be the way to confine the solid particles in the gas disk longer. In a protoplanetary disk, this occurs in the case of strong coupling (small particles) and high particle densities.

In this section we explore this idea, looking at the evolution of a peculiar distribution of the solid particles. The initial condition is an annular density enhancement inside a gas disk in nearly Keplerian rotation around the star. This initial density distribution of the solid particles could be due to a migration mechanism similar to that investigated by Youdin & Shu (2002). The above density distribution mimics the one obtained by these two authors. In this case the time evolution of the spatial distribution of the gas and solid phases cannot be found analytically but we can rely on the numerical code developed and tested in the previous sections. For the numerical computations we adopted a disk with the same size as in Sect. 4 (inner radius at 5 AU and outer radius at 10 AU) but used a higher resolution $(n_r,n_\theta)=(400,200)$ in order to be confident of computing the drag force with high accuracy. The initial distribution of the gas component is assumed the same as in Sect. 4, given by Eqs. (12) and (17). Further, the over-density of solid particles is assumed to be in a ring centered at 7.5 AU from the star with an initial distribution expressed as a sum of two functions: (i) a power law surface density 100 times smaller than the standard gas density; and (ii) a Gaussian bump given by:

 \begin{displaymath}\displaystyle \sigma_{\mbox{\scriptsize d}} = \sigma_{\mbox{\... ...mes 10^{-2} + \exp \left( -9 \times (r-7.5)^2 \right) \right). \end{displaymath} (26)

The solid component is assumed to be in the form of decimeter-size particles whose "stopping'' time is of the order of the orbital period and permits exchange of energy and momentum in a reasonable time scale. The time evolution of the two components, both gas and solid, was computed with our code over a duration of 100 yr (approximately 10 rotations).

Figure 11 shows in some snapshots the time evolution obtained for the surface densities and the radial velocities of the two phases. Radial motions are very different following the value of the solid to gas density ratio (nearly equals to 1 inside the bump and to 0.01 outside). At t=0, the outward migration of the gas is very small except in the bump while the inward drift of the particles is large (increasing outward) except in the bump where it is four times smaller. Such differences come from the exchanges of energy and angular momentum between the two phases and were computed using the analytical formulation of Sect. 2.2 (Eq. (16) in the case $\alpha = 1$). At t=100 yr, the evolution mainly concerns the bump region in which the gas surface density become flatter. Consequently, due to the equation of state Eq. (3), the radial pressure gradient is weaker and the particle drift velocity becomes four times smaller than that of the initial condition, as shown in Fig. 11.

The time evolution of a two-phase disk is, then, significantly different from a single-phase one in which the solid particles are only considered as a passive component with no back reaction onto the gas. Particularly, when the gas disk contains an over-dense bump of the solid component, the two-phase approach shows that the inward drift of particles in the bump is reduced by more than an order of magnitude. This indicates that the exchange of energy and angular momentum can play an important role in the evolution of protoplanetary disks by reducing the mass loss rate of particles in high-concentration regions. This mechanism opens some interesting perspectives but cannot solve the lifetime problem encountered in the formation of the planetesimals.

7 Conclusions

A high accuracy hydrodynamical code has been developed and tested to simulate the evolution of a protoplanetary disk with a two-phase approach. The code accuracy has been tested in three different cases: (i) a steady state single phase disk; (ii) a two-phase disk in which the back reaction of the particles onto the gas is neglected; (iii) a two-phase disk in which the gas and the particles are fully coupled to one another. In the second case the total energy and angular momentum in the system are not conserved but the test is relevant since analytical solutions are known. We checked that our code has the required accuracy to describe the evolution of a two-phase disk. A reasonable grid resolution was found to be 200 cells in the r-direction and 100 cells in the $\theta $-direction.

The disk evolves through the transfer of energy and angular momentum from the particles to the gas. If the back reaction of the gas is completely neglected, solid particles lose energy and angular momentum and drift inward toward the star on a time scale of the order of 103 yr. If the back reaction of the gas is taken into account in a semi-analytical way, under a quasi-steady state assumption, the residence time of the particles in the disk is increased by a factor 4. Finally if the coupling of the two phases is fully taken into account, computations are completely self consistent and our simulations show that particle lifetime can also depend on the solid phase concentration. It is found that in over-dense annular regions, particles can survive 16 times longer than in the uncoupled case.

The present work shows that the lifetime of solid particles in a protoplanetary disk must be determined using a self-consistent evolution of the disk in which the gaseous and the solid phases are coupled. It opens some interesting perspectives but does not presently allow a solution of the lifetime problem in the continuous growth scenario for the formation of planetesimals.

Another way to solve the problem is to trap the particles in anti-cyclonic gaseous structures before they fall into the star (Barge & Sommeria 1995). Johansen et al. (2004) performed hydrodynamical simulations of compressible coupled solid particles-gas equations under a shearing box approximation. They showed that an anticyclone survives for at least one orbital rotation and can accumulate solid particles with a density enhancement of 17. This possibility will be explored in a second paper. Of course, to correctly address the problem would also require one to describe how solid particles sink toward the equatorial plane of the disk. This is planned in a future work which will include the third dimension in the problem with the goal to develop a 3D two-phase code to study protoplanetary disks.

Acknowledgements
We express our gratitude to H. Li for giving us the data of the eigenfunctions for Rossby wave instability. S. Inaba received financial support from Research Fellowships of the Japan Society for Promotion of Science for Young Scientists (No. 06423). This work has received support from INSU under grant AS-N°22531 (Oct. 2000: "Formation Planétaire et Milieux Diphasiques'').

References

 
Copyright ESO 2005