EDP Sciences
Highlight
Free Access
Issue
A&A
Volume 501, Number 3, July III 2009
Page(s) 1139 - 1160
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/200810301
Published online 27 May 2009

Extrasolar planet population synthesis

I. Method, formation tracks, and mass-distance distribution

C. Mordasini1,[*] - Y. Alibert1,2 - W. Benz1

1 - Physikalisches Institut, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland
2 - Institut UTINAM, CNRS-UMR 6213, Observatoire de Besançon, BP 1615, 25010 Besançon Cedex, France

Received 1 June 2008 / Accepted 15 April 2009

Abstract
Context. With the high number of extrasolar planets discovered by now, it has become possible to use the properties of this planetary population to constrain theoretical formation models in a statistical sense. This paper is the first in a series in which we carry out a large number of planet population synthesis calculations within the framework of the core accretion scenario. We begin the series with a paper mainly dedicated to the presentation of our approach, but also the discussion of a representative synthetic planetary population of solar like stars. In the second paper we statistically compare the subset of detectable planets to the actual extrasolar planets. In subsequent papers, we shall extend the range of stellar masses and the properties of protoplanetary disks.
Aims. The last decade has seen a large observational progress in characterizing both protoplanetary disks, and extrasolar planets. Concurrently, progress was made in developing complex theoretical formation models. The combination of these three developments allows a new kind of study: the synthesis of a population of planets from a model, which is compared with the actual population. Our aim is to obtain a general overview of the population, to check if we quantitatively reproduce the most important observed properties and correlations, and to make predictions about the planets that are not yet observable.
Methods. Based as tightly as possible on observational data, we have derived probability distributions for the most important initial conditions for the planetary formation process. We then draw sets of initial conditions from these distributions and obtain the corresponding synthetic planets with our formation model. By repeating this step many times, we synthesize the populations.
Results. Although the main purpose of this paper is the description of our methods, we present some key results: we find that the variation of the initial conditions in the limits occurring in nature leads to the formation of planets of wide diversity. This formation process is best visualized in planetary formation tracks in the mass-semimajor axis diagram, where different phases of concurrent growth and migration can be identified. These phases lead to the emergence of sub-populations of planets distinguishable in a mass-semimajor axis diagram. The most important ones are the ``failed cores'', a vast group of core-dominated low mass planets, the ``horizontal branch'', a sub-population of Neptune mass planets extending out to 6 AU, and the ``main clump'', a concentration of giant gaseous planets at around 0.3-2 AU.

Key words: stars: planetary systems - stars: planetary systems: formation - stars: planetary systems: protoplanetary disks - planets and satellites: formation - solar system: formation - methods: numerical

1 Introduction

As of spring 2009, more than 300 extrasolar planets have been discovered (Schneider's Extrasolar Planet Encyclopedia at http://exoplanet.eu). The richness and diversity of the characteristics of these exoplanets, like their mass or semimajor axis, is impressive, and was not necessarily expected from the single example - our own solar system - that was available to us before the discovery of 51 Peg b (HD 217014b) by Mayor & Queloz (1995).

Since then, the observational field of extrasolar planet searches has seen a rapid evolution leading to numerous additional discoveries of planets orbiting other stars. These discoveries have also triggered numerous theoretical studies about the formation and evolution of these planets. Key physical processes in planet formation and evolution could be identified whose importance was not fully realized in previous works based on the solar system alone.

Some of these discovered planets, and multiple planetary systems, are sufficiently interesting by themselves to warrant individual theoretical studies. Examples are the extrasolar planetary system with three Neptune-mass planets around HD 69830 (Lovis et al. 2006; Alibert et al. 2006), or the transiting Neptune mass planet GJ 436b (Butler et al. 2004; Gillon et al. 2007; Figueira et al. 2009). Of course, the giant planets in our own solar system provide a much larger and more detailed set of constraints than any known extrasolar planet. Therefore, each formation model applied to discuss extrasolar planet formation should also be put to the test to reproduce the characteristics of our own giant planets (Pollack et al. 1996; Alibert et al. 2005b; Hubickyj et al. 2005; Benvenuto & Brunini 2005).

The modeling of the formation of such single systems, while a necessary condition to validate formation models, is not satisfactory by itself. Indeed, the number of model parameters is generally large while the number of constraints deriving from a single system is small, and not strong enough to completely constrain any formation model.

Thanks to the rapid growth of the number of known extrasolar planets, the situation has dramatically changed: instead of having only a single object or a single system to study, we now begin to be able to describe an entire population of extrasolar planets orbiting FGK stars in the solar neighborhood. While this population is still smaller than one would ideally like, it nevertheless already allows us to extract statistically a wealth of information (e.g. Udry & Santos 2007; Cumming et al. 2008) to constrain formation models that exceeds by far what one extrasolar planet can do. This is especially true since most of the extrasolar planets have been discovered by radial velocity measurements so that only a few orbital elements and a minimum mass are known for one individual object. For the growing number of transiting planets more physical properties can be derived and compared with internal structure models (Baraffe et al. 2008; Figueira et al. 2009). Unfortunately, transiting planets known so far are all in close proximity to their host star. Hence it is sometimes unclear to what extend their characteristics are related to their formation or rather to subsequent evolution (e.g. evaporation).

Parallel to the discovery of more and more end-products of the planetary formation process i.e. planets, much observational progress (e.g. Meyer et al. 2006) has also been made in characterizing the initial conditions for this process, i.e. the protoplanetary disks. Thanks to these observations, we begin to be able to determine the probability of occurrence of any particular initial condition for planetary formation, like disk metallicity, mass or lifetime.

With these two sets of observational data at hand, a new interesting class of theoretical planet formation studies has become possible, where a theoretical model serves as the link between these two groups of observations: the synthesis of populations of planets by Monte Carlo methods. In this approach the observed distributions of disk properties are used as varying initial conditions for the model. The final characteristics of the synthetic planets that form in the model can then be compared statistically to those of the actual observed populations. This addresses the question of whether the observed diversity of extrasolar planets is simply the consequence of the diversity of disk properties.

As we shall show, such studies have proven to be very fruitful, as they not only allow us to reproduce observations but also show the links and correlations between the different initial conditions and the characteristics of the resulting planets. Thereby they provide insights into the formation mechanism. Such an approach, by predicting the actually existing planet population as opposed to the actually detected one, allows to optimize future searches and instruments when coupled to a synthetic detection bias for a particular detection method.

Compared to similar studies e.g. the pioneering work of Ida & Lin (2004a,b, 2005, 2008), or the studies of Kornet & Wolf (2006), Robinson et al. (2006) and Thommes et al. (2008), we have attached particular importance to three distinct areas: first, we use the detailed extended core accretion formation model of Alibert et al. (2005a) that has been successfully applied to quantitatively explain the many observed constraints of the giant planets of our own solar system (Alibert et al. 2005b). Second, we stick as tightly as possible to probability distributions derived from observations, and third (cf. the companion paper Mordasini et al. 2009, hereafter Paper II), we use quantitative statistical methods to compare model outcomes and observations and require that as many different observational constraints as possible be satisfied at the same time. In this way, we can check which observed properties can be reproduced by our formation model. But also, discrepancies between the synthetic and the actual population provide new insights, allowing us to improve the models.

In this first paper, we present the methods we use to generate the synthetic population, in particular the formation model and the probability distributions for the Monte Carlo variables. We then show the resulting planetary formation tracks in the semi-major vs. mass plane. These tracks are of fundamental importance to understand the characteristics of the resulting synthetic population, as they illustrate how and why planets reach their final position in the mass-distance diagram. This a-M distribution at the end of the formation phase is characterized by a number of structures (clumps, concentrations and depletions). Some particular regions in this diagram are identified and discussed. In the companion Paper II, we use a synthetic detection bias for the radial velocity method to identify the subset of detectable synthetic planets. We then compare this sub-population with statistical, quantitative methods to the actual extrasolar planet population. In these first two papers, we assume a mass of the host star of 1 $M_\odot $, as most known extrasolar planets are found around solar type stars. In later papers in this series we shall study the influence of different host star masses as well as of different formation environments (i.e. disk properties). Some of these have already been observationally identified, such as the well known correlation between the stellar metallicity and the detection probability of giant planets.

The outline of the paper is as follows: in Sect. 2 we give an overview of our formation model, with a focus on the modifications and necessary simplifications[*] relative to Alibert et al. (2005a). In Sect. 3 the Monte Carlo approach is described, and we determine the probability distributions for the initial conditions in Sect. 4. Section 5 illustrates the numerical population synthesis process with example formation tracks in the distance-mass plane and discusses the properties of the planet population. The conclusions are drawn in Sect. 6.

2 Giant planet formation model

The link between the initial conditions i.e. the properties of the protoplanetary disk, and the final outcome i.e. the planets, can only be given by a theoretical formation model. This link is usually very complicated, involving many feed-back mechanisms and nonlinearities. For the simulations presented in this paper, we calculate in a consistent way the formation of a protoplanet, its migration, and the structure and evolution of the protoplanetary disk.

2.1 Disk structure and evolution

The structure and evolution of the protoplanetary disk is modeled as a non-irradiated, 1+1D $\alpha $-disk (Shakura & Sunyaev 1973), following the method originally presented in Papaloizou & Terquem (1999). We thus solve the diffusion equation (effective viscosity  $\tilde{\nu}$) describing the evolution of the gas surface density $\Sigma$ as a function of time t and distance a to the star:

\begin{displaymath}%
{{\rm d} \Sigma \over {\rm d} t} = {3 \over a} {\partial \o...
...l a} ( \tilde{\nu} \Sigma a^{1/2}) \right] + \dot{\Sigma}_w(a)
\end{displaymath} (1)

The photo-evaporation term $\dot{\Sigma}_w$ is given by (Veras & Armitage 2003):

\begin{displaymath}%
\dot{\Sigma}_w =\left\{ \begin{array}{ll}
0 & \textrm{for}\...
...a_{\rm max}-R_{g}) a} & \textrm{otherwise}
\end{array} \right.
\end{displaymath} (2)

where Rg is taken to be 5 AU, $a_{\rm max}$ is the size of the disk, and the total mass loss $\dot{M}_{w}$ due to photo-evaporation is an input parameter which together with the $\alpha $ parameter determines the lifetime of the disk.

For simplicity, we adopt an initial profile of the gas disk surface density according to the phenomenological model of Hayashi (1981), $\Sigma(a,t=0)=\Sigma_{0}\left(a/a_0\right)^{-3/2}$ where $\Sigma _{0}$ is the surface density at our reference distance (a0 = 5.2 AU), and the computational disk extends from $a_{\rm min}$ = 0.1 AU to $a_{\rm max}$ = 30 AU. The initial total gas disk mass in the computational disk is then $4 \pi \Sigma_{0}a_0^{3/2}(a_{\rm max}^{1/2}-a_{\rm min}^{1/2})$. For the initial profile with $\Sigma\propto a^{-3/2}$ the accretion rate decreases from the inner to the outer parts of the disk. As shown in e.g. Papaloizou & Terquem (1999), the inner parts of the disk evolve rapidly toward a state of constant accretion rate  $\dot{M}_{\rm *}$. Therefore, the inner initial gas disk profile is truncated in order to obtain an accretion rate lower than a constant value of order 3 $\times $ $10^{-7}~M_\odot/$yr. This allows us to speed up the calculation of the disk evolution.

The initial solid surface density is given by $\Sigma_{\rm D}= f_{\rm D/G}f_{\rm R/I}\Sigma_{0}\left(a/a_0\right)^{-3/2}$ (Hayashi 1981; Weidenschilling et al. 1997) where $f_{\rm D/G}$ is the dust-to-gas ratio of the disk, and $f_{\rm R/I}$ is a factor describing the degree of condensation of ices. Its value is set to 1/4 in the regions of the disk for which the initial mid-plane temperature exceeds the sublimation of water ice ( $T_{\rm mid} > 170$ K), and 1 otherwise. The semimajor axes  $a_{\rm ice}$ where this temperature is reached as a function of initial gas surface density  $\Sigma _{0}$ is plotted in Fig. 1. For a minimum mass solar nebula (MMSN) like $\Sigma _{0}$ (100-200 g/cm2), the iceline is as expected found between 2 and 4 AU (Hayashi 1981). Note that in the active disk model we use, the effect of stellar irradiation on the temperature structure of the disk is not included.

 \begin{figure}
\par\includegraphics[width=9cm,clip]{10301fg1.eps}
\end{figure} Figure 1:

Position of the iceline $a_{\rm ice}$ as a function of the initial gas surface density  $\Sigma _{0}$ at 5.2 AU (upper three lines). It corresponds to an initial  $T_{\rm mid}$ of 170 K. The iceline is plotted for three values of $\alpha $: 0.01 (dashed line), 0.007 (solid line) and 0.001 (dotted line). The lower three lines correspond to an initial  $T_{\rm mid}$ of 1600 K, roughly the evaporation temperature of rock. The rockline  $a_{\rm rock}$ is however not taken into account in the nominal model, due to the difficulty in defining its relevant location, as disk evolution is very rapid close-in and irradiation effects might be important (cf. Paper II).

Open with DEXTER

2.2 Migration rate

The migration of the protoplanet occurs in two main regimes depending upon its mass. Low mass planets undergo type I migration (Ward 1997; Tanaka et al. 2002) which depends linearly on the body's mass. The prevalence of extrasolar planets has led us to suspect that the actual type I migration rate is probably significantly lower than currently estimated (Menou & Goodman 2003; Nelson & Papaloizou 2004). For this reason, we allow for a arbitrary reduction of the type I migration rate as calculated in Tanaka et al. (2002) by a constant efficiency factor $\f1$.

The migration type changes from type I to type II when the planet becomes massive enough to open a gap in the disk. We assume that this happens when the Hill radius of the planet becomes greater than the density scale height $\tilde H$ of the disk (Lin & Papaloizou 1986). Planetary masses where the migration regime changes can be low with such a thermal criterion only, as found also by Papaloizou & Terquem (1999) who use a similar condition. This is especially the case as due to disk evolution, the disk scale height $\tilde H$ decreases with time, so that the minimal mass needed to open a gap decreases. This effect is emphasized by the fact that our disk model currently does not include irradiation, so that especially towards the end of disk evolution, $\tilde H$ becomes smaller than in a disk including it, and smaller planets can open a gap (Edgar et al. 2007). The order of magnitude we obtain is however consistent with the one derived from Armitage & Rice (2005), since they give a gap opening condition (including the effect of viscosity) of $ M_{\rm planet}/ M_*\ga \alpha^{1/2}(\tilde H /a_{\rm planet})^2$. In our simulation, the transition typically occurs when the aspect ratio of the disk has become tiny, between 2 and 3%, meaning a transition at tens of Earth masses. We note that Crida et al. (2006) have derived a new criterion for gap opening which depends on both the disk aspect ratio and the Reynolds number. Using such a modified transition mass has some influence on the planetary formation tracks (see Sect. 5.3.5).

