1 - Physikalisches Institut, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland
2 - Theoretical Astrophysics Group, University of Leicester, Leicester LE1 7RH, UK

Abstract
We present a new model of giant planet formation that extends the core-accretion model of Pollack et al. (1996, Icarus, 124, 62) to include migration, disc evolution and gap formation. We show that taking these effects into account can lead to much more rapid formation of giant planets, making it compatible with the typical disc lifetimes inferred from observations of young circumstellar discs. This speed up is due to the fact that migration prevents the severe depletion of the feeding zone as observed in in situ calculations. Hence, the growing planet is never isolated and it can reach cross-over mass on a much shorter timescale. To illustrate the range of planets that can form in our model, we describe a set of simulations in which we have varied some of the initial parameters and compare the final masses and semi-major axes with those inferred from observed extra-solar planets.

Key words: stars: planetary systems - stars: planetary systems: formation - solar system: formation

The standard giant planet formation scenario is the so-called core-accretion model. In this model, a solid core is formed first by accretion of solid planetesimals which themselves were formed by sedimentation and coagulation of small dust grains (Wetherill & Steward 1989; Lissauer 1993). As the core grows, it eventually becomes massive enough to gravitationally bind some of the nebular gas, thus surrounding itself by a tenuous envelope. The subsequent evolution of this core/envelope structure has been studied first by Perri & Cameron (1974) and subsequently in great detail by Bodenheimer & Pollack (1986, hereafter referred to as BP86) and Pollack et al. (1996, hereafter referred to as P96). These authors have found that both the solid core and the gaseous envelope subsequently grow in mass, the envelope remaining in quasi-static and thermal equilibrium. During this phase, the energy radiated by the gas is supplied by the energy released from the accretion of planetesimals. As the planet reaches a high enough mass (of the order of $20{-}30 \;M_\oplus$ at 5 AU, but depending on different physical parameters such as the solid accretion rate), radiative losses become so large that they can no longer be offset by planetesimal accretion alone and the envelope starts to contract. At this stage, the mass is nearly equally distributed between the mass of accreted gas and the mass of accreted planetesimals. The contraction of the envelope increases the gas accretion rate which in turn increases the energy losses through radiation and the process runs away, very rapidly building up a massive envelope.

Detailed numerical calculations, in particular by P96, have shown that this core-accretion formation model can be subdivided in three distinct phases. During phase 1, the core is formed by accretion of planetesimals located inside the growing planet's feeding zone, which extends to a few Hills radii. The accretion rate of solids during this phase is typically of the order of $10^{-5} \;M_\oplus/$ yr while the gas accretion rate is several orders of magnitude lower. Phase 1 ends when the feeding zone of the planet becomes severely depleted, which generally occurs before the planet has reached a mass high enough to accrete a large amount of gas. During phase 2, the mass increase is essentially due to the slow accretion of gas in the envelope. Note that if the envelope mass increases, the feeding zone of the planet also increases which in turns allows the accretion of more planetesimals. Both accretion rates (solid and gas) turn out to be relatively constant during this phase with the gas accretion rate exceeding by a fraction of an order of magnitude the solid accretion rate. This phase lasts until the planet has reached the cross-over mass (mass of accreted planetesimals equals mass of accreted gas), at which time the system enters phase 3. At this point, evolution proceeds extremely fast by runaway gas accretion as the envelope is no longer able to maintain quasi-static equilibrium. The mass of the planet increases correspondingly. The timescale for the formation of a giant planet in this core-accretion scenario is almost entirely determined by phase 2, which, as shown by P96, turns out to be an extremely sensitive function of the disc surface density. P96 found for their preferred model for the formation of Jupiter (model J1, the surface density being equal to $\sim$ 4 times the value in the minimum mass solar nebula) a formation timescale close to 8 Myr while reducing the surface density to 75% of this value leads to a formation time of nearly 50 Myr.

The core-accretion scenario was motivated by the apparent existence of a solid core, long estimated to be of the order of the mass of accreted planetesimals when the planet reaches the cross-over mass, in all the giant planets of the solar system. This scenario has also been supported by the enrichment in heavy elements (heavier than He) of both Jupiter and Saturn compared to the solar value, deduced from interior structure models and various remote/in situ measurements (radius, mass, surface abundance, gravitational moments, see Guillot et al. 2004).

However, the major difficulty affecting the core-accretion scenario, which has been pointed out repeatedly, is related to the timescale required to form a giant planet. Based on astronomical observations, protoplanetary discs are believed to transport mass inward at a rate of $10^{-8 \pm 1}\; M_\odot/$ yr (Hartmann et al. 1998) while their total mass has been estimated to lie between 10^-3 and $10^{-1}\; M_\odot$ (Beckwith & Sargent 1996). From these numbers as well as from the observations that about half the stars in young clusters loose their discs within 3 Myr (Haisch et al. 2001) one infers a typical circumstellar disc lifetime of 1-10 Myr which is of the same order as (if not smaller than) the giant planet formation timescale. To circumvent this problem, Boss (2002) proposed that giant planets form directly from the gravitational fragmentation and collapse of a protoplanetary disc. However, there are still a number of open issues about this scenario such as the formation and survival of bound structures, since most calculations so far have used an isothermal equation of state and/or too low resolution. Furthermore, while solving the timescale problem, it is not clear at all if the peculiar composition and structure of Jupiter and Saturn can be explained within this model. More work is definitively required to investigate this scenario further.

