Consequences of the simultaneous formation of giant planets by the core accretion mechanism

O. M. Guilera; A. Brunini; O. G. Benvenuto

doi:10.1051/0004-6361/201014365

Home

All issues

Volume 521 (October 2010)

A&A, 521 (2010) A50

Full HTML

Free Access

Issue		A&A Volume 521, October 2010


Article Number		A50
Number of page(s)		14
Section		Planets and planetary systems
DOI		https://doi.org/10.1051/0004-6361/201014365
Published online		20 October 2010

A&A 521, A50 (2010)

Consequences of the simultaneous formation of giant planets by the core accretion mechanism

O. M. Guilera^1,2, - A. Brunini^1,2, - O. G. Benvenuto^1,2,

1 - Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque s/n, B1900FWA La Plata, Argentina
2 - Instituto de Astrofísica de La Plata, IALP, CCT-CONICET-UNLP, Argentina

Received 5 March 2010 / Accepted 17 May 2010

Abstract
Context. The core accretion mechanism is presently the most widely accepted cause of the formation of giant planets. For simplicity, most models presently assume that the growth of planetary embryos occurs in isolation.
Aims. We explore how the simultaneous growth of two embryos at the present locations of Jupiter and Saturn affects the outcome of planetary formation.
Methods. We model planet formation on the basis of the core accretion scenario and include several key physical ingredients. We consider a protoplanetary gas disk that exponentially decays with time. For planetesimals, we allow for a distribution of sizes from 100 m to 100 km with most of the mass in the smaller objects. We include planetesimal migration as well as different profiles for the surface density $\Sigma$ of the disk. The core growth is computed in the framework of the oligarchic growth regime and includes the viscous enhancement of the planetesimal capture cross-section. Planet migration is ignored.
Results. By comparing calculations assuming formation of embryos in isolation to calculations with simultaneous embryo growth, we find that the growth of one embryo generally significantly affects the other. This occurs in spite of the feeding zones of each planet never overlapping. The results may be classified as a function of the gas surface density profile $\Sigma$ : if ${\Sigma \propto r^{-3/2}}$ and the protoplanetary disk is rather massive, Jupiter's formation inhibits the growth of Saturn. If ${\Sigma \propto r^{-1}}$ isolated and simultaneous formation lead to very similar outcomes; in the the case of ${\Sigma \propto r^{-1/2}}$ Saturn grows faster and induces a density wave that later accelerates the formation of Jupiter.
Conclusions. Our results indicate that the simultaneous growth of several embryos impacts the final outcome and should be taken into account by planet formation models.

Key words: planets and satellites: formation - planet-disk interactions - methods: numerical

1 Introduction

The core instability mechanism is usually considered the way in which giant planets formation proceeds. This mechanism considers planetary formation as a consequence of a two-step process. A solid material seed with a mass of the order of that of the Moon is immersed in a protoplanetary disk. This object will be the core of the planet. The protoplanetary disk has a population of planetesimals coexisting with a gaseous component. The Moon-sized object begins to accrete solid material from its feeding zone. The feeding zone is often assumed to be an annulus extended across few times the Hill radius $R_{\rm H}$ (defined to be $R_{\rm H}= a ( M / 3 M_* )^{1/3}$ , where a is the radius of the orbit of the planet of mass M, while M_* is the mass of the central star of the system) on both sides of the orbit of the protoplanet. The core gravitationally binds a tiny amount of gas and the whole protoplanet remains in hydrostatic and thermal equilibrium. The gas accretion rate is initially far lower than that of solids. This is true until the core reaches a mass of ${\approx}10~M_{\oplus}$ , where $M_{\oplus}$ is the Earth mass. At these stages, the amount of gas bound to the protoplanet becomes comparable to the core mass. It is then no longer possible to be in thermal equilibrium conditions and the protoplanetary envelope begins to contract. This corresponds to the onset of the runaway gas accretion. On a very short timescale, the planet accretes a large amount of gas and reaches its final mass. This sequence of events was envisaged by Mizuno (1980) by employing static models and later with evolutionary models by Bodenheimer & Pollack (1986) and Pollack et al. (1996). Core instability calculations of giant planet formation have been carried out by many groups, e.g., Alibert et al. (2005a,b), Hubickyj et al. (2005), and Dodson-Robinson et al. (2008).

Fortier et al. (2007, 2009) were the first to consider the oligarchic growth regime for the accretion of planetesimals. This regime is known to provide a good description of planetesimal dynamics when the embryos present in the nebula have masses much lower than the Earth mass. Thus, this is the correct description of the accretion regime for almost the entire formation process. Oligarchic growth takes into account the effects of the gravitational perturbations of growing planets. In this regime, it is known that protoplanets grow on a longer timescale that in other accretion regimes such as those considered by, e.g., Pollack et al. (1996). One of the most remarkable characteristics of the core instability mechanism is that runaway gas accretion occurs at a core mass value largely independent of the final mass of the planet. This is in close agreement with the estimate of the core masses of the giant planets of our Solar System (see Guillot 2005; Militzer et al. 2008).

In spite of being considered the most likely mechanism of planetary formation, core instability may have some serious difficulties. One of the main drawbacks was considered to be the long timescale for planetary formation. In many simulations, planets form on a time interval longer than (or need timescales uncomfortably close to) the dissipation timescale of the protoplanetary nebula, usually considered to be of the order of 10 Myr (Hillenbrand 2005). In any case, we should recall the number of simplifying assumptions considered when constructing the models that lead to this apparent paradox. For example, most of the calculations available in the literature considered a single size distribution of planetesimals. Relaxing this assumption in favor of a size distribution for which most of the mass belongs to small planetesimals, Benvenuto et al. (2009) demonstrated that the four giant planets of the Solar System can be formed in a timely way with core masses in close agreement with current estimate.

Another usual assumption in detailed simulations of planetary growth is that each planet grows alone in the disk. This would be correct if the population of planetesimals to be accreted by one planet were not appreciably perturbed by the presence of another embryo. At first sight, it may be understood that this is the case if the feeding zone of each planet does not overlap the one corresponding to any other planet. However, this is not the case if we include planetesimal migration. This process leads to a net inward motion of planetesimals. A planet will perturb the swarm of planetesimals that may be later accreted by another planet moving along an inner orbit. Moreover, as we show below, even the presence of an inner planet will be able to affect the accretion process of an outer object. To our knowledge, no calculation has been performed to quantitatively analyze these effects. It is the main aim of the present paper to perform a first step towards filling this gap. In doing so, we consider the simultaneous growth of two planets at distances from the central star as those of Jupiter and Saturn. We consider a distribution of sizes of planetesimals and also different profiles for the surface density of solids and gas corresponding to the protoplanetary disk. We also consider two values of the gas to solid ratio.

Oligarchic growth predicts the simultaneous formation of many embryos on orbits separated by about 10 $R_{\rm H}$ from each other. Thus, starting with this initial configuration would be more realistic than the one adopted here. However, to perform a fully detailed simulation in this context we would have to be able to compute the merging of embryos, which in the presence of a gaseous envelope for each protoplanet is a very complex process. This is beyond the scope of the present paper. In any case, even with the adopted initial configuration we expect to realistically estimate the effect of the simultaneous growth of more than one planet. In this work, we do not consider planetary migration.

This paper is organized as follows: in Sect. 2, we describe the protoplanetary disk inside which giant planets grow. Section 3 is devoted to describing the oligarchic growth regime. Section 4 describes how we compute the growth of each protoplanetary embryo. With all the ingredients of our model quantitatively defined, in Sect. 5 we present our results corresponding to the formation of Jupiter and Saturn, comparing it with the cases of isolated formation. In Sect. 6, we explore different profiles for the surface density of the nebula, and in Sect. 7 we apply our model to the standard Hayashi nebula. Finally, in Sect. 8 we present a discussion of our results and some concluding remarks.

2 A brief description of the protoplanetary disk model

We consider an axisymmetric disk with inner and outer radii of $R_{\rm min}=0.4$ AU and $R_{\rm max}=30$ AU, respectively. The disk has both gaseous and planetesimal components. For the surface density profile of both of them, we consider a power law distribution of the form

	$\displaystyle \Sigma_{\rm g} \propto a^{-p}$	(1)
	$\displaystyle \Sigma_{\rm s} \propto a^{-p},$	(2)

with a temperature profile given by