Type II migration (Ward 1997) itself comes in two forms: As long as the local disk mass is large compared to the planet's mass  $M_{\rm planet}$ (called ``disk dominated'' migration in Armitage 2007), the planet is coupled to the viscous evolution of the disk and its migration rate is independent of its mass. The planetary migration timescale is then the same as the gas viscous timescale (e.g. Ida & Lin 2004a). Once the local disk mass and the planet's mass become comparable, migration slows down (Lin & Papaloizou 1986) and eventually stops. Due to the inertia of the planet the disk can no longer deliver the amount of angular momentum necessary to force the planet to migrate at the gas' radial speed (e.g. Trilling et al. 1998; called ``planet dominated'' migration in Armitage 2007).

As Armitage (2007) and Thommes et al. (2008), we have found that this braking phase plays a key role in determining the final semi-major axis of massive planets. The reason for this can be seen in Fig. 2 where $2 \Sigma a^2$ is plotted as a function of time and semimajor axis for an example disk evolving under the influence of viscosity and photo-evaporation. The quantity  $2 \Sigma a^2$ serves as the measure of the local disk mass to which the planet's mass is compared (Lin & Papaloizou 1986; Syer & Clarke 1995; Armitage 2007).

 \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10301fg2.eps}
\end{figure} Figure 2:

Example of the evolution of $2 \Sigma a^2$. Seven different moments in time are plotted (time in units of Myr). At 4.6 Myr the disk has nearly vanished, and the calculations are stopped.

Open with DEXTER

The plot shows that except at the very end, $2 \Sigma a^2$ always increases with a. Initially, in the region where most giant planets begin their formation ($\sim$5-10 AU), a mass of at least a few Jupiter masses is needed to enter into the braking phase. However, after 1-2 Myr, which is the typical timescale to build protoplanets that have a sufficient mass to migrate in type II migration, a mass of the order of ${\sim}10~M_\oplus$ at ${\sim}1$ AU is already sufficient to enter into the slower planet dominated type II mode. Combined with Eq. (3), Fig. 2 also shows that once the braking starts, it steadily increases as $2 \Sigma a^2$ decreases with decreasing a and with time, while the planet's mass can continue to grow.

A consequence of the temporal evolution of the gaseous disk is the slowing down of the type II migration rate, thus providing a natural mechanism to halt planets at intermediate distances. At very small distances from the star (${\la}0.1$ AU), other, special stopping mechanisms additionally might be at work (Lin et al. 1996).

As in Alibert et al. (2005a), the migration rate during the braking phase ( $M_{\rm planet}>2 \Sigma(a_{\rm planet}) a_{\rm planet}^2$) is calculated as

