A&A 378, 180-191 (2001)
DOI: 10.1051/0004-6361:20011183

Diversity of planetary systems from evolution of solids in protoplanetary disks

K. Kornet1 - T. F. Stepinski2 - M. Rózyczka1


1 - Nicolaus Copernicus Astronomical Center , Bartycka 18 , Warsaw, 00-716, Poland
2 - Lunar and Planetary Institute, 3600 Bay Area Blvd., Houston, TX 77058, USA

Received 1 June 2001 / Accepted 6 August 2001

Abstract
We have developed and applied a model designed to track simultaneously the evolution of gas and solids in protoplanetary disks from an early stage, when all solids are in the dust form, to the stage when most solids are in the form of a planetesimal swarm. The model is computationally efficient and allows for a global, comprehensive approach to the evolution of solid particles due to gas-solid coupling, coagulation, sedimentation, and evaporation/condensation. The co-evolution of gas and solids is calculated for 107 yr for several evolution regimes and starting from a comprehensive domain of initial conditions. The output of a single evolutionary run is a spatial distribution of mass locked in a planetesimal swarm. Because swarm's mass distribution is related to the architecture of a nascent planetary system, diversity of swarms is taken as a proxy for a diversity of planetary systems. We have found that disks with low values of specific angular momentum are bled out of solids and do not form planetary systems. Disks with high and intermediate values of specific angular momentum form diverse planetary systems. Solar-like planetary systems form from disks with initial masses $\le$0.02 $M_{\odot}$ and angular momenta $\le$ $3\times 10^{52}$ gcm2s-1. Planets more massive than Jupiter can form at locations as close as  1 AU from the central star according to our model.

Key words: accretion disks - solar system: formation


1 Introduction

The architecture of a planetary system results from a chain of processes that start at the formation of a protoplanetary disk (PPD) around a nascent star. A protoplanetary disk is a mixture of gas and solids and is made up from a material that, during the formation process, did not become directly part of the star, but, due to its excess angular momentum, remained in orbit around the newborn star. Global properties of protoplanetary disks, such as their masses and sizes are known in abundance from astronomical observations. Typical disk masses are between 0.001 and 0.1 $M_{\odot}$ and typical disk sizes are between a few and a few hundred AU (see Beckwith & Sargent 1993; Beckwith 1994 and references therein). The ranges of about two orders of magnitude in observed masses and sizes of protoplanetary disks are partially due to different stages at which different disks are observed (Stepinski 1998a) and partially due to an intrinsic scatter in initial conditions. On the other hand, the structure of a planetary system is known only for a single object - the Solar System. High accuracy spectroscopy has revealed the existence of over 60 solar-type stars with unresolved companions that could plausibly be planets (cfa-www.harvard.edu/planets/). If these companions are indeed planets formed in the same way as planets in the Solar System, a broad distribution of their orbital parameters suggests a large diversity in possible configuration of planetary systems. We postulate that indeed the dispersion in protoplanetary disks initial conditions leads to architecture diversity among planetary systems.

A protoplanetary disk is an evolving object; its properties such as mass and size are functions of time. The evolutions of gaseous and solids components of a protoplanetary disk accompany each other but are not identical. One important distinction between the evolution of the two components is that the gas maintains its form but the solids do not. Initially the solids are all in the form of a fine dust, but, given enough time, they convert into roughly 1-10 km sized solid bodies (planetesimals). It is currently thought (for a review see Weidenschilling & Cuzzi 1993) that this is achieved via a buildup of progressively more massive particles by the process of coagulation culminating in planetesimals. In general, the mass distribution of the solids evolves due to gas-solid coupling, coagulation, sedimentation, and evaporation/condensation. Once planetesimals are formed, their further aggregation leads to planets. This final accumulation stage received a fair amount of attention (Safronov 1968; Greenberg et al. 1978; Kaula 1979; Wetherill 1980; Nakagawa et al. 1983) because its physics is fairly well understood, inasmuch as it can be formulated as the N-body problem with gravitational forces. The earlier accumulation stages, from dust to planetesimals, are governed by a complex combination of the processes mentioned above. Although these processes have been discussed and modeled individually (for a review see Weidenschilling & Cuzzi 1993), little work has been done to model the evolution of the solids component of the protoplanetary disk associated with their collective action.

Cassen (1996) constructed a simple model of disk evolution coupled to an equally simple prescription for the rate at which condensible material is decoupled from the gas to become a part of a surviving planetary system. In his model evaporation/condensation is the only process affecting the coupling of solids to the gas (coagulation, sedimentation and aerodynamic drag are not considered). Cassen's model was designed to reproduce the observed abundances of elements in chondritic meteorites and not to address the issue of diversity of planetary systems. Schmitt et al. (1997) constructed a model of a protoplanetary disk in which dust evolution is directly coupled to the evolution of the gas. Their model includes all relevant processes without any further simplification. However it was designed to study the effects of coagulation on gas opacity, and the computation followed the dust evolution for only 100 years. Stepinski & Valageas (1996, 1997) (hereafter referred to as SV96 and SV97) developed a model capable of simultaneously following the evolution of gas and solid particles due to gas-solid coupling, coagulation, sedimentation and evaporation/condensation for up to 107 years. Their model was designed specifically to obtain the radial distribution of solid material circumnavigating a star in the form of the planetesimal swarm. Despite far-reaching simplifications that made those calculations possible, the model required significant computational resources and only a handful of evolutionary runs were calculated.

Our current work is based on the ideas developed in SV96 and SV97. Our first objective is to develop and utilize a computationally efficient model that tracks the evolution of solids from an early stage, when they are in the dust form, to the stage when most solids are in the form of planetesimals. In order to achieve a necessary computational efficiency, the model developed here introduces further simplifications on top of those used by SV97. Our second objective is to use the model to establish the dependence between a distribution of matter at the early evolutionary stage of a PPD (initial conditions), and a distribution of matter at the final evolutionary stage (the planetesimal swarm). Our model is not designed to track the further evolution of planetesimals into planets, thus, in this paper, the planetesimal swarm is a surrogate for a planetary system. The computational efficiency of the model is necessary in order to evolve an extensive set of initial conditions into an extensive set of swarms.