$\displaystyle T \propto a^{-1/2}.$

(3)

It can be shown that the volumetric density distribution at the disk plane is

$\displaystyle \rho_{\rm g} \propto a^{-p - 5/4}.$

(4)

Here we consider a density distribution for the size of planetesimals. We study 31 different sizes with radii between 100 m and 100 km in steps selected such that the quotient of masses of consecutive sizes is a factor of two, in a similar way to Brunini & Benvenuto (2008). We assume that the material that composes the planetesimals has a density of $\rho_{\rm p}= 1.5~{\rm g~cm}^{-3}$ .

We consider a number of planetesimals per unit of mass distribution given by ${\rm d}N/{\rm d}m \propto m^{-5/2}$ (Kokubo & Ida 1998, 2000; Ormel et al. 2010), for which most of the mass of solids is in the smallest planetesimals ( $M_{\rm T} \propto m^{-1/2}$ ). Numerical simulations indicate that the mass distribution of planetesimals may be represented by a single or piecewise power law ${\rm d}N/{\rm d}m \propto m^{-q}$ . For a constant value of q, the total mass of solids is $M_{\rm T} \propto m^{2-q}$ , so if q<2, the total mass of solids is contained in the biggest planetesimals; in contrast, if q>2 the total mass of solids is contained in the smallest planetesimals. Wetherill & Stewart (1993) studied the evolution of a planetesimal system considering an initial population of planetesimals with radius of $\sim$ 10 km that evolved only by means of collisions and fragmentation. They found that a planetesimal size distribution relaxes to a piecewise power law: a population of small planetesimals produced by fragmentation (with $q\sim1.7$ ) and a population of large planetesimals that follow an accretive regime (with $q\sim2.5$ ). Kokubo & Ida (2000) studied the evolution of planetesimal size using N-body simulations finding that, in the oligarchic regime, large planetesimals follow a continuous power-law distribution with 2<q<3. Ormel et al. (2010) performed statistical simulations that include several physical processes such as dynamical friction, viscous stirring, gas drag, and fragmentation, finding that the transition between the runaway growth and oligarchic growth is characterized by a power-law size distribution of mass index $q\sim2.5$ . On the other hand, planetesimal formation continues to be studied and to date it is not well understood how meter to kilometer-sized planetesimals appear in the disk. Hence, the primordial size distribution of these bodies has not yet been established. In the paper of Wetherill & Stewart (1993), where the system evolves only through collisions and fragmentation, meter-sized planetesimals can only appear in the disk as fragments of larger bodies. Ida et al. (2008) performed some calculations showing that the magneto - rotational instability (Balbus & Hawley 1991; hereafter MRI) turbulence is a serious problem for planetesimal formation. They show that it is very difficult to form kilometer-sized planetesimals from centimeter-meter sized particles because of the predominance of an erosive rather than accretive regime. Under these conditions, it does not seem possible that kilometer-sized planetesimals formed by the accretion of smaller ones. However, they also suggest that, if MRI were inefficient, planetesimals may have formed in ``dead zones'', where the interaction between solids and gas is weak (Gammie 1996; Sano et al. 2000). This would favor the survival of centimeter and meter-sized bodies. Kilometer-sized planetesimals may have formed by the accretion of smaller planetesimals, and the coexistence of meter and kilometer sized planetesimals may have been possible. This possibility is important for the formation of the giant planets in the solar system, since smaller planetesimals favor the rapid formation of the solid embryos. Moreover, if our giant planets formed in an environment where the MRI was effective they should have migrated towards the Sun becoming hot giant planets (Matsumura et al. 2009). On the other hand, laboratory experiments show that reaccumulation of fragmentation debris can lead to the formation of planetesimals (Teiser & Wurm 2009). Owing to the present status of the theory of planetesimal formation, especially with regard to the period when they grow from meter to kilometer size, we consider that the size and mass distribution of the planetesimals adopted here is a valid hypothesis.

Planetesimals are affected by gaseous friction, which causes an inwards decay and can alter the solid distribution in the disk. As demonstrated by Thommes et al. (2003) and Chambers (2006), this effect has a strong influence on the timescale of accretion and on the final mass values reached by the planets at different positions in the protoplanetary disk. The mean orbital evolution of planetesimals is given by (Adachi et al. 1976)

$\displaystyle \frac{{\rm d}a}{{\rm d}t} =$

-