\begin{displaymath}\frac{{\rm d} a_{\rm planet}}{{\rm d}t}=-\frac{3 \nu}{a_{\rm ...
...ac{\Sigma(a_{\rm planet},t) a_{\rm planet}^2}{M_{\rm planet}},
\end{displaymath} (3)

where $a_{\rm planet}$ is the semimajor axis and $\nu$ the viscosity. This is a modification of Eq. (64) in Ida & Lin (2004a) in the sense that as Edgar (2007) we use the disk properties like $\Sigma$ directly at the planet's position and not at the radius of maximum viscous coupling, assuming that for the rather high $\alpha $ we are using, wave dissipation and thus angular momentum exchange will occur essentially in the proximity of the planet (Lin & Papaloizou 1984).

A slowing down of ${\rm d}a/{\rm d}t \propto\Sigma a_{\rm planet}^2 / M_{\rm planet}$ explains naturally why larger planets should stop further out, provided that $\Sigma a^2$ increases with a. This behavior is indicated by observations (Zucker & Mazeh 2002)

2.3 Protoplanet structure and evolution

The structure of the forming planetary envelope is calculated by solving the standard equations of planet evolution as in Alibert et al. (2005a), but assuming that the luminosity of the envelope L is uniform, and equal to the accretion luminosity of planetesimals:

\begin{displaymath}%
L ={ G M_{\rm core}\dot{M}_{\rm core}\over R_{\rm core}}\cdot
\end{displaymath} (4)

In this equation, $\dot{M}_{\rm core}$ is the accretion rate of planetesimals, $M_{\rm core}$ is the mass of the core of the planet, and $R_{\rm core}$ its radius. The accretion rate is calculated according to Greenzweig & Lissauer (1992), using the same prescription for the planetesimal random velocities  $v_{\rm disp}$ as in Pollack et al. (1996). We assume an $\dot{M}_{\rm core}$ that is independent of migration, for the reasons given in Ida & Lin (2008). We therefore do not explicitly compare the timescales of planetesimal random velocity excitation with the migration timescale to see whether the protoplanet acts as a ``predator'' or ``shepherd'' (Tanaka & Ida 1999). If the type I efficiency factors $\f1$ are understood as a consequence of a ``random walk'' type migration (e.g. Nelson & Papaloizou 2004), where the single modifications of the orbit occur on a short timescale, then we are in a regime that remains to be explored in detail (Daisaka et al. 2006). In the slow, planet dominated type II regime, where planets are massive, accretion of planetesimals is usually no longer important, as ejection dominates (Sect. 5.1.4). When calculating $\dot{M}_{\rm core}$, we take into account the focusing effect of the planetary envelope (Inaba & Ikoma 2003). The iterative procedure to obtain the effective capture radius of the planet  $R_{\rm capt}$ is the same as in Pollack et al. (1996). In these calculations the effects of ablation are included (Benvenuto & Brunini 2008). We have found that ignoring the focussing effect of the envelope leads to a planetary population with similar general properties, but where roughly only half as many giant planets can form.

The assumption that L is uniform and due to planetesimal accretion only, constitutes a major difference between the models in this work and the ones of Alibert et al. (2005a): we do not take into account the exact location of the energy deposition of infalling planetesimals in the envelope, nor the energy released by the contraction of the envelope. Tests have shown that these assumptions do not strongly affect the formation, as also shown by Rice & Armitage (2003).

Solving the structure equations using the local disk temperature and pressure as external boundary conditions gives us the gas accretion rate of the planet as well as the critical mass for gas runaway accretion. Note that Miguel & Brunini (2008) have recently shown that the large uncertainties affecting the constants used in parameterized gas accretion laws as in Ida & Lin (2004a) lead to large variations of the predicted final planetary mass distribution.

The aforementioned method is valid as long as the disk can supply enough mass to keep the outer radius equal to the Hill (or the accretion) radius, i.e. if the gas accretion rate deduced from the envelope structure calculation is below the maximum rate at which gas can be delivered by the disk onto the planet.

This latter quantity can be influenced by the response of the disk on the planet's tides. Indeed, hydrodynamical simulations (Lubow et al. 1999; D'Angelo et al. 2002) have shown that, when the planet opens a gap in the disk, the accretion rate of gas is highly reduced. However, it has also been shown (see Kley & Dirksen 2006) that when the mass of the planet becomes of the order of 3-5  $M_{\textrm{\jupiter }}$ ( $M_{\textrm{\jupiter }}$ is the mass of Jupiter ${\approx}318$ $M_\oplus $), the disk-planet system can undergo a dynamic instability, leading to a substantial increase of the accretion rate of gas. For that reason, we assume in this paper that the planetary gas accretion rate in the disk limited case is simply equal to the accretion rate in the disk, namely

\begin{displaymath}%
\frac{{\rm d}M_{\rm planet}}{{\rm d}t} = \dot{M}_{\rm disk} = 3 \pi \tilde{\nu} \Sigma.
\end{displaymath} (5)

This setting constitutes another difference to the models in Alibert et al. (2005a), where we had limited the planet's accretion rate across a gap according to Veras & Armitage (2003).

As mentioned in Alibert et al. (2005a), since we are primarily interested in the mass and semi-major axis evolution of forming planets, the planet internal structure is no longer calculated once the limiting accretion rate  $\dot{M}_{\rm disk}$ is reached. Therefore, in this second phase, we can no longer explicitly compute the capture radius  $R_{\rm capt}$ of the planet. Rather, it is simply assumed to scale with the core radius, i.e. the ratio  $R_{\rm capt}/R_{\rm core}$ is kept constant. As a consequence, the amount of solids accreted after the limiting accretion rate as been reached is uncertain. This affects however mainly large, gas-dominated mass planets, and not low mass ones (like ``Hot Neptunes'').

2.4 Limitations of the model

From its conception, our model is well suited to describe the formation of giant gaseous planets, but only to a lesser extend the formation of very low mass (terrestrial) planets. This is mainly due to three assumptions that are made. Additionally, our model is a formation, not an evolutionary model (Sect. 2.4.4).

2.4.1 Initial embryo mass

First, we always assume an initial seed embryo mass of $M_{\rm emb,0}$ = 0.6 $M_\oplus $, similar to Pollack et al. (1996) or Bodenheimer & Pollack (1986), as at such a mass the opacity in the envelope becomes sufficient to justify the diffusion approximation for the radiative flux (Bodenheimer & Pollack 1986). This assumption implies that our models are only valid for planets with final masses exceeding this value by some significant margin. More quantitatively, the assumption of $M_{\rm emb,0}$ = 0.6 $M_\oplus $ is reasonable when the local isolation mass  $M_{\rm iso}$ (cf. Sect. 4.4) is significantly larger than 0.6 $M_\oplus $. For disks similar to the MMSN, $M_{\rm iso}$ is larger than 0.6 $M_\oplus $ only beyond the iceline  $a_{\rm ice}$ (e.g. Lissauer & Stewart 1993).

2.4.2 Growth after disk dispersal

Second, we stop the calculations when the gas disk has disappeared (or when the planet has migrated close to the sun, cf. below), and ignore all processes taking place later on. While for giant gaseous planets this assumption is reasonable (except for effects due to the concurrent growth of many planets, see below), inside the iceline, for disks similar to the MMSN, growth from $M_{\rm iso}\sim0.01{-}0.1$ $M_\oplus $ to the final masses occurs through giant impacts on timescales that exceed typical gas disk lifetimes by roughly one order of magnitude (Goldreich et al. 2004a). Thus, for terrestrial mass planets, growth after the dispersal of the gas disk is of greater importance. Planetary accretion proceeds at a slower pace at larger distances (Sect. 4.5), so that our assumption of an essentially completed formation of the planets at the time of disk dispersion is also not fulfilled for the formation of gas-free ice giants at large distances ( $a \ga 10{-}30$ AU, Ida & Lin 2004a).

We therefore caution that our synthetic planetary populations are incomplete for masses of less than a few earth masses for $a<a_{\rm ice}$ and less than a few 10 $M_\oplus $ for $a>a_{\rm ice}$. We stress however, that the fact that the vast majority of seed embryos do not become giant planets (cf. below) is not an artifact of the model but a consequence of the fact that the most common protoplanetary disks allow only the formation of relatively low mass planets.

2.4.3 One embryo per disk approach

An important limitation of the model that is linked to the last point is the fact that we follow, as in previous similarly detailed giant planet formation models (Pollack et al. 1996; Alibert et al. 2005a), the growth of only one embryo per disk, revolving on a circular orbit.

In reality it is clear that more than one embryo will emerge in the same protoplanetary disk, which typically begin to form in rapid succession and close proximity (Thommes et al. 2008). This multiplicity can have several kinds of effects during formation, as shown by Thommes et al. (2008): gravitational interactions between forming planets can modify their migration rate during formation, in particular by locking into resonances. Moreover, similar interactions can lead to modification of the semi-major and eccentricity distributions, also after the formation process itself (see also Adams & Laughlin 2003). In addition, planets forming in the same protoplanetary disk act as competitors for planetesimals and gas accretion, where e.g. one planet can cut off the gas supply of the other ones. The competition for planetesimal accretion was in particular addressed by Alibert et al. (2005b) in the case of the Solar System formation, and of the HD 69830 system (Alibert et al. 2006). In these two models, the internal structure of forming planets was calculated (as in the present models), but not the gravitational interactions between forming planets. In a different, and complementary approach, Thommes et al. (2008) have calculated a large set of multi-planet formation models (following a much larger number of embryos compared to the two afore-mentioned studies but without determining the internal structure of planets) accounting for both competition for gas and solids accretion, and gravitational interactions between planets. Including several embryos while keeping the detailed physics describing one single embryo is a difficult, but important step to be taken in future models.

2.4.4 Planets very close to the star

Other complications arise when planets migrate very close to their host star (${\la}0.1$ AU). This close to the star, the disk structure is more complex than further out, due for example to magnetic field effects (Lin et al. 1996) or tidal interactions (Trilling et al. 2002), which influence the formation of planets entering this zone, by altering the accretion or the migration rate (Papaloizou & Terquem 1999). Also after formation, very close-in planets can be subject to mass loss by evaporation (Vidal-Madjar et al. 2003; Baraffe 2004). Planets in close proximity to the star thus can only be described if a more detailed disk model is used, if additional stopping mechanisms for migration are taken into account and if a subsequent evolutionary model for evaporation (Baraffe et al. 2006) is included. These complications also illustrate the importance of discovering more (small) planets at ``safe'' distances from the star, as they are a more direct constraint on formation models.

Our model currently does not include any of the aforementioned effects. We therefore simply stop the calculations when a planet of mass  $M_{\rm planet}$ has migrated to  $a_{\rm touch}$, which is defined as the semimajor axis where the inner boundary of the planet's feeding zone touches the inner boundary of our computational disk at $a_{\rm min}$ = 0.1 AU, (``the feeding limit'') i.e. at $a_{\rm touch}=a_{\rm min}/(1{-}4(M_{\rm planet}/(3~M_*))^{1/3})$. If a planet has migrated to  $a_{\rm touch}$, all we can state is that its final semimajor axis would be ${\leq}a_{\rm touch}$ (it is also possible that it eventually would have fallen into the host star), and what its mass at $a_{\rm touch}$ was.

3 Monte Carlo method

The basic idea of using a Monte Carlo method to synthesize planetary populations is to sample all possible combinations of initial conditions (protoplanetary disk mass, metallicity, etc.) with a realistic probability of occurrence. This leads to all possible final outcomes of the formation process (i.e. planets) also occurring with their relative probabilities. We first explain the general six step procedure that we used.

In the first step, we identified four crucial initial conditions, and studied the domain of possible values they can take (Sect. 3.1). Some other initial conditions had to be kept constant during the synthesis of one population, for simplicity or computational time restrictions (Sect. 3.2). In the second step, we derived probability distributions for each of the four Monte Carlo variables (Sect. 4). In the third step, we draw in a Monte Carlo fashion large numbers of sets of initial conditions. The forth step consists of using the formation model for each set of initial conditions, giving the temporal evolution of the planet (formation tracks, Sect. 5.1) as well as its final properties (mass, semimajor axis, composition etc., Sect. 5.2).

Many of these synthetic planets would remain undetected by current observational techniques. So, to be able to compare the synthetic planet population with the observed one, we apply in the fifth step a detailed synthetic detection bias (Paper II). In this way, we obtain the sub-population of observable synthetic planets. Ultimately, in the sixth step, we performed quantitative statistical tests (Paper II) to compare the properties of this observable synthetic exoplanet sub-population with a comparison sample of real extrasolar planets.

3.1 Monte Carlo variables

We use four Monte Carlo variables to describe the varying initial conditions for the planetary formation process. Three describe the protoplanetary disk and one the seed embryo.

1.
The dust-to-gas ratio in the protoplanetary disk  $f_{\rm D/G}$ determines (together with  $\Sigma _{0}$) the solid surface density. Models with $f_{\rm D/G}$ between 0.013 and 0.13 were computed. Combined with the domain of $\Sigma _{0}$, this corresponds to initial solid surface densities at a0=5.2 AU of between 0.65 and 130 g/cm2. For comparison, the MMSN has a value of approximately 2.5 g/cm2 (Hayashi 1981).

2.
The initial gas surface density $\Sigma _{0}$ at 5.2 AU gives the amount of gas available. Values between between 50 and 1000 g/cm2 were used. The MMSN is estimated to have had a value of about 100-200 g/cm2 (Hayashi 1981).

3.
The last variable that characterizes a disk is the rate at which it loses mass due to photoevaporation  $\dot{M}_{w}$. For the population presented below, it was allowed to vary between 5 $\times $ $10^{-10}~M_\odot$/yr and 3 $\times $ $10^{-8}~M_\odot$/yr.

4.
The initial semimajor axis of the seed embryo within the disk, $a_{\rm start}$, is the fourth variable. It can take values of 0.1  $\leq a_{\rm start}\leq$ 20 AU.

3.2 Parameters

Some other initial conditions of the model were kept constant for all planets of a given population. We mention only the most important parameters here. More details can be found in Alibert et al. (2005a). For the nominal population discussed in Sect. 5, we use a viscosity parameter $\alpha $ for the disk model of 0.007 and an efficiency factor for type I migration $\f1$ of 0.001. The influence of these two important parameters is briefly discussed in Sect. 5.3.3, and will be further considered in forthcoming publications. In this and the companion paper the mass of the central star M* is kept constant at 1 $M_\odot $.

4 Probability distributions

In the next step we determine the probability of occurrence of a certain combination of initial conditions. In the ideal case, the probability distributions for all our variables would be derived directly from observations. Unfortunately, in reality, this is not possible either because in some cases observations do not exist or, even if they exist, a certain amount of modeling is necessary to extract the distributions from the observations.

4.1 Dust to gas ratio fD/G - [Fe/H]

To establish a link between the dust-to-gas ratio $f_{\rm D/G}$, which is the computational variable required by our model, and the corresponding observable, the stellar metallicity [Fe/H], we assume: (1) the stellar content in heavy elements is a good measure of the overall abundance of heavy elements in the disk during formation time. Support for this assumption comes from the small differences between solar photospheric and meteoritic abundances (Asplund et al. 2005); (2) a scaled solar composition and (3) a negligibly small influence of the change of the relative heavy element content on the relative hydrogen content in the comparatively small [Fe/H] domain of interest for planet formation in the solar neighborhood ($-0.5\leq$ [Fe/H] $\leq$ 0.5). Then, similar to Murray et al. (2001), we can write

\begin{displaymath}%
\frac{f_{\rm D/G}}{f_{\rm D/G,\odot}}=10^{[\rm Fe/H]}
\end{displaymath} (6)

where $f_{\rm D/G,\odot}$ is the dust to gas ratio corresponding to [Fe/H] = 0. This formula implies that we assume that iron is a good tracer of the relevant overall amount of solids available for planet formation. Robinson et al. (2006) have found that at a given iron abundance, planet host stars are enriched in silicon and nickel over stars without planets, indicating that the above relation is a simplification.

Measurements of the heavy element abundance in the Sun yield the amount (for complete condensation) of high Z material that existed initially in the form of uniformly mixed fine dust grains. However, what is relevant for our simulations is the concentration of solids in the innermost 20 AU of the disk at a later stage, namely when the dust has evolved into the 100 km planetesimals used in our model.

As has been shown by Kornet et al. (2001), the transition from the very early dust phase to the later planetesimal phase involves a number of coupled mechanisms of dust-dust and dust-gas interactions like dust settling to the midplane, dust growth by coagulation and radial drift. This leads to a redistribution of the solids within the disk, which can in turn have important effects on planetary formation (Kornet et al. 2005). The key point is that these processes lead to an increase of the solid to gas ratio in the inner ( ${\la}10{-}20$ AU) planet forming regions of the disk by advection of solids from the outer parts where a lot of mass resides (Kornet et al. 2004). The factor of increase from the initial ``(dust-)'' $f_{\rm D/G}$ to the ``(planetesimal-)'' $f_{\rm D/G}$ varies depending on the initial conditions and is not completely uniform across the inner disk, but Kornet et al. (2004) typically find values of 2 to 4.

To take this effect into account at least to first order, we set the effective planetesimal  $f_{\rm D/G}$ in the inner disk to a value about 3 times higher than the value that is inferred from the solar photosphere. Unfortunately, the latter value has been debated recently: Anders & Grevesse (1989) found a value for the protosolar Z of 0.0189, whereas a more recent study (Lodders 2003) indicates lower values (Z = 0.0149). Multiplying this value by roughly 3 gives a  $f_{\rm D/G,\odot}$ of ${\approx}0.04$, the value that we regard as nominal in our simulations. Note that if dust redistribution mechanisms are at work with similar consequences in all disks, we can still use the observed [Fe/H] distribution to scale  $f_{\rm D/G}$ for other stars.

Next, we determine the probability distribution for [Fe/H]. The distribution of the metallicity of solar like stars in the solar neighborhood has been well studied (e.g. Nordström et al. 2004), and can be well approximated by a Gaussian distribution (mean $\mu $ and dispersion $\sigma $), i.e.

\begin{displaymath}%
p([{\rm Fe/H}])=\frac{1}{\sigma \sqrt{2 \pi}} {\rm e}^{-\frac{([{\rm Fe/H}]-\mu)^2}{2\sigma^2}}.
\end{displaymath} (7)

The probability density for $f_{\rm D/G}$ is then found by using Eq. (6) and the fundamental transformation law of probabilities, $\vert p([{\rm Fe/H}]) d[{\rm Fe/H}]\vert=\vert p(f_{\rm D/G})df_{\rm D/G}\vert$ leading to

\begin{displaymath}%
p(f_{\rm D/G})=\frac{\log(e)}{f_{\rm D/G}\sigma \sqrt{2 \pi...
...{\rm e}^{-\frac{(\log(f_{\rm D/G})-\tilde{\mu})^2}{2\sigma^2}}
\end{displaymath} (8)

where $\tilde{\mu}=\log(f_{\rm D/G,\odot})+\mu$. This is simply a lognormal distribution.

Since we aim to quantitatively compare observations and theoretical calculations and compare, for example, real and theoretical detection probabilities, we need the metallicity distribution of a sample of stars that are actually included in planet searches. Such samples can have a different metallicity distribution than volume limited samples as shown by Fischer & Valenti (2005) who find that their planet search sample is shifted by ${\sim}0.1$ dex towards higher metallicities relative to a volume limited subset. The CORALIE planet search sample (Udry et al. 2000) consists of about 1650 G and K dwarfs within 50 pc, and is volume limited. In Fig. 3 the [Fe/H] distribution (obtained from calibrations to the CORALIE cross-correlation function) for about 1000 non-binary, slow rotating stars within the CORALIE search sample is plotted (Santos et al. 2003).

To derive an analytical [Fe/H] probability distribution, we have assumed that the CORALIE data is Gaussian, and performed a non-linear least square fit to obtain the parameters describing this distribution, $\mu $ and $\sigma $. The results are $\mu=-0.02\approx0.0$ and $\sigma=0.22$ (Table 1) which we have used in the calculations. The fit is also plotted in Fig. 3. It approximates the data very well between -0.45 $\leq$ [Fe/H] $\leq$ 0.4. We have repeated the same procedure also for the planet search sample, and the volume limited subsample described in Fischer & Valenti (2005). The results are also given in Table 1. For these two data sets, the Gaussian fit is however somewhat less good, especially for the volume limited sample.

 \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg3.eps}
\end{figure} Figure 3:

Solid lines: histogram and Gaussian fit to the [Fe/H] distribution of the CORALIE planet search sample. Dotted line: for comparison, distribution of photometric metallicities [Me/H] of Nordström et al. (2004). Probability densities are given, i.e. the area under the curves has been normalized to unity in each case.

Open with DEXTER

Table 1:   Parameters $\mu $ and $\sigma $ describing the Gaussian distribution of [Fe/H], for different observation samples.

The distribution derived from the CORALIE survey combined with our fiducial value for $f_{\rm D/G,\odot}$, and the range of $f_{\rm D/G}$, result in the following range of stellar metallicities: -0.49< [Fe/H] < 0.51. This range includes almost all known planet-hosting main sequence stars in the solar neighborhood. Finally, we note that despite the complications of linking  $f_{\rm D/G}$ to [Fe/H] mentioned here, one can reasonably assume that it is the variable that is most directly constrained by observational data.

4.2 Gas surface density $\Sigma _{0}$ - disk mass Mdisk

The probability distribution of gas surface densities  $\Sigma _{0}$ can also be at least partially inferred from observations of star forming regions. The flux density of thermal continuum emission originating from cold dust orbiting young stellar objects (YSO) at millimeter and submillimeter wavelengths allows an estimate of the total dust disk mass (Beckwith et al. 1990). By adding a gas content typical of the interstellar medium (Beckwith & Sargent 1996; Andrews & Williams 2005), total disk masses  $M_{\rm disk}$ are found. The study of many YSO forming concurrently in a given star formation region then gives distributions for  $M_{\rm disk}$.

While the details of these distributions differ across known star formation regions, most have roughly the shape of a lognormal distribution for $M_{\rm disk}$, i.e.  $\log(M_{\rm disk})$ is Gaussian distributed with a mean $\mu $ and a standard deviation $\sigma $ (Andrews & Williams 2005; Robinson et al. 2006). To compute the parameters $\mu $ and $\sigma $ we use Beckwith & Sargent (1996), who give in their Fig. 5 histograms for the distribution of total disk gas masses in the Taurus-Auriga and Ophiuchus star formation region, and perform as for [Fe/H] a nonlinear least square fit to obtain the two sets of parameters.

 \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg4.eps}
\end{figure} Figure 4:

Histograms and fits to circumstellar disk masses  $M_{\rm disk}$ for Ophiuchus (solid lines) and Taurus-Auriga (dotted lines). The actual gas masses in the computational domain out to 30 AU are between 0.004 and 0.09 $M_\odot $.

Open with DEXTER

Table 2:   Parameters describing the lognormal distribution of circumstellar disk masses $\mu $ and $\sigma $, from different sources.

The observational data, as well as the fits are plotted in Fig. 4, while the values for $\mu $ and $\sigma $ are given in Table 2. For both star formation regions, the disk mass distribution can be reasonably well fitted by Gaussians. The resulting distributions are however very broad which leads to significant probabilities for disks with masses in excess of 0.1 $M_\odot $ (or even 0.3 $M_\odot $). For typical thermodynamic conditions, such disks are self-gravitationally unstable in the outer regions (e.g. Ida & Lin 2004a). It is likely that this result is due to the presence of the remains of the envelope from which the star formed, contributing to the observed flux (Andrews & Williams 2005).

The large spread in masses is also due to the spread in age of the observed YSOs in one cluster (Robinson et al. 2006) and thus an evolutionary effect. May be the spread additionally enhanced by the spread of stellar masses, if the the mass of stars and disks are correlated.

Since our planet formation model includes disk evolution, we require a mass distribution taken as early as possible, before strong evolutionary effects have taken place. Therefore the distribution of Ophiuchus, which is about 2-3 times younger than Taurus-Auriga (White & Hillenbrand 2004), is more appropriate for our goal and is used as the nominal distribution for the simulations.

When converting the observed quantity (disk masses) to the corresponding numerical quantity in the model (the initial gas surface density  $\Sigma _{0}$) we are confronted with a situation similar to the [Fe/H] to $f_{\rm D/G}$ conversion (Matsuo et al. 2007): Observed protoplanetary disks have physical radii  $a_{\rm phys}$ typically of a few hundred AU (e.g. Beckwith & Sargent 1996). For computational reasons however we simulate only the inner part of the disk out to $a_{\rm max}=30$ AU. Putting the observed disk masses into such a small disk would lead to unrealistically high surface densities. Therefore, when converting the observed disk masses to  $\Sigma _{0}$, we assume that the mass is contained in a disk with a physical radius $a_{\rm phys}=300$ AU of which we simulate only the innermost 30 AU. The outcomes of our simulations are not very sensitive to the exact value of the assumed physical outer radius, as the disk mass scales only with  $\sqrt{a_{\rm phys}}$. Thus, the probability density for  $\Sigma _{0}$ is

\begin{displaymath}%
p(\Sigma_{0})=\frac{\log(e)}{\Sigma_{0}\sigma \sqrt{2 \pi} } {\rm e}^{-\frac{(\log(\Sigma_{0})-\tilde{\mu})^2}{2\sigma^2}}
\end{displaymath} (9)

where $\tilde{\mu}=\mu-\log\left(4 \pi
a_0^{3/2}(\sqrt{a_{\rm phys}}-\sqrt{a_{\rm min}} ~)\right)$, which is of a lognormal type (Bronstein et al. 1999).

With the boundaries of the computational disk at $a_{\rm min}=0.1$ AU and $a_{\rm max}=30$ AU and given the range of $\Sigma _{0}$ of 50-1000 g/cm2, we cover a range of disk gas masses in the computational domain between 0.004 and 0.09 $M_\odot $. These masses should be stable against self-gravitional collapse (Mayer et al. 2004).

  4.3 Photo-evaporation rate $\dot{M}_{w}$ - disk lifetime tdisk

The values for the photo-evaporation rate $\dot{M}_{w}$ determine, together with the viscosity parameter $\alpha $, the timescale  $t_{\rm disk}$ on which the gas disk disappears. Constraints on its possible values can therefore also be inferred from observation.

Haisch et al. (2001) have shown by near-infrared observations of hot dust that the circumstellar disk fraction in young stellar clusters is a roughly linearly decreasing function of age with essentially all disks having disappeared after $\sim$6 Myr. They argue that this timescale should be a tracer of the evolution of the bulk of the disk material as well, especially also of the gas disk. A very similar time dependence is also found by Hillenbrand (2005).

Assuming a uniform distribution in log for $\dot{M}_{w}$, we have determined the bounds of the $\dot{M}_{w}$ distribution by requiring that the lifetime of our synthetic disks for a fixed value of the viscosity parameter $\alpha $ and for the Ophiuchus distribution of initial disk masses follows the observed disk lifetime distribution found by Haisch et al. (2001). To determine the bounds, we first generated a distribution of disk masses as described in Sect. 4.2. We then set some trial boundaries for the distribution of  $\dot{M}_{w}$ and drew values uniformly in log inside these limits. We can then calculate the resulting lifetime $t_{\rm disk}(\alpha, \Sigma_{0}, \dot{M}_{w})$ for each synthetic disk, so that we get a distribution of synthetic  $t_{\rm disk}$. The bounds were then adjusted in an iterative fashion until the observed and the synthetic disk lifetime distribution are similar.

It is found that a reasonable fit is obtained by allowing a range between 5 $\times $ 10-10 and 3 $\times $ $10^{-8}~M_\odot$/yr for $\alpha=7$ $\times $ 10-3. For other $\alpha $, the same iterative procedure was repeated. As one expects, the higher $\alpha $, the lower the $\dot{M}_{w}$ boundaries necessary in order to reproduce the Haisch et al. (2001) values. We have for example found that for $\alpha=10^{-2}$, bounds of 1 $\times $ 10-10 and 1.5 $\times $  $10^{-8}~M_\odot$/yr are appropriate. Such values are compatible with the ones found elsewhere (Armitage et al. 2003).

In Fig. 5 the fraction of stars in the model with remaining disks using the nominal $\dot{M}_{w}$ distribution and two comparison cases are plotted, together with the fit of Haisch et al. (2001). The model and the observed distribution have a similar decrease of the remaining disk fraction if the systematic errors of the ages in the observational data are taken into account.

 \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg5.eps}