Section 2 contains the description of our model. In Sect. 3 we explain the numerical methods employed to implement the model. The results are presented in Sect. 4, and the fine details of how exactly gas and solid particles evolve in a particular disk for 107 yr are discussed in Sect. 5. Finally, in Sect. 6 we present conclusions and discussion.

2 Model of a protoplanetary disk

We model a protoplanetary disk as a two-component turbulent fluid. The dynamically dominant component is a gaseous protoplanetary disk (GPD). The collection of all solids in the GPD constitutes the second component, which we call solid protoplanetary disk (SPD). SPDs evolves together with, but not identically to GPDs. The evolution of each component is governed by equations of continuity and momentum conservation. Both, the GPD and the SPD are considered to be perfect fluids. The SPD is also a pressureless fluid, subject to an additional force arising from the friction between solid particles and the gas. This force is neglected in calculating the evolution of the GPD. Initial conditions are imposed at the beginning of what is frequently referred to as a dissipative stage (Cameron 1985) of the protoplanetary disk, when the GPD evolves due to an anomalous viscosity caused by (unspecified) disk instability. It is during this stage that the major part of the transition from dust to planetesimals occurs. The disk is assumed to maintain an axial symmetry ( $\partial/\partial \phi =0$) and to be geometrically thin ($z/r\ll1$). We analyze the disk in the cylindrical polar coordinates $(r,\ \phi,\ z)$ with the star at the origin and the disk's midplane located in the plane z=0.

2.1 Evolution of the GPD

The gas is explicitly assumed to be turbulent. The precise nature of the turbulence is not identified, but it is expected to be a mechanism intrinsic to the GPD such as shear-induced instability (Dubrulle 1993), magnetorotational instability (Blaes & Balbus 1994) or thermal convection (Lin & Papaloizou 1980), and not a process induced by the presence of solid particles (Weidenschilling 1980). Thus, in our model, the evolution of the GPD is not influenced by the presence of the SPD.

The average velocity field ${\vec V}(t,r)=\{V_{r},V_{\phi}, V_{z}\}$ and the average surface density, $\Sigma$ of the GPD can be calculated under a particular model of velocity correlations (Reynolds stress). The details of Reynolds averaging in this context are given in SV97, here it is sufficient to note that the turbulent viscosity is modeled by

 \begin{displaymath}\nu_{\rm t} = {1 \over 3} \alpha {C_{\rm s}} H \end{displaymath} (1)

where $\alpha<1$ is a free parameter introduced by Shakura and Sunyaev (1978), $C_{\rm s}$ is the speed of sound, and H is disk's half-thickness. For a given disk model $\alpha $ is assumed to be constant and uniform.

The average surface density of the GPD, $\Sigma$, is given by the following evolution equation,

 \begin{displaymath}\frac{\partial \Sigma}{\partial t} - \frac{3}{r} \frac{\parti... ...c{\partial}{\partial r} (r^{1/2} \nu_{\rm t} \Sigma)\right]=0. \end{displaymath} (2)

It can be shown (see, for example, Frank et al. 1985) that $\nu_{\rm t}$ is not an explicit function of time, but instead depends only on the local disk's quantities and can be expressed as $\nu_{\rm t}=\nu_{\rm t}(\Sigma,r)$. The form of this function depends on the opacity law, in our calculations we used the Rossenland mean opacity given by Ruden & Pollack (1991). Thus, (2) is closed and can be solved given an initial condition.

We obtain solutions to the evolution Eq. (2) analytically, following the method developed by Stepinski (1998b). In Stepinski's method initial conditions (the functions $\Sigma(t=0, r)$) are restricted to those having a self-similar form parameterized by the following three parameters: j0, an angular momentum of the GPD at t=0 (in units of $10^{52}\, {\rm g}\ {\rm cm}^2\ {\rm s}^{-1}$), m0, the mass of the GPD at t=0 (in units of $M_{\odot}$), and $m_{\star}$, the mass of the central star (in units of $M_{\odot}$). It is assumed that $m_{\star}$ stays constant during the evolution despite the fact that there is accretion from the disk onto the star. The initial outer radius of the GPD, r0 (in units of AU) is given by

 \begin{displaymath}r_0=0.021\ m^{-2}_{0}\ j^2_{0}\ m^{-1}_{\star}. \end{displaymath} (3)

Thus, if we consider only disks around stars with $m_{\star}=1$, the initial conditions are parameterized conveniently by only two numbers, j0 and m0. The analytical method yields an explicit formula $\Sigma (t, r; \alpha, j_{\rm0}, m_{\rm0})$. Other GPD quantities, including the average velocity field ${\vec V}$ and the midplane temperature T are also given by explicit formulas in the same format.

2.2 Evolution of the SPD

In our model the SPD is a collection of solid particles of different sizes entrained in the GPD. The crucial approximation is that the size distribution of particles at any given radial location of a disk is narrowly peaked about a mean value particular for this location and time instant. Such an approximation was first proposed by Morfill (1985) who noted that the resulting particle growth rate is not sensitively dependent on the size distribution. The practical implementation of this assumption means that a particle size, a, can be expressed as a single value function of time and position, a=a(t,r). Collectively, all particles forming the SPD can be regarded as a turbulent fluid characterized by the average surface density $\Sigma_{\rm d}(t,r)$ and the average velocity field ${\vec V}_{\rm d}(t,r)=\{V_{{\rm d}r},V_{{\rm d}\phi},V_{{\rm d}z}\}$. The evolution of the SPD is governed by the set of two equations.

 \begin{displaymath}\frac{\partial \Sigma_{\rm d}}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}(r V_{{\rm d}r} \Sigma_{\rm d})=0 \end{displaymath} (4)

and

 \begin{displaymath}\frac{\partial \Sigma_{a}}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}(r V_{{\rm d}r} \Sigma_{a}) = f \Sigma_{\rm d} \end{displaymath} (5)