$$\displaystyle \frac{2a}{T_{\rm fric}}\bigg(\eta^2 + \frac{5}{8}e^2 +$

(5)

where

$\displaystyle T_{\rm fric}= \frac{8\rho_{\rm p} r_{\rm p}}{3C_{\rm D}\rho_{\rm g}v_k}\cdot$

(6)

Here $\rho_{\rm p}$ is the density of planetesimals, $r_{\rm p}$ is the radius of the planetesimal, $\eta=(v_k - v_{\rm gas})/v_k$ , and $\alpha$ is the exponent of the power-law density of the gas in the disk mid plane ( $\rho_{\rm g} \propto a^{-\alpha}$ ).

The evolution of planetesimal disks follows the equation of continuity

$\displaystyle \frac{\partial \Sigma_{\rm s}}{\partial t} - \frac{1}{a}\frac{\partial}{\partial a} \left(a\frac{{\rm d}a}{{\rm d}t}\Sigma_{\rm s}\right) = F(a),$

(7)

where F(a) describes the sinks of disk material (accretion by the forming planets and solid sublimation across the ice line, as done by Brunini & Benvenuto 2008). For simplicity, we also assume that the gaseous component dissipates following an exponential decay of its density as

$\displaystyle \rho_{\rm g}(a,t)= \rho_{\rm g}(a,0){\rm e}^{-t/\tau}.$

(8)

We set $\tau=6$ Myr as the characteristic timescale for the dissipation of the protoplanetary nebula (Haisch et al. 2001).

3 Accretion onto planets: oligarchic growth regime

The process of the accretion of solids is described by the ``particle in a box'' approximation (Inaba et al. 2001)

$\displaystyle \frac{{\rm d}M_{\rm C}}{{\rm d}t}= \frac{2\pi \Sigma(a_P)R_{\rm H}^2}{P}P_{\rm coll},$

(9)

where $M_{\rm C}$ is the mass of the core, $\Sigma(a_{\rm p})$ is the surface density of solids at the location of the planet, $R_{\rm H}$ is the Hill radius, P is the orbital period, and $P_{\rm coll}$ is the collision probability, which is a function of the core radius, the Hill radius of the planet, and the relative velocity of planetesimals $P_{\rm coll}=P_{\rm coll}(R_{\rm C},R_{\rm H}, v_{\rm rel})$ . For regimes of high $(\hat{e},\hat{i}>2)$ , medium $(0.2<\hat{e},\hat{i}<2)$ and low velocities $(\hat{e},\hat{i}<0.2)$ , $P_{\rm coll}$ is given by

		$$\displaystyle P_{\rm coll\_high} = \frac{(R_{\rm C}+r_{\rm p})^2}{2\pi R_{\rm H...$	(10)
		$\displaystyle P_{\rm coll\_med} = \frac{(R_{\rm C}+r_{\rm p})^2}{4\pi R_{\rm H}^2 \hat{i}}\bigg[17.3 + \frac{232~R_{\rm H}}{(R_{\rm C}+r_{\rm p})}\bigg],$	(11)
		$\displaystyle P_{\rm coll\_low} = 11.3\bigg[\frac{R_{\rm C}+r_{\rm p}}{R_{\rm H}}\bigg]^{1/2},$	(12)

where $R_{\rm C}$ is the radius of the core, $\hat{e}$ and $\hat{i}$ are the reduced eccentricity and inclination, defined by $\hat{e}= ea_P/R_{\rm H}$ and $\hat{i}= ia_P/R_{\rm H}$ , where $\beta=\hat{i}/\hat{e}$ , e and i are the quadratic mean values of the eccentricity and inclination of planetesimals at the feeding zone of the planet, and $I_{\rm F}(\beta)$ and $I_{\rm G}(\beta)$ are functions that can be expressed in terms of elliptic integrals that, in the interval $0 < \beta \le 1$ can also be approximated by (Chambers 2006)

		$\displaystyle I_{\rm F}(\beta) = \frac{1+0.95925 \beta + 0.77251 \beta^2}{\beta(0.13142 + 0.12295\beta)},$	(13)
		$\displaystyle I_{\rm G}(\beta) = \frac{1 + 0.39960\beta}{\beta(0.0369 + 0.048333\beta + 0.006874 \beta^2)}\cdot$	(14)

In contrast to the works of Brunini & Benvenuto (2008) and Chambers (2006), we adopt the results of Inaba et al. (2001), who assume that $P_{\rm coll}$ is given in the whole range of velocities for $(\hat{e},\hat{i})$ as

$\displaystyle P_{\rm coll}= \mbox{min}\big[P_{\rm coll\_med}, (P_{\rm coll\_low}^{-2} + P_{\rm coll\_high}^{-2})^{-1/2}\big].$

(15)

We consider the drag force that planetesimals experience on entering the planetary envelope, which largely increases their capture cross-section. Inaba & Ikoma (2003) found an approximate solution to the equations of motion, which allows a rapid estimation of the critical radius for capture $r_{\rm p}$ as a function of the radius of the captured planetesimal, the density of the gaseous envelope $\rho$ , and the enhanced radius $\tilde{R}_{\rm C}$

$\begin{displaymath}r_{\rm p}= \frac{3 \rho(\tilde{R}_{\rm C}) \tilde{R}_{\rm C}}... ...{v_{\infty}^2 + 2GM_P(\tilde{R}_{\rm C}) / R_{\rm H}}\right), \end{displaymath}$

(16)

where $v_{\infty}$ is the relative velocity of the planet and planetesimal when the two are far apart and $M_P(\tilde{R}_{\rm C})$ is the total mass of the planet contained whitin $\tilde{R}_{\rm C}$ . Inaba & Ikoma (2003) propose replacing $\tilde{R}_{\rm C}$ for $R_{\rm C}$ in the expressions of collision probability

		$$\displaystyle P_{\rm coll\_high} = \frac{(\tilde{R}_{\rm C}+r_{\rm p})^2}{2\pi ...$	(17)
		$\displaystyle P_{\rm coll\_med} = \frac{(\tilde{R}_{\rm C}+r_{\rm p})^2}{4\pi R... ...\hat{i}}\bigg[17.3 + \frac{232~R_{\rm H}}{(\tilde{R}_{\rm C}+r_{\rm p})}\bigg],$	(18)
		$\displaystyle P_{\rm coll\_low} = 11.3\bigg[\frac{\tilde{R}_{\rm C}+r_{\rm p}}{R_{\rm H}}\bigg]^{1/2},$	(19)

and again

$\displaystyle P_{\rm coll}= \mbox{min}\big[P_{\rm coll\_med}, (P_{\rm coll\_low}^{-2} + P_{\rm coll\_high}^{-2})^{-1/2}\big],$

(20)

which gives $P_{\rm coll}$ = $P_{\rm coll}(\tilde{R}_{\rm C},R_{\rm H}, v_{\rm rel})$ . Finally, we derived a generalized version of Eq. (9) following Brunini & Benvenuto (2008)

$\displaystyle \frac{{\rm d}M_{\rm C}}{{\rm d}t}$	=	$\displaystyle \int_{DS} {\rm d}m \int_{FZ} 2\pi \psi(a, R_{\rm H}, a_P)$
		$\displaystyle \times ~\frac{2\pi\Sigma(a,m)R_{\rm H}^2}{P}~P_{\rm coll}(a,m)~a~{\rm d}a,$	(21)

where DS represents the integration over the distribution of sizes of planetesimals and FZ indicates that the integration extends over the feeding zone. The functional form of $\psi$ is the same as used by Brunini & Benvenuto (2008).

3.1 Relative velocities out of equilibrium

The relative velocity, $v_{\rm rel}$ , between a planetesimal and the protoplanet may be described by

$\displaystyle v_{\rm rel}= \sqrt{\frac{5}{8}e^2 + \frac{1}{2}i^2} v_k,$

(22)

where v_k is the Keplerian velocity at the distance a_P. The relative velocity is governed by gravitational stirring caused by the protoplanets and damping caused by gas drag. This can be modeled following Ohtsuki et al. (2002) as

		$\displaystyle \bigg(\frac{{\rm d}\langle e^2 \rangle}{{\rm d}t}\bigg)_{\rm stirr} = \bigg(\frac{M_P}{3bM_{\star}P}\bigg) P_{VS},$	(23)
		$\displaystyle \bigg(\frac{{\rm d}\langle i^2 \rangle}{{\rm d}t}\bigg)_{\rm stirr} = \bigg(\frac{M_P}{3bM_{\star}P}\bigg) Q_{VS},$	(24)

where

		$\displaystyle P_{VS} = \bigg[\frac{73\hat{e}^2}{10\Lambda^2}\bigg] \ln \bigg(1 ... ...g) + \bigg[\frac{72I_{PVS}(\beta)}{\pi\hat{e}\hat{i} }\bigg] \ln (1+\Lambda^2),$	(25)
		$\displaystyle Q_{VS} = \bigg[\frac{4\hat{i}^2 + 0.2 \hat{i} \hat{e}^3}{10\Lambda^2 \hat{e}}\bigg] \ln (1 + 10\Lambda^2)$
		$\displaystyle \qquad\quad +~ \bigg[\frac{72I_{QVS}(\beta)}{\pi\hat{e}\hat{i} }\bigg] \ln(1+\Lambda^2),$	(26)

where $\Lambda^2= \hat{i}(\hat{i}^2 + \hat{e}^2)/12$ , and $I_{PVS}(\beta)$ and $Q_{PVS}(\beta)$ are given in terms of elliptic integrals than can be approximated by (Chambers 2006)

		$\displaystyle I_{PVS}(\beta) = \frac{\beta - 0.36251}{0.061547+0.16112\beta +0.054473\beta^2},$	(27)
		$\displaystyle Q_{PVS}(\beta) = \frac{0.71946-\beta}{0.21239+0.49764\beta +0.14369\beta^2}\cdot$	(28)

Nevertheless, these velocities decline appreciably when we move far away from the protoplanet. Hasegawa & Nakazawa (1990) demonstrated that when the distance from the protoplanet is larger than 3.5-4 times its Hill radius, the excitation of the relative velocities of planetesimals weakens significantly. Thus, we consider that

		$\displaystyle \bigg(\frac{{\rm d}\langle e^2 \rangle}{{\rm d}t}\bigg)_{\rm stir... ...= f(\Delta)\bigg(\frac{{\rm d}\langle e^2 \rangle}{{\rm d}t}\bigg)_{\rm stirr},$	(29)
		$\displaystyle \bigg(\frac{{\rm d}\langle i^2 \rangle}{{\rm d}t}\bigg)_{\rm stir... ...= f(\Delta)\bigg(\frac{{\rm d}\langle i^2 \rangle}{{\rm d}t}\bigg)_{\rm stirr},$	(30)

where

$\displaystyle f(\Delta) = \Bigg(1 + \bigg\vert\bigg\vert \frac{\Delta}{4~R_{\rm H}}\bigg\vert\bigg\vert^5\Bigg)^{-1},$

(31)

where $\Delta$ represents the distance from the protoplanet, and $f(\Delta)$ guarantees that the velocity profile of planetesimals is smooth along the entire disk and the planetary excitation on planetesimals is restricted to the protoplanetary neighborhood. This is important for an adequate solution of the solid distribution.

The friction caused by the gaseous component of the protoplanetary disk decreases the orbital eccentricities and inclinations of planetesimals at a rate given by (Adachi et al. 1976)

		$\displaystyle \bigg(\frac{{\rm d}e}{{\rm d}t}\bigg)_{\rm gas} = \frac{\pi e r_{\rm p}^2C_{\rm D}\rho_{\rm g}v_k}{2m_{\rm p}}\bigg( \eta^2 + v_{\rm rel}^2\bigg),$	(32)
		$\displaystyle \bigg(\frac{{\rm d}i}{{\rm d}t}\bigg)_{\rm gas} = \frac{\pi i r_{\rm p}^2C_{\rm D}\rho_{\rm g}v_k}{4m_{\rm p}}\bigg( \eta^2 + v_{\rm rel}^2\bigg),$	(33)

where $C_{\rm D}$ is a dimensionless coefficient that describes the gaseous friction (it is ${\approx}1$ for the case of spherical bodies), $m_{\rm p}$ is the mass of planetesimals, $\rho_{\rm g}$ is the density of nebular gas at a distance a from the central star, and $\eta$ is the ratio of the gas velocity to the local Keplerian velocity, given by

$\displaystyle \eta= \frac{v_k - v_{\rm gas}}{v_k}= \frac{\pi}{16}(\alpha + \beta)\bigg(\frac{c_{\rm s}}{v_k}\bigg)^2,$

(34)

where $\beta$ is the exponent of the power-law profile of temperature of the nebula ( $T(a) \propto a^{-\beta}$ ) and $c_{\rm s}$ is the local velocity of sound.

4 Models of giant planet formation

As previously stated, the formation of planets was computed in the framework of the core instability scenario, in which three aspects are of key importance:

a.: The solid accretion rate that increases the protoplanetary core mass. This is fundamental because it largely determines its formation timescale. The rate of solid accretion also provides the energy release necessary to support the gaseous envelope against gravitational contraction. This rate was provided by the particle in a box approximation (Eq. (21)) as described above.
b.: The rate of gas accretion and the evolution of the envelope. The calculation of the structure of the gaseous protoplanetary envelope was performed by solving the standard equations of stellar structure (see below for additional details).
c.: The interaction between planetesimals and the gaseous envelope. For simplicity, we assumed that they do not disaggregate on entering the protoplanetary envelope and are deposited in the nucleus. The gravitational energy release due to planetesimal accretion was incorporated at the bottom of the gaseous envelope.

The equations that govern the evolution of the protoplanetary gaseous envelope are those of stellar evolution theory, namely

			(35)
			(36)
		$% latex2html id marker 4950 $\displaystyle \frac{\partial L_r}{\partial m_r} = \... ...m pl}-T\frac{\partial S}{\partial t}\qquad \mbox{equation of energetic balance}$$	(37)
		$% latex2html id marker 4951 $\displaystyle \frac{\partial T}{\partial m_r} = - \frac{G m_r T}{4\pi r^4 P} \nabla\qquad \mbox{equation of energy transport,}$$	(38)

where $\rho$ is the density of the envelope, G is the universal gravitational constant, $\epsilon_{\rm pl}$ is the energy release rate due to the accretion of planetesimals, S is the entropy per unit mass, and $\nabla \equiv \frac{{\rm d} \ln T}{{\rm d} \ln P}$ is the dimensionless temperature gradient, which depends on the type of energy transport.

We assumed that the outer boundary conditions change in full accordance with the evolution of the protoplanetary disk, i.e., are function of time. At the outermost point of the protoplanet model, we define $\rho= \rho_{\rm neb}(t)$ . However, we assumed the temperature to be fixed. The remaining ingredients of the calculations presented below are as in Fortier et al. (2007, 2009).

4.1 Methodology of the code

The goal of this work was to construct self-consistent models of giant planet formation inside a realistic protoplanetary disk that evolves with time and, in particular, to analyze how the simultaneous growth of more than one object affects the growth of each individual object. For this purpose, we separated the calculation into several steps. We generalized the code that computes the in-situ formation of one giant planet to compute the simultaneous formation of N planets. Then, we coupled this code with the one that solves the disk evolution, which, in turn, has the role of being the main program.

Once the parameters model have been defined, the calculation starts with the initial disk model. The Hill's radius of each planet was then defined and the migration velocities computed for each planetesimal size. These velocities are employed to solve the surface density of the disk. Then we computed the evolution of the gaseous component of the disk. We then computed a growth step for each planet. To do so, we had to consider the contributions due to each planetesimal size inside the planetary feeding zone and also the collision probabilities for each planet. For each planet, we computed the accretion rate of solid and the new structure for the gaseous envelope. We adopted the new masses for each planet and computed the sinks of planetesimals, necessary for solving the continuity equation. This procedure took computational cycle, and when necessary the whole procedure was repeated. To achieve a stable and accurate sequence of models, for the next time step we adopted the minimum timescale of all the ingredients considered in the calculation.

We considered the simulation to have ended when all planets had reached the adopted mass values. In the case that one or more planets had reached their final masses while some others had not, we computed the evolution of the whole system. In this case, we did not consider any additional growth of the planets that reached their final masses but still took them into account when computing the evolution of the disk, perturbing the population of planetesimals.

5 Application to the simultaneous formation of Jupiter and Saturn

We apply the model described in Sects. 2-4 to quantitatively analyze the effects of the simultaneous formation of Jupiter and Saturn comparing with the results corresponding to the case of isolated formation.

5.1 A standard nebula with profile ${\Sigma \propto r^{-3/2}}$

We first considered the standard model of the Solar nebula proposed by Hayashi (1981). It states that

		$\displaystyle \Sigma_{\rm s}(a) = \left\lbrace \begin{array}{ll} 7.1 \left(\dis... ...}}\right)^{-3/2}~{\rm g~cm}^{-2}, & \mbox{$a>2.7$ ~{\rm AU}} \end{array}\right.$	(39)
		$\displaystyle \Sigma_{\rm g}(a) = 1700 \left(\frac{a}{1~{\rm AU}} \right)^{-3/2}~{\rm g~cm}^{-2}$	(40)
		$\displaystyle T(a) = 280 \left(\frac{a}{1~{\rm AU}}\right)^{-1/2}~{\rm K}$	(41)
		$\displaystyle \rho_{\rm g}(a) = 1.4\times10^{-9} \left(\frac{a}{1~{\rm AU}} \right)^{-11/4}~{\rm g~cm}^{-3}.$	(42)

The discontinuity at 2.7 AU in the surface densities of solids is caused by the condensation of volatiles (the ``snow line''). Hayashi (1981) employed a gas/solid ratio of 240 inside the snow line. Nevertheless, we employ the same model but with a gas/solid ratio of 100 as in Mordasini et al. (2009). In that work, Mordasini et al. (2009) adopted a value of z= 0.0149 (Lodders 2003) for the abundance of heavy elements in the Sun; but this abundance value is higher by a factor of 2 to 4 in the internal regions of the Solar disk because of a redistribution of solids (Kornet et al. 2004). Mordasini et al. (2009) adopted a factor of 3 considering a value of z= 0.04, which corresponds to a gas/sold ratio of 100 inside the snow line. Thus, our standard Solar nebula is defined by

		$\displaystyle \Sigma_{\rm s}(a) = \left\lbrace \begin{array}{ll} 7.1 \left(\dis... ...}}\right)^{-3/2}~{\rm g~cm}^{-2}, & \mbox{$a>2.7$ ~{\rm AU}} \end{array}\right.$	(43)
		$\displaystyle \Sigma_{\rm g}(a) = 710 \left(\frac{a}{1~{\rm AU}} \right)^{-3/2}~{\rm g~cm}^{-2}$	(44)
		$\displaystyle T(a) = 280 \left(\frac{a}{1~{\rm AU}}\right)^{-1/2}~{\rm K}$	(45)
		$\displaystyle \rho_{\rm g}(a)= 5.92 \times10^{-10} \left(\frac{a}{1~{\rm AU}} \right)^{-11/4}~{\rm g~cm}^{-3}.$	(46)