\end{figure} Figure 5:

Fraction of stars in the model possessing a gaseous disk as a function of time for a uniform distribution in log of $\dot{M}_{w}$ between 5 $\times $ 10-10 and 3 $\times $ $10^{-8}~M_\odot$/yr (thin solid line). The thick solid line is the fit of Haisch et al. (2001) to the observed JHKL excess/disk fraction as a function of mean cluster age. The error bar at the top right is the overall systematic uncertainty in age of these observations, also from Haisch et al. (2001). For comparison, the fraction is also given if two different distributions of  $\dot{M}_{w}$ are used in the model: uniformly in log between 1 $\times $ 10-8-1 $\times $ $10^{-7}~M_\odot$/yr (dash-dotted line), leading to disk lifetimes that are clearly too short, or uniform in log in 1 $\times $ 10-10-1 $\times $ $10^{-9}~M_\odot$/yr, leading to disk lifetimes that are too long (dotted line).

Open with DEXTER

It is clear that the distribution of $\dot{M}_{w}$ could in reality be much more complicated and that also the assumption of a temporally constant UV flux might not be justified. One should for example, when external UV sources are considered, take into account the stellar environment, as O or B stars in rich clusters can greatly enhance the far-ultraviolet flux (Adams et al. 2004; Armitage 2000).

The boundary conditions of our disk model (Sect. 2.1) are similar to the ones presented in Alibert et al. (2005a). In particular, at the outer disk boundary, we assume no mass influx from the outer parts of the disk. This assumption could lead to an underestimation of the disk dispersal time, since the outer parts of the disk could act as a mass reservoir. However, as explained above, we adjust the photoevaporation rate in order to obtain disk lifetimes similar, by construction, to observed ones. With the initial power law surface density, the disk accretion rate is found to decrease towards the outer parts of the disk, reaching values below $10^{-9}~M_\odot$/yr, which is of the order of, or lower than, the photoevaporation rate. Therefore, a large part of the mass present in the outer disk is photoevaporated rather than accreted onto the inner parts. This means that if we would include inflow, we would have to increase the required photoevaporation rates to obtain disk lifetimes consistent with observations, but otherwise our results would remain similar.

4.4 Embryo start position a start

The starting position of the planetary embryos $a_{\rm start}$ inside the disk cannot be constrained by observations. Therefore, we derive a probability distribution using theoretical arguments only.

The derivation can be made relatively easily if one assumes that during the early phases of planetary accretion, radial motions can be neglected. In this case, it follows (e.g. Lissauer & Stewart 1993) from the restricted 3-body problem that an embryo can accrete background planetesimals only within its feeding zone, which has a half width of BL times ( $B_L\approx3{-}5$) its Hill sphere radius $R_{\rm H}$. The Hill sphere radius of the planetary embryo at a semimajor axis a and a mass  $M_{\rm emb,0}$ is given by $R_{\rm H}=\left(M_{\rm emb,0}/(3~M_*0\right))^{1/3} a$, which is, for fixed $M_{\rm emb,0}$, $\propto$a. Thus, as confirmed by many different numerical simulations (e.g. Weidenschilling et al. 1997; Kokubo & Ida 2000), runaway bodies emerge with relative separations of their semimajor axis $\Delta$ proportional to their semimajor axes a, so $\Delta /a$ = const. This situation also prevails later during oligarchic growth stages, although the numerical value of BL increases (Ida & Lin 2004a).

The probability p(a) of starting in an interval ${\rm d}a$ is inversely proportional to the spacing $\Delta$, so

\begin{displaymath}%
p(a) {\rm d}a \propto \frac{{\rm d}a}{\Delta} \propto \frac{{\rm d}a}{a} = {\rm d} \log (a) \propto {\rm const.},
\end{displaymath} (10)

which means that the probability distribution for the starting position is uniform in log, as assumed also by Ida & Lin (2004a).

While the distribution of the initial position is clear, there remains the issue that we start our simulations with a seed embryo of mass $M_{\rm emb,0}=0.6$ $M_\oplus $. In order to be self-consistent, we have to require that this amount of mass is available at the starting position. In other words, we have to ensure that for a set of initial values ( $f_{\rm D/G}$, $\Sigma _{0}$ and  $a_{\rm start}$) there is indeed enough mass available to build up such a seed. For this, we consider the amount of heavy elements in the planet's feeding zone $2\pi a_{\rm start}2 B_L R_{\rm H}\Sigma_{\rm D}$ to calculate the isolation mass  $M_{\rm iso}$ (Lissauer & Stewart 1993):

\begin{displaymath}%
M_{\rm iso}= \frac{(4 \pi B_L a_{\rm start}^2 \Sigma_{\rm D})^{3/2}}{(3 M_*)^{1/2}}\cdot
\end{displaymath} (11)

For the disk profile we use this becomes

\begin{displaymath}%
M_{\rm iso}=\frac{(4 \pi B_L f_{\rm D/G}f_{\rm R/I}\Sigma_{0}a_0^{3/2} a_{\rm start}^{1/2})^{3/2}}{(3 M_*)^{1/2}},
\end{displaymath} (12)

where $f_{\rm R/I}$ is itself a function of $\Sigma _{0}$ (and $\alpha $, Fig. 1).

If the isolation mass calculated from Eq. (12) with $B_L=2\sqrt{3}$ (Lissauer & Stewart 1993) for a set of $f_{\rm D/G}$, $\Sigma _{0}$ and $a_{\rm start}$ exceeds 0.6 $M_\oplus $, we conclude that these are self-consistent initial conditions for the formation model.

As $M_{\rm iso}$ $\propto$ $a_{\rm start}^{3/4}$, the isolation mass criterion sets an inner border for the possible domain of $a_{\rm start}$. We also require that the starting time of the embryo  $t_{\rm start}$ is smaller than the disk lifetime  $t_{\rm disk}$ (cf. the next Sect. 4.5). In practice, we first calculate for a given $f_{\rm D/G}$, $\Sigma _{0}$ and $\dot{M}_{w}$ the range of $a_{\rm start}$ in which $M_{\rm iso}\geq M_{\rm emb,0}$ and $t_{\rm disk}\geq t_{\rm start}$, and then draw  $a_{\rm start}$ from that range.

 \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg6.eps}
\end{figure} Figure 6:

Probability distribution for $a_{\rm start}$. There is a marked peak around 3 AU. Outside this distance, the probability distribution falls off as expected for a distribution that is uniform in log. The fraction inside about 2 AU is due to disks with the ability to form seeds sufficiently massive ( ${\geq }M_{\rm emb,0}$) inside the iceline as well.

Open with DEXTER

In Fig. 6 the resulting distribution of $a_{\rm start}$ is plotted. Its particular shape is the result of the combination of the probability distribution for  $a_{\rm start}$ (uniform in log) and the criteria that $t_{\rm disk}\geq t_{\rm start}$ and $M_{\rm iso}\geq M_{\rm emb,0}$. The fact that the latter is fulfilled only beyond the iceline for most disks, and the property of the uniform-in-log distribution to emphasize the smallest possible values, lead to the peak near 3 AU. Further out, the probability decreases, a consequence of the uniform-in-log distribution and the increasing $t_{\rm start}$. Only disks with a high surface density (i.e. a concurrently high  $\Sigma _{0}$ and $f_{\rm D/G}$) which can form sufficiently massive seeds ( ${\geq }M_{\rm emb,0}$) inside the iceline contribute to the distribution inside ${\sim}2$ AU. The transition to higher probabilities is not sharp, as the iceline is itself a function of $\Sigma _{0}$ (Fig. 1). Note that via the $M_{\rm iso}$ and $t_{\rm disk}$ criteria, the shape of the $a_{\rm start}$ distribution depends on the $f_{\rm D/G}$, $\Sigma _{0}$ and $\dot{M}_{w}$ distributions.

4.5 Embryo start time t start

Depending upon initial conditions, the time needed to form just the initial seed embryo of 0.6 $M_\oplus $ is not short compared to the overall disk lifetime. We take this effect into account by introducing a time delay $t_{\rm start}$ between the beginning of the disk evolution and the time at which the embryo is placed in the disk. This time delay is calculated in a deterministic way as a function of $f_{\rm D/G}$, $\Sigma _{0}$ and $a_{\rm start}$. It is therefore not an independent Monte Carlo variable. Similarly to  $a_{\rm start}$, the time delay cannot be constrained by observations. We therefore again rely on theoretical arguments.

During the early stages of planetary accretion, the process of embryo growth proceeds mainly by accumulation of background planetesimals and can be well described in the two-body approximation (Kokubo & Ida 2002). The accretion rate of an embryo growing from background planetesimals is thus (e.g. Lissauer & Stewart 1993)

\begin{displaymath}%
\frac{{\rm d} M_{\rm emb}}{{\rm d}t}\simeq \pi R_{\rm emb}^...
...m D}\left( 1+\frac{v_{\rm esc}^2}{v_{\rm disp}^2} \right)\cdot
\end{displaymath} (13)

Here, $R_{\rm emb}$ is the radius of the embryo, $v_{\rm esc}$ the escape velocity from it and $v_{\rm disp}$ the typical random velocity of the planetesimals.

As the embryo grows, the solid surface density of planetesimals $\Sigma_{\rm D}$ must decrease correspondingly (Thommes et al. 2003):

\begin{displaymath}%
\frac{{\rm d} \Sigma_{\rm D}}{{\rm d}t} = - \frac{(3 M_*)^{...
..._L M_{\rm emb}^{1/3}}\frac{{\rm d} M_{\rm emb}}{{\rm d}t}\cdot
\end{displaymath} (14)

Looking at the dependence of the embryo growth time on semimajor axis, it is found (Steward & Ida 2000; Weidenschilling & Davis 2001) that $t_{\rm start}$ increases significantly with distance to the star. The reason is that both the surface density of planetesimals and the frequency of collisions decrease with distance, as the pace on which the later occur scales with $\Omega$ (Kokubo & Ida 2002). This increase of the growth timescale with distance is well reproduced in numerical simulations where a ``wave'' of embryo growth is seen to propagate from the inner parts to the outer parts of the disk (Weidenschilling & Davis 2001; Thommes et al. 2003).

Thus, building up a body of fixed mass, in our case the initial seed, takes longer for larger  $a_{\rm start}$, a trend that holds everywhere but at the iceline where $\Sigma_{\rm D}$ increases suddenly. To derive the corresponding $t_{\rm start}$, we integrate numerically Eqs. (13) and (14) until the mass has reached $M_{\rm emb,0}=0.6$ $M_\oplus $. As in the formation model, we use the procedure of Pollack et al. (1996) to calculate the planetesimal random velocity $v_{\rm disp}$ and assume a constant planetesimal size of 100 km. As initial values at t=0, we use a mass of $0.66m_{\rm pla}^{3/5}M_{\rm iso}^{2/5}$ (Chambers 2006).

Compared to the formation model which calculates the evolution of the protoplanet starting with $M_{\rm emb,0}=0.6$ $M_\oplus $, we neglect several factors in the calculation of  $t_{\rm start}$: first, migration is neglected. Tests have shown that with the type I migration efficiency factor $\f1$ we use for the nominal case, migration of bodies this small is indeed negligible. Second, we use just the simple particle in a box law for the gravitational focusing from Eq. (13) instead of the full three body results from Greenzweig & Lissauer (1992), and third, we neglect the increase of the capture radius  $R_{\rm capt}$ due to the gaseous envelope. The latter assumption is appropriate given the small seed masses and the assumed planetesimal size as shown by our own tests but also by Kornet & Wolf (2006) or Fortier et al. (2007).

 \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg7.eps}\hspace*{7mm}
\includegraphics[width=8.5cm,clip]{10301fg8.eps}
\end{figure} Figure 7:

Snapshots of the embryo mass (solid line) as a function of semimajor axis at four moments in time for two different solid surface densities. The dashed line is the isolation mass. The dotted line is $M_{\rm emb,0}=0.6$ $M_\oplus $. The initial solid surface density at 1 AU is 7 g/cm2 ( left panel) and 35 g/cm2 ( right panel). It should be kept in mind that this kind of calculation is needed to generate the start time  $t_{\rm start}$ when the embryo is put into the formation model. The real evolution of the solid core for M>0.6 $M_\oplus $ is in general much more complex than plotted here. In this figure, we have continued the calculations up to the isolation mass to allow comparison with other models.

Open with DEXTER

Figure 7 illustrates the result of the numerical integrations to obtain  $t_{\rm start}$. For a solid surface density at 1 AU of 7 g/cm2 ($\sim$MMSN) and for 35 g/cm2 ($\sim$$\times $ MMSN) four snapshots in time of the embryo mass as a function of semimajor axis are plotted.