where $\Sigma_{a}=a \Sigma_{\rm d}$, and f is the source function for particle growth by coagulation (see Sect. 2.4). The first of these two equations is simply a vertically integrated mass continuity equation for solid particles. The distribution of mass in the SPD is driven by the radial component of its average velocity field that depends on frictional coupling between particles and the gas. This coupling depends on local particle size and local gas properties. The particle size changes due to coagulation, whose rate in turn depends on the density of the SPD. Thus (4) is coupled to both the gas evolution and the particle size evolution.

Equation (5) represents the particle size evolution and encapsulates the coagulation process under assumptions of our model. It is convenient to convert the usual coagulation equation ${\rm d} a/{\rm d} t = f$ (see Sect. 2.4) into the conservative form (5) that can be viewed as the continuity equation for "size-weighted'' surface density $\Sigma_{a}$. The particle size evolution Eq. (5) is coupled to the evolution of the gas via $V_{{\rm d}r}$ and to the evolution of the SPD via the source term and the definition of $\Sigma_{a}$.

The evolution of the SPD is determined by solving the set of Eqs. (4)-(5) starting from prescribed initial conditions. In order to close the set, the average radial velocity, ${\vec V}_{{\rm d}r}$, and the source function for particle coagulation, f, must be expressed in terms of gas parameters, particle size, and $\Sigma _{\rm d}$(4)-(5).

2.3 The average radial velocity of the SPD

Because of the coupling to the turbulent GPD, the velocity field of the SPD consists of average and fluctuating parts. The average velocity field, ${\vec V}_{\rm d}$ can be calculated under a particular Reynolds stress model for the SPD. The details of such calculations are given in SV97. The SV97 method of obtaining the components of ${\vec V}_{\rm d}$ leads to the system of equations that is nonlinear and not separable, and thus needs to be solved numerically. Moreover, in general, the system depends on particle densities, requiring a simultaneous solution of velocities and particles densities. This high computational load makes this approach impractical for our current purpose. Therefore, in this paper we calculate ${\vec V}_{\rm d}$ in a simpler, more approximate fashion that leads to an explicit formula for ${\vec V}_{{\rm d}r}$. For the purpose of calculating ${\vec V}_{\rm d}$ we assume that solid particles are coupled to the average flow of the gas. Under such an assumption and because of the previous assumption about particle sizes (Sect. 2.2), the velocity field of the SPD is tantamount to velocities of specified individual particles. The radial component of ${\vec V}_{\rm d}$ can be expressed as

 \begin{displaymath}V_{{\rm d}r}= V_{\rm k} { {\left [\Delta \Lambda + {V_r\over V_{\rm k}}\right ]}\over {1+\Lambda^2} }\cdot \end{displaymath} (6)

In (6) $V_{\rm k}=\Omega_{\rm k} r$ is the Keplerian velocity, $\Delta=(1/\rho)(\partial P/\partial r)(r/V_{\rm k}^2)$ is the ratio of the residual gravity to the central gravity, and $\Lambda=\Omega_{\rm k} t_{\rm s}$, where $t_{\rm s}$ is the so-called stopping time. The residual gravity is the name given to the support against the central gravity provided by the pressure gradient in the gaseous disk. The stopping time is a convenient measure of the strength of the frictional force. The value of $t_{\rm s}$ depends on particle size, the density of the gas and the speed of the particle relative to the gas, v. The dimensionless quantity $\Lambda$ determines efficiency of gas-solids coupling. If $\Lambda \rightarrow 0$, the stopping time is small compared to the orbital period, and solids are strongly coupled to the gas. If $\Lambda \rightarrow \infty$, the stopping time is very long compared to the orbital period and solids are decoupled from the gas.

Although (6) is written in the format suggesting that  $V_{{\rm d}r}$ is given explicitly, this is not necessarily the case because the quantity $\Lambda$ depends, in general, on the relative velocity v and thus on $V_{{\rm d}r}$. Recall that the stopping time is defined as $t_{\rm s}=m_{\rm p} v/F_{\rm D}$, where $m_{\rm p}$ is a mass of a solid particle and $F_{\rm D}$ is the drag force. The drag law takes different forms in different physical regimes. SV96 gives a complete description of $F_{\rm D}$ in all regimes. Here it is sufficient to note that when the mean free path of gas molecules is larger than the size of the particle, $F_{\rm D}$ is given by the Epstein law, otherwise it is given by the Stokes law. Furthermore, the Stokes law takes three different forms depending on the value of the Reynolds number. Taking these laws into consideration we find out that $\Lambda$ does not depend on v in the Epstein regime and in the first Stokes regime. In the second Stokes regime $\Lambda \sim v^{-0.4}$, and in the third Stokes regime $\Lambda \sim v^{-1}$. These two regimes apply when solids are large. Relative velocity between a large solid and the gas is dominated by its transverse component which is of the order of $\Delta V_{\rm k} \sim C_{\rm s}^2/V_{\rm k}$. For the sake of simplicity we assume that in these regimes $v=\epsilon \ C_{\rm s}^2/V_{\rm k}$, where $\epsilon$ is a free constant parameter ( $\epsilon=0.5$ is used in our calculations). With such an approximation the formula (6) becomes indeed explicit and the radial velocity of a solid particle and thus the radial component of ${V}_{\rm d}$ is given in terms of particle size and the gas properties.

The presence of two terms in the numerator of (6) reflects the existence of two mechanisms driving the radial motion of solid particles. The term $\Delta \Lambda$ corresponds to engendering radial motion of particles by virtue of the difference between tangential velocities of the gas and the particles and its magnitude is of the order of $\Lambda C_{\rm s}^2/V_{\rm k}^2$. The term $V_r/V_{\rm k}$corresponds to radial motion of particles induced by the radial flow of the gas, and its magnitude is of the order of $\alpha C_{\rm s}^2/V_{\rm k}^2$. The relative importance of these two terms is $\Lambda/\alpha$. Thus, only small particles (those that are strongly coupled to the gas) are carried radially by the gas accretion flow, whereas larger particles moves radially because their transverse (orbital) velocities are decreased by an aerodynamic drag.