For numerical reasons, we followed Thommes et al. (2003) spreading the snow line in a region of about 1 AU with a smooth function, so the surface densities of solids are described by

$\displaystyle \Sigma_{\rm s}(a)$	=	$\displaystyle \left\lbrace 7.1 + \left( 30-7.1\right) \left[\frac{1}{2}\tanh\left(\frac{a-2.7}{0.5}\right) + \frac{1}{2}\right]\right\rbrace$
		$\displaystyle \times ~\left( \frac{a}{1~{\rm AU}}\right) ^{-3/2}~{\rm g~cm}^{-2}.$	(47)

This model evidently, corresponds to a nebula less massive than that of Hayashi (1981). In our model, the mass contained between 0.4 and 30 AU

$\displaystyle 0.005~M_{\hbox{$\odot$ }}\leq M_{\rm gas} \leq 0.05~M_{\hbox{$\odot$ }},$

(48)

which corresponds to nebulae of between 1 and 10 MMSN (where we assume that the minimum mass solar nebula corresponds to a disk of 0.005 solar masses).

5.1.1 Jupiter and Saturn: isolated formation

After defining our disk model, we computed the isolated formation of Jupiter and Saturn. When we refer to isolated formation, we mean that we consider only one planet forms in the Solar nebula while the disk evolves. The initial conditions and final masses we considered here are given in Table 1.