To allow comparison with other models (Thommes et al. 2003; Ida & Lin 2004a; Chambers 2006), we have continued the integration up to  $M_{\rm iso}$ for this plot instead of stopping at  $M_{\rm emb,0}$ as we do when generating the initial conditions. The dotted line in Fig. 7 is  $M_{\rm iso}$ and shows that for the MMSN, the earliest possible time to start is somewhat more than 1 Myr, just outside the iceline. Then, the domain of possible start positions grows only slowly to larger semimajor axes. After 10 Myr, when the gas disks will have disappeared (cf. Fig. 5), the largest possible start position is still only 6-7 AU. In contrast, in the 5 $\times $ MMSN case, $M_{\rm iso}$ is larger than $M_{\rm emb,0}$ inside the iceline too, and the earliest embryo can start at around 0.8 AU, slightly before 0.1 Myr. At later times, embryos can start from all semimajor axes.

Compared to the results of Thommes et al. (2003), our calculations show a faster growth of the cores, especially at large distances, which is due to the different way of computing planetesimal eccentricities and inclinations in their model, as illustrated by Fortier et al. (2007). Compared to Ida & Lin (2004a), the results are quite similar, even if core growth proceeds at large orbital distances somewhat faster in our model. Compared to Chambers (2006) one finds that core growth in our model is faster than in his simple equilibrium model, but slower than in his complete model that is considerably more complex, including e.g. planetesimal fragmentation.

As mentioned in Sect. 4.4, we only start embryos in that part of the disk where $M_{\rm iso}\geq M_{\rm emb,0}$ and $t_{\rm disk}\geq t_{\rm start}$. The latter condition gives an outer bound for possible starting positions. The reasoning behind it is that if one of the numerous planetary seeds can form while the disk is still present, it would have done so, and that it is a candidate to eventually become a giant planet observable today. In other parts of the disk, seed embryos also form, but they remain very small during the presence of the gaseous disk. Thus, we aim at minimizing the negative side effects of having only one seed per disk on the population of giant planets, but at the same time make our populations incomplete at small masses (cf. Sect. 2.4).

For a significant fraction (${\sim}28$%) of the sets of initial conditions we draw, one or both of the two aforementioned conditions cannot be fulfilled anywhere in the disk, namely when $f_{\rm D/G}$ and/or $\Sigma _{0}$ come from the low tail of their distributions, while  $\dot{M}_{w}$ is high. In such cases, no calculations were made, but we keep the record of the corresponding initial conditions where the formation of sizable planets is not possible and correct for them when calculating for example overall detection probabilities (Paper II).

5 Results

Once all Monte Carlo variables have been drawn, the next step consists of computing the formation of the planet corresponding to these initial conditions. This process can be illustrated by means of formation tracks in the mass-distance plane. Except where otherwise stated, all results are obtained for a population with $\alpha=$ 0.007 and $\f1=0.001$. The reason for this choice is that the resulting sub-population of observable synthetic planets reasonably well reproduces the observed population (Paper II).

5.1 Planetary formation tracks

Figure 8 shows formation tracks of about 1500 randomly chosen synthetic planets. The tracks lead from the initial position at $a(t=0)=a_{\rm start}$ and the fixed $M(t=0)=M_{\rm emb,0}$ to the final position marked by a large black symbol when planet growth and migration stops. The color of the track indicates the migration mode: Red for type I migration, blue for ordinary (disk dominated) type II migration and green for the braking phase. In this phase, planet dominated type II migration occurs (Eq. (3)) and the planetary gas accretion rate is given by the rate at which the disk can supply gas (Eq. (5)).

 \begin{figure}
\par\includegraphics[width=16cm,clip]{10301fg9.eps}
\end{figure} Figure 8:

Planetary formation tracks in the mass-distance plane. The large black symbols show the final position of a planet. The shape of the symbols is explained in the text. Planets reaching the feeding limit at  $a_{\rm touch}$ (indicated by the long dashed line) have arbitrarily been set to 0.1 AU. The short dashed lines have a slope of $-\pi $ (discussion in Sect. 5.1.3). Each track is color-coded according to the migration mode, and small black dots are plotted on the tracks every 0.2 Myr to indicate the temporal evolution of a planet.

Open with DEXTER

Even if the tracks show a great diversity, one can distinguish groups of planets with similar tracks. These groups are due to different formation stages that planets might undergo. In the next sections, we study representative tracks of four such groups.

5.1.1 Tracks of ``failed cores''

During the first stage of formation at low masses, type I migration (red) occurs. Since for this example population type I migration is very slow ($\f1=0.001$), the tracks are almost vertical. Planets that have migrated as type I only are represented by filled circles in Fig. 8.

For most embryos, this first stage is also the final one. Their evolution stops at low masses because most initial conditions do not allow the formation of more massive planets during the lifetime of the disk. Therefore, most seeds (50-75%, see Paper II) contribute to building up a large population of ``failed cores'' with $M \sim 1{-}10$ $M_\oplus $ which, from the point of view of giant planet formation, failed to accrete a significant amount of gas. The population synthesis calculations of Ida & Lin (2004a,2008) also contain a large sub-population of low mass planets. This is compatible with the non-detection of giant planets around 90 to 95% of nearby solar like stars.

In Fig. 9, left panel, exemplary formation tracks for a number of such planets are plotted. As expected from Eq. (12) for  $M_{\rm iso}$, ``failed cores'' can reach larger masses at larger distances. The right panel of Fig. 9 shows the temporal evolution of the mass and semimajor axis of one typical ``failed core''. This seed starts at $a_{\rm start}=3.7$ AU in a disk with $f_{\rm D/G}=0.028$ ([Fe/H] = -0.15) and $\Sigma_{0}=165$ g/cm2. This initial position is situated not far outside the iceline. For such a solid surface density, forming the initial seed takes a significant amount of time (cf. Fig. 7), namely about 1.1 Myr.

 \begin{figure}
\par\mbox{ \resizebox{9cm}{8.7cm}{\includegraphics{10301f10.eps}}...
...ox{9cm}{!}{\includegraphics[angle=270,origin=rb]{10301f11.eps}} }
\end{figure} Figure 9:

Planetary formation tracks in the mass-distance plane for ``failed cores'' ( left panel). The black points represent the final position of the planets. The thick line starting at 3.7 AU is the track of one prototypical example, for which the right panels shows the temporal evolution. Its final position is represented by a large square. In the right panel, the total mass M (solid line), the mass of accreted solids MZ (dashed line) and the mass of the envelope  $M_{\rm env}$ multiplied by a factor of 10 for better visibility (dotted line) are plotted (scale on the left) as a function of time t. The temporal evolution of the planet's semimajor axis a is also plotted (dash-dotted line, scale on the right).

Open with DEXTER

As is shown in the right panel of Fig. 9, the core then quickly accretes all planetesimals in its reach. Gas accretion is of negligible importance. At about 1.2 Myr, the mass of the core approaches the local isolation mass[*]. For the remaining 0.2 Myr of evolution, the core grows only very slowly. The envelope now becomes more massive, due to the reduced luminosity of the core. The evolution of this planet corresponds to the two first phases described by Pollack et al. (1996), with the difference that further evolution is inhibited by the dispersion of the protoplanetary nebula after 1.45 Myr. At this time, we are left with a ``failed core'', consisting of about 3.6 $M_\oplus $ of heavy elements, and ${\sim}0.1$ $M_\oplus $ of gas. The extent over which migration occurred is tiny because of $\f1=0.001$, roughly 0.004 AU, much less than the extent of the planet's Hills radius. The fact that further growth is inhibited by the disappearance of the gaseous disk is characteristic for this type of planet.

The vast sub-population of ``failed cores'' is not identical to the final terrestrial planet population, expected to be located in a similar a-M region. Rather, they represent an earlier moment in evolution. ``Failed cores'' are formed from one large embryo accreting small field planetesimals while the gas disk is still present. Terrestrial planets on the other hand derive their final properties from giant impacts between bodies of a similar size (several ``failed cores'') on much longer timescales, a phase missing in our model. We expect that after disk dispersal, all the ``failed cores'' of one disk would start to interact gravitationally, leading to scattering, ejections and collisions, until the remaining planets have settled into stable orbits (e.g. Ford & Chiang 2007; Thommes et al. 2008).

  5.1.2 Tracks of ``horizontal branch'' planets

In some other cases the core grows so large (and does so sufficiently quickly) that the planet can open a gap in the gas disk long before the latter disappears. At this point, the migration mode changes to the disk dominated type II migration (blue lines in Fig. 8), which is the second phase.

After a short transitional phase, planets starting inside about 4-6 AU begin to move, in disk dominated type II migration, along nearly, but not completely horizontal tracks at M $\sim$ 7-30 $M_\oplus $. These nearly horizontal tracks are clearly seen in Fig. 8, forming a ``horizontal branch'' of planets. We identify the planets having passed through such a phase during their formation a posteriori by the condition that while in type II migration, one finds $\frac{{\rm d}\log M}{{\rm d}\log a}<0.1$. Physically this means that migration occurs on a significantly shorter timescale than accretion. Planets having passed through the ``horizontal branch'' have their final position marked by triangles in Fig. 8.

Figure 10, left panel, shows some exemplary formation tracks of planets that stay in the ``horizontal branch'' until the gas disk has disappeared (so that they end up at intermediate distances), or until they reach the feeding limit (so that they have a final position ${\la}0.1$ AU).

 \begin{figure}
\par\mbox{\resizebox{9cm}{8.55cm}{\includegraphics{10301f12.eps}} \resizebox{9cm}{!}{\includegraphics{10301f13.eps}} }\end{figure} Figure 10:

As Fig. 9, but for planets of the ``horizontal branch''. In this case, the prototypical track (thick line, and large square at the final position) starts at $a_{\rm start}=4.7$ AU. It ends as a ``Hot Neptune'' planet at the feeding limit at 0.1 AU. Its temporal evolution is shown in the right panel.

Open with DEXTER

The prototypical example of the right panel of Fig. 10 is formed in a disk whose mass is similar to the mean values adopted in this work ( $\Sigma_{0}=270$ g/cm2). The dust to gas ratio is with $f_{\rm D/G}=0.03$ somewhat smaller than the mean value. This leads to a formation time of the initial seed of about 1 Myr. For this disk, the iceline is located at 4.2 AU, therefore the planet accretes icy planetesimals at the beginning of its formation. During the first 1.6 Myr, the formation is similar to the one described in classical core accretion papers, namely a rapid core formation (up to about 8 $M_\oplus $, the isolation mass, in 0.15 Myr) followed by a phase of low mass growth, similar to phase 2 of Pollack et al. (1996). Just before t=1.7 Myr, however, migration switches to type II, due to the concurrent growth of the planet and the decrease of the disk scale height with time, so that a planet of a relatively low mass also can open a gap in the disk.

In the population presented here, we have reduced type I migration by a large factor. Therefore changing from the (strongly reduced) type I to (normal) type II results in a net increase of the migration rate and the planet moves into regions of the disks that have not yet been depleted in planetesimals (Alibert et al. 2004). This significantly increases the core growth and and hence its luminosity. The latter translates to a slight decrease of the envelope mass at 1.7 Myr.

Shortly after switching migration type, the planet crosses the iceline, which reduces the solid accretion rate (see the small kink in the mass lines just after 1.7 Myr). During the remaining accretion inside the iceline, another ${\sim}17$ $M_\oplus $ of rocky planetesimals are collected. Comparing this to the total amount of heavy elements initially inside the iceline (${\sim}20$ $M_\oplus $) shows that the planet is quite efficient in emptying the planetesimal disk. As a consequence, the final core mass is approximately the sum of the mass collected while passing through the ``horizontal branch'', plus the mass of the icy planetesimals accreted during the initial in-situ growth phase.

At roughly 1.9 Myr, the local disk mass becomes comparable to the planet's mass so that the migration rate starts to slow down (Fig. 2) and the planets starts to accrete gas. However, migration continues (albeit at a reduced rate) and at t=2.2 Myr the planet enters the feeding limit at roughly 0.1 AU, where we stop the calculations.

Contrary to the two remaining types of representative tracks described below, the gas accretion rate of ``horizontal branch'' planets is governed during their entire formation by the planet i.e. the maximum accretion rate limited by the disk is never reached. Since the gas accretion rate is moderate (due to the fact that the core luminosity is most of the time quite large), the planet at the end of its formation is made of a large core (${\sim}26$ $M_\oplus $) and a significant, but still much smaller envelope (${\sim}6$ $M_\oplus $). We find that the sub-population of close-in, low mass planets ( $M \la 30{-}40$ $M_\oplus $) is characterized by a ratio  $M_{\rm env}/M_{\rm core}$ that varies between $\sim$0.02 (at $M\sim10~M_\oplus$) and $\sim$0.3 (at $M\sim40$ $M_\oplus $), i.e. these planets have a structure roughly comparable to Neptune. This could be different for close-in Neptune mass planets formed through giant impacts between initially smaller bodies (Ida & Lin 2008), especially as the impacts could remove their tenuous gaseous envelopes. Transiting Neptune mass planets around solar like stars can serve as a diagnostic to distinguish these different formation channels. The recently detected transiting ``Hot Neptune'' HAT-P-11b (Bakos et al. 2009) has a radius compatible with a rock/ice core with a 10% H/He envelope, whereas a pure rock/ice planet (as well as a miniature gaseous planet) are excluded (Bakos et al. 2009). This suggest that this planet was formed in a similar way as described here.

The ``horizontal branch'' is thus the ``conveyor belt'' by which Neptune-like planets are transported close to the star. These ``Hot Neptunes'' are a sub-population of planets that high precision radial velocity surveys (using e.g. HARPS, see Lovis et al. 2006) now find in increasing numbers. Note that subsequent evolutionary effects, namely evaporation, can sometimes significantly modify the structure of these bodies (Baraffe et al. 2006).

5.1.3 Tracks of ``main clump'' planets

The fact that the tracks in the ``horizontal branch'' are not completely horizontal, i.e. that growth in mass continues in this phase, is important for the further evolution of the third group of planets, the planets of the ``main clump''.

These are planets with final masses mostly between the mass of Saturn and three Jupiter masses, and semimajor axes mainly between ${\sim}0.3$ and 2 AU (see Fig. 13). For these planets, the core grows to a size that triggers runaway gas accretion while the planet collects solids as it passes through the ``horizontal branch''. Once runaway is triggered, the gas accretion rate increases very rapidly. With its rapidly growing mass, the planet soon exceeds the local disk mass and migration enters the planet dominated type II regime. The planet leaves the ``horizontal branch'' upwards in mass, thereby starting its third phase of formation (plotted in green in Fig. 8).

 \begin{figure}
\par\mbox{\includegraphics[width=9cm,height=8.8cm]{10301f14.eps} \includegraphics[width=9cm]{10301f15.eps} }
\end{figure} Figure 11:

As Fig. 9, but for planets passing through the ``horizontal branch'' to become members of the ``main clump''. Here, the prototypical embryo (thick line, square at the final position, temporal evolution in the right panel) eventually leads to a Saturnian planet situated at 0.9 AU.

Open with DEXTER