The most important feature of the radial velocity field of the SPD is its dependence on $\Lambda$. Within our approximations, (6) is an explicit formula, and the character of the function $V_{{\rm d}r}(\Lambda)$ can easily be deduced. The $V_{{\rm d}r}(\Lambda)$ has a maximum at $\Lambda \sim 1$. Thus, particles with $t_{\rm s} \sim \Omega_{\rm k}^{-1}$ acquire the highest inward velocity, $V_{{\rm d}r}^{\rm max}(t,r)$. Because the stopping time depends on both particle size and gas properties, that velocity is not reached by particles of only one, well defined size. Instead, particles of different sizes reach $V_{{\rm d}r}^{\rm max}$depending on location in the disk and the time. From (6) it is clear that very small particles, those with $\Lambda \rightarrow 0$, move radially with $V_{{\rm d}r}=V_r$, and very large particles, those with $\Lambda \gg 1$, have $V_{{\rm d}r}=0$.

2.4 Coagulation

We assume that solid particles are spheres with bulk density  $\rho_{\rm s}$ and mass $m_{\rm d}$ each, where $m_{\rm d}(t,r)=(4/3)\pi a^3(t,r) \rho_{\rm s}$. Because of our assumption about the size distribution of particles (Sect. 2.2) only particles of the same size are present at a particular location, at any given time. The density of matter concentrated into solid particles is $\rho_{\rm d}(r)$ and the particles' number density is $n_{\rm d}(r)$. The geometrical cross section for collision between two such particles is $\sigma=\pi {(a+a)}^2$ and the mean time between collisions is $\tau_{\rm coll}=1/(n_{\rm d} \sigma V_{\rm rel})$ where $V_{\rm rel}$ is the mean relative speed between particles. Assuming that particles always stick to each other upon collision, the growth of particle mass $m_{\rm d}$ in unit time can be expressed as follow

 \begin{displaymath}\frac{{\rm d}m_{\rm d}}{{\rm d}t}\approx \frac{m_{\rm d}}{\ta... ...}\ n_{\rm d}\ m_{\rm d} = 4\pi a^2\ V_{\rm rel}\ \rho_{\rm d}. \end{displaymath} (7)

The magnitude of $V_{\rm rel}$ has been derived in SV97 to be

 \begin{displaymath}V_{\rm rel}^2 = 2^{3/2} R_{\rm o} \alpha C_{\rm s}^2 \frac{{Sc}-1}{{Sc}^2} \end{displaymath} (8)

where $R_{\rm o}$ is the Rossby number for turbulent motion and Sc denotes the Schmidt number. In the present calculations we assume $R_{\rm o}=3$ and ${Sc}=1+\Lambda$. Substituting (8) into (7) we obtain the size evolution equation

 \begin{displaymath}\frac{{\rm d}a}{{\rm d}t}= \frac{\rho_{\rm d}}{\rho_{\rm s}}\... ... d}}{\rho_{\rm s} H_{\rm d}}\ \frac{\sqrt{\Lambda}}{1+\Lambda} \end{displaymath} (9)

where we used the relation $\rho_{\rm d}=\Sigma_{\rm d}/2 H_{\rm d}$. An approximate formula for the half-thickness of the SPD, $H_{\rm d}$, has been derived by SV97 and can be expressed as
 
$\displaystyle {\left(\frac{H_{\rm d}}{H}\right)}^2 = \frac{\alpha \sqrt{2}}{2} ... ...+8\Lambda (\Lambda+\sqrt{2}\alpha)}} {(1+\Lambda)(\alpha\sqrt{2}+\Lambda)}\cdot$     (10)

The right hand side of (9) is the source function for particle growth by coagulation, f, present in Eq. (5). From (9), (10), and the fact that $\Lambda$ is given in terms of a particle size and the gas properties it is clear that f is given in terms of only a, $\Sigma _{\rm d}$ and the gas properties. Thus substituting f and  $V_{{\rm d}r}$ into the the system of Eqs. (4)-(5) yields the closed system with two unknowns, $\Sigma _{\rm d}$ and  $\Sigma_{a}$.

2.5 Evaporation

The midplane temperature, T, of the gas in the disk is calculated from our GPD model (Sect. 2.1). For any particular time the temperature decreases with distance from the star. The particle traveling inward will evaporate when it finds itself at the location where the ambient gas is sufficiently hot. Such a location defines an evaporation radius, $r_{\rm evap}$, that depends on the composition of the particle. The evaporation radius decreases with time, as the entire disk cools down during its evolution. We consider, in separate calculations, three distinct, idealized forms of solid particles, water ice, low-temperature silicates, and high-temperature silicates to account for solids species with different evaporation temperatures and bulk densities. The corresponding evaporation temperatures, $T_{\rm evap}$, are 150 K, 400 K and 1350 K, and the corresponding bulk densities, $\rho_{\rm s}$, are 1, 2.6 and 3.3 gcm-3. In our calculations we treat vapor as particles with $a \rightarrow 0$ that are completely coupled to the gas. The model permits the outward movement of vapor particles into the $T<T_{\rm evap}$ environment where they condense. However, in our calculations, the radial velocity of solids at the location of the evaporation radius was always directed inward and such condensation has not occurred.

3 Numerical methods

Equations (4) and (5) are solved numerically to obtain spatio-temporal distributions of mass in the SPD. The SPD expands and/or shrinks during its evolution, and its maximum extent is a priori unknown. Because of this feature we solve our equations on a moving grid, whose outer boundary follows the motion of the outer edge of the SPD. We define a new independent variable, $x=r/s_{\rm d}$ where $s_{\rm d}(t)$ is the instantaneous outer radius of the SPD. In order to preserve the forms of Eqs. (4) and (5) we also introduce new dependent variables $\Sigma_{\rm d}'$ and $\Sigma_{a}'$ and rescale the SPD's radial velocity
 
$\displaystyle \Sigma_{\rm d}'$ = $\displaystyle \Sigma_{\rm d} s_{\rm d}^2$ (11)
$\displaystyle \Sigma_{a}'$ = $\displaystyle \Sigma_{a} s_{\rm d}^2$ (12)
$\displaystyle V_{{\rm d}r}'$ = $\displaystyle \frac{(V_{{\rm d}r}-\dot{s}_{\rm d})}{s_{\rm d}}\cdot$ (13)