Table 1: Parameters for the planets considered in the calculations.

The results we obtained are summarized in Table 2. We found that for 4 to 10 MMSN, Jupiter is formed in less than 6 Myr. Saturn is formed in less than 10 Myr for 6 to 10 MMSN. Although theoretical models find small cores for Jupiter (0-12 $M_{\oplus}$ Guillot 2005; 14-18 $M_{\oplus}$ Militzer et al. 2008) and Saturn (9-22 $M_{\oplus}$ Guillot 2005), the models also predict $10{-}40~M_{\oplus}$ of heavy elements in Jupiter's envelope and $20{-}30~M_{\oplus}$ of heavy elements in Saturn's envelope (Guillot 2005; Guillot & Gautier 2009). Because we assumed that all the falling planetesimals reach the core, $M_{\rm c}$ corresponds to the total heavy element mass in the core and envelope.

Table 2: Isolated formation of Jupiter and Saturn for a disk with surfaces densities of solids and gas ${\propto }a^{-3/2}$ .

Table 3: Same as Table 2 but for the simultaneous formation of Jupiter and Saturn.

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg1.eps} \end{figure}$	Figure 1: Comparison of cumulative masses as a function of time for the simultaneous formation of Jupiter and Saturn and isolated formation of Saturn for a 6 MMSN disk with power index p= 3/2. Clearly, when Jupiter reaches its gaseous runaway, Saturn formation is inhibited in comparison with the case of its isolated formation.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg2.eps} \end{figure}$	Figure 2: Same as Fig. 1 but for a 8 MMSN. In this case, Jupiter does not inhibit Saturn formation but largely delays it.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg3.eps} \end{figure}$	Figure 3: The evolution of planetesimal migration mean velocities at Saturn's feeding zone for a 6 MMSN disk with power index p= 3/2. The increase in the migration mean velocities, for the case of Saturn's simultaneous formation, is due to Jupiter's rapid formation as it reaches gaseous runaway.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg4.eps} \end{figure}$	Figure 4: Same as Fig. 3 but for the case of a 8 MMSN nebula.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg5.eps} \end{figure}$	Figure 5: A comparison of the evolution of the solid accretion rate for the simultaneous formation of Jupiter and Saturn and isolated formation of Saturn for a 6 MMSN disk with power index p= 3/2. We see that Jupiter's formation cuts off the solid accretion rate in the case of simultaneous formation with Saturn.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg6.eps} \end{figure}$	Figure 6: Same as Fig. 5 but for a 8 MMSN disk. The formation of Jupiter significantly modifies the solid accretion rate of Saturn, in contrast to the isolated formation, which delays Saturn's.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg7.eps} \end{figure}$	Figure 7: Planetesimals's eccentricity and inclination modulation function due to Jupiter gravitational stirring (Eq. (31)). When Jupiter achieves its final mass, the planetesimals's eccentricity and inclination excitation caused by Jupiter at Saturn's location, is less than 1% of its value at Jupiter's orbit.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{14365fg8.eps} \end{figure}$

Figure 8:

Time evolution of the eccentricities along the disk for the case of the isolated formation of Jupiter. Curves correspond to 0, 1 $\times$ 10^-5, 1 $\times$ 10^-4, 1 $\times$ 10^-3, 0.01, 0.05, 1, 1.73, and 1.76 Myr. The last two curves correspond to times before and after Jupiter gaseous runaway. At these moments, the mass of the planet was of $52~M_{\oplus}$ and $318~M_{\oplus}$ , respectively.

Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{14365fg9.eps} \end{figure}$	Figure 9: Same as Fig. 8, but for the time evolution of the inclinations along the disk. It is clear that inclination excitations are smaller than eccentricity excitations.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{14365fg10.eps} \end{figure}$	Figure 10: Same as Figs. 8 and 9, but for the time evolution of the planetesimal migration velocities along the disk. The planetesimal migration velocities, at 9.5 AU, are increased by Jupiter, only when it achieves its final mass.
Open with DEXTER

5.1.2 Jupiter and Saturn: simultaneous formation

We next assumed that the planets form simultaneously; in this case, the main results we obtained are given in Table 3. For 6 and 7 MMSN, we see that Jupiter formation completely inhibits Saturn formation (simulations are halted when the time of formation exceeds 10 Myr). The Saturn formation timescale is about 10 times longer than the Jupiter formation timescale for the remaining cases.

Comparing the results given in Tables 2 and 3, we see that Saturn has almost no effect on the formation of Jupiter. However, the opposite is not true: the formation of Jupiter, clearly inhibits the formation of Saturn in some cases (see Fig. 1) and largely increases the formation time of Saturn in others (see Fig. 2). This is due to an increment in the migration velocities of planetesimals in Saturn's neighborhood when Jupiter reaches gaseous runaway that, in turn, accelerates the migration of planetesimals to the Saturn's feeding zone (Figs. 3 and 4). This increment in the migration velocity of planetesimals causes the solid accretion timescale to become longer than planetesimal migration timescales, and the solid accretion rate of Saturn (when it is formed simultaneously with Jupiter) becomes less efficient than for the isolated Saturn formation (Figs. 5 and 6). Remarkably, Jupiter's and Saturn's feeding zones do not overlap at any time. Although Saturn's feeding zone is well beyond 4 Hill radii from Jupiter, the increase in the migration velocity of planetesimals at this location, is caused by an eccentricity and inclination excitation of the planetesimals related to Jupiter's perturbations. When Jupiter reaches its final mass ( $318~M_{\oplus}$ ), its corresponding Hill radius is ${\sim}0.35$ AU. At 9.5 AU, Saturn is a distance of ${\sim}12~R_{\rm H}$ from Jupiter. As we can see in Fig. 7, when Jupiter reaches its final mass, the tail of the modulation function (see Eq. (30)) produces an excitation (in planetesimal eccentricity and inclination) lower than 1% at 9.5 AU compared to that produced at Jupiter's location. It is this excitation that causes the increment in the migration velocity of planetesimals at the Saturn's neighborhood when both planets are formed simultaneously. In Figs. 8 and 9, we show the time evolution of eccentricity and inclination along the disk for the case of Jupiter formed in isolation. The increases in eccentricity and inclination when Jupiter reaches gaseous runaway, and reaches its final mass, are responsible for the increment in the migration velocities of planetesimals at 9.5 AU (Fig. 10). We note that the choice of the modulation function is rather arbitrary. However, it should satisfy some key properties: it should produce a planetary excitation on planetesimals restricted to the protoplanetary neighborhood; it must have a tail that is capable of exciting the surrounding planetesimals to be accreted; and, from a numerical point of view, it should guarantee a smooth velocity profile along the entire disk. Moreover, it ought to have an appreciable tail at distances well beyond the assumed function for representing the feeding zone ( $\psi \propto {\rm e}^{- (\frac{\Delta}{4~R_{\rm H}})^2}$ ). If this were not the case, there would be no planetary formation, even for the case of isolated embryos. The choice of another possible modulation function that satisfies the aforementioned conditions, but has a narrower tail would produce even smaller perturbations caused by Jupiter (in eccentricities and inclinations) at Saturn's location. However, if Jupiter and Saturn formed closer to each other, as predicted by the ``Nice Model'' of Tsiganis et al. (2005) (Jupiter $\sim$ 5.5 AU, Saturn $\sim$ 8.3 AU), an inhibition in the formation of Saturn should still occur. To our knowledge, a detailed study of this modulation function is not yet available and its calculation is a difficult task, beyond the scope of this paper. In any case, it seems an unavoidable conclusion that the modulation function should be significantly higher than zero on a long tail, sufficient to allow the excitation of planetesimals by an already formed Jupiter at Saturn's orbit.

When the formation timescales for Jupiter are much shorter than those corresponding to Saturn, the formation of the latter is inhibited. One possible way out of this effect is to consider a smoother profile for the nebular surface density. This point is investigated in the next section.

6 Exploring different disk profiles

Models of steady accretion $\alpha$ -disks predict profiles of ${\Sigma \propto r^{-1}}$ (Shakura & Sunyaev 1973; Pringle 1981). We used this profile and the profile ${\Sigma \propto r^{-1/2}}$ , previously proposed by Lissauer (1987). We repeated the calculations reported in Sect. 5, but using instead these profiles for the disk density.