Since the discovery by Mayor & Queloz (1995) of the first extrasolar giant planet at a short distance to its star, there is mounting evidence that planets might have formed at locations which do not necessarily correspond to those where they are observed today. Gravitational interactions between the growing planet and the gaseous disc lead to angular momentum transfer resulting in a net inward migration of the planet and possibly gap formation (Lin & Papaloizou 1986; Lin et al. 1996; Ward 1997; Tanaka et al. 2002). For low mass planets, the migration rate is linear with mass (type I migration, Ward 1997). Higher mass planets open a gap and the migration rate is set by the viscosity independently of planetary mass (type II migration, Ward 1997). While the general physical understanding of the origin of migration is clear, the actual migration rates obtained, especially for type I migration, are so short that they are inconsistent with the number of extrasolar planets actually detected. Type II migration timescales are found to lie between 0.1 to 10 Myr, a timescale comparable to (if not shorter than) the typical lifetime of the disc as well as the planet formation timescale in the core-accretion scenario.

Since all relevant timescales (planet formation, disc evolution, and migration) are of the same order of magnitude, it appears difficult to obtain a self-consistent picture while omitting anyone of these processes. For this purpose, we have developed a new code, structured around three modules. These modules calculate the disc structure and evolution (using the simple formalism of $\alpha$ viscosity), including planet migration, the interaction of planetesimals with the planet atmosphere, and the planet structure and evolution. In Sect. 2 we will describe these modules in more detail and present the numerical tests we have performed to validate them. Section 3 will be devoted to the effect of migration and disc evolution on formation timescales, and to formation models of extrasolar planets. Finally, we will discuss these results and conclude in Sect. 4.

2 Equations and assumptions

Our code to compute the formation of giant planets consists of three modules. The first module follows the method given by Papaloizou & Terquem (1999, referred to PT99 in the following) to calculate the disc structure and its time evolution. The second calculates the interaction between the infalling planetesimals and the atmosphere of the growing planet; the last is a planetary structure and evolution code, written especially for this project. We now describe these modules, and present some of the tests we have performed.

2.1 Disc structure and evolution

2.1.1 Vertical structure

The first module aims at determining the structure (both vertical and radial) of a protoplanetary disc. The method is identical to the one used by PT99 and therefore we only briefly recall the main points. Cylindrical symmetry is assumed and therefore cylindrical coordinates $(r,z,\theta)$ are a natural choice. The disc is assumed to be thin. For each distance rto the star, the vertical structure is calculated by solving the equation of hydrostatic equilibrium,

$\begin{displaymath}{1 \over \rho} {{\partial {P}} \over {\partial {z}}} = - \Omega^2 z , \end{displaymath}$

(1)

$\begin{displaymath}{{\partial {F}} \over {\partial {z}}} = {9 \over 4} \rho \nu \Omega^2 , \end{displaymath}$

(2)

$\begin{displaymath}F = {- 16 \pi \sigma T^3 \over 3 \kappa \rho} {{\partial {T}} \over {\partial {z}}} , \end{displaymath}$

(3)

The boundary conditions for this part of the calculation are the same as in PT99, formally,

$\begin{displaymath}T(z=H) = T(\tau_{\rm ab},T_{\rm b},r,\dot{M_{\rm st}},\alpha) , \end{displaymath}$

(4)

$\begin{displaymath}P(z=H) = {\Omega^2 H \tau_{\rm ab} \over \kappa(T(z=H),P(z=H))} , \end{displaymath}$

(5)

This system of 3 equations with 4 boundary conditions has in general no solution, except for a certain value of H. This value is found iteratively: Eqs. (1)-(3) are numerically integrated from z=H to z=0, using a fifth-order Runge-Kutta method with adaptive step length (Press et al. 1992) until F(z=0)=0 to a given accuracy.

Using this procedure, we calculate, for each distance to the star r and each value of the equilibrium accretion rate $\dot{M_{\rm st}}$ , the distributions of pressure, temperature and density $T(z;r,\dot{M_{\rm st}})$ , $P(z;r,\dot{M_{\rm st}})$ , $\rho(z;r,\dot{M_{\rm st}})$ .

Using these distributions, we finally calculate the mid-plane temperature ( $T_{\rm mid}$ ) and pressure ( $P_{\rm mid}$ ), as well as the effective viscosity $\tilde{\nu}(r,\dot{M_{\rm st}})$ , the disc density scale height $\tilde{H} (r,\dot{M_{\rm st}})$ defined by $\rho(z=\tilde{H}) = e^{-1/2} \rho(z=0)$ . The surface density $\Sigma(r,\dot{M_{\rm st}})$ is also given as a function of $\dot{M_{\rm st}}$ (for each radius). By inverting this former relation, we finally obtain relations $T_{\rm mid}(r,\Sigma)$ , $P_{\rm mid}(r,\Sigma)$ , $\tilde{\nu}(r,\Sigma)$ and $\tilde{H}(r,\Sigma)$ for each value of r (and each value of the other parameters $\alpha$ , $\tau_{\rm ab}$ and $T_{\rm b}$ ).

2.1.2 Evolution of the surface density

The time evolution of the disc is governed by the well-known diffusion equation (Lynden-Bell & Pringle 1974):

$\begin{displaymath}{{\rm d} \Sigma \over {\rm d} t} = {3 \over r} {\partial \ove... ... \over r} {\partial \over \partial r } \left( r J(r) \right) , \end{displaymath}$

(8)

$\begin{displaymath}{{\rm d} \Sigma \over {\rm d} t} = {3 \over r} {\partial \ove... ...da (r) \right] + \dot{\Sigma}_w(r) + \dot{Q}_{\rm planet}(r) . \end{displaymath}$

(9)

$\begin{displaymath}\Lambda(r) ={ f_{\Lambda}\over 2 r}\sqrt{GM_{\rm star}}\left(... ... {r \over {\rm max}( \vert r-a\vert , \tilde{H} )} \right)^4 , \end{displaymath}$

(10)