During this final braking phase, the planets migrate at the reduced type II rate (Eq. (3)), and accrete gas at a rate given by the disk's evolution (Eq. (5)). This has the interesting consequence that if we combine these equations, we find that planets move on formation tracks which are straight lines in the $\log(a)-\log(m)$ plane with a slope ${\rm d} \log M/{\rm d} \log a= -\pi$. This behavior is clearly visible in the formation tracks. In fact, we have used the criterion $\vert\ {\rm d}\log M/{\rm d}\log a+\pi\vert<0.1$ to identify the braking phase a posteriori. The evolution along these straight lines slows down in time, as can be seen by the increasing number of small black ticks on the track near the final position of the planet (Fig. 8). This is a consequence of the concurrent decrease of the gas accretion rate due to disk evolution, and the slowing down of migration. At the end, the planet has accreted all the outer disk gas that has not been photo-evaporated.

Examples of planets that undergo such a three staged evolution are plotted in Fig. 11. Two specific planets have reached the feeding limit near 0.1 AU, illustrating how Pegasi planets form, even though our model does not further treat the formation of ``Hot'' planets once they have migrated to  $a_{\rm touch}$.

The prototypical planet of the ``main clump'' in the right panel of this figure is formed in a disk with $f_{\rm D/G}$ = 0.02, i.e. [Fe/H] = -0.3, and $\Sigma_{0}=280$ g/cm2. The planet starts its formation beyond the iceline ( $a_{\rm start}=6$, $a_{\rm ice}=4.3$ AU), and the beginning of its formation (up to 4 Myr) is similar to that of the ``horizontal branch'' planet: the planets empties its feeding zone, reaching an isolation mass of about 6 $M_\oplus $. Just before t=4 Myr, migration switches to type II. Shortly thereafter at 4.1 Myr the planet crosses the iceline (see the prominent kink in the solid line in Fig. 11). Finally, the braking phase starts at about 4.4 Myr. Switching from type I to type II migration results in an increase of both the solid accretion rate and the core luminosity with a corresponding loss of envelope. The difference to the previous case originates from the longer living disk and the fact that the planet started at a larger distance from the central star. The solid mass available while migrating through the ``horizontal branch'' is therefore larger (it scales with $r^2 \Sigma_{\rm D}\propto r^{0.5}$) and the core reaches a mass large enough (21 $M_\oplus $) to trigger a runaway accretion of gas around 4.6 Myr, significantly before the gas disk disappears (or the planet reaches the feeding limit, as for the ``horizontal branch'' case). The gas accretion rate then increases rapidly, but soon (at 4.7 Myr) reaches the maximum rate allowed by the disk, which is a rather moderate 1.6 $\times $ 10-4 $M_\oplus $/yr. The disk limited accretion rate then decreases slowly with time as the disk evolves. From 4.7 Myr onwards the planet is in the braking phase, with the characteristic $-\pi $ slope in the a-M diagram until the disk disappears at 6 Myr. At this point, a Saturnian planet has formed at 0.9 AU with a total mass of about 130 $M_\oplus $ and a  $M_{Z}\approx 25~M_\oplus$. The prototypical planet also illustrated how the combination of effects that are per se straightforward (changing the migration regime, crossing the iceline, disk limited gas accretion rate) can lead to a complex formation history.

It is interesting to note that the ``main clump'' region of the a-M iagram is somewhat overpopulated in the observed population as well. Examples are (among many others) HD 100777b (Naef et al. 2007) or HD 142b (Tinney et al. 2002).

5.1.4 Tracks of ``outer group'' planets

For starting positions larger than $\sim$4-7 AU, the formation tracks are different. In the case of a high  $f_{\rm D/G}$ and $\Sigma _{0}$ and a low  $\dot{M}_{w}$ drawn together, the core growth timescale is short compared to the disk depletion timescale, even at these large distances. As the amount of solid material available is large, embryos can then grow to a supercritical mass almost in-situ (nearly vertical tracks up to several tens of $M_\oplus $ in Fig. 8), without the need to collect solid material by migration. Significant migration can nevertheless occur, but only in the planet dominated type II mode and when the gas accretion rate is regulated by the disk. Such planets, which do not have a ``horizontal branch'' phase, are represented in Fig. 8 by filled squares. The final semimajor axes of the in-situ supercritical planets are outside $\sim$0.4-1 AU, but overlap with ``main clump'' planets. Sometimes ``outer group'' planets become extremely massive ``Super Jupiters'', with masses more than one order of magnitude larger than that of Jupiter.

 \begin{figure}
\par\mbox{\includegraphics[width=9cm,height=8.3cm]{10301f16.eps} \includegraphics[width=9cm]{10301f17.eps} }
\end{figure} Figure 12:

As Fig. 9, but for in-situ critical cores which end up in the ``outer group''. The prototypical embryo starts in a massive, metal rich, long lived disk which allows the formation of a very massive ${\sim}11~M_{\textrm{\jupiter }}$ ``Super Jupiter'' ultimately located at 4.5 AU.

Open with DEXTER

Examples of such tracks are plotted in Fig. 12. Comparing these tracks with the ones of the ``main clump'' shows that there is a continuous transition between the two types. The prototypical ``outer group'' planet in the right panel of Fig. 12 is formed in a massive, metal-rich ( $\Sigma_{0}=500$ g/cm2, [Fe/H] = 0.3) and long lived disk. Due to the large  $a_{\rm start}$ (10 AU), the isolation mass at the initial location is then very large (around 150 $M_\oplus $). The seed starts at 1.1 Myr, and rapidly switches to type II migration at 1.2 Myr. The critical mass (around 25 $M_\oplus $) is already attained at 1.24 Myr at a position very close to  $a_{\rm start}$, triggering a rapid accretion of gas. At 1.25 Myr, the gas accretion rate reaches the value limited by the disk of initially 2.5 $\times $ 10-3 $M_\oplus $/yr.

The final amount of heavy elements in this planet is quite uncertain. Indeed, ``outer group'' planets undergo a large part of their formation history in the disk limited gas accretion phase. As mentioned above, the internal structure of planets during this phase is no longer calculated, thus their solid accretion rate (via  $R_{\rm capt}$) is uncertain. As a consequence, the final MZ found here (around 250 $M_\oplus $) is unsure, and only a lower boundary (${\sim}25$ $M_\oplus $) is actually well determined. Internal structure models (Baraffe et al. 2008) do however require comparable amounts of heavy elements to reproduce the observed mass-radius relation of the transiting ``Hot Super Jupiter'' HD 147506b (Bakos et al. 2007). Thus, even if this planet is in terms of semimajor axis very different to the planet discussed here, it still illustrates that objects with very large amounts of solids in their interior may exist. Figure 12 shows that the forming planet is eventually strongly dominated by gas (about 10  $M_{\textrm{\jupiter }}$ are accreted). Therefore, the uncertainty on the final total mass (which is of primary interest in this study) is low.

The figure also shows that MZ (with all the uncertainties discussed above) reaches its final value between 8 and 9 AU. The planet could in principle continue to accrete all the planetesimals down to its final location at 4.5 AU. However, the planet is so massive at this stage that nearly all planetesimals are rather ejected and not accreted (Ida & Lin 2004a). We, however, find that the large size of forming planets makes them less effective in ejecting planetesimals than old, compact objects of the same mass.

5.2 Mass-semimajor axis diagram

The planetary formation tracks illustrate how planets grow in mass and migrate to their final position. In this section, we discuss some aspects of the distribution of final masses and positions.

5.2.1 Diversity of planets

Figure 13 shows the mass and the semi-major axis of $N_{\rm synt}\approx50~000$ synthetic planets. The first striking result is that core accretion allows for a very diverse planet population. The final mass of the planets varies between the smallest possible mass (0.6 $M_\oplus $) and a few extremely massive ${\sim}40$  $M_{\textrm{\jupiter }}$ planets. The final position varies between the innermost possible radius ( ${\approx}0.1$ AU) and about 18 AU. We conclude that the observed diversity of extrasolar planets is, within the core accretion paradigm, a natural consequence of the observed diversity of the properties of the protoplanetary disks, which we have varied within the observed range. It is therefore obvious that a better understanding of the properties of the protoplanetary disks, especially the innermost region, has important implications for this planet formation theory. In our population synthesis, a number of parameters are kept fixed at all times. We can speculate that in nature, most of these quantities also fluctuate to some degree. The core accretion mechanism should therefore be able to produce planets of an even greater diversity than found here. On the other hand, the question of whether the core accretion can explain all the planets (cf. Matsuo et al. 2007) is a much more difficult one. The models are probably not yet developed enough and observations are still too incomplete to allow a definitive statement.

Figure 13 reveals that the synthetic planets are not randomly distributed inside the mass-distance plane. Instead, various concentrations, bars and depleted regions can be distinguished. One can also study the form of the quite well defined envelope filled by the synthetic population, and why certain regions remain empty. When comparing Fig. 13 with actual discoveries, one should, however, bear in mind the incompleteness of the model, in particular affecting low mass planets (Sect. 2.4).

 \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301f18.eps}
\end{figure} Figure 13:

Final mass M versus final distance a of $N_{\rm synt}\approx50~000$ synthetic planets of the nominal planetary population. The feeding limit at  $a_{\rm touch}$ is plotted as a dashed line. Planets migrating into the feeding limit have been set to 0.1 AU. As  $a_{\rm touch}$ becomes very large for $M \protect\ga 20~M_{\textrm{\jupiter }}$, a few extremely massive planets are also in the feeding limit, which should, however, be regarded as a simulation artifact because our simplification of putting planets that reach the feeding limit to 0.1 AU ceases to be justified.

Open with DEXTER

5.2.2 The limiting envelope

At large masses, the planetary population is bounded by the largest overall mass for a given semimajor axis, $M_{\rm max}(a)$. Outside $a\sim 3$ AU, $M_{\rm max}$ is a decreasing function of a, falling to about 100 $M_\oplus $ at ${\sim}20$ AU. The reason is that the time needed to build up a massive core at large distances becomes comparable to and ultimately longer than the gas disk dispersion timescale (Ida & Lin 2004a). See Sect. 5.3.4 for processes that could modify this behavior.

Inside 3 AU, the behavior of $M_{\rm max}$ is inverted, i.e.  $M_{\rm max}$ is an increasing function of a. This means that there is an absence of very massive planets at small orbital radii. This is what has been discovered in the observed extrasolar planet population if only planets around single host stars are considered (Udry et al. 2003; Zucker & Mazeh 2002). With the initial surface density profile used here ( $\Sigma\propto a^{-3/2}$), we find that inside 3 AU, $M_{\rm max}$ scales approximately as a3/4 (as  $M_{\rm iso}$), provided that type I migration is slow, see Sect. 5.3.3. For populations obtained with higher type I migration rates, $M_{\rm max}$ is flat inside ${\sim}3$ AU. Future discoveries of a very large number of single giant planets out to several AU (Ge et al. 2007) around single stars will help to better define the exact shape of  $M_{\rm max}(a)$.

5.2.3 Structures in the a - M diagram

The different phases of planet formation and migration that were identified in the formation tracks leave traces also in the final a-M of the planets. One can distinguish the ``failed cores'', the ``horizontal branch'', the ``main clump'', and the ``outer group'' planets. As a new feature, Fig. 13 also shows a depletion of planets with masses between 30 to 100 $M_\oplus $. This is the analogue of the ``planetary desert'' first discussed by Ida & Lin (2004a). Compared to their results, the depletion is much less severe in our simulations.

The reason for this difference is difficult to pinpoint, as both formation models differ in many aspects, but it is at least partially due to the way the maximal gas accretion rate of the planets is calculated. Both models use the criterion that the gas accretion rate given by the planet's Kelvin-Helmholtz timescale (which we implicitly obtain by the structure calculations) must be smaller than the mass transfer rate in the disk. Ida & Lin (2004a) however use this criterion only if the mass of the planet is additionally larger than the local gas isolation mass calculated from the unperturbed disk profile. This quantity is usually clearly larger (Ida & Lin 2004a) than the minimal planet mass required to fulfill the viscous and the thermal condition (Lin & Papaloizou 1985). At these masses, the planet becomes able to tidally open a (partial) gap in the protoplanetary disk (e.g. Ida & Lin 2008). This could reduce the amount of gas that is effectively in the planet's direct gravitational reach, i.e. its gas isolation mass. Therefore, it seems possible that the planet mass where the accretion rate in the disk becomes the limiting factor might be significantly smaller than the gas isolation mass as assumed in Ida & Lin (2004a), so that we disregard the second additional criterion, and always limit the planet's gas accretion rate to  $\dot{M}_{\rm disk}$.

As illustrated by the prototypical planet of the ``main clump'', the time at which planets start runaway gas accretion occurs for typical initial conditions generally at quite advanced stages of disk evolution. By then, the disk has undergone significant mass loss, and  $\dot{M}_{\rm disk}$ is low. Therefore, disk limited planetary accretion rates are usually down to a few $10^{-4}~M_\oplus$/yr so that growing a Jupiter-mass object requires an amount of time comparable to the remaining disk lifetime (see Sect. 5.3.2 for a population with very long, unrealistic disk lifetimes). Hence, the probability that the disk disappears while the planet is at an intermediate mass is not vanishingly small. This populates the ``planetary desert'' with intermediate mass objects, so that we predict only a moderate decrease of the relative number of planets with masses between 30 and 100 $M_\oplus $ (about a factor of 2-3 compared to Jovian planets, see the planetary initial mass function in Paper II). We have investigated this point by synthesizing a test population where we limit the gas accretion rate of a planet by  $\dot{M}_{\rm disk}$ only if its envelope mass is larger than the gas isolation mass, as in Ida & Lin (2004a), and found that this results in a population that has indeed a stronger depletion of intermediate mass planets, which are then about 5 to 6 times less frequent than Jovian planets.

Long baseline, high precision RV surveys will test for the existence of the ``horizontal branch'' and the ``planetary desert'', and determine their extension and intensity. The upper mass limit of the branch constrains the minimal mass needed to go into runaway gas accretion, while the lower mass boundary of the branch constrains the speed of type I migration. The degree of depletion of intermediate mass planets constrains planetary accretion rates during runaway gas accretion.

In the past 12 years, a significant number of extrasolar planets with masses much larger than one Jupiter mass were discovered. We also find such planets ($M \ga 10$  $M_{\textrm{\jupiter }}$) in our model, mainly in the ``outer group'', as shown in Fig. 13. Their formation is the consequences of the fact that we do not reduce gas accretion due to gap formation. About 0.4% of all initial conditions lead to the formation of planets with a mass exceeding 12-13  $M_{\textrm{\jupiter }}$, which is commonly regarded as the brown dwarf limit (Chabrier et al. 2000). A handful of planets (out of the $\sim$50 000) even reach a mass between 20 to 40  $M_{\textrm{\jupiter }}$. As it was shown recently that such objects burn deuterium in the layers above the core (Baraffe et al. 2008), they form an interesting population of deuterium burning planets.