We decided to normalize these profiles at Jupiter's location. In this way, for the smoother slopes of the nebular surface densities, Saturn will find more material inside its corresponding feeding zone and, in turn, the formation timescales will be shorter. We note that in this case the mass of the disks between 0.4-30 AU never exceeds 0.1 $M_{\odot}$ , so these disks should be stable against self-gravitational collapse (Mayer et al. 2004).

We have shown that for our model of Solar nebula, the surface density is given by Eq. (43). It can be rewritten, scaling to Jupiter's location as

$\displaystyle \Sigma_{\rm s}(a)$

=

$\displaystyle \left\lbrace \begin{array}{ll} 0.6 \left(\displaystyle{\frac{a}{5... ...}\right)^{-3/2}~{\rm g~cm}^{-2}, & \mbox{$a>2.7$ ~{\rm AU}}. \end{array}\right.$

(49)

Thus, we employ a solid profile

$\displaystyle \Sigma_{\rm s}(a)$

=

$\displaystyle \left\lbrace \begin{array}{ll} 0.6 \left(\displaystyle{\frac{a}{5... ...}\right)^{-3/2}~{\rm g~cm}^{-2}, & \mbox{$a>2.7$ ~{\rm AU}}. \end{array}\right.$

(50)

Using the same profile for the temperature and the same gas/solid ratio for the different values of p (-3/2, -1, -1/2) at Jupiter's location (5.2 AU), the value of the surface densities of solids and gas and the volumetric gas density in the mid plane of the disk were found to be the same regardless of the value of p.

A nebula with profile ${\Sigma \propto r^{-1}}$

For example, we considered p=1 in the profile given by Eq. (50)

$\displaystyle \Sigma_{\rm s}(a)$

=

$\displaystyle \left\lbrace \begin{array}{ll} 0.6 \left(\displaystyle{\frac{a}{5... ...}\right)^{-3/2}~{\rm g~cm}^{-2}, & \mbox{$a>2.7$ ~{\rm AU}}. \end{array}\right.$

(51)

Rescaling to 1 AU and using the same gas to solid ratio, we found that the disk becomes defined as

		$\displaystyle \Sigma_{\rm s}(a) = \left\lbrace \begin{array}{ll} 3.13 \left(\di... ...}}\right)^{-1}~{\rm g~cm}^{-2}, \;& \mbox{$a>2.7$ ~{\rm AU}} \end{array}\right.$	(52)
		$\displaystyle \Sigma_{\rm g}(a) = 313 \left(\frac{a}{1~{\rm AU}} \right)^{-1}~{\rm g~cm}^{-2}$	(53)
		$\displaystyle T(a) = 280 \left(\frac{a}{1~{\rm AU}}\right)^{-1/2}~{\rm K}$	(54)
		$\displaystyle \rho_{\rm g}(a) = 2.61 \times10^{-10} \left(\frac{a}{1~{\rm AU}} \right)^{-9/4}~{\rm g~cm}^{-3}.$	(55)

Spreading the snow line in a region of about 1 AU with a smooth function

$\displaystyle \Sigma_{\rm s}(a)$	=	$\displaystyle \left\lbrace 3.13 + \left(13.15-3.13\right) \left[ \frac{1}{2}\tanh\left(\frac{a-2.7}{0.5}\right) + \frac{1}{2}\right]\right\rbrace$
		$\displaystyle \times ~\left( \frac{a}{1~{\rm AU}}\right)^{-1}~{\rm g~cm}^{-2},$	(56)

the disk mass enclosed between 0.4 and 30 AU was found to be

$\displaystyle 0.0065~M_{\hbox{$\odot$ }}\leq M_{\rm gas} \leq 0.065~M_{\hbox{$\odot$ }},$

(57)

which corresponds to values for nebulae between 1 and 10 MMSN (now the MMSM correspond to a disk of 0.0065 solar masses).

After defining the disk, we computed planetary formation.

6.1.1 Jupiter and Saturn: isolated formation

The results we obtained are presented in Table 4. For 3 to 10 MMSN, Jupiter's formation time was found to be short enough not to conflict with observational estimations. However, the final core mass of Jupiter is higher than expected. Similar results were obtained for the isolated formation of Saturn. We found that the timescales for the isolated formation of Jupiter and Saturn (except for the case of 4 MMSN) are in closer agreement (in all cases, Jupiter is formed before Saturn). We next considered the case of simultaneous formation.

6.1.2 Jupiter and Saturn: simultaneous formation

Comparing Tables 4 and 5, we see that, except for the case of 4 MMSN, for which the formation timescales are very different, in the remaining cases the simultaneous formation of Jupiter and Saturn is analogous to the case of isolated growth.

Although in all cases, the final core masses are higher than observational estimations, we note that ablation of accreted planetesimals was not considered here. Small planetesimals could completely disintegrate before they reach the core, reducing the final core masses of Jupiter and Saturn.

Table 4: Isolated formation of Jupiter and Saturn for a disk with power index p= 1.

Table 5: Same as Table 4 but for the simultaneous formation of Jupiter and Saturn.

6.2 A nebula with profile ${\Sigma \propto r^{-1/2}}$

For the profile given by Eq. (50) with power of p=1/2, rescaling at 1 AU, employing the same gas/solid ratio, and spreading the snow line in a region of about 1 AU with a smooth function we found that the disk is characterized by

		$$\displaystyle \Sigma_{\rm s}(a) = \left\lbrace 1.36 + \left(5.77-1.36\right) \l...$
		$\displaystyle \qquad\quad\; \times ~\left( \frac{a}{1~{\rm AU}}\right) ^{-1/2}~{\rm g/cm}^{-2},$	(58)
		$\displaystyle \Sigma_{\rm g}(a) = 136 \left(\frac{a}{1~{\rm AU}} \right)^{-1/2}~{\rm g~cm}^{-2},$	(59)
		$\displaystyle T(a) = 280 \left(\frac{a}{1~{\rm AU}}\right)^{-1/2}~{\rm K},$	(60)
		$\displaystyle \rho_{\rm g}(a) = 1.13 \times10^{-10} \left(\frac{a}{1~{\rm AU}} \right)^{-7/4}~{\rm g~cm}^{-3}.$	(61)

The mass of gas between 0.4 AU and 30 AU is

$\displaystyle 0.01~M_{\hbox{$\odot$ }}\leq M_{\rm gas} \leq 0.1~M_{\hbox{$\odot$ }},$

(62)

which correspond to nebulae of masses between 1 and 10 MMSN (where the MMSM correspond to a disk of 0.01 solar masses).

Again, we first compute isolated formation for Jupiter and Saturn.

6.2.1 Jupiter and Saturn: isolated formation

Table 6: Isolated formation of Jupiter and Saturn for a disk with power index p= 1/2.

Table 7: Same as Table 6 but for the simultaneous formation of Jupiter and Saturn.

For each planet, formation timescales are again, shorter than estimate of disk lifetimes, but the final core masses are higher than observational estimate.

Except for the case of 3 MMSN, Saturn formation occurs remarkably on a timescale shorter than that corresponding to Jupiter. This result may seem rather paradoxical because the solid accretion rate is lower at larger distances from the central star (Eq. (9)), even for this moderate profile. However, it has a simple explanation. For a solid density profile $\Sigma_{\rm s} \propto r^{-p}$ , if we consider an annulus of the disk of thickness ${\rm d}r$ , the mass of the annulus is ${\rm d}m= 2\pi r\Sigma_{\rm s}(r){\rm d}r$ so, ${\rm d}m \propto r^{1-p}$ . If p>1, the mass increases inwards and due to the planetesimal migration, the incoming flux of mass towards the annulus is lower than the outcoming one. On the other hand, if p<1, the mass grows outwards, and the incoming flux of mass is greater than the outcoming. The evolution of the solid density profile is governed by Eq. (7), where sinks, caused by planetary accretion, are also considered. Figure 13 clearly shows that initially the mean density of solids in Jupiter's feeding zone is greater than the mean density of solids in Saturn's feeding zone, but as time advances, the mean density of solids in Saturn's feeding zone become greater than that corresponding to Jupiter. Because of this, Saturn's formation occurs before that of Jupiter's. It is clear that the migration of planetesimals plays a very important role in planetary formation, especially for p <1.

We now consider simultaneous formation.

6.2.2 Jupiter and Saturn: simultaneous formation

We see that for this nebular profile the temporal order of the phenomena is inverted. Jupiter formation does not affect the formation of Saturn (the latter forms first except in the case of 3 MMSN), but Saturn formation shortens the Jupiter formation timescale and largely increases their final cores (Fig. 11).

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg11.eps} \vspace*{6mm} \end{figure}$	Figure 11: Comparison of cumulative masses as a function of time for the simultaneous formation of Jupiter and Saturn and isolated formation of Jupiter for a 6 MMSN disk with power index p= 1/2. The rapid formation of Saturn significantly favors the formation of Jupiter.
Open with DEXTER