$\displaystyle \left\{ \begin{array}{l} \dot{\Sigma}_w = 0 ~~~~~~~~~~~ {\rm for}... ...t{\Sigma}_w \propto R^{-1} ~~~ {\rm for} ~~~ R > R_{\rm g} , \end{array}\right.$

(11)

To solve the diffusion Eq. (9) we need to specify two boundary conditions. The first one is given at the outer radius of the disc (in our simulations this radius is usually taken at 50 AU). At this radius, one can either give the surface density $\Sigma$ or its temporal derivative. Since the characteristic evolution time of the disc is the diffusion timescale

$\begin{displaymath}T_{\nu} \propto {r^2 \over \tilde{\nu}} \propto {1 \over \alpha \Omega } \left( {r \over H} \right)^2 , \end{displaymath}$

(12)

The second condition is specified at the inner radius where we have used the following condition:

$\begin{displaymath}r {{\partial { \tilde{\nu} \Sigma}} \over {\partial {r}}} \bigg\vert _{\rm inner ~~~ radius} = 0 . \end{displaymath}$

(13)

$\begin{displaymath}\Phi(r) \equiv 2 \pi r J(r) = 3 \pi \tilde{\nu} \Sigma + 6 \pi r {{\partial {\tilde{\nu} \Sigma}} \over {\partial {r}}} , \end{displaymath}$

(14)

$\begin{displaymath}\Phi(r) \Big\vert _{\rm inner ~~~ radius} = 3 \pi \tilde{\nu} \Sigma = \dot{M_{\rm st}} , \end{displaymath}$

(15)

2.2 Migration rate

Dynamical tidal interactions of the growing protoplanet with the disc lead to two phenomena: inward migration and gap formation (Lin & Papaloizou 1979, Ward 1997, Tanaka et al. 2002). For low mass planets, the tidal interaction is linear, and migration is of type I (Ward 1997), whereas higher mass planets open a gap, leading to a reduction of the inward migration (referred to type II migration).

Analytical models of type I migration have been computed by Ward (1997). The resulting migration timescales are much shorter than both the disc lifetime and the planet growth timescale, making survival of forming planets unlikely: the planet is accreted onto the central star. Migration could be stopped if there is an inner cavity in the disc, but planets at larger distances remain difficult to explain. Tanaka et al. (2002) performed new analytical calculations of type I migration, in two- or three-dimensional discs and found longer migration timescales but still too short to ensure survival. Their migration rate is nevertheless confirmed by recent three-dimensional numerical calculations of disc structure and planet migration (Bate et al. 2003).

On the other hand, suggestions of increased type I migration timescales can be found in Nelson & Papaloizou (2004). As shown by these authors, the torques exerted on at least low mass planets ( $M_{\rm planet}< 30 \;M_\oplus$ ) embedded in turbulent MHD discs are strongly fluctuating, resulting in a slowing down of the net inward motion. Contrary to laminar discs (as considered by Tanaka et al. 2002; and Bate et al. 2003) the migration proceeds as a random walk, and the mean value of the migration velocity seems to be highly reduced, compared to the laminar case. Moreover, as shown by Menou & Goodman (2004), type I migration of low-mass planets can be slowed down by nearly one order of magnitude in regions of opacity transitions.

These considerations seem to indicate that the actual type I migration timescale may in fact be considerably longer than the one originally estimated by Ward (1997) or even by Tanaka et al. (2002). For these reasons, and for lack of better knowledge, we actually use for type I migration the formula derived by Tanaka et al. (2002) reduced by an arbitrary numerical factor $f_{\rm I}$ chosen between 1/10 and 1/100. Tests have shown that provided this factor is small enough to allow planet survival, its exact value does not change the formation timescale but just the extent of the migration (see Sect. 3.1).

$\begin{displaymath}{{{\rm d} {a_{\rm planet}}} \over {{\rm d} {t}}} = - 2 f_{\rm I} a_{\rm planet}{ \Gamma \over L_{\rm planet}} , \end{displaymath}$

(16)

$\begin{displaymath}\Gamma = (1.364 + 0.541 \alpha_{\Sigma,{P}} )\left( {M_{\rm p... ...\rm s},P}} \right)^2 \Sigma_{\rm P} r_P ^4 \Omega_{\rm p} ^2 , \end{displaymath}$

(17)

For type II migration, two cases have to be considered. For low mass planets (when their mass is negligible compared that of the disc) the inward velocity is given by the viscosity of the disc. As the mass of the planet grows and becomes comparable that of the disc, migration slows down and eventually stops. In the latter case, the variation of the planet orbital angular momentum is equal to the angular momentum transport rate in the gaseous disc (Lin et al. 1996; Ida & Lin 2004):

$\begin{displaymath}{{\rm d} \over {\rm d} t} \left[ M_{\rm planet}a_{\rm planet}^2 \Omega \right] = {3 \over 2} \Sigma \tilde \nu \Omega r^2 . \end{displaymath}$

(18)