The outer edge of the disk is now always at $x_{\rm max}=1$, whereas the inner edge has to be arbitrarily assigned a small number $x_{\rm min}=\beta$. The value of $\beta=10^{-3}$ was used in our calculations. We use a grid of 50 points equidistant in $\log x$, and distributed between $\log x_{\rm min}$ and $\log x_{\rm max}$. The transformed equations are solved using standard methods of numerical hydrodynamics. We employ a second-order code originally described by Rózyczka (1985), and based on the same philosophy as the popular ZEUS code (Stone & Norman 1992). However, we use a monotonization procedure proposed by van Albada et al. (1982), instead of the standard monotonization method (van Leer 1977) used in Rózyczka's code. Because $V_{{\rm d}r}'(t, x_{\rm max})=0$ no boundary condition at $x=x_{\rm max}$ is required. The value of $V_{{\rm d}r}'(t, x_{\rm min})< 0$ because the inner edge is always in the accretion zone of the disk. Therefore a free outflow boundary conditions is applied at  $x_{\rm min}$.

Several tests of the code were performed (analytical solutions of simple evolutionary cases were recovered, and in a few more complicated cases full agreement was reached with numerical solutions of the original set of equations). We also checked that the results were not sensitive to the choice of $\beta$ (provided its value was not too large).

4 Results

We have performed a large suite of calculations to answer the following question: What is the dependence of the global properties of the nascent planetary system (represented here by the swarm of planetesimals at t=107 yr) on initial conditions of the GPD? Each individual calculation tracks the evolution of the GPD and the SPD from a particular initial condition (at the time arbitrarily set to t=0) to a configuration reached t=107 years later when all solids are either in the form of planetesimals or have been accreted onto the star or evaporated. Each evolutionary run is labeled by values of four parameters $(j_0, m_0, \alpha, {\rm species})$. The first two labels, total angular momentum and total mass of the GPD at t=0, characterize the initial conditions (see Sect. 2.1). To cover the range of masses and angular momenta indicated by observations, we use 14 values of j0 between 0.5 and 50 and 11 values of m0 between 0.02 and 0.2 (i.e. we consider 154 different initial conditions). The initial radial distribution of solids is $\Sigma_{\rm d}(t=0,r)=\delta \Sigma (t=0,r)$, where $\delta=10^{-2}$for ice and $\delta = 6 \times 10^{-3}$ for silicates to account for cosmic abundance. The third label, the Shakura-Sunayev parameter $\alpha $, characterizes the turbulence in the GPD and influences the rate of the disk's evolution. We use four values of $\alpha $, 10-4, 10-3, 10-2 and 10-1. Finally, the fourth label describes the composition of solids, their evaporation temperature and the bulk density. We consider three different species of solids (see Sect. 2.5). Altogether 1848 evolutionary runs were calculated.

The global properties of the swarm of planetesimals are represented by two numbers, the total mass of the SPD at t=107 yr, $m_{\rm d}$, and the radius of the SPD at t=107 yr, $r_{\rm d}$. Thus, we use the results of our calculations to construct a transformation $(j_0, m_0) \mapsto (m_{\rm d}, r_{\rm d})$.

  \begin{figure} \par\includegraphics[width=14cm,clip]{KSM1537f1.eps} \end{figure} Figure 1: The radius (the bottom face of each cube) and the total mass (the top face of each cube) of icy planetesimal swarm at t=107 years as functions of initial conditions. The functions are shown in the form of contour plots. The contours' labels for $r_{\rm d}$ are in units of AU and the contours' labels for $m_{\rm d}$ are in units of  $M_{\oplus }$. Cubes are marked by corresponding values of $\alpha $.
Open with DEXTER

Figure 1 shows the results of our calculations for solids composed of water ice. The figure has four panels, each for a different value of $\alpha $. For a particular value of $\alpha $ the plot shows a "cube'' with $r_{\rm d}(j_0,m_0)$ displayed on the bottom face and $m_{\rm d}(j_0,m_0)$ displayed on the top face. These two functions are represented as contour plots. Each labelled curve on the bottom face of the cube represents a set of initial conditions that result in a planetesimal swarm with the radius indicated by the label. Similarly, each labelled curve on the top face of the cube represents a set of initial conditions that results in a planetesimal swarm with the total mass indicated by the label. The functions $r_{\rm d}(j_0,m_0)$ and $m_{\rm d}(j_0,m_0)$ are obtained via interpolation from the 154 values actually calculated. A wavy character of some contours lines is due to limited coverage of calculated values and would disappear for a dense enough grid of (j0, m0). The cube is constructed from the two contour plots in order to better illustrate the full transformation $(j_0, m_0) \mapsto (m_{\rm d}, r_{\rm d})$. One can choose a particular initial condition, (j0, m0), on, say the bottom face of the cube, find a corresponding value of $r_{\rm d}$ and then move straight up to the top face of the cube to find a corresponding value of $m_{\rm d}$.

Analyzing Fig. 1 we find that the domain of initial conditions (j0, m0) can be divided into two portions by contours labeled "0''. These contours, that coincide on the bottom and the top face of each cube (less than perfect coincidence is due to imperfections of interpolation), are the demarcation lines between initial conditions that lead to formation of a planetesimal swarm and those that do not. Disks evolving from initial conditions located below a 0-labeled contour are, over time, bled from all solids leaving no material behind to form planetesimals. Note that such, no-swarm, initial conditions are those that correspond to the GPD that is initially heavy and small (small values of j0 correspond to small values of r0) and thus hot. The no-swarm initial conditions are, in the first approximation, those for which $r_0<r_{\rm evp}(t=0)$. In such a disk all solids are initially in the vapor form. As the GPD evolves it expands and cools, its radius increases and the value of  $r_{\rm evp}$ decreases, leading, eventually, to the emergence of the region of the disk where solids may exist. However, in the process of disk expansion the major portion of gas and vapor moves inward and is accreted onto the star. Only a small portion of gas and vapor moves outward, increasing the size of the disk. This small amount of the vapor condenses into small solids right outside  $r_{\rm evp}$. On a relatively short timescale these solid particles coagulate, gain inward radial velocities and move back to the location $r<r_{\rm evp}$ where they evaporate again. Eventually, all vapor is accreted onto the star and the disk is bled from the solid component. The region of no-swarm initial conditions is the largest for $\alpha =10^{-1}$ because, everything else being equal, disks with higher values of $\alpha $ are hotter. Decreasing $\alpha $ results in progressively smaller set of no-swarm initial conditions. Finally, for $\alpha=10^{-4}$ all initial conditions considered here lead to the formation of a swarm.