The rapid formation of Saturn produces a density wave of solids, responsible for the acceleration of Jupiter's formation (see Fig. 12). This wave increases the mean density of solids in Jupiter's feeding zone (Fig. 13) increasing the accretion rate of solids (Fig. 14). In our model, we do not consider some potentially relevant phenomena such as planetesimals resonant trapping, planetesimals shepherding, and gap opening in the planetesimal disk (Tanaka & Ida 1997, 1999; Zhou & Lin 2007; Shiraishi & Ida 2008). These phenomena could change the process of planetary formation slowing down or even inhibiting solid accretion. Unfortunately, they are very difficult to incorporate in a semi-analytical model such as the one presented here.

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg12.eps} \end{figure}$

Figure 12:

Time evolution of the planetesimal density profiles for a 6 MMSN disk with power index p= 1/2 for the simultaneous formation of Jupiter and Saturn. The density wave of solids, produced by the rapid formation of Saturn is a very dispersive phenomena and strongly depends on planetesimal size. As time advances, density profiles decrease at the planets's location. Curves correspond to 0, 1 $\times$ 10^-7, 1 $\times$ 10^-6, 1 $\times$ 10^-5, 1 $\times$ 10^-4, 1 $\times$ 10^-3, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, and 0.5 Myr.

Open with DEXTER

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg13.eps} \vspace*{6mm} \end{figure}$	Figure 13: Evolution of the mean density of solids in the Jupiter and Saturn feeding zones for a 6 MMSN disk with power index p= 1/2. It is clear that the presence of Saturn increases the mean density of solids in Jupiter's feeding zone, because of an inner density wave of planetesimals produced by the rapid formation of Saturn (see Fig. 12).
Open with DEXTER

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg14.eps} \end{figure}$	Figure 14: A comparison of the evolution of the solid accretion rate for the simultaneous formation of Jupiter and Saturn and isolated formation of Jupiter for a 6 MMSN disk with power index p= 1/2. We see that Saturn's rapid formation significantly increases the solid accretion rate for the simultaneous formation of Jupiter.
Open with DEXTER

7 The case of the Hayashi nebula

Finally, for the sake of completeness, we repeated the calculations for the case of the standard model of Hayashi. Since this model only differs from our own model for p=3/2 in terms of the amount of gas, similar results are expected because a large amount of gas can reduce the relative velocities of the planetesimals, and make the accretion of solids more efficient. The corresponding results are given below.

7.1 Jupiter and Saturn: isolated formation

Comparing the two considered values of the gas/solid ratio (Tables 2 and 8) we find the expected results. The higher the gas/solid ratio, the lower the formation timescale and the higher the core mass. The effects are clearly much larger for the case of Saturn. We now consider whether qualitatively, the main characteristics of simultaneous formation remain the same.

7.1.1 Jupiter and Saturn: simultaneous formation

We see (Table 9) that, except for the case of 5 MMSN, the simultaneous formation changes quantitatively and even qualitatively relative to the results for the case of a ratio gas/solid = 100 (Table 3). For the cases of 6, 7, and 8 MMSN, the formation times of Saturn shorten, while those of Jupiter become longer and the core masses are also lowered considerably. This is because, while an increment in the amount of gas improves the efficiency of the solid accretion rates, it also accelerates the planetesimal migration.

Table 8: Isolated formation of Jupiter and Saturn for the Hayashi nebula.

Table 9: Same as Table 8 but for the simultaneous formation of Jupiter and Saturn.

8 Discussion and conclusions

In the framework of the core instability hypothesis, we have considered the in situ formation of Jupiter and Saturn occurring simultaneously in a protoplanetary disk populated by planetesimals with a size distribution for which most of the mass is in small objects. We considered protoplanetary nebulae with power-law surface densities for gas and solids, and assumed that the gas component of the disk dissipates following an exponential law with a characteristic timescale of 6 Myr.

In the first instance, we calculated the isolated and simultaneous formation of Jupiter and Saturn for a standard nebula with profile ${\Sigma \propto r^{-3/2}}$ but with a lower gas/solid ratio than that adopted by Hayashi (1981), following the work by Mordasini et al. (2009). We quantitatively analyzed how the isolated formation of Jupiter (Saturn) is affected when it occurs simultaneously with Saturn (Jupiter). For the isolated formation, we have found that Jupiter and Saturn achieved its final masses on a timescale in agreement with observational estimate for 6 to 10 MMSN (the corresponding disk masses are between 0.003 $M_{\odot}$ and 0.05 $M_{\odot}$ ). However, the final core masses are higher than estimates in most of the cases (Table 2). The most important result is that Jupiter's formation timescales are somewhat shorter than the corresponding ones for Saturn. We have found that the rapid formation of Jupiter inhibits - or largely increases- the timescale of Saturn's formation when they grow simultaneously inside the disk (Table 3). Jupiter's formation increases the migration mean velocities of the planetesimals at Saturn's feeding zone (Figs. 3 and 4). This phenomena causes the timescale of solids accretion rate to becomes longer than the planetesimal migration timescale such that, Saturn's solid accretion rate becomes less efficient (Figs. 5 and 6). In Fig. 15, we show that the presence of Jupiter significantly modifies the time evolution of solids in Saturn's neighborhood and, in turn, modifies Saturn formation.

$\begin{figure} \par\includegraphics[width=6.5cm,clip]{14365fg15.eps} \end{figure}$

Figure 15:

Comparison of the time evolution of planetesimal density profiles for a 8 MMSN disk with power index p= 3/2 in Saturn's neighborhood for its isolated formation, and also its simultaneous formation with Jupiter. As time advances, the density profiles decrease in both cases; however, profiles evolution are very different. Planetesimals of 1, 10, 100 km in size produce a density wave for Saturn's isolated formation. This density wave is inhibited when Saturn forms simultaneously with Jupiter: for all sizes of planetesimals, the evolution of density profiles decreases more quickly when Saturn forms in the presence of Jupiter.

Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{14365fg16.eps} \end{figure}$	Figure 16: Comparison of the simultaneous formation of Jupiter and Saturn for a 6 MMSN disk with different power indices ( p= 3/2;1;1/2). We see that simultaneous formation of Jupiter and Saturn changes as the power index decreases.
Open with DEXTER

We then explored smoother profiles in trying to solve the problem of the different formation timescales. Accretion $\alpha$ disks predict profiles of ${\Sigma \propto r^{-1}}$ . We used this profile and another proposed by Lissauer (1987) - ${\Sigma \propto r^{-1/2}}$ - that could solve formation timescales of the outer Solar System giant planets. We normalized these profiles at Jupiter's position (5.2 AU) so that, at this position, the surface densities of solids and gas, and the volumetric density of gas in the mid plane disk are the same for all considered profiles ( p=3/2,1,1/2). The smoother the profiles, the higher the density of solids at Saturn's position, and the closer the similarity in the formation timescales.

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{14365fg17.eps} \vspace*{-2.5mm} \end{figure}$	Figure 17: Comparison of the simultaneous formation of Jupiter and Saturn for a 6 MMSN disk for the Hayashi nebula ( $\Sigma _{\rm g}/\Sigma _{\rm s}= 240$ ) and our model ( $\Sigma _{\rm g}/\Sigma _{\rm s}= 100$ ). The formation processes are very different in each case.
Open with DEXTER

For the profile ${\Sigma \propto r^{-1}}$ , we found that the formation timescales of Jupiter and Saturn are more similar to those for isolated formation. In this case, the simultaneous formation of Jupiter and Saturn is practically does not unchanged with respect to the corresponding isolated formation of each planet. However, for the profile ${\Sigma \propto r^{-1/2}}$ the temporal order of the phenomena is inverted. In most cases, Saturn forms before Jupiter. The rapid formation of Saturn induces a density wave of solids (Fig. 12) that significantly favors Jupiter's formation (Fig. 11). This wave increases the amount of solids in Jupiter's feeding zone (Fig. 13), such that the solid accretion rate becomes more efficient (Fig. 14).