The synthetic population does not contain analogs of Uranus and Neptune both in term of mass and semimajor axis. Indeed, no planets at all are found outside ${\sim}20$ AU. Partially, this is simply an artifact of the lack of planet formation after the dispersion of the gas disk in our model (Sect. 2.4.2). Uranus and Neptune are however challenging because these planets accreted significant hydrogen-helium envelopes i.e. because they were apparently formed while the gas disk was still present (Goldreich et al. 2004a; Chambers 2006).

It is well known that core accretion requires for the formation of Uranus and Neptune at their current position formation timescales which exceed typical gas disk lifetimes by a large factor (Pollack et al. 1996; Thommes et al. 2003). Hence, it was proposed that for Uranus and Neptune, one must give up the principle of in-situ formation and take into account the possibility of an ejection because of N-body interactions with the giant planets (Thommes et al. 2002; Tsiganis et al. 2005). The fact that our population contains in the ``horizontal branch'' a significant number of planets with the same mass and composition as the ice giants, but at smaller orbital distances, supports this second interpretation, and the general idea that planetary systems start with more crowded and compact configurations (e.g. Ford & Chiang 2007).

5.3 Non-nominal populations

The representative tracks we have identified in the previous sections and the final a-M were all calculated for one particular population, i.e. a particular choice of underlying assumptions and parameters of the model, like $\alpha $ or the type I migration efficiency factor. Here we discuss the effect of changing these parameters.

 \begin{figure}
\par\mbox{\resizebox{9cm}{!}{\includegraphics{10301f19.eps}}\resi...
...10301f20.eps}}\resizebox{9cm}{!}{\includegraphics{10301f22.eps}} }\end{figure} Figure 14:

Planetary formation tracks as in Fig. 8 but for non-nominal populations. Top left: for planetesimal eccentricities and inclinations as in Thommes et al. (2003). Top right: for $\dot{M}_{w}=0$ (no photoevaporation). Bottom left: for a type I migration efficiency factor $\f1$ = 0.1. Bottom right: for a type II migration rate as in Ida & Lin (2004a).

Open with DEXTER

5.3.1 Core growth regime

The accretion rate of planetesimals is crucial not only for the growth of the cores themselves, but also indirectly for the accretion of the envelopes. The solid accretion rate is a function of the velocity dispersion  $v_{\rm disp}$ of the planetesimals. We use the prescription of Pollack et al. (1996) for $v_{\rm disp}$ corresponding to a situation between the shear and the dispersion dominated regime. Different results have been obtained concerning the importance of these two regimes (Ida & Makino 1993; Rafikov 2003; Goldreich et al. 2004b). In the dispersion dominated picture of Thommes et al. (2003), planetesimal random velocities are higher than in the Pollack et al. (1996) description. Fortier et al. (2007) have studied the effects of these high  $v_{\rm disp}$ on giant planet formation and found an increase of the formation time of Jupiter by around one order of magnitude relative to Pollack et al. (1996), depending on the mass of the disk.

We have synthesized a population where we follow Fortier et al. (2007) in calculating the planetesimal eccentricity and inclination as in Thommes et al. (2003, their Eq. (10)). The top left panel of Fig. 14 shows the resulting formation tracks and makes it clear that there is a very strong effect, especially at large distances. Most seeds do not even grow from 0.6 to 1 $M_\oplus $. Only an extremely small fraction of initial conditions (extremely high $f_{\rm D/G}, \Sigma_{0}, t_{\rm disk}$) allows the formation of giant planets which are all inside 0.5 AU and have low masses.

It is obvious that Kolmogorov-Smirnov tests (Paper II) indicate a negligible probability that this synthetic population and the observed one come from the same parent distribution. We therefore conclude that in order to reproduce the observed extrasolar planets by the core accretion mechanism, solid accretion must occur on timescales clearly shorter than predicted by the model of Thommes et al. (2003). Indeed, including additional effects like planetesimal fragmentation (Chambers 2006; Kenyon & Bromley 2009) leads to significantly reduced core growth and planet formation timescales (see e.g. Thommes et al. 2008). These effects will have, at least qualitatively, similar consequences as the low  $v_{\rm disp}$ we assume.

5.3.2 Photoevaporation

Another population was synthesized with photoevaporation switched off ( $\dot{M}_{w}=0$), so that the disks evolve only due to viscosity. As expected, disk lifetimes are clearly increased for such a case. The mean disk lifetime is now about three times larger than in the nominal population. This is clearly incompatible with the observed distribution of disk lifetimes (Sect. 4.3). The top right panel of Fig. 14 shows that the final positions of the planets in the $\dot{M}_{w}=0$ population differ significantly from the nominal case, even if the tracks leading there are similar.

As expected, migration becomes more important, and a new group of massive planets inside about 0.5 AU is formed. The fraction of embryos that reach the inner border of the computational disk increases by about a factor of three compared to the nominal case. The fraction of initial conditions that leads to giant planets increases too, as planets have more time to grow, which also explains why a much emptier ``planetary desert'' is seen than in the nominal case.

Statistical tests (Paper II) show that this synthetic population is not compatible with the observed one: the KS tests for the $M\sin i$ and the semimajor axis distributions indicate a probability of only a few percent that the two samples come from the same parent population; the synthetic planets are too massive and too close-in. We therefore conclude that a distribution of disk lifetimes that is incompatible with the observational data (Haisch et al. 2001) leads to a synthetic population that is incompatible with the observed exoplanets, too.

5.3.3 Type I migration

Motivated by the short migration timescales found by Tanaka et al. (2002), which indicate that cores are lost to the star before they can grow large enough to trigger runaway gas accretion, several authors have been re-examining type I migration rates and found considerable cause for uncertainties (e.g. Nelson & Papaloizou 2004; Paardekooper & Mellema 2006). Even if it is clear that the simple multiplication of the migration rates of Tanaka et al. (2002) by a constant factor $\f1$ is only a first order approximation at best of the true type I migration rate (Kley & Crida 2008), the Monte Carlo simulations represent a suitable way to study the effects of various magnitudes of type I migration on the population as a whole, thereby possibly ruling out some values for $\f1$ (Ida & Lin 2008).

We have therefore synthesized populations using $\f1$ = 0.01, 0.1 and 1.0, as well as the value of 0.001 used above. In Fig. 14, bottom left, planetary formation track are plotted for a population using $\f1=0.1$ and otherwise nominal settings. For the high migration rates, one should keep in mind that migration is neglected before the seed embryo is put into the disk, therefore the true migration is partially underestimated (Ida & Lin 2008).

At $\f1=0.1$, migration causes many ``failed cores'' to migrate into the feeding limit, especially at small distances, as found also by Ida & Lin (2008). A significant number of ``failed cores'' nevertheless remains at intermediate distances, due to the disappearance of the nebula. Varying the migration rate has also related, important consequence for the occurrence of close-in (${\la}0.3$ AU) planets with masses $3 \la M \la 10$ $M_\oplus $. Such planets do not exist for a low migration rate (cf. Figs. 8 and 13) but become increasingly more numerous and smaller with increasing type I rate. This absence is the result of the combination of two effects. First, our model does not include the in-situ formation of terrestrial planets after disk dissipation. Hence, the formation of Mercury type planets and larger counterparts, or other mechanisms that lead to the formation of close-in terrestrial planets (see Raymond et al. 2007, for an overview) are not taken into account. This in-situ accretion would fill the empty area ``from below'' (in mass). Due to the limited amount of solids very close to the star, the mass to which planets can grow in-situ is however limited. Second, in most disks seed embryos can only start beyond the iceline (cf. the distribution of  $a_{\rm start}$, Sect. 4.4). Thereafter, they must be brought in by migration without growing too much in order to have a small final mass. This is only possible for fast migration (compared to accretion), and sets a minimal mass of planets that fill the empty region ``from outside'' (in semimajor axis).

For $\f1=0.1$ and 1.0, the mass of planets that migrate to $a_{\rm touch}$ is usually larger than about 6 and 2 $M_\oplus $, respectively. For all $\f1$, a few planets with even smaller masses are also found inside the feeding limit which come from very small initial positions  $a_{\rm start}$. It should further be noted that this result depends on the assumed radial surface density profile of solids. The $\Sigma_{\rm D}\propto a^{-3/2}$ we use for simplicity down to 0.1 AU leads to high solid surface densities near the star. For a truncated disk profile used in tests, which might be more appropriate, and $\f1=0.1$, even planets as small as ${\sim}2$ $M_\oplus $ have reached the feeding limit, and the region inside ${\sim}0.3$ AU which is empty at low $\f1$ is then filled.

Observations able to detect planets in the range 1 to 10 Earth masses close to the star will show if there is a minimal mass for planets inside a few tens of an AU, and whether the mass spectrum there is continuous, or contains a gap between planets formed in-situ and cores brought in from further out by migration. Hence, such objects will potentially provide us with an indicator of the efficiency of type I migration.

Considering the very massive planets (${\ga}10$  $M_{\textrm{\jupiter }}$), it is found that at low type I migration rates, no such planets are found within about 0.4 AU of the star. At high rates, some very massive planets have, in contrast, migrated to the feeding limit. The type I migration, which directly affects small planets, also has a significant and measurable effect on the most massive planets, as it allows planets forming at small distances to have a more rapid access to larger amounts of solids than the locally available  $M_{\rm iso}$ for $\f1=0.001$.

5.3.4 Type II migration

The formation tracks of giant planets of both of the ``main clump'' and the ``outer group'' show the importance of the planet dominated migration phase, as well as the disk limited gas accretion rate (which itself results from the process of gas transport inside the protoplanetary disk) in the braking phase. Using planet formation models following the growth and migration of a large number of embryos, Thommes et al. (2008) came to the same conclusion regarding the importance of these two processes.