It is interesting to understand how the size and the mass of the swarm depend on initial conditions. Consider a family of initial conditions defined by $m_0={\rm const}$ line on any frame in Fig. 1. Such a family can be sorted by an increasing value of j0 starting from  $j_{0{\rm min}}$ that corresponds to an initial condition located just on the 0-labeled contour. The masses of swarms formed in disks evolving from succeeding initial conditions increase quickly at first, but then stay constant after reaching a certain "saturation'' mass. The radii of swarms increase without saturation, but at a decreasing rate. This pattern is independent of the values of m0 and/or $\alpha $, although the value of the saturation mass depends on these parameters.

These results can be explained as follows. As discussed above, disks with $j_0 \le j_{0{\rm min}}(m_0,\alpha)$, accrete all solids onto the star. For disks with j0 just larger than  $j_{0{\rm min}}$, a small outer portion of the disk has a temperature smaller than  $T_{\rm evp}$and small solid particles exist there at t=0. They coagulate on a relatively short timescale to achieve the maximum-speed size. At that stage $V_{{\rm d}r} \sim V_{{\rm d}r}^{\rm max}$ and most of particles are lost to the evaporation zone. Some particles survive because they manage to grow to sizes bigger than the maximum-speed size and slow down before falling into the evaporation zone (see Sect. 2.3). These particles form a residual swarm, whose mass is much smaller than the original mass of the SPD. The larger the value of j0, the larger is the initial disk, the more particles avoid falling into the evaporation zone, and the more massive swarm is formed. For large enough values of j0, the great majority of particles have enough time to grow to large sizes that prevent them from migration into the evaporation zone. The mass of the resulting swarm is about equal to the initial mass of the SPD. In such cases the solid material is reshuffled within the disk, but its total mass is basically preserved.

  \begin{figure} \par\includegraphics[width=14cm,clip]{KSM1537f2.eps} \end{figure} Figure 2: The radius (the bottom face of each cube) and the total mass (the top face of each cube) of high-temperature silicates planetesimal swarm at t=107 years as functions of initial conditions. The functions are shown in the form of contour plots. The contours' labels for $r_{\rm d}$ are in units of AU and the contours' labels for $m_{\rm d}$ are in units of  $M_{\oplus }$. Cubes are marked by corresponding values of $\alpha $.
Open with DEXTER

Figure 2 shows the results of our calculations for solids composed of high-temperature silicates. The transformation $(j_0, m_0) \mapsto (m_{\rm d}, r_{\rm d})$ for the SPD consisting of silicates is qualitatively similar to the transformation for the SPD consisting of ice. This is to be expected considering that mechanism of such a transformation depends on particle material via just one parameter ( $T_{\rm evap}$). The key differences between silicates and the ice are as follows. The no-swarm portion of the (j0, m0) plane is slightly smaller for silicates than it is for the ice because $r_{\rm evap}({\rm silicates})<r_{\rm evap}({\rm ice})$ so, other things being equal, $j_{0{\rm min}}$ is smaller for the silicates. For the same, swarm-forming initial conditions, the size of the ice swarm is larger than the size of the silicate swarm, mostly because the inner radius of the silicate swarm is smaller than the inner radius of the ice swarm. Note that, although the original protoplanetary disk contains much more ices than silicates, the final silicates swarm is not necessarily less massive than the final ice swarm. Results for low-temperature silicates (not shown here) fall between those for ice and high-temperature silicates.


  \begin{figure} \par\includegraphics[width=14cm,clip]{KSM1537f3.eps} \end{figure} Figure 3: Summary of the evolution of protoplanetary disk with j0=10.4, m0=0.182 and $\alpha =10^{-1}$. The bottom face of each cube is a contour plot depicting the surface density as function of t and r. The contours' labels are in units of g/cm2. The top face of a gas cube shows T(t,r), contours' labels are in K. The top faces of solids cubes show a(t,s), contours' labels are in cm. Note that cubes have different radial ranges.
Open with DEXTER

5 Details of disk evolution

Although the focus of this paper is on the $(j_0, m_0) \mapsto (m_{\rm d}, r_{\rm d})$ transformation (see Sect. 4), it is also interesting to track the evolution of an example protoplanetary disk in detail. In this section we present results of our calculations pertaining to the evolution of a disk characterized by $\alpha =10^{-1}$ and starting from initial conditions parameterized by j0=10.4 and m0=0.182. We have chosen this particular case because it provides an example of rather dramatic evolutionary changes of gas and solids components, and it illustrates all relevant physical processes. In other, possibly more realistic, scenarios not all evolutionary features may be observed in a single run.

At t=0 the GPD extends to 68 AU. Ice SPD, low-temperature silicates SPD and high-temperature silicates SPD extend from 12.2 AU to 68 AU, 3.4 AU to 68 AU and 1.1 AU to 68 AU, respectively. Their masses are, respectively, equal to 606, 364, 364  $M_{\oplus }$. The mass of the SPD is entirely in the form of small particles (assumed size is 10-3 cm). At t=107 yr the GPD has lost most of its mass, but it extends to a very large radius (7300 AU). The ice SPD extends from 0.64 AU to 11.6 AU and has a mass of 56  $M_{\oplus }$. The low-temperature silicates SPD extends from 0.17 AU to 5.4 AU and has a mass of 38  $M_{\oplus }$. The high-temperature silicates SPD extends from 0.08 AU to 5.5 AU and its mass is 84  $M_{\oplus }$. Thus, during the 107 yr of evolution the disk lost more than 90% of ice and low-temperature silicates and about 80% of high-temperature silicates. All solids have been redistributed, their final distributions bearing no resemblance to the original ones. Most mass of the final SPD is locked in planetesimals with a > 104 cm. Note that the sizes of swarms given above are rather formal values, as the bulk of the mass resides in a much narrower ring (see below).