In Fig. 16, we compare the simultaneous formations of Jupiter and Saturn for the case of a 6 MMSN disk with different power indices ( p= 3/2;1;1/2). For the case of p=3/2, Jupiter inhibits Saturn's formation. For smoother profiles (p= 1;1/2), the process becomes more efficient, for various reasons. First, at Saturn's position there is a higher density of solids and gas. In addition, as power indices decrease the mass of solids grows outward, and due to planetesimal migration the replenishment of solids in the feedings zones is more efficient.

Finally, we repeated the calculations for the standard model of Hayashi. Because Hayashi's nebula is more massive than that of our model with p= 3/2 - only differing in the amount of gas - our results were expected to be qualitatively similar. This was found to be true for the isolated formation of Jupiter and Saturn (Table 8), where formation timescales decrease and final core masses increase. This is because a greater amount of gas leads to the accretion rate of solids becoming more efficient. However, the properties of the simultaneous formation of Jupiter and Saturn differ qualitatively and quantitatively in the case of our model with p= 3/2. While for our model Jupiter inhibits Saturn's formation, for the Hayashi's nebula the formation times of Saturn are decreased by the presence of Jupiter. In addition, Jupiter final core masses are significantly decreased by the presence of Saturn.

In Fig. 17, we show the differences between the simultaneous formations of Jupiter and Saturn for a 6 MMSN disk between the Hayashi nebula and our model with p= 3/2. The gas/solid ratio clearly has an important role in the formation process. Different gas/solid ratios may change qualitatively and quantitatively the whole formation of a planetary system.

The primary hypothesis of our model is that most of the planetesimal mass is in small objects. Most of the works that simulate giant planet formation do not consider oligarchic growth regime. In these works, a much faster time-dependent regime is adopted for the growth of the core (Pollack et al. 1996; Alibert et al. 2005a,b; Hubickyj et al. 2005; Dodson-Robinson et al. 2008). With this solid accretion rate, the final core mass of the giant planet is quickly achieved. Relative velocities are slower and the solid accretion rate becomes more efficient, so planetesimals of 100 km in size are usually used. While some studies predict this type of planetesimals (Johansen et al. 2007; Morbidelli et al. 2009), N-body simulations predict that when the embryos reach a mass similar to the Moon, the runaway regime switches to oligarchic regime (Ida & Makino 1993; Kokubo & Ida 1998, 2000, 2002). We note that if our primary hypothesis is relaxed, and most of the planetesimal mass does not reside in small objects - or larger planetesimals are considered - formation timescales become longer than disk lifetimes.

Acknowledgements

We are especially grateful to A. Fortier for the close collaboration and the lively discussions.

References

Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Prog. Theor. Phys., 56, 1756 [NASA ADS] [CrossRef] [Google Scholar]
Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005a, A&A, 434, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Alibert, Y., Mousis, O., Mordasini, C., & Benz, W. 2005b, ApJ, 626, L57 [NASA ADS] [CrossRef] [Google Scholar]
Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
Benvenuto, O. G., Fortier, A., & Brunini, A. 2009, Icarus, 204, 752 [NASA ADS] [CrossRef] [Google Scholar]
Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391 [NASA ADS] [CrossRef] [Google Scholar]
Brunini, A., & Benvenuto, O. G. 2008, Icarus, 194, 800 [NASA ADS] [CrossRef] [Google Scholar]
Chambers, J. 2006, Icarus, 180, 496 [NASA ADS] [CrossRef] [Google Scholar]
Dodson-Robinson, S. E., Bodenheimer, P., Laughlin, G., et al. 2008, ApJ, 688, L99 [NASA ADS] [CrossRef] [Google Scholar]
Fortier, A., Benvenuto, O. G., & Brunini, A. 2007, A&A, 473, 311 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fortier, A., Benvenuto, O. G., & Brunini, A. 2009, A&A, 500, 1249 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
Guillot, T. 2005, Annu. Rev. Earth Planet. Sci., 33, 493 [Google Scholar]
Guillot, T., & Gautier, D. 2009, Treatise of Geophysics, Planets and Moons, ed. G. Schubert, & T. Spohn (2007), 10, 439 [arXiv:0912.2019] [Google Scholar]
Haisch, K. E., Jr., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
Hasegawa, M., & Nakazawa, K. 1990, A&A, 227, 619 [NASA ADS] [Google Scholar]
Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Hillenbrand, L. A. 2005, in A Decade of Discovery: Planets Around Other Stars STScI Symp. Ser., 19, ed. M. Livio, in press [arXiv:astro-ph/0511083] [Google Scholar]
Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415 [NASA ADS] [CrossRef] [Google Scholar]
Ida, S., & Makino, J. 1993, Icarus, 106, 210 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Ida, S., Guillot, T., & Morbidelli, A. 2008, ApJ, 686, 1292 [NASA ADS] [CrossRef] [Google Scholar]
Inaba, S., & Ikoma, M. 2003, A&A, 410, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Inaba, S., Tanaka, H., Nakazawa, K., Wetherill, G. W., & Kokubo, E. 2001, Icarus, 149, 235 [NASA ADS] [CrossRef] [Google Scholar]
Johansen, A., Oishi, J. S., Low, M.-M. M., et al. 2007, Nature, 448, 1022 [Google Scholar]
Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [NASA ADS] [CrossRef] [Google Scholar]
Kokubo, E., & Ida, S. 2000, Icarus, 143, 15 [NASA ADS] [CrossRef] [Google Scholar]
Kokubo, E., & Ida, S. 2002, ApJ, 581, 666 [NASA ADS] [CrossRef] [Google Scholar]
Kornet, K., Rózyczka, M., & Stepinski, T. F. 2004, A&A, 417, 151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lissauer, J. J. 1987, Icarus, 69, 249 [NASA ADS] [CrossRef] [Google Scholar]
Lodders, K. 2003, ApJ, 591, 1220 [NASA ADS] [CrossRef] [Google Scholar]
Matsumura, S., Pudritz, R. E., & Thommes, E. W. 2009, ApJ, 691, 1764 [NASA ADS] [CrossRef] [Google Scholar]
Mayer, L., Quinn, T., Wadsley, J., & Stadel, J. 2004, ApJ, 609, 1045 [NASA ADS] [CrossRef] [Google Scholar]
Militzer, B., Hubbard, W. B., Vorberger, J., Tamblyn, I., & Bonev, S. A. 2008, ApJ, 688, L45 [NASA ADS] [CrossRef] [Google Scholar]
Mizuno, H. 1980, Prog. Theor. Phys., 64, 544 [Google Scholar]
Morbidelli, A., Bottke, W. F., Nesvorný, D., & Levison, H. F. 2009, Icarus, 204, 558 [NASA ADS] [CrossRef] [Google Scholar]
Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ohtsuki, K., Stewart, G. R., & Ida, S. 2002, Icarus, 155, 436 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Ormel, C. W., Dullemond, C. P., & Spaans, M. 2010, ApJ, 714, L103 [NASA ADS] [CrossRef] [Google Scholar]
Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 [NASA ADS] [CrossRef] [Google Scholar]
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
Shiraishi, M., & Ida, S. 2008, ApJ, 684, 1416 [NASA ADS] [CrossRef] [Google Scholar]
Tanaka, H., & Ida, S. 1997, Icarus, 125, 302 [NASA ADS] [CrossRef] [Google Scholar]
Tanaka, H., & Ida, S. 1999, Icarus, 139, 350 [NASA ADS] [CrossRef] [Google Scholar]
Teiser, J., & Wurm, G. 2009, MNRAS, 393, 1584 [NASA ADS] [CrossRef] [Google Scholar]
Thommes, E. W., Duncan, M. J., & Levison, H. F. 2003, Icarus, 161, 431 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Wetherill, G. W., & Stewart, G. R. 1993, Icarus, 106, 190 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Zhou, J.-L., & Lin, D. N. C. 2007, ApJ, 666, 447 [NASA ADS] [CrossRef] [Google Scholar]

Footnotes

...: Fellow of the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Argentina
...: Member of the Carrera del Investigador Científico, CONICET, Argentina.
...: Member of the Carrera del Investigador Científico, Comisión de Investigaciones Científicas de la Provincia de Buenos Aires, Argentina.

All Tables