$\begin{displaymath}{{\rm d} a_{\rm planet} \over {\rm d} t} = -{3 \nu \over 2 a_... ...s {\rm Min} (1,{2 \Sigma a_{\rm planet}^2 / M_{\rm planet}}) . \end{displaymath}$

(19)

2.3 The planetesimals

2.3.1 Surface density and physical properties

The initial amount of heavy elements in the disc is a poorly constrained quantity. For this reason, the dust-to-gas ratio is varied in our simulations, and takes two values depending on the mid-plane temperature of the disc: f_D/G for temperatures below 150 K and 1/4 f_D/G for higher temperatures. In principle, the position of the iceline should evolve because of the viscous evolution of the disc. However, since our treatment of the planetesimals disc is very simple, we do not take into account this evolution.

We assume that due to the scattering effect of the planet, the surface density of planetesimals is constant within the current feeding zone but decreases with time proportionally to the mass accreted (and/or ejected from the disc) by the planet. The feeding zone is assumed to extend to a distance of 4 $R_{\rm H}$ on either side of the planetary orbit, where $R_{\rm H}\equiv \left({ M_{\rm planet}\over 3 M_*} \right)^{1 / 3} a_{\rm planet}$ is the Hills radius of the planet. For the inclinations and eccentricities of the planetesimals, we use the following prescription (P96):

$\begin{displaymath}i = {1 \over a_{\rm planet}} \sqrt{2 G M_{\rm planetesimal} \over r_{\rm planetesimal}} {1 \over \sqrt{3} \Omega} , \end{displaymath}$

(20)

$\begin{displaymath}e = {\rm max}\left(2 i, 2 {R_{\rm H}\over a_{\rm planet}}\right) \cdot \end{displaymath}$

(21)

$\begin{displaymath}{{\rm accretion ~ rate } \over {\rm ejection ~ rate }} = \left( { V_{\rm esc,disk} \over V_{\rm surf,planet}} \right)^4 , \end{displaymath}$

(22)

Note that our model for the evolution of the disc of planetesimals remains a very simple one in which a number of effects are neglected. For example, we omit the effect of gas drag on the planetesimals which, given their assumed size of 100 km, is reasonable. Using the gas density obtained in our model we can calculate the typical timescale for radial drift. For regions above $\sim$ 4 AU we find values of order of $\sim$ 10⁶ to 10⁷ the orbital time. In regions closer to the sun and at the begining of our simulation, this timescale may be much lower owing to the higher gas density. However, these inner regions evolve rapidly to lower densities, and gas drag becomes quickly negligible.

We also neglect, apart for ejection and accretion onto the planet, the perturbations created by the growing planet in the planetesimal disc such as heating, gap formation and shepherding (Tanaka & Ida 1997; Tanaka & Ida 1999; Thommes et al. 2003). For the latter effect, we have calculated the critical migration rate below which shepherding occurs (Tanaka & Ida 1999) and found that in all cases the migration rate exceeds this critical value. The planet acts as a "predator'' and not as a "shepherd'' (see Tanaka & Ida 1999).

Apart from their radius, the planetesimals are characterized by a number of bulk properties like their density $\rho_{\rm b}$ , their tensile strength $\sigma_{\rm T}$ , heat of ablation (vaporization or melting) $Q_{\rm abl}$ and some more material-dependent parameters. This allows the simulation of the fate of different planetesimal types in the atmosphere of the growing protoplanet, as described in the next part.

2.3.2 Interaction with the growing atmosphere

Given a core and the structure of the surrounding atmospheric envelope (pressure P, temperature T, density $\rho$ etc. as a function of the distance from the planetary center), the second module computes the trajectory and destruction of planetesimals inside this region by integrating a system of ordinary differential equations using a fourth-order Runge-Kutta method with automatic step size control (Press et al. 1992). This approach is similar to the one described by Podolak et al. (1988, hereafter referred to as PPR88). We will present this module and its results in more details in an oncoming separate paper (see Mordasini et al. 2005), and restrict ourselves here to a description of the physical aspects we take into account. There are four main mechanisms controlling the fate of a planetesimal:

Gravity: The gravitational attraction is calculated in the two-body approximation (planet - planetesimal) assuming a spherical mass distribution. This is justified because inside $R_{\rm H}$ , where the calculation takes place, the effect of the planet generally largely predominates over the third body effect.

Aerodynamic drag force: Apart from gravity, the aerodynamic drag force $F_{\rm D}$ is the second force that defines the trajectory of the planetesimal. Using the standard formula for the drag force at high Reynolds numbers (see e.g. Landau & Lifshitz 1959) we have:

Thermal mass loss: The basic mechanism of thermal ablation is very simple: the drag force leads to a dissipation of kinetic energy, some fraction of which is used to heat the gas, while the remainder goes into heating up the body. When the heat flux is sufficiently high, the surface of the body becomes so hot that melting or vaporization starts. From these considerations, the simple classical ablation equation known from the study of terrestrial meteors is (Öpik 1958; Bronsthen 1983):

$\begin{displaymath}\frac{{\rm d}M}{{\rm d}t}=-\frac{1}{2}C_{\rm H}\rho v^3 S\frac{1}{Q_{\rm abl}} \cdot \end{displaymath}$

(24)

One of the problems with the classical equation - apart from the neglect of re-radiation from the hot surface and heat conduction into the interior - is that the rate of thermal mass loss is heavily dependent upon the value of $C_{\rm H}$ which can vary by several orders of magnitude depending upon the size of the infalling body (from 10^-5 to $\approx$ 0.5, Svetsov et al. 1995).

For terrestrial, cm-sized meteorites a mean heat transfer coefficient of about 0.1 can reasonably explain the observations (Bronsthen 1983). Whether this value can be extrapolated to impactors of tens or even hundreds of kilometers in size remained for a long time unclear (PPR88). Fortunately, the large number of studies of the impact of comet Shoemaker-Levy 9 (SL9) onto Jupiter have now shown that modelling such an impact with the above equation setting $C_{\rm H}\approx 0.1$ greatly overestimates the mass loss and predicts too small a penetration depth in comparison to detailed hydrodynamic simulations (Sekanina 1993; Field & Ferrara 1994; Ahrens et al. 1994). One can conclude from the SL9 event that for large, km-sized bodies, thermal ablation is of minor importance compared to mechanical destruction by aerodynamical forces (Svetsov 1995, see below) and that $C_{\rm H}$ is small for such large objects, of the order of 10^-4 - 10^-3 (see e.g. Field & Ferrara 1995, Zahnle & Mac Low 1994). Anyway, it is clear that thermal ablation must be considered in all cases even when fragmentation occurs, as it controls the final dissolution of the impactor (Korycansky et al. 2000), especially if one is interested not only in the energy deposition profiles (which was of major interest in the SL9 case) but also in the mass deposition (as is the case here, since we aim to calculate the fraction of mass dissolved inside the atmosphere and the actual mass of the solid core).

In order to obtain a realistic model of the energy input, we follow Zahnle (1992), replacing the conductive heat influx term per surface unit $\propto$ $\frac{1}{2}\rho v^3$ in the supersonic continuum flow by $C_{\rm H,rad}(T_{\rm ps})\sigma T_{\rm ps}^4$ whenever the latter is smaller, a way to take into account that in this regime the energy input is proportional to the radiative flux from the shock front (Öpik 1958). In this equation $T_{\rm ps}$ is the post-shock temperature, which is directly computed as a function of the local atmospheric properties and the impactor velocity by solving numerically the normal shock wave jump conditions for a real gas (Landau & Lifshitz 1959), using the equation of state of Saumon et al. (1995). Our results for $T_{\rm ps}$ are very similar to the ones obtained by Chevalier & Sarazin (1994). $C_{\rm H,rad}(T_{\rm ps})$ is a coefficient which denominates the fraction ( $\leq$ 1) of radiation reaching the planetesimal surface due the screening effect of ablated material. We use here the results of Bibermann et al. (1980).

We also include, where appropriate, radiation from the undisturbed atmosphere (PPR88) and a convective heat input mechanism which is proportional to the temperature difference between the surrounding gas and the impactor surface temperature. We use the expressions presented by Sibulkin (1952), corrected - if necessary - for a turbulent boundary layer. A similar expression can also be found in Svetsov (1995).

The actual surface temperature of the planetesimal and the amount of ablated mass are obtained by equating the energy influx to the amount of re-radiation and energy used for ablation. For the reasons presented in PPR88 we neglect the heat conduction towards the interior of the planetesimal. For vaporization, ${\rm d}M/{\rm d}t$ is linked to the wall temperature via the Knudsen-Langmuir formula (Bronsthen 1983, PPR88).

The thermal mass loss rates found for large, non-fragmented bolides are in our simulation considerably lower than in a model using the simple ablation equation and $C_{\rm H}=0.1$ .

Mechanical destruction: In deeper layers of the atmosphere the aerodynamic pressure acting on the front of the impactor can overcome the internal strength of the body (either tensile strength or self-gravity, see PPR88). The bolide will then approximately act as a fluid and undergo deformation, fragmentation and mechanical ablation (Svetsov et al. 1995). Numerical simulations of the SL9 event have proven the prevalent importance of these effects for understanding the atmospheric destruction of a large impactor (see e.g. Zahnle & Mac Low 1994). These effects occur in two different regimes. In the first one (static regime, when the aerodynamic pressure loading builds up on a time-scale larger than the time required for a sound wave to travel through the body, see Svetsov et al. 1995; Korykansky et al. 2000), we use a model which combines the effects of lateral spreading (using the so called "pancake'' model of Zahnle 1992) and the growth of Rayleigh-Taylor instabilities (Sharp 1984; Roulsten & Ahrens 1997; Korycansky et al. 2000) at the front of the fluidized impactor to obtain a multi-staged fragmentation model. In the second (dynamical regime, when the pressure loading is dynamical), where no lateral spreading is observed (Svetsov et al. 1995; Korycansky et al. 2000) but where instead mass is removed at the front of the bolide, we use a simple model for mechanical ablation which is also based on the growth of RT instabilities. For large impactors which are found in this regime, mass loss by mechanical ablation overcomes thermal ablation usually by a large factor.

After the calculation, the three dimensional energy and mass deposition profiles are converted back into one dimension to make them usable by our third module described in the next section. The planetesimal impact code is also used to calculate the capture radius, a calculation performed in a same way as described in P96.

2.4 Protoplanet structure and evolution

2.4.1 Equations of the internal structure

The third module calculates the internal structure of the planet including a growing core (including the energy deposited by accreted planetesimals) and the gas accretion due to both the contraction of the envelope and the increase of the outer radius of the planet. The standard equations of planet evolution are solved:

$\displaystyle {{\rm d} r^3 \over {\rm d} m}$	=	$\displaystyle {3 \over 4 \pi \rho},$	(25)
$\displaystyle {{\rm d} L \over {\rm d} m}$	=	$\displaystyle -{{\rm d} U \over {\rm d} t} + {P \over \rho^2} {{\rm d} \rho \over {\rm d} t} + \epsilon_{\rm acc},$	(26)
$\displaystyle {{\rm d} P \over {\rm d} m}$	=	$\displaystyle {- G ( m + M_{\rm core}) \over 4 \pi r^4},$	(27)
$\displaystyle {{\rm d} T \over {\rm d} P}$	=	$\displaystyle \nabla_{\rm ad} ~~~ {\rm or} ~~~ \nabla_{\rm rad},$	(28)

Two approximations can be used regarding the fate of the matter of planetesimals, after their destruction inside the envelope. In the "sinking'' approximation (see P96), this mass is assumed to slowly sink toward the core, leading to an extra term in the core luminosity. In the "no sinking'' approximation, the matter is assumed to remain inside the envelope. Note that, as in P96, the matter is added to the core mass in Eq. (27) in both cases. Under this aspect, the sinking approximation appears to be more self-consistent.

Two inner boundary conditions are necessary to solve Eqs. (25)-(28): the core radius is set to $R_{\rm core} =\left( { M_{\rm core} \over 4/3 \pi \rho_{\rm core} } \right)^{1/3}$ and the core luminosity $L_{\rm core}$ is given by the sum of the remaining kinetic energy plus the corresponding potential energy of planetesimals, after having passed through the atmosphere. The density of the core is fixed to 3.2 g/cm³ as in P96.

The mass of the envelope is given by requiring that the outer radius of the planet is equal to the minimum of the Hill radius and the accretion radius $R_{\rm accr} \equiv { G M_{\rm planet}\over C_{\rm s}^2 }$ where $C_{\rm s}$ is the sound velocity inside the disc at the location of the planet (see P96).

Finally, two outer boundary conditions are necessary. They are given by requiring that the disc and the planet join smoothly at the outer radius:

This condition is valid as long as the disc can supply enough mass to keep the outer radius equal to the Hill (or the accretion) radius, i.e. if the gas accretion rate is below the maximum accretion rate (see Sect. 2.4.2).

2.4.2 Solid accretion rate

We use the expression of Greenzweig & Lissauer (1992) to calculate the solid accretion rate:

$\begin{displaymath}{{{\rm d} {M_{\rm solids}}} \over {{\rm d} {t}}} = \pi R_{\rm c}^2 \Sigma_{\rm P}\Omega F_{\rm g} , \end{displaymath}$

(31)

2.4.3 Gas accretion rate

The condition $R(m=M_{\rm planet}) = {\rm min}(R_{\rm H},R_{\rm accr})$ can be fulfilled only if the disc can supply enough gas to the planet. Once a gap opens in the disc, the maximum gas accretion rate is set to the rate given by (Veras & Armitage 2004):

$\begin{displaymath}{\dot{M}_{\rm gas,max} \over \dot{M}_{\rm disc} }= 1.668 \lef... ...ight)^{1/3} {\rm e}^{- {M_{\rm planet}\over 1.5 M_J}} + 0.04 , \end{displaymath}$

(33)

2.5 Initial conditions and physical assumptions

We start our calculations with a core of $0.6\; M_\oplus$ , at a distance $a_{\rm start}$ from the star. The initial surface density of the disc is usually taken as a power law, the total mass being given (for fixed boundaries) by the normalisation at 5.2 AU. Its lifetime is given by the value of $\alpha$ and the rate of photo-evaporation. Finally, the solid surface density is fully characterized by the dust-to-gas ratio f_D/G.

In addition, two numerical parameters have to be specified, $f_{\rm I}$ (reduction of type I migration) and $f_{\Lambda}$ , the numerical factor in the expression of the momentum transfer between the planet and the disc. As shown in Alibert et al. (2004), modifications in the outer planet boundary conditions due to the gap formation have a small effect. Therefore, the gap essentially limits the gas accretion onto the planet, and we will restrict ourselves to the case $f_\Lambda = 0$ .

Some other parameters of the model, which are not changed in the calculations presented here are summarized in Table 1.

2.6 Numerical tests

Various tests have been performed to validate our entire model. Concerning the disc part of our code, we have compared the resulting functions $\tilde{\nu}(r,\Sigma)$ for different values of $\alpha$ to the analytical fits given in PT99. Figure 1 shows a comparison for $\alpha = 10^{-3}$ . We recall that PT99 mention a difference between the numerical results and their fits of the order of 50% at most. This is also comparable to the maximum difference between our results and the fits.

$\begin{figure} \par\includegraphics[height=7.7cm,width=7.5cm,clip]{2032fig1.ps}\end{figure}$	Figure 1: $\dot{M_{\rm st}}(r,\Sigma)$ given by our calculations (solid lines), and the analytical fits by PT99 (dotted lines), for $\alpha = 10^{-3}$ . The values of the radii (in AU) at which these rates are computed are indicated on the plot.
Open with DEXTER

The interior structure module has been used to reproduce the results of BP86. Figure 2 shows the evolution of a forming Jupiter, for a planet below the cross-over mass (until the mass of accreted planetesimals and the mass of accreted gas are equal). In this calculation, the accretion rate of planetesimals is constant, and there is no interaction between the envelope and the planetesimals. The boundary conditions are those of case 1, 6 and 7 of BP86. Figure 2 shows the central temperature and density for these models, which is in good agreement with their results (see Fig. 1 in BP86). We conclude that our internal structure module works properly.

Our planetesimal accretion module has been tested extensively by comparing results of simulations of impacts into terrestrial and Venusian atmosphere (Revelle 1979; Hills & Goda 1993; Zahnle 1992), as well as into Jupiter's atmosphere (Zahnle & Mac Low 1994). Figure 3 shows one of these tests: energy deposition profiles are plotted for compact ice comets of different initial radii hitting the Jovian atmosphere. The initial conditions are chosen to match the SL9 event. As expected, larger bolides penetrate more deeply and have a larger peak energy deposition. Our results are virtually identical to those of the analytical model of Zahnle & Mac Low (1994), which is not surprising, as the pancake model is used in both cases. The edge at an altitude of about 90 km comes from the radiative limit for thermal ablation and appears at the same point as in Zahnle & Mac Low (1994) and Crawford (1996). The thick line schematically represents an energy release profile obtained from a high-resolution hydrodynamic simulation by Mac Low & Zahnle (1994), predicting a very similar peak energy deposition altitude and value. A model using the classical ablation equation and a high heat transfer coefficient ( $C_{\rm H} \sim 0.1$ ), but no mechanical effects, cannot reproduce this profile.

$\begin{figure} \par\includegraphics[height=6.0cm,width=7.5cm,clip]{2032fig3.eps}\end{figure}$	Figure 3: Energy deposition profiles of SL9-type impactors with different initial radii in the Jovian atmosphere. For comparison, a schematic representation of the energy deposition profile of a 0.5 km impactor in the high-resolution hydrodynamic simulations of Mac Low & Zahnle (1994) is given (thick line).
Open with DEXTER

Finally, the entire code has been tested using the same initial conditions as P96 (case J1) turning disc evolution and migration off. In this case, we obtain a cross-over time of $\sim$ 8 Myr, in close agreement with their result (compare Fig. 4 with Fig. 1 of P96). We conclude that our code properly reproduces giant planet formation in the absence of migration and disc evolution.

$\begin{figure} \par\includegraphics[height=7.7cm,width=6.3cm,clip]{2032fig4.ps}\end{figure}$	Figure 4: Mass of accreted planetesimals ( $M_{\rm heavy}$ , solid line), mass of accreted gas (dotted line), and total mass of the planet (dashed line) as a function of time, for conditions similar to case J1 of P96.
Open with DEXTER

3 Results

3.1 Formation timescales

As shown in BP86 and P96, the formation timescale is essentially given by the time needed to reach the cross-over mass. To quantify the effects of migration and disc evolution on planet formation we compare different models, all reaching the cross-over phase at the same location in the disc, namely 5.5 AU.

For comparison purpose with P96, we consider an initial disc density profile given by $\Sigma \propto r^{-2}$ , like in P96, again for comparison purposes. The constant is chosen to yield $\Sigma = 525$ g/cm² at 5.2 AU; this corresponds to $\sim$ 3 times the surface density in the minimum mass solar nebula at the position of Jupiter. This surface density profile yields isolation masses that do not depend on the distance to the sun. We do not take into account photo-evaporation, and we start with an embryo of 0.6 $M_\oplus$ initially at 7, 8 or 15 AU, depending upon the choice of the parameter $f_{\rm I}$ (0.01, 0.03 and 0.1 respectively). The viscosity parameter $\alpha$ is set to $2\times10^{-3}$ and the dust-to-gas ratio is equal to 1/70 for disc mid-plane temperature below 150 K, and 1/280 in the opposite case. Finally, we use for these simulations the no-sinking approximation (see P96).

$\begin{figure} \par\includegraphics[height=7.7cm,width=6.5cm,clip]{2032fig5.ps}\end{figure}$

Figure 5: Mass of accreted planetesimals (lines starting at $0.6\; M_\oplus$ ), and mass of accreted gas as a function of time, until the cross-over mass is reached. Solid lines: starting at 8 AU (with $f_{\rm I} = 0.03$ ), dotted lines: starting at 15 AU (with $f_{\rm I} = 0.1$ ) and dashed lines: in situ model. The model starting at 7 AU is very close to the one starting at 8 AU and is not shown.

Open with DEXTER

Figure 5 shows the mass of accreted planetesimals and gas as a function of time for the models considered here. Note that the mass of accreted planetesimals does not correspond to the core mass since some fraction of them are being destroyed while traversing the envelope and never reach the core. For the in situ case (without disc evolution), the time to reach the cross-over mass is around 30 Myr, much longer than the typical disc lifetimes (Haisch et al. 2001). For the simulations with migration and disc evolution, the corresponding timescales are respectively $\sim$ 1, $\sim$ 1 and $\sim$ 2.2 Myr for an embryo starting at 7, 8 or 15 AU respectively.

Allowing for migration and disc evolution, we obtain a time to reach the cross-over mass of about 1-3 Myr, i.e. more than ten times faster than in our identical model in which migration and disc evolution have been switched off. As explained in Alibert et al. (2004), the main reason for this speed-up is found to be migration which prevents the severe depletion of the feeding zone as observed in the in situ formation model. Migration therefore suppresses the need to accrete a gas envelope in order to extend the feeding zone (phase 2 of P96) and the planet reaches cross-over mass much faster.

$\begin{figure} \par\includegraphics[height=7.6cm,width=6.8cm,clip]{2032fig6.ps}\end{figure}$	Figure 6: Distance to the central star for the same models as in Fig. 5. The kinks around 0.8 and 2 Myr signal the change from type I to type II migrations. Solid line: starting at 8 AU (with $f_{\rm I} = 0.03$ ), dotted line: starting at 15 AU (with $f_{\rm I} = 0.1$ ), and dashed line: in situ model.
Open with DEXTER

As stated in Sect. 2.3.2, the effect of ablation is found to be negligible, disruption of planetesimals occurring very deep in the envelope. The two cases, sinking and no-sinking, considered by P96 differ therefore much less in our calculations than in P96. For the simulation started at 8 AU, the time to reach the cross over mass is found to be of order $\sim$ 1.2 Myr using the sinking approximation, compared to $\sim$ 1 Myr in the no-sinking case. For that reason, and since it is more self-consistent, we will use in the following the sinking approximation, as in P96's standard calculations.

3.2 Formation of extrasolar planets

Depending on the initial parameters used, very different planets can form. Figure 7 shows the results of a set of 1000 simulations, in a mass vs. semi-major axis diagram. The initial parameters for these simulations are $a_{\rm start}= 3$ , 4, 5, 6, 7.5 10, 15 and 20 AU, f_D/G = 1/25, 1/50, 1/70, 1/100 and 1/200, and a photoevaporation rate of $10^{-9}\; M_\odot/$ yr, $2 \times 10^{-9} \;M_\odot/$ yr, $4 \times 10^{-9}\; M_\odot/$ yr, $10^{-8} \;M_\odot/$ yr and $2 \times 10^{-8}\; M_\odot/$ yr. The initial disc is given by $\Sigma \propto r^{-3/2}$ , the normalization at 5.2 AU being between 100, 300, 500, 700 and 900 ${\rm g/cm^2}$ . This corresponds to disc masses of $\sim$ 0.01 to $\sim$ $0.1 \;M_\odot$ between 0.25 and 50 AU. The viscosity parameter is fixed to $2\times10^{-3}$ , and we use $f_{\rm I} = 0.03$ for the whole set of simulations. The edges of the disc are at 0.25 and 50 AU, so planets whose inner radius of the feeding zone is below 0.25 AU are considered to have migrated inside the star, and are represented by dots on the vertical axis.

One must keep in mind that the diagram shown here does not take into account any statistical weighting of the different initial parameters, so one should be careful when comparing the obtained distribution with the observed one. We present it simply to illustrate the possible diversity of planets that can form in our model. In a forthcoming paper we will discuss the probability of occurrence of these planets.

$\begin{figure} \par\includegraphics[height=6cm,width=6cm,clip]{2032fig7.ps}\end{figure}$

Figure 7: Final mass and semi-major axis (dots), for a set of simulations, using different initial parameters (see text for details). The crosses are the parameters of the observed exoplanets as of November 2003. The big circle represents the final parameters of the model detailed in Sect. 3.2. The dashed line gives the inner limit of our simulations (when the inner radius of the planet feeding zone is below 0.25 AU).

Open with DEXTER

$\begin{figure} \par\includegraphics[height=7.6cm,width=7.8cm,clip]{2032fig8.ps}\end{figure}$	Figure 8: Mass of the different components for the model of Sect. 3.2. Solid line: mass of accreted planetesimals, dashed line: mass of H/He, and heavy solid line: mass of the disc.
Open with DEXTER

$\begin{figure} \par\includegraphics[height=7.6cm,width=7.3cm,clip]{2032fig9.ps}\end{figure}$	Figure 9: Distance to the Sun for the model of Sect. 3.2. The kink around 1 Myr corresponds to the change from type I to type II migration.
Open with DEXTER

As an illustration taken from this set of simulations, we provide some details on the formation process of one particular planet. This object evolved from the following: $a_{\rm start}= 15$ AU, a disc photo-evaporation rate $4 \times 10^{-9}\; M_\odot/$ yr, an initial normalisation of the disc at 5.2 AU equal to 500 g/cm², and f_D/G = 1/70. With a final mass of $\sim$ $3.5\; M_{\rm J}$ , and a final distance $\sim$ 2.5 AU from its star, this planet is indicated by the circle in Fig. 7. Figure 8 shows the evolution of the mass of accreted gas, the mass of accreted planetesimals, as well as the mass of the disc, and Fig. 9 gives the distance to the central star as a function of time. The cross-over mass is reached after $\sim$ 1.6 Myr, and shortly after, due to gap formation, the accretion rate of gas is limited to its maximum value, which decreases with decreasing disc mass. The formation process ends after 5.5 Myr when the disc has essentially disappeared. Note however, that $90 \%$ of the final mass of the planet has already been reached after $\sim$ 4 Myr. At the end of the simulation, the planet has accreted $\sim$ $30 \;M_\oplus$ of planetesimals. As stated in Sect. 3.1, this mass does not correspond to the mass of the core, since some of the accreted planetesimals may have been destroyed during the crossing of the envelope. Our model provides an estimate of the core mass by tracking the location where the planetesimals are destroyed (second module). For this simulation, $\sim$ $6.4 \;M_\oplus$ of accreted planetesimals reach the core without being destroyed in the envelope. Ignoring further processes, this provides the mass of the core which is similar to the one inferred from Jupiter's internal structure models (Guillot et al. 2004). The total content in heavy elements will finally depend on the metallicity of the accreted gas. For solar composition (prior to the formation of planetesimals), the accreted gas adds an additional $\sim$ $7\; M_\oplus$ of heavy elements.

The migration of the planet can be divided into three phases. Before $\sim$ 1 Myr the planet undergoes type I migration at which time a gap opens and migration switches to type II. Shortly after $\sim$ 2 Myr, the mass of the planet becomes non-negligible compared with the disc mass and migration slows down and eventually stops at the time the disc disappears.

4 Summary and discussion

We have presented a new code devoted to the calculation of giant planet formation models, taking into account the protoplanetary disc structure, as well as the migration of the forming planet. These calculations show that the formation of giant planets, at least in the first phase until runaway gas accretion, can be greatly sped-up if one takes into account the effect of migration. This is mainly due to the suppression of phase 2 as described in P96. Using an initial disc model similar to the one of P96, we obtain a time to reach the cross-over mass of the order of 1 Myr, significantly shorter than the typical disc lifetimes. Moreover, this speed-up due to migration has been found to be robust against reasonable changes in various parameters (Alibert et al. 2004). Therefore, migration not only accounts for the orbital distribution of extra-solar planets, but also considerably shortens the formation timescale in the core-accretion scenario. The formation of giant planets in this scenario is therefore in excellent agreement with inferred disc lifetimes without having to consider discs significantly more massive than the minimum mass solar nebula.

We note that the speed-up due to migration does not preclude other effects to further reduce this timescale. For example, reducing the dust opacity would decrease the critical mass (Ikoma et al. 2001), thus leading to another reduction in the formation timescale. Furthermore, such a reduced opacity could account for the existence of giant planets with small central cores.

Using different initial conditions can lead to the formation of a wide variety of planets. However, the comparison of our results with observations of extrasolar planets needs to take into account the probability distribution of various initial conditions. Work is under progress to obtain, using a Monte-Carlo approach, synthetic distributions which can then be compared in a statistical way with the observed distributions.

These calculations are however subject to some uncertainties, among them the simplified treatment of the disc of planetesimals, or the calculation of the ejection rate, which directly determines the final heavy elements content. Work is in progress to improve these aspects.

Finally, it seems to be very difficult to form a planet, and to prevent it from spiraling into the sun if the amount of type I migration as computed today is not reduced by a factor of at least 10.

Models of giant planet formation with migration and disc evolution

1 Introduction