Figure 3 shows the summary of the disk evolution. The figure has four panels, one for the gas, and the rest for the three species of solids. Each panel is in the form of a "cube'', in the same sense as in Figs. 1-2, to represent a pair of functions of two independent variables, t and r. For GPD, these functions are T(t,r) (shown on the top face of the cube), and $\Sigma(t,r)$(shown on the bottom face of the cube). For SPD the functions are a(t,r) (shown on the top face of the cube), and $\Sigma_{\rm d}(t,r)$ (shown on the bottom face of the cube). All functions are represented as contour plots obtained via interpolation from a finite number of values that actually were computed. Note that the computed values are not given on a regular grid of (t,r), and a triangular interpolation has to be used.

Functions T(t,r) and $\Sigma(t,r)$ describe the evolution of the GPD. There are six constant temperature contours, three of them chosen to coincide with  $T_{\rm evap}$ for the three species of solids. The curvatures of the temperature contours point to the overall cooling of the disk. In particular the evaporation radii for all species of solids decrease with time. There are 10 constant surface density contours to illustrate the evolution of the gaseous component. These contours show that the mass of the GPD decreases rapidly. They also illustrate the spreading of the GPD (note that contours corresponding to $\Sigma \le 30$ gcm-2 start at progressively later times).

Functions a(t,r) and $\Sigma_{\rm d}(t,r)$ indicate the evolution of the SPD. The evolution of all solids species are qualitatively similar. For t > 10 yr the contour a=10-2 cm coincides with the  $T_{\rm evap}$ contour and demarcates disk locations where only vapor exists from locations where solid bodies may exist. The "mohawk-like'' pattern formed by the collection of constant a contours sheds light on how the evolution of solids proceeds. A t=t0 line on the contour plot depicting the function a(t,r) gives the radial distribution of particle size at t0, whereas the r=r0 line gives the history of particle sizes at r0. The inward migration of the outer edge of the SPD can be deduced from the fact that contours corresponding to progressively larger sizes end at progressively smaller radii. At earlier times the contours of constant a are almost parallel to the r axis, indicating small variability of particle sizes along the radius of the disk. For later times the bending of the contours increases indicating more particle size variability.

To fully understand the evolution of the SPD, we need to inspect the function $\Sigma_{\rm d}(t,r)$ in addition to the function a(t,r). The most important, additional information from studying $\Sigma_{\rm d}(t,r)$ is that it falls off rapidly at a certain radius that is almost time-independent and is close to $r_{\rm evap}(t=0)$. Thus, although some solids exist between $r_{\rm evap}(t)$ and $r_{\rm evap}(t=0)$ (as shown by the particle size plot), their surface density there is negligible and we can take $r=r_{\rm evap}(t=0)$ as an effective inner radius of the SPD. The explanation of this phenomenon is as follows. The viscous timescale of this model is long, and up to t=104 yr the location of $r_{\rm evap}$ does not change much. However, the coagulation timescale is shorter, and by that time solids, initially in the form of small particles, have coagulated into relatively large bodies. From t=104 yr to t=107 yr $r_{\rm evap}$ decreases steadily and the solids are free, in principle, to migrate inward. However, again the coagulation timescale is faster than the viscous timescale and few solid bodies actually migrate inward, whereas the majority manage to stop their migration by increasing their size. This phenomenon was already seen in SV97 calculations and its detailed explanation is given there.

Figure 4 shows, on the linear scale, functions $\Sigma_{\rm d}(t=10^7,r)$ for all species of solids. This figure carries the same information, albeit in a more readable form, as the information along t=107 yr lines on contour plots of Fig. 3. In our calculations we consider the swarm to be fully assembled at t=107yr, and thus Fig. 4 illustrates the final distribution of mass locked in planetesimals. In this particular swarm the high and low temperature silicates coexist between 2 AU and 5.5 AU, but only high-temperature silicates are found between 1 AU and 2 AU. The amount of solids at radii smaller than 1 AU is negligible. Icy planetesimals do exist in a significant amount between 6 AU and 11.6 AU, i.e. they are completely separated from silicates. Most of the ice resides at about 10. There is enough material in the swarm to form a planetary system composed of an icy planet (or planets) with the total mass of $M\leq0.14\,M_{\rm J}$ located at about 10 AU and a silicate planet (or planets) with the total mass of $M\leq0.18\,M_{\rm J}$ located between 1 AU and 5 AU.

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{KSM1537f4.eps} \end{figure} Figure 4: Radial distribution of the surface density of solids in a disk characterized by $\alpha =10^{-1}$, j0=10.4 and m0=0.182 at t=107 yr. The values of  $\Sigma _{\rm d}$ were computed at finite number of locations indicated by circles. Light gray circles are for the ice, dark gray circles are for the low-temperature silicates and black circles are for high-temperature silicates. The size of the circle is proportional (nonlinearly) to the logarithm of the solids size. The dashed line indicates the surface density of the gas.
Open with DEXTER

6 Discussion and conclusions

We have developed a formalism that allows us to calculate globally the evolution of solids in a protoplanetary disk. The formalism is based on the method developed by SV97, but it relies on two additional approximations. First, we use an analytical model of GPD in place of the numerical model used in SV97. Second, for gas-solids interactions we assume a coupling between solid particles and average flow of the gas instead of a coupling between particles and instantaneous, turbulent flow as considered in SV97. Because of these changes the computational efficiency of our new formalism is significantly higher. Both formalisms produce qualitatively similar results, as can be assessed by comparing details of disk evolution (Sect. 5) with analogous results in Sects. 4 and 5 of SV97. Thus, as far as the evolution of a particular disk is considered, we confirm the earlier results discussed in detail in Sect. 6 of SV97.

The efficiency of the new formalism allows for a broad search through plausible initial conditions in order to establish the relation between those conditions and the properties of emerging planetesimal swarms. To investigate the diversity of the swarms, we performed a systematic search through initial conditions by sampling a 2-D parameter space confined by the minimum and maximum values of j0 and m0 as indicated by observations.