It is clear that our description of the planetary ${\rm d}m/{\rm d}t$ and  ${\rm d}a/{\rm d}t$ in this phase remains a first approximation that could be significantly improved by a more sophisticated description of planet-disk interactions. For example, the eccentric instability that justifies the assumption that the planetary gas accretion rate is the same as the disk accretion rate occurs only for planets larger than a certain mass (Kley & Dirksen 2006). At smaller masses, the accretion rate could be smaller (D'Angelo et al. 2002). Also planet-planet interactions of several giant planets forming concurrently can have important consequences for the migration and accretion rate. This was shown by Thommes et al. (2008), where a chain of migrating giant planets, locked together by mean motion resonances, is a frequent configuration. Then, only the outermost planet directly exchanges torques with the gas disk, and accretes gas from it. While with our one embryo per disk approach, we obviously cannot model such a behavior, we note that this configuration typically leads to an even earlier transition to the planet dominated regime, as multiple planets are better at holding off the outer disk than one (Thommes et al. 2008).

The ${\rm d}a/{\rm d}t$ of a satellite that is massive compared to the local disk mass was studied by Syer & Clarke (1995) and Ivanov et al. (1999). They both agree that in this situation, the migration timescale becomes longer than the viscous timescale, but quantitatively their results differ. The differences are linked to the question of what part of the angular momentum flux is effective in moving the planet (Ida & Lin 2005), and whether gas overflows the gap formed by the planet (Armitage 2007). Therefore, there is some uncertainty on the exact migration rate in the case of a massive planet.

To have an estimate of the sensitivity of our results on the migration rate in the braking phase, we have synthesized a test population replacing Eq. (3) by the prescription of Ida & Lin (2004a). The bottom right panel of Fig. 14 shows the resulting formation tracks. One sees that there is a visible modification of the tracks, but that the basic types remain. The comparison of the tracks with the dashed lines (slope equals $-\pi $) shows that with this alternative prescription, type II migration occurs on a shorter timescale. The reason is that the ${\rm d}a/{\rm d}t$ is now proportional to  $\Sigma(R_{m})R_{m}^2$ instead of $\Sigma(a_{\rm planet})a_{\rm planet}^2$, where $R_{\rm m}$ is the radius of the maximum viscous couple, which moves outwards very quickly, and is in almost all cases larger than  $a_{\rm planet}$ once planets are large enough to migrate by type II migration. Therefore the first quantity is usually larger than the second one (Fig. 2).

There are also two examples of tracks where outward migration occurs. As expected from the rapid outward movement of $R_{\rm m}$, only a very small group of about 10 seeds (out of ${\sim}10~000$ initial conditions) has undergone significant outward migration[*]. The maximal planetary mass  $M_{\rm max}$ at ${\sim}20$ AU is therefore increased relative to the nominal model (Sect. 5.2.2) by about one order of magnitude.

The discovery of many more giant planets outside ${\sim}20$ AU by techniques like astrometry or direct imaging (Kalas et al. 2008; Marois et al. 2008) would indicate that besides outward migration, one or several of the following mechanisms not included in the nominal model is important: Scattering of planetary seeds to large  $a_{\rm start}$ early during formation (``monarchical growth'', Weidenschilling 2005, 2008), a completely different formation mechanism (direct gravitational collapse, e.g. Boss 2001; Mayer et al. 2004) or scattering of planets in initially crowded multiple planetary systems after formation (e.g. Veras & Armitage 2003). Note, however, for the last possibility that the simulations of Thommes et al. (2008), which combine in a self-consistent way the formation and subsequent N-body interactions, do not usually lead to the formation of giant planets further out than ${\sim}20$ AU. In another non-nominal population, we have addressed the early embryo ejection and used a modified prescription to calculate the starting time  $t_{\rm start}$ that approximately mimics a ``monarchical growth'' mode (Weidenschilling 2005, 2008). In this case, while leaving the population in the inner system relatively untouched, a population of very massive planets outside 3 AU comes into existence so that the maximal planetary mass  $M_{\rm max}$ does not decrease any longer outside ${\sim}3$ AU, but remains constant out to  ${\sim}10{-}20$ AU.

5.3.5 Transition from type I to type II migration

As mentioned at the beginning of the paper, the transition from type I to type II occurs when the planet can open a gap in the protoplanetary disk. In our nominal model, we assume that this occurs when the disk scale height becomes smaller than the planet's Hill radius (which is the relevant criterion in the limit of very large Reynolds numbers). In order to infer the influence of this part of our model, we have synthesized non-nominal populations using the type I/type II transition condition derived by Crida et al. (2006). Using the low type I migration of the nominal model, only very few ``Hot'' planets are then formed, since the transition to (faster) type II occurs at higher masses. However, when $\f1$ is increased, we obtain similar formation tracks and final sub-populations. Note that in this case, very few planet migrate in the disk dominated type II regime, and the majority of them passes directly from type I to the planet dominated type II regime. As outlined by e.g. Armitage and Rice (2005), the transition, and especially the migration rate during the transition, is not fully understood and could result in special migration regimes not considered here (e.g. Papaloizou & Terquem 2006). This incomplete understanding of migration in general remains one of the major uncertainties in planet formation theories of today.

6 Summary and conclusion

We have presented our extrasolar planet population synthesis calculations. As the formation model we use a slightly simplified version of the extended core accretion model presented in Alibert et al. (2005a) which has been shown to reproduce many observed properties of the giant planets of our own solar system (Alibert et al. 2005b).

We use four random variables to describe the possible initial conditions for planet formation: the dust-to-gas ratio  $f_{\rm D/G}$, the initial gas surface density  $\Sigma _{0}$, the photoevaporation rate  $\dot{M}_{w}$ and the starting position of the embryo in the disk, $a_{\rm start}$. The distributions for the first three Monte Carlo variables can be derived from observed properties of stars of the solar neighborhood, or protoplanetary disks. For the dust-to-gas ratio, we use the distribution of [Fe/H] of the stars in the CORALIE planet search sample (Santos et al. 2003). For the gas surface densities, we use the disk mass distribution in $\rho$ Ophiuchi (Beckwith & Sargent 1996). For the photoevaporation rate, we use the distribution of disk lifetimes by Haisch et al. (2001). We have also discussed the complications that arise from the conversion of such observed quantities into figures that can be used in numerical simulations. The last random variable, $a_{\rm start}$, cannot be derived from observations. Here we follow Ida & Lin (2004a) and use a distribution that is uniform in $\log(a)$.

With the adopted random distributions, we find that our formation model predicts a population of synthetic planets that is very diverse, not unlike the actually observed population. The final semimajor axis varies by two orders of magnitude, and the mass by four orders. Similarly, the internal bulk composition is very different, and covers gas giants and also ``Super Earth'' mass planets with a envelope to core mass ratio similar to Venus. We conclude that this diversity illustrates the ability of the underlying concepts of core formation to serve as a unified formation model for planets of masses ranging from a few times the Earth mass to beyond the brown dwarf limit. The observed diversity of extrasolar planet is a natural consequence of the different properties of protoplanetary disks.

Planet formation is a process of concurrent mass accretion and migration which can be well represented by planetary formation tracks in the mass-distance plane. These tracks show a number of distinct phases, which are visible in the final a-M distribution of the planets, too.

In the first phase, at small masses, planets migrate in type I migration (which must however be slow to be compatible with observations, cf. Paper II), and accrete mostly solids. Planets which remain in this stage until the moment when the disk disappears form a vast sub-population of low mass, core-dominated planets, ``failed cores''. The distributions of protoplanetary disks properties are such that this is the most likely outcome of the formation process. This is consistent with the observation that 90 to 95% of FGK stars in the solar neighborhood apparently remained without giant planets.

If the disk properties are instead such that the planet grows massive enough to open a gap in the disk, the second phase starts. Then, planets migrate inwards in type II migration, collecting solids on their way in. In this phase, migration occurs on a shorter timescale than accretion, so that the tracks of the planets show a ``horizontal branch'' phase in the a-M plot. Planets in this phase have masses between 10 and 30 $M_\oplus $, a static gaseous envelope and an internal gas to solid ratio similar to Neptune.

Some planets on the ``horizontal branch'' collect enough solids to become supercritical for gas runaway accretion. Their gas accretion rate then quickly reaches the disk limited value, and their mass becomes large compared to the local disk mass. Therefore, giant gaseous planets reach their final mass and semimajor axis in the third phase, when they accrete gas at the rate controlled by the disk and migrate in the planet dominated, slower type II regime. This leads to the formation of a concentration of giant planets between 0.3 to 2 AU in the ``main clump'', with masses between ${\sim}100$ $M_\oplus $ and 3  $M_{\textrm{\jupiter }}$.

Embryos starting far from the parent star in metal rich disk can grow supercritical for gas runaway accretion in-situ, without passing through the ``horizontal branch''. This leads to the formation of an ``outer group'' of giant planets.

Apart from these four sub-populations, we find in the final nominal a-M diagram (1) an absence of massive ( ${\ga}10~M_{\textrm{\jupiter }}$) planets both close (${\la}0.5$ AU) to the star and very far (${\ga}10$ AU) from it; (2) a certain depletion of planets with masses between 30 and 100 $M_\oplus $, in analogy to the ``planetary desert'' (Ida & Lin 2004a) which is however not very strong (a factor of 2-3 relative to giant planets), as at the time planets start runaway gas accretion, the disk limited accretion rate, which is decisive in this regime, is usually already quite low; (3) a handful very massive ( ${\ga}20~M_{\textrm{\jupiter }}$) deuterium burning planets (Baraffe et al. 2008), showing that this mass domain can be populated by different formation channels; (4) an absence of Neptune analogs in terms of both mass and semimajor axis, but many such planets at smaller distances.

Population synthesis is a valuable tool to asses the global consequences of model settings. From the synthesis of a number of non-nominal populations, we find that (1) cores must form quickly (Thommes et al. 2003; Chamber 2006); (2) disks with lifetimes too long compared to observations lead to too massive planets too close-in; (3) fast type I migration brings low mass planets close to the parent star, so that the depletion of planets with $3\la M/M_\oplus\la 10$ inside 0.3 AU in the nominal population vanishes; (4) the early ejection of seeds, and outward type II migration could result in very massive planets out to about 20 AU, but only in small numbers; (5) the transition regime between type I and type II migration has important implications for the final population.

With improving detection methods leading to a rapidly increasing number of known extrasolar planets, the actual distributions of masses and semi-major axis can be determined with growing statistical confidence and compared to model predictions. Hence, planet populations studies are the ideal tools to extract from these distributions constraints on various aspects of the formation models, thereby reaping the scientific benefits of the large observational efforts invested. In the companion paper (Paper II), we carry out a direct statistical comparison between a carefully chosen sample of actual exoplanets and detectable synthetic planets and use this comparison to extract relevant constraints on the formation mechanism. Also, population synthesis can itself be used to guide and perhaps optimize the next generation of detection instruments.

Acknowledgements
We thank Stephane Udry and Shigeru Ida for useful discussions. We are grateful for detailed and interesting comments by two anonymous referees. This work was supported in part by the Swiss National Science Foundation. Computations were made on the ISIS, ISIS2 and UBELIX clusters at the University of Bern, and on the cluster of the Observatoire de Besançon funded by the Conseil Général de Franche-Comté.

References

 

Footnotes

...[*]
Current address: Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
... simplifications[*]
The synthesis of a population of ${\sim}30~000$ planets takes several days on a 50 CPU cluster. Most of the time is spent in solving the planetary envelope structure equations.
... mass[*]
From Eq. (11) one would calculate a  $M_{\rm iso}$ of about $2.8~M_\oplus$, using BL=4. However, as we do not reduce the initial solid surface density by the amount of material already in the initial seed, a value larger for the mass by about 3/2 $\times $ $M_{\rm emb,0}$ is obtained.
... migration[*]
Note that we calculate $R_{\rm m}$ in this test population as in Ida & Lin (2004a), and not in a self-consistent way from the disk model, where photoevaporation can influence the evolution of $R_{\rm m}$ at late times and hence the outward migration (Veras & Armitage 2003).

All Tables

Table 1:   Parameters $\mu $ and $\sigma $ describing the Gaussian distribution of [Fe/H], for different observation samples.

Table 2:   Parameters describing the lognormal distribution of circumstellar disk masses $\mu $ and $\sigma $, from different sources.

All Figures

  \begin{figure}
\par\includegraphics[width=9cm,clip]{10301fg1.eps}
\end{figure} Figure 1:

Position of the iceline $a_{\rm ice}$ as a function of the initial gas surface density  $\Sigma _{0}$ at 5.2 AU (upper three lines). It corresponds to an initial  $T_{\rm mid}$ of 170 K. The iceline is plotted for three values of $\alpha $: 0.01 (dashed line), 0.007 (solid line) and 0.001 (dotted line). The lower three lines correspond to an initial  $T_{\rm mid}$ of 1600 K, roughly the evaporation temperature of rock. The rockline  $a_{\rm rock}$ is however not taken into account in the nominal model, due to the difficulty in defining its relevant location, as disk evolution is very rapid close-in and irradiation effects might be important (cf. Paper II).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.8cm,clip]{10301fg2.eps}
\end{figure} Figure 2:

Example of the evolution of $2 \Sigma a^2$. Seven different moments in time are plotted (time in units of Myr). At 4.6 Myr the disk has nearly vanished, and the calculations are stopped.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg3.eps}
\end{figure} Figure 3:

Solid lines: histogram and Gaussian fit to the [Fe/H] distribution of the CORALIE planet search sample. Dotted line: for comparison, distribution of photometric metallicities [Me/H] of Nordström et al. (2004). Probability densities are given, i.e. the area under the curves has been normalized to unity in each case.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg4.eps}
\end{figure} Figure 4:

Histograms and fits to circumstellar disk masses  $M_{\rm disk}$ for Ophiuchus (solid lines) and Taurus-Auriga (dotted lines). The actual gas masses in the computational domain out to 30 AU are between 0.004 and 0.09 $M_\odot $.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg5.eps}
\end{figure} Figure 5:

Fraction of stars in the model possessing a gaseous disk as a function of time for a uniform distribution in log of $\dot{M}_{w}$ between 5 $\times $ 10-10 and 3 $\times $ $10^{-8}~M_\odot$/yr (thin solid line). The thick solid line is the fit of Haisch et al. (2001) to the observed JHKL excess/disk fraction as a function of mean cluster age. The error bar at the top right is the overall systematic uncertainty in age of these observations, also from Haisch et al. (2001). For comparison, the fraction is also given if two different distributions of  $\dot{M}_{w}$ are used in the model: uniformly in log between 1 $\times $ 10-8-1 $\times $ $10^{-7}~M_\odot$/yr (dash-dotted line), leading to disk lifetimes that are clearly too short, or uniform in log in 1 $\times $ 10-10-1 $\times $ $10^{-9}~M_\odot$/yr, leading to disk lifetimes that are too long (dotted line).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg6.eps}
\end{figure} Figure 6:

Probability distribution for $a_{\rm start}$. There is a marked peak around 3 AU. Outside this distance, the probability distribution falls off as expected for a distribution that is uniform in log. The fraction inside about 2 AU is due to disks with the ability to form seeds sufficiently massive ( ${\geq }M_{\rm emb,0}$) inside the iceline as well.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301fg7.eps}\hspace*{7mm}
\includegraphics[width=8.5cm,clip]{10301fg8.eps}
\end{figure} Figure 7:

Snapshots of the embryo mass (solid line) as a function of semimajor axis at four moments in time for two different solid surface densities. The dashed line is the isolation mass. The dotted line is $M_{\rm emb,0}=0.6$ $M_\oplus $. The initial solid surface density at 1 AU is 7 g/cm2 ( left panel) and 35 g/cm2 ( right panel). It should be kept in mind that this kind of calculation is needed to generate the start time  $t_{\rm start}$ when the embryo is put into the formation model. The real evolution of the solid core for M>0.6 $M_\oplus $ is in general much more complex than plotted here. In this figure, we have continued the calculations up to the isolation mass to allow comparison with other models.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=16cm,clip]{10301fg9.eps}
\end{figure} Figure 8:

Planetary formation tracks in the mass-distance plane. The large black symbols show the final position of a planet. The shape of the symbols is explained in the text. Planets reaching the feeding limit at  $a_{\rm touch}$ (indicated by the long dashed line) have arbitrarily been set to 0.1 AU. The short dashed lines have a slope of $-\pi $ (discussion in Sect. 5.1.3). Each track is color-coded according to the migration mode, and small black dots are plotted on the tracks every 0.2 Myr to indicate the temporal evolution of a planet.

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{ \resizebox{9cm}{8.7cm}{\includegraphics{10301f10.eps}}...
...ox{9cm}{!}{\includegraphics[angle=270,origin=rb]{10301f11.eps}} }
\end{figure} Figure 9:

Planetary formation tracks in the mass-distance plane for ``failed cores'' ( left panel). The black points represent the final position of the planets. The thick line starting at 3.7 AU is the track of one prototypical example, for which the right panels shows the temporal evolution. Its final position is represented by a large square. In the right panel, the total mass M (solid line), the mass of accreted solids MZ (dashed line) and the mass of the envelope  $M_{\rm env}$ multiplied by a factor of 10 for better visibility (dotted line) are plotted (scale on the left) as a function of time t. The temporal evolution of the planet's semimajor axis a is also plotted (dash-dotted line, scale on the right).

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{\resizebox{9cm}{8.55cm}{\includegraphics{10301f12.eps}} \resizebox{9cm}{!}{\includegraphics{10301f13.eps}} }\end{figure} Figure 10:

As Fig. 9, but for planets of the ``horizontal branch''. In this case, the prototypical track (thick line, and large square at the final position) starts at $a_{\rm start}=4.7$ AU. It ends as a ``Hot Neptune'' planet at the feeding limit at 0.1 AU. Its temporal evolution is shown in the right panel.

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{\includegraphics[width=9cm,height=8.8cm]{10301f14.eps} \includegraphics[width=9cm]{10301f15.eps} }
\end{figure} Figure 11:

As Fig. 9, but for planets passing through the ``horizontal branch'' to become members of the ``main clump''. Here, the prototypical embryo (thick line, square at the final position, temporal evolution in the right panel) eventually leads to a Saturnian planet situated at 0.9 AU.

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{\includegraphics[width=9cm,height=8.3cm]{10301f16.eps} \includegraphics[width=9cm]{10301f17.eps} }
\end{figure} Figure 12:

As Fig. 9, but for in-situ critical cores which end up in the ``outer group''. The prototypical embryo starts in a massive, metal rich, long lived disk which allows the formation of a very massive ${\sim}11~M_{\textrm{\jupiter }}$ ``Super Jupiter'' ultimately located at 4.5 AU.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{10301f18.eps}
\end{figure} Figure 13:

Final mass M versus final distance a of $N_{\rm synt}\approx50~000$ synthetic planets of the nominal planetary population. The feeding limit at  $a_{\rm touch}$ is plotted as a dashed line. Planets migrating into the feeding limit have been set to 0.1 AU. As  $a_{\rm touch}$ becomes very large for $M \protect\ga 20~M_{\textrm{\jupiter }}$, a few extremely massive planets are also in the feeding limit, which should, however, be regarded as a simulation artifact because our simplification of putting planets that reach the feeding limit to 0.1 AU ceases to be justified.

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{\resizebox{9cm}{!}{\includegraphics{10301f19.eps}}\resi...
...10301f20.eps}}\resizebox{9cm}{!}{\includegraphics{10301f22.eps}} }\end{figure} Figure 14:

Planetary formation tracks as in Fig. 8 but for non-nominal populations. Top left: for planetesimal eccentricities and inclinations as in Thommes et al. (2003). Top right: for $\dot{M}_{w}=0$ (no photoevaporation). Bottom left: for a type I migration efficiency factor $\f1$ = 0.1. Bottom right: for a type II migration rate as in Ida & Lin (2004a).

Open with DEXTER
In the text


Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.