Our calculations produce planetesimal swarms, which are taken here as surrogates for planetary systems. In reality, further redistribution of mass may occur during the planetesimals to planets accumulation, as well as during subsequent gravitational interactions between just formed planets and a residual swarm (Del Popolo et al. 2001) or a gaseous disk (Ward 1997). These processes introduce further diversity of planetary systems that is not addressed by our calculations. It should be also mentioned that a profound redistribution of various PPD components, accompanied by a significant mass loss from the disk, may result from eruptive events observed in those objects (e.q. FU Ori outbursts). Such events, likely to occur prior to a dissipative stage of the disk, are beyond the scope of the present paper.

Our most important finding is that the diversity of planetary systems emerges naturally from the process of the evolution of solids in protoplanetary disks of diverse initial configurations. In other words, if protoplanetary disks form in variety of masses and sizes, this variation alone leads to formation of extremely different planetary systems. In particular, we have found that for certain initial configurations of protoplanetary disks, a planetesimal swarm, and thus the planetary system, does not form at all. These configurations are characterized, in general, by small specific angular momentum (small initial sizes relative to their initial masses). In these disks solids migrate into the evaporation zone where they are destroyed and accreted onto the star before they have time to grow and stop their inward motion. The size of the no-swarm domain of initial conditions depends on the species of solids and on the timescale of disk evolution (related to the viscosity parameter $\alpha $). Configurations with large values of specific angular momentum (large initial sizes with respect to their masses) tend to redistribute solids during their evolution without losing them to the star. The residual swarms have relatively large masses. Finally, configurations with intermediate values of specific angular momentum lose a fraction of solids during their evolutions, producing less massive swarms. Statistical analysis of astronomical observations may, in principle, establish relative frequencies of different types of initial conditions. Such analysis, which is beyond the scope of this paper, coupled to our model may predict relative frequencies of different sorts of potential planetary systems.

It is interesting to see which initial conditions lead to a planetary system that resembles our Solar System. We can approximate the size of the planetesimal swarm from which the Solar System formed by the current radius of the Kuiper Belt (30-50 AU). For the mass of the swarm we can take the mass of solids in the present Solar System, <100  $M_{\oplus }$. Figure 1 shows that regardless of the value of $\alpha $, swarms with such properties can originate only from a protoplanetary disks with $m_0 \le 0.02$ and $j_0 \le 3$. Specifically, assuming $\alpha=10^{-2}$, only two models among those actually calculated fulfill roughly the above conditions. The first model, with j0=3 and m0=0.02, originally extends to 440 AU and has the mass of 67  $M_{\oplus }$ in icy particles and 40  $M_{\oplus }$ in high-temperature silicates. The final icy swarm extends from 1 AU to 54 AU and has the mass of 66  $M_{\oplus }$, while the final silicates swarm extends from 0.3 AU to 37 AU and has the mass of 39  $M_{\oplus }$. The second model, with j0=1 and m0=0.02, originally extends to 50 AU and has the mass of 67  $M_{\oplus }$ in icy particles and 40  $M_{\oplus }$ in high-temperature particles. The final icy swarm extends from 2.25 AU to 18 AU and has the mass of 64  $M_{\oplus }$, while the final silicates swarm extends from 0.2 AU to 12 AU and has the mass of 39  $M_{\oplus }$.

Neither of these two models fits perfectly with the Solar System. The j0=3 model leaves too much ice in terrestrial planets zone, while the j0=1 model produces a swarm that is smaller than the present-day Solar System. Nevertheless, both models produce swarms that, at least in the global, qualitative sense, can lead to the formation of the planetary system somewhat similar to our own. Interestingly, disks that produce solar-like planetary systems have initial global properties similar to that postulated in the so-called "minimum-mass'' solar nebula model (Weidenschilling 1977). This similarity is coincidental inasmuch as in our model there is no a priori requirement of preserving the total mass of the solids. In other words, the fact that disks leading to the solar-like planetary system preserve the original mass of the solid component is a result (and perhaps an unexpected one) from our model, and not a build-in property.

More massive swarms are possible if initial disks are more massive and have enough angular momentum. The amount of gas left in the disk after 107 years depends on the value of $\alpha $ (the smaller the $\alpha $, the more gas is left). Portions of this gas can be incorporated into planetary masses, while the remaining gas interacts with newly formed planets, causing further changes to the configuration of a system. Within the limits of the approximations adopted here, the in situ formation of massive planets at 1-5 AU from a star can be explained in terms of SPD evolution alone (i.e. without planetary migration). However, in situ formation of such planets at 0.1 AU (the so-called "hot Jupiters'', of which the companion to 51 Peg is a primary example) is rather unlikely, as reasonable initial conditions do not produce swarms with adequate surface density distributions.

Finally, it is important to remember that our model is built on the premise that transformation of solids from dust to planetesimals occurs through the process of hierarchical coagulation. If coagulation is not the major factor in such transformation, our model is invalid. One cannot definitively exclude a possibility that a collective process, such as, for example, the gravitational instability of a dust layer (Goldreich & Ward 1973) transforms small particles rapidly into planetesimals. However, the formation of a midplane layer of dust, thin enough to be susceptible to the gravitational instability is unlikely in a turbulent disk. It has been proposed (Barge & Sommeria 1995; Tanga et al. 1996) that the gravitational instability may work for the solids component in vortices appearing in the disk, but neither the existence of such vortices nor their ability to form planetesimals was demonstrated. Thus, coagulation remains the most likely mechanism of a planetesimal formation. The "single-size'' approximation used here for the description of coagulation is obviously rather crude. However, we believe that, despite its limitations, our model captures the essence of processes leading to formation of a planetesimal swarm, and it provides a good qualitative illustration of the evolution of solids in a protoplanetary disk.

Acknowledgements

This research was conducted at the Lunar and Planetary Institute, which is operated by the Universities Space Research Association under contract No. NASW-4574 with the National Aeronautics and Space Administration, and at the Nicolaus Copernicus Astronomical Center. This is Lunar and Planetary Institute Contribution No. 1096. K. K. and M. R. acknowledge support from the Committee for Scientific Research through the grant 2P03D01419. The authors thank the referee, Dr. J. Lunine, for constructive comments.

References

 
Copyright ESO 2001