© ESO, 2012

1. Introduction

Since the discovery of the first extrasolar planet around a solar-type star (Mayor & Queloz 1995) more than 800 extrasolar planets have been identified. Observations indicate that planets are abundant in the universe. Planets orbiting stars show a great variety of semi-major axis (from less than 0.01 AU to more than hundreds of AU) and masses (from less than an Earth mass to several Jupiter masses), as can be found in The Extrasolar Planets Encyclopaedia¹. Planet formation models should be able to explain this observed diversity. The sequential accretion model, also called core-nucleated accretion model, is currently the most accepted scenario for planetary formation (e.g. Mizuno 1980; Pollack et al. 1996; Alibert et al. 2005a,b, among others), because it can account naturally for the formation of planets in all mass ranges². It proposes that planetary growth occurs mainly in two stages. In the first stage, the formation of planets is dominated by the accretion of solids. If the protoplanet is able to grow massive enough (~10 M_⊕) while the gas component of the protoplanetary disc is still present, it can gravitationally bind some of the surrounding gas, giving birth to a gas giant planet. The accretion of gas is slow at the beginning: the planet growth is dominated by the accretion of solids and the energy released from the accreted planetesimals slows down the accretion of the envelope. When the accretion of planetesimals declines (generally because the protoplanet has emptied its feeding zone), the accretion of gas is triggered, and the planet can accrete hundreds of Earth masses in a very short time. The runaway accretion of gas characterises the second stage of the sequential accretion model, where the growth of the planet is dominated by the accretion of gas.

Since it was first proposed (Mizuno 1980), the sequential accretion model has been extensively studied and improved, trying to include the many fundamental processes that occur simultaneously with the growth of the planet, and that impact directly on it. Constructing a complete model that accounts for all these processes in a reasonable way is a hard task. Among its main ingredients, it has to include a realistic model for an evolving protoplanetary disc, a model for the accretion of solids and gas to form the planets (which itself requires knowledge of the internal structure of the planet), a model for the interactions between the planets and the disc, and a model for the interactions between the forming planets. Each of these topics itself represents an independent, ongoing area of research. Alibert et al. (2005, A05 from now on) included some of these processes in a single-planet formation model. Given the complexity of the problem or the unknowns related to some of these processes, many simplifications have to be assumed to keep the problem tractable from the physical and computational point of view. This is particularly important because we aim to compute thousands of simulations to account for the wide range of possible initial conditions (see Mordasini et al. 2009, M09 from now on, for more details). Therefore, our models represent a compromise between accuracy and simplicity in their physical description.

The first stage of planetary formation corresponds to the growth of the solid embryo, which is dominated by the accretion of planetesimals. The growth of an embryo proceeds in two different regimes (Ida & Makino 1993; Ormel et al. 2010). At the beginning, big planetesimals, which have larger cross-sections, are favoured to grow even bigger by accreting planetesimals that they encounter on their way. Being more massive, in turn, enlarges the gravitational focusing, which leads to accretion in a runaway fashion. However, at some point, these runaway embryos become massive enough to stir the planetesimals around them. This results in an increase of the relative velocities and the corresponding reduction of the gravitational focusing. Growth among small planetesimals is stalled and only big embryos have the possibility to continue accreting, although at a slower pace. This second regime in the growth of solid embryos is known as oligarchic growth, as only the larger planetesimals or embryos (the oligarchs) are able to keep on growing. One important aspect is that the transition between runaway and oligarchic growth occurs for very small embryos. As shown in Ormel et al. (2010), the actual mass for this transition depends upon many factors: the size of the accreted planetesimals, the location in the disc, and the surface density of solids, among others. In most of the cases, an embryo of  ~0.01 M_⊕, or even smaller, is already growing in the oligarchic regime.

Our model builds up on the model of A05 and M09. Our primary aim is to study the formation of planets of different sizes, and in particular the cases for which the accretion of gas is important. Therefore, in the computations presented here, we focus on the first phase of planetary formation, when the gas component of the disc is still present (in the case of small rocky planets, collisions between embryos after the dissipation of the disc should be included to calculate their final masses). For the formation of giant planets, the growth of the solid core is dominated by the oligarchic growth. One of the weak points of the majority of previous giant planet formation models (e.g. Pollack et al. 1996; Hubickyj et al. 2005; A05; M09; Lissauer 2009; Mordasini et al. 2012a, hereafter M12) is the description of the solid disc, in particular of the interactions between forming planets and planetesimals. These simplified models lead to an overestimation of the solid accretion rate which, in turn, results in an underestimation of the formation time of the whole planet.

Indeed, in those works, the model for the accretion rate of solids is oversimplified: the whole formation of giant planet cores is assumed to proceed very fast, underestimating the excitation that planetesimals suffer due to the presence of the embryos. When oligarchic growth is adopted as the dominant growth model, giant planet formation turns out to be more difficult. Formation times become much longer than the typical lifetime of the protoplanetary disc. Fortier et al. (2007, 2009) and Benvenuto et al. (2009) studied the formation of giant planets adopting the oligarchic growth for the core. Assuming in situ formation for the planets and a simple, unevolving protoplanetary disc, these authors showed that the formation of giant planets is unlikely if the planetesimals that populate the disc are big (more than a few kilometres in size). However, formation could be accelerated if most of the accreted mass is bound in small planetesimals (less than 0.1 km). Guilera et al. (2010, 2011), also considering in situ models, studied the simultaneous formation of several planets where planetesimal drifting is included. These authors considered different density profiles for the disc and found that only for massive discs the formation of the giant planets as in the solar system is possible only if planetesimal radii are smaller than 1 km. These models, however, do not take into account that planets would likely migrate during their formation.

We include in our planet formation model a more realistic description of the accretion of solids. In Sect. 2, we review the basics of the A05 formation model, presenting some improvements in the computation of the disc structure, the internal structure, and migration. In Sect. 3 we describe the new treatment of the planetesimal accretion. In Sect. 4 we present the results obtained for the formation of isolated planets (the formation of planetary systems is described in Alibert et al. 2012; and Carron et al. 2012). In Sect. 5 we discuss our results and put them in context. Finally, in Sect. 6 we summarise our results and underline the main conclusions.

2. Formation model

The model and the numerical code used to calculate the formation of planets is in essence the same as in A05. In what follows, we summarise the most relevant aspects of the model and the improvements that have been introduced since that work. In the next section, we focus on the accretion rate of solids and describe in detail the adopted model for the protoplanet-planetesimal interactions.

2.1. Protoplanetary disc: gas phase

The structure and evolution of the protoplanetary disc is computed by first determining the vertical structure of the disc for each distance to the central star, and second, computing the radial evolution due to viscosity, photoevaporation, and mass accretion by forming planets.

2.1.1. Vertical structure

The vertical disc structure is computed by solving the following equations: $\begin{array}{ccc}& & \\ & & \end{array}$\arraycolsep1.75ptand $\mathit{F}\mathrm{=}\frac{\mathrm{-}\mathrm{16}\mathit{\pi }{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}^{\mathrm{3}}}{\mathrm{3}\mathit{\kappa }{\mathit{\rho }}_{\mathrm{gas}}}\frac{\mathit{\partial T}}{\mathit{\partial z}}\mathrm{·}$(3)They reflect the hydrostatic equilibrium, the energy conservation, and the diffusion for the radiative flux. In these equations, z is the vertical coordinate, ρ_gas the gas density, P the pressure, T the temperature, ν the macroscopic viscosity, F the radiative flux, κ is the opacity (Bell & Lin 1994), and σ_SB is the Stefan-Boltzmann constant. The Keplerian frequency, Ω, is given by ${\mathrm{\Omega }}^{\mathrm{2}}\mathrm{=}\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathit{/}{\mathit{a}}^{\mathrm{3}}\mathit{,}$(4)Gbeing the gravitational constant, M _⋆  the mass of the central star and a the distance to the star³.

Equations (1) − (3) are solved with four boundary conditions. The first three are the temperature, the pressure, and the energy flux at the surface. The surface of the disc is defined as the place where the vertical optical depth (between the surface and infinity) is equal to 0.01. The fourth boundary condition is that the energy flux equals 0 in the midplane (see A05 for details). The three differential equations, together with the four boundary conditions, have a solution only for one value of the disc thickness H, which gives the location of the disc surface. The macroscopic viscosity ν is calculated using the standard Shakura & Sunyaev (1973) α-parametrization, $\mathit{\nu }\mathrm{=}\mathit{\alpha }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{\Omega }$. The speed of sound c_s is determined from the equation of state (Saumon et al. 1995). The temperature at the surface T_surf is computed as in A05. In the models presented in this paper, α is set to 7 × 10^-3. This value of the alpha parameter has to be taken as an example.

In the calculations of this work we neglected irradiation and the possible presence of a dead zone. These effects will be included in future works.

2.1.2. Evolution

The evolution of the gas disc surface density ($\mathrm{\Sigma }\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{-}\mathit{H}}^{\mathit{H}}{\mathit{\rho }}_{\mathrm{gas}}\mathrm{d}\mathit{z}$) is computed by solving the diffusion equation: $\frac{\mathrm{d\Sigma }}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{3}}{\mathit{a}}\frac{\mathit{\partial }}{\mathit{\partial a}}\left[{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{\mathit{\partial }}{\mathit{\partial a}}\begin{array}{}\mathrm{˜}\\ \mathit{\nu }\end{array}\mathrm{\Sigma }{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{2}}\right]\mathrm{+}\mathrm{\Sigma ̇}\mathit{w}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{+}\mathit{Q̇}\mathrm{planet}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{,}$(5)where $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\nu }\end{array}$ is the effective viscosity, $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\nu }\end{array}\mathrm{\equiv }\frac{\mathrm{1}}{\mathrm{\Sigma }}{{}^{\mathrm{\int }}}_{\mathrm{-}\mathit{H}}^{\mathit{H}}\mathit{\nu }{\mathit{\rho }}_{\mathrm{gas}}\mathrm{d}\mathit{z}$. Photoevaporation is included using the model of Veras & Armitage (2004), $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(6)where R_g = 5 AU. The total mass loss due to photoevaporation is a free parameter. The sink term is equal to the gas mass accreted by the forming planets. For every forming planet, mass is removed from the protoplanetary disc in an annulus centred on the planet, with a width equal to the planet’s Hill radius $\begin{array}{ccc}{\mathit{R}}_{\mathrm{H}}\mathrm{=}{\mathit{a}}_{\mathit{M}}{\left(\frac{\mathit{M}}{\mathrm{3} {\mathit{M}}_{\mathit{\star }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}& & \end{array}$(7)where M is the total mass of the planet and a_M is the location of the planet.

Equation (5) is solved on a grid that extends from the innermost radius of the disc to 1000 AU. At these two points, the surface density is constantly equal to 0. The innermost radius of the disc is of the order of 0.1 AU.

Figure 1 presents a typical evolution of a disc, whose parameters correspond to the first row of Table 1, where the curves are plotted every 10⁵ years. In this model, the photoevaporation term is adjusted to obtain a disc lifetime equal to 3 Myr.

Fig. 1 Disc model. Gas surface density as a function of radius at different times. For the initial model (red line) the parameters are those of the first row of Table 1, the disc lifetime being 3 Myr. Each line corresponds to a different time (from top to bottom 0 to 3 Myr, every 0.1 Myr).

The characteristics of the protoplanetary disc are chosen to match the observations as closely as possible. The initial disc density profiles we consider are given by, $\mathrm{\Sigma }\mathrm{=}\mathrm{(}\mathrm{2}\mathrm{-}\mathit{\gamma }\mathrm{)}\frac{{\mathit{M}}_{\mathrm{disc}}}{\mathrm{2}\mathit{\pi }{\mathit{a}}_{\mathrm{C}}^{\mathrm{2}\mathrm{-}\mathit{\gamma }}{\mathit{a}}_{\mathrm{0}}^{\mathit{\gamma }}}{\left(\frac{\mathit{a}}{{\mathit{a}}_{\mathrm{0}}}\right)}^{\mathrm{-}\mathit{\gamma }}\exp \left[\mathrm{-}{\left(\frac{\mathit{a}}{{\mathit{a}}_{\mathrm{C}}}\right)}^{\mathrm{2}\mathrm{-}\mathit{\gamma }}\right]\mathit{,}$(8)where a₀ is equal to 5.2 AU. The mass of the disc (M_disc), the characteristic scaling radius (a_C) and the power index (γ) are derived from the observations of Andrews et al. (2010). Adopting this kind of initial density profile is a difference from previous works (A05 and M09). For numerical reasons, the innermost disc radius, a_inner, is always greater than or equal to 0.1 AU, and differs in some cases from the one cited in Andrews et al. (2010). The parameters used to generate the initial disc’s profile are listed in table 1. Andrews et al. (2010) also derived a value for the viscosity parameter α. In contrast, and for simplicity, we assume here that the viscosity parameter is the same for all protoplanetary discs we consider (α = 7 × 10^-3). Using a different α parameter will be the subject of future work. In the observations of Andrews et al. (2010), the mass of the central star ranges from 0.3 to 1.3 M_⊙. However, we assume here that these disc profiles are all suitable for protoplanetary discs around solar mass stars. Future disc observations will help to improve this part of our models.

As in A05, the planetesimal-to-gas ratio is assumed to scale with the metallicity of the central star. For every protoplanetary disc we consider, we select the metallicity of a star at random from a list of  ~1000 CORALIE targets (Santos, priv. comm.).

Finally, following Mamajek (2009), we assume that the cumulative distribution of disc lifetimes decays exponentially with a characteristic time of 2.5 Myr. When a lifetime T_disc is selected, we adjust the photoevaporation rate so that the protoplanetary disc mass reaches 10^-5   M_⊙ at the time t = T_disc, when we stop the calculation.

2.2. Protoplanetary disc: solid phase

We consider that the planetesimal disc is composed of rocky and icy planetesimals. Here we assume a mean density of 3.2 g/cm³ for the rocky and 1 g/cm³ for the icy planetesimals. The rocky planetesimals are located between the innermost point of the disc (given by the fourth column of Table 1), and the initial location of the ice line, whereas the disc of icy planetesimals extends from the ice line to the outermost point in the simulation disc.

The location of the ice line is computed from the initial gas disc model, using the central temperature and pressure. The ice sublimation temperature we use depends upon the pressure. In our model the location of the ice line does not evolve with time. In particular, no condensation of moist gas or sublimation of icy planetesimals is taken into account. Moreover, the location of the ice line is based on the central pressure and temperature, meaning that the ice line is taken to be independent of the height in the disc. In reality the ice line is likely to be an ice surface whose location depends upon the height inside the disc (see Min et al. 2011).

For the models presented here, we assume that all planetesimals have the same radius. The planetesimal mass is calculated assuming that they are spherical and have constant density (which depends on their location in the disc), and it does not evolve with time. The extension of our calculations towards non-uniform and time-evolving planetesimal mass function is something we are working on and will be included in the next paper.

The surface density of planetesimals, Σ_m, is assumed to be proportional to the initial surface density of gas Σ₀. This means that $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathit{m}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{D}\mathit{/}\mathrm{G}}{\mathit{f}}_{\mathrm{R}\mathit{/}\mathrm{I}}\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{,}& & \end{array}$(9)where f_D/G is the dust-to-gas ratio of the disc, and it scales with the metallicity of the central star (the lowest value being 0.003 and the highest value 0.125), and f_R/I takes into account the degree of condensation of heavy elements. As in M09, we consider f_R/I = 1/4 inside the ice line, and f_R/I = 1 beyond it, for rocky and icy planetesimals.

The surface density of planetesimals evolves as a result of accretion and ejection by the forming planets. The same procedure as in A05 is adopted. Planetesimals that can be accreted by a growing planet are those within the planet’s feeding zone, here assumed to be an annulus of 5   R_H at each side of the orbit. The solids surface density inside the feeding zone is considered to be constant, i.e. the protoplanet instantaneously homogenises it as a result of the scattering it produces among planetesimals. Ejected planetesimals are considered to be lost.

We stress that this work is the first one of a series of papers. Here, effects such as planetesimal drifting due to gas drag, fragmentation, and planetesimal size distribution are neglected, but will be included in future works.

2.3. Gas accretion: attached phase

Equations

Planetary growth proceeds through solids and gas accretion. Gas accretion is the result of a planet’s contraction and is computed by solving the standard internal structure equations $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{r}}^{\mathrm{3}}}{\mathrm{d}{\mathit{M}}_{\mathit{r}}}& \mathrm{=}& \\ \frac{\mathrm{d}\mathit{P}}{\mathrm{d}{\mathit{M}}_{\mathit{r}}}& \mathrm{=}& \\ \frac{\mathrm{d}\mathit{T}}{\mathrm{d}\mathit{P}}& \mathrm{=}& \end{array}$where r,P,T are the radius, the pressure and the temperature inside the envelope. These three quantities depend upon the gas mass, M_r, included in a sphere of radius r. The temperature gradient is given by the adiabatic (∇_ad) or by the radiative gradient (∇_rad), depending upon the stability of the zone against convection, which we check using the Schwarzschild criterion. These equations are solved using the equation of state (EOS) of Saumon et al. (1995). The opacity, which enters in the radiative gradient ∇_rad, is computed according to Bell & Lin (1994). In this work we assume that the grain opacity is the full interstellar opacity. However, Podolak (2003) and Movshovitz & Podolak (2008) showed that the grain opacity in the envelope of forming planets should be much lower than the interstellar one. Reducing the grain opacity accelerates the formation of giant planets (Pollack at al. 1996; Hubickyj et al. 2005) because it allows the runaway of gas to start at smaller core masses. Since the objective of the present paper is to explore the consequences of the planet-planetesimal interactions, we only consider a single opacity, corresponding to the full interstellar opacity.

We omitted the energy equation that yields the luminosity, which in turn is part of the radiative gradient (in the parts of the planet that are stable against convection), and of the determination of the planet stability to convective motions. Including the energy equation in its standard form $\frac{\mathrm{d}\mathit{L}}{\mathrm{d}{\mathit{M}}_{\mathit{r}}}\mathrm{=}{\mathit{ϵ}}_{\mathrm{m}}\mathrm{-}\mathit{T}\frac{\mathrm{d}\mathit{S}}{\mathrm{d}\mathit{t}}$ (the first term results from the accretion of planetesimals, the second from the contraction of the planet) usually entails numerical difficulties. Here we follow Mordasini et al. (2012b) to calculate the luminosity in an easier way. Note, however, that we improved the approach of these authors to take the energy of the core analytically into account.

The total luminosity is given by L = L_cont + L_m, L_cont being the contraction luminosity and L_m the accretion luminosity. We assume L_cont to be constant in the whole planetary envelope. L_cont is computed as the result of the change in total energy of the planet between two time steps t and t + dt, ${\mathit{L}}_{\mathrm{cont}}\mathrm{=}\mathrm{-}\frac{{\mathit{E}}_{\mathrm{tot}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{d}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{E}}_{\mathrm{tot}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{\mathit{E}}_{\mathrm{gas}\mathit{,}\mathrm{acc}}}{\mathrm{d}\mathit{t}}\mathit{,}$(13)where E_tot is the total planetary energy and ${\mathit{E}}_{\mathrm{gas}\mathit{,}\mathrm{acc}}\mathrm{=}\mathrm{d}\mathit{t}\mathit{Ṁ}{}_{\mathrm{gas}}\mathit{u}_{\mathrm{int}}$ is the energy advected during gas accretion (u_int being the internal specific energy). E_gas,acc is negligible compared to the other terms. The luminosity due to accretion of planetesimals is ${\mathit{L}}_{\mathrm{m}}\mathrm{=}\mathit{G}\frac{\mathit{Ṁ}\mathrm{core}{\mathit{M}}_{\mathrm{core}}}{{\mathit{R}}_{\mathrm{core}}}\mathrm{·}$(14)However, the energy at the time t + dt is not known before computing the internal structure at this given time. To circumvent this problem, we use the following approach: the energy is split into two parts, one related to the core, one to the envelope.

The core energy is given by ${\mathit{E}}_{\mathrm{core}}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{3}\mathit{/}\mathrm{5}\mathrm{)}\hspace{0.17em}\mathit{G}{\mathit{M}}_{\mathrm{core}}^{\mathrm{2}}\mathit{/}{\mathit{R}}_{\mathrm{core}}$, the core density is assumed to be uniform. The envelope energy is assumed to follow a similar functional form: E_env =  −k_envM_envg, where g is a mean gravity, taken to be G(M_core/R_core + M_tot/R_planet). This last formula defines k_env, in which all our ignorance of the internal structure is hidden. To calculate the envelope energy at time t + dt, the value of k_env is first taken to be the value resulting from the structure at time t. Then, iteration on k_env is performed until convergence is reached. In general, a first order correction is sufficient to reach a satisfactory solution.

Boundary conditions

The internal structure equations are solved for M_r varying between the core mass M_core and the total planetary mass. Four boundary conditions are given, namely the core radius, the total planetary radius, and T_surf,planet and P_surf,planet, the temperature and pressure at this point. Given the boundary conditions, the differential equations have only one solution for a given total planetary mass.

The core radius is given as a function of the core luminosity and the pressure at the core surface by the following formula, which constitute a fit to the results of Valencia et al. (2010), $\begin{array}{ccc}\frac{{\mathit{R}}_{\mathrm{core}}}{\mathrm{9800} \mathrm{km}}& \mathrm{=}& {\left(\frac{{\mathit{M}}_{\mathrm{core}}}{\mathrm{5} {\mathit{M}}_{\mathrm{\oplus }}}\right)}^{\mathrm{0.28}\mathrm{+}\mathrm{0.02}\sqrt{\frac{{\mathit{M}}_{\mathrm{core}}}{\mathrm{5}{\mathit{M}}_{\mathrm{\oplus }}}}}\\ & & \mathrm{\times }{\mathrm{10}}^{\mathrm{-}{\left[{\log }_{\mathrm{10}}\left(\mathrm{1}\mathrm{+}{\mathit{P}}_{\mathrm{core}}\mathit{/}\sqrt{\frac{{\mathit{M}}_{\mathrm{core}}}{\mathrm{5} {\mathit{M}}_{\mathrm{\oplus }}}}\mathrm{)}\right)\mathit{/}\mathrm{7}\right]}^{\mathrm{3}}}\mathit{,}\end{array}$(15)where P_core is the pressure at the core-envelope interface, expressed in GPa.

The radius of the planet, R_M, is given by (Lissauer et al. 2009), ${\mathit{R}}_{\mathit{M}}\mathrm{=}\frac{\mathit{GM}}{\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{k}}_{\mathrm{1}}}\mathrm{+}\frac{\mathit{GM}}{{\mathit{k}}_{\mathrm{2}}{\mathit{R}}_{\mathrm{H}}}}\mathit{,}$(16)where ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ is the square of the sound velocity in the disc midplane at the planet’s location, k₁ = 1 and k₂ = 1/4. At the planet’s surface, the temperature and pressure are given by, ${\mathit{T}}_{\mathrm{surf}\mathit{,}\mathrm{planet}}\mathrm{=}{\left({\mathit{T}}_{\mathrm{disc}}^{\mathrm{4}}\mathrm{+}\frac{\mathrm{3}\mathit{\tau L}}{\mathrm{16}\mathit{\pi }{\mathit{\sigma }}_{\mathit{P}}{\mathit{R}}_{\mathit{M}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}$(17)and ${\mathit{P}}_{\mathrm{surf}\mathit{,}\mathrm{planet}}\mathrm{=}{\mathit{P}}_{\mathrm{disc}}\mathit{,}$(18)with τ = κ(T_disc,ρ_disc)ρ_discR_M, L is the planet luminosity, and T_disc, ρ_disc, P_disc are the temperature, density and pressure in the disc midplane at the location of the planet.

2.3.1. Gas accretion: detached phase

By solving the differential Eqs. (10) to (12) with the boundary conditions mentioned above, one can derive the planetary envelope mass as a function of time, and therefore the gas accretion rate Ṁ_gas. However, the rate of gas accretion that can be sustained by the protoplanetary disc is not arbitrary, and is in particular limited by the viscosity. When the gas accretion rate required by the forming planet is higher than the one that can be delivered by the disc, Ṁ_gas,max, the planet goes into the detached phase.

In the detached phase, the planetary growth rate by gas accretion does not depend upon its internal structure, but is rather given by the structure and evolution of the disc. During this phase, the internal structure is given by solving the same Eqs. (10) to (12), this time  for a  mass  M_r ranging  from M_core

to M_planet (which is known). The boundary conditions are the same, except for two:

• The pressure, which includes the dynamical pres-sure from gas free-falling from the disc to the planet,${\mathit{P}}_{\mathrm{surf}\mathit{,}\mathrm{planet}}\mathrm{=}{\mathit{P}}_{\mathrm{disc}}\mathrm{+}\frac{\mathit{Ṁ}\mathrm{gas}}{\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathit{M}}^{\mathrm{2}}}{\mathit{v}}_{\mathrm{ff}}\mathit{.}$(19)In this equation, v_ff is the free-falling velocity from the Hill radius to the planetary radius, ${\mathit{v}}_{\mathrm{ff}}\mathrm{=}\mathrm{-}\sqrt{\mathrm{2}\mathit{GM}\mathrm{\times }{}^{\mathrm{(}}\mathrm{1}\mathit{/}{\mathit{R}}_{\mathit{M}}\mathrm{-}\mathrm{1}\mathit{/}{\mathit{R}}_{\mathrm{H}}{}^{\mathrm{)}}}$. The planetary radius is not known a priori, but computed as a result of integrating Eqs. (10) to (12).

• The maximum accretion rate, Ṁ_gas,max, which is equal to $\mathit{Ṁ}\mathrm{gas}\mathit{,}\max \mathrm{=}\begin{array}{}\max \end{array}\left[\mathit{F}\mathrm{(}{\mathit{a}}_{\mathit{M}}\mathrm{+}{\mathit{R}}_{\mathrm{H}}\mathrm{)}\mathit{,}\mathrm{0}\right]\mathrm{+}\begin{array}{}\min \end{array}\left[\mathit{F}\mathrm{(}{\mathit{a}}_{\mathit{M}}\mathrm{-}{\mathit{R}}_{\mathrm{H}}\mathrm{)}\mathit{,}\mathrm{0}\right]\mathit{,}$(20)where $\mathit{F}\mathrm{=}\mathrm{3}\mathit{\pi \nu }\mathrm{\Sigma }\mathrm{+}\mathrm{6}\mathit{\pi r}\frac{\mathit{\partial }\mathit{\nu }\mathrm{\Sigma }}{\mathit{\partial a}}$ is the mass flux in the disc. Geometrically, the maximum accretion rate that can be provided by the disc is equal to the mass flux entering the planet’s gas feeding zone. The gas can enter either from the outer parts of the disc (which is the general case), or from the inner part of the disc (which can be the case in the outer part of the disc).

Therefore, during a time step, a planet has access to the mass delivered at its location by the disc (Ṁ_gas,max × dt), and a mass reservoir made of the gas mass already present in the planet’s gas feeding zone (see also Ida & Lin 2004). This reservoir of gas is assumed to be empty when the planet is massive enough to open a gap (which coincides with the transition to type II migration, see next section). However, the feeding zone continues to receive gas due to viscosity at the local accretion rate (Eq. (20)).

2.4. Orbital evolution: disc-planet interaction

Disc-planet interaction leads to planet migration, which can occur in different regimes. For low-mass planets, which are not massive enough to open a gap in the protoplanetary disc, migration occurs in type I (Ward 1997; Tanaka et al. 2002; Paardekooper et al. 2010, 2011). For higher mass planets, migration is again subdivided into two modes: disc-dominated type II migration, when the local disc mass is larger than the planetary mass (the migration rate is then simply given by the viscous evolution of the protoplanetary disc), and planet-dominated type II migration in the opposite case (see M09). The transition between type I and type II migration occurs when $\frac{\mathrm{3}}{\mathrm{4}}\frac{{\mathit{H}}_{\mathrm{disc}}}{{\mathit{R}}_{\mathrm{H}}}\mathrm{+}\frac{\mathrm{50} {\mathit{M}}_{\mathit{\star }}}{\mathit{M}\mathrm{Re}}\mathrm{=}\mathrm{1}\mathit{,}$(21)(Crida et al. 2006), where H_disc is the disc scale-height at the location of the planet, and $\mathrm{Re}\mathrm{=}\frac{{\mathit{a}}_{\mathit{M}}^{\mathrm{2}}\mathrm{\Omega }}{\mathit{\nu }}$ is the macroscopic Reynolds number at the location of the planet (ν is the same as the one used for the disc evolution).

First models of type I migration (Ward 1997; Tanaka et al. 2002) predicted such rapid migration rates that it was necessary to arbitrarily reduce the migration rate by a constant factor, named f_I in A05 and M09, to reproduce observations. Since these first calculations, type I migration has been studied in great detail, and new formulations for type I migration rates are now available (Paardekooper et al. 2010, 2011). We use in our model an analytic description of type I migration that reproduces the results of Paardekooper et al. (2011). A detailed description of this model is presented in Dittkrist et al. (in prep.), and preliminary results have been presented in Mordasini et al. (2010).

3. The accretion rate of solids

The growth of the solid component of a protoplanet, M_core, is assumed to be caused by to the accretion of planetesimals. Adopting the particle-in-a-box approximation, its growth rate can be calculated with $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{M}}_{\mathrm{core}}}{\mathrm{d}\mathit{t}}\mathrm{=}\left(\frac{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathit{m}}{\mathit{R}}_{\mathrm{H}}^{\mathrm{2}}}{{\mathit{P}}_{\mathrm{orbital}}}\right) {\mathit{P}}_{\mathrm{coll}}\mathit{,}& & \end{array}$(22)(Chambers 2006), where Σ_m is the solids surface density at the location of the protoplanet and P_orbital is its orbital period. The collision rate, P_coll, is the probability that a planetesimal is accreted by the protoplanet. This probability depends upon the relative velocity between planetesimals and the protoplanet which, in turn, depends upon the planetesimal’s eccentricities and inclinations. We represent by e (i) the root mean square of the eccentricity (inclination) of planetesimals. Planetesimals are found to be in different velocity regimes depending on their random velocities. These regimes are known as high-, medium-, and low-velocity regime. Each regime is characterised by a range of values of the planetesimal’s reduced eccentricities ($\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathrm{=}\mathit{ae}\mathit{/}{\mathit{R}}_{\mathrm{H}}$) and inclinations ($\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}\mathrm{=}\mathit{ai}\mathit{/}{\mathit{R}}_{\mathrm{H}}$): the high-velocity regime is defined by $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathit{,}\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}\gtrsim \mathrm{2}$, the medium-velocity regime by $\mathrm{2}\gtrsim \begin{array}{}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathit{,}\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}\gtrsim \mathrm{0.2}$ and the low-velocity regime by $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathit{,}\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}\lesssim \mathrm{2}$. This leads to different collision rates $\begin{array}{ccc}{\mathit{P}}_{\mathrm{high}}& \mathrm{=}& \\ {\mathit{P}}_{\mathrm{med}}& \mathrm{=}& \\ {\mathit{P}}_{\mathrm{low}}& \mathrm{=}& \end{array}$(see Inaba et al. 2001, and references therein), where R is the radius of the protoplanet (in the case of a solid body without a gaseous envelope R is its geometrical radius), r_m is the radius of the planetesimals, $\mathit{\beta }\mathrm{=}\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}\mathit{/}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}$, and the functions I_F and I_G are well-approximated by $\begin{array}{ccc}{\mathit{I}}_{\mathrm{F}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}& \mathrm{\simeq }& \\ {\mathit{I}}_{\mathrm{G}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}& \mathrm{\simeq }& \end{array}$for 0 < β ≤ 1, which is the range of interest for this work (Chambers 2006).

According to Inaba et al. (2001), the mean collision rate can be approximated by ${\mathit{P}}_{\mathrm{coll}}\mathrm{=}\min \mathrm{\left(}{\mathit{P}}_{\mathrm{med}}\mathit{,}\mathrm{\right(}{\mathit{P}}_{\mathrm{high}}^{-2}\mathrm{+}{\mathit{P}}_{\mathrm{low}}^{-2}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}\mathit{.}$(28)When an embryo is able to gravitationally bind gas from its surroundings, it becomes more difficult to define its radius, which is not just the core radius. For the purpose of the collision rate, the capture radius of the protoplanet should depend upon the mass of the protoplanet, upon the planetesimals’ velocity with respect to the protoplanet, upon the density profile of the envelope, ρ(r), and upon the size of the accreted planetesimals (smaller planetesimals are more affected by the gas drag of the envelope and therefore are easier to capture). As in Guilera et al. (2010), here we adopt the prescription of Inaba & Ikoma (2003) where the capture radius R can be obtained by solving the following equation ${\mathit{r}}_{\mathit{m}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\frac{\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{R}}{{\mathit{\rho }}_{\mathit{m}}}\left(\frac{{\mathit{v}}_{\mathrm{rel}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{GM}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{/}\mathit{R}}{\hspace{0.17em}{\mathit{v}}_{\mathrm{rel}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{GM}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{/}{\mathit{R}}_{\mathrm{H}}}\right)\mathit{,}$(29)where ρ_m is the planetesimals’ bulk density, G is the gravitational constant, and the relative velocity v_rel is given by ${\mathit{v}}_{\mathrm{rel}}\mathrm{=}{\mathit{v}}_{\mathrm{k}}\sqrt{\mathrm{5}\mathit{/}\mathrm{8}\hspace{0.17em}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}\hspace{0.17em}{\mathit{i}}^{\mathrm{2}}}\mathit{,}$(30)with v_k the Keplerian velocity (v_k = Ωa). This simple formula for the capture radius approximates more complex models well (as the one described in A05) with the advantage that it reduces the computational time.

It is clear from these equations that the accretion rate of solids depends upon the eccentricities and inclinations of planetesimals, which define their relative velocities with respect to the embryo: the higher the relative velocity, the less likely it is that planetesimals are captured by the embryo. The eccentricities and inclinations of planetesimals are affected by the damping produced by the nebular gas drag, by the gravitational stirring of the protoplanet (protoplanet-planetesimal interactions) and, to a lesser extent, by their mutual gravitational interactions (planetesimal-planetesimal interactions): $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{e}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}& \mathrm{=}& \\ \frac{\mathrm{d}{\mathit{i}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}& \mathrm{=}& \end{array}$The first term represents the effect of the nebular gas drag, the second term the viscous stirring produced by an embryo of mass M, and the third term represents the planetesimal-planetesimal viscous stirring.

The drag force experienced by a spherical body depends upon its relative velocity with respect to the gas. If we consider that the protoplanetary nebula is mainly composed of H₂ molecules, the mean free path of a molecule of gas is $\mathit{\lambda }\mathrm{=}\mathrm{(}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{\sigma }}_{{\mathrm{H}}_{\mathrm{2}}}{\mathrm{)}}^{-1}\mathit{,}$(33)where n_H2 is the number density of H₂ molecules and σ_H2 is the collision cross-section of an H₂ molecule. Depending upon the ratio between the planetesimal’s radius and the mean free path of the molecules, three drag regimes can be defined (Rafikov 2004, and references therein). The first two drag regimes are for planetesimals whose radii are larger than the mean free path, r_m ≳ λ. These are the quadratic and the Stokes regime. To distinguish these regimes we adopt the criterion proposed by Rafikov (2004) in terms of the molecular Reynolds number Re_mol ≡ v_relr_m/ν_mol, where ν_mol is the molecular viscosity, ν_mol = λc_s/3. If Re_mol ≳ 20, we assume that the gas drag is in the quadratic regime, and the differential equations for the evolution of the eccentricity and inclination are given by $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{e}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{drag}}& \mathrm{=}& \\ \frac{\mathrm{d}{\mathit{i}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{drag}}& \mathrm{=}& \end{array}$(Adachi et al. 1976 corrected by Inaba et al. 2001)⁴, with ξ ≃ 1.211. The value of η depends upon the distance to the star, on the gas density and on the pressure gradient, dP/da, $\mathit{\eta }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathrm{\Omega }}^{\mathrm{2}}\mathit{a}{\mathit{\rho }}_{\mathrm{gas}}}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{a}}\mathit{,}$(36)where ρ_gas and dP/da are derived from the disc model.

The gas drag timescale is ${\mathit{\tau }}_{\mathrm{drag}}\mathrm{=}\frac{\mathrm{8}{\mathit{\rho }}_{\mathit{m}}{\mathit{r}}_{\mathit{m}}}{\mathrm{3}{\mathit{C}}_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{gas}}{\mathit{v}}_{\mathrm{k}}}\mathit{,}$(37)where C_D is the drag coefficient, which is of the order of unity.

The Stokes regime occurs for r_m ≳ λ and Re_mol < 20, and the equations for the eccentricity and inclinations of planetesimals are $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{e}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{drag}}& \mathrm{=}& \\ \frac{\mathrm{d}{\mathit{i}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{drag}}& \mathrm{=}& \end{array}$(Adachi et al. 1976; Rafikov 2004).

In this paper, as for example in Stepinski & Valageas (1996), we defined two different Reynolds numbers, Re and Re_mol, and two different viscosities, ν and ν_mol. The macroscopic quantities (Re and ν) are a measure of the fluid dynamics of the disc in a global scale (to compute the evolution of the disc) while the microscopic quantities (Re_mol and ν_mol) characterise the local state of the gas and are used to calculate the eccentricities and inclinations of planetesimals.

When r_m ≲ λ, the third regime, Epstein regime, takes place, and the evolution of eccentricities and inclinations follows the equations $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{e}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{drag}}& \mathrm{=}& \\ \frac{\mathrm{d}{\mathit{i}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{drag}}& \mathrm{=}& \end{array}$(Adachi et al. 1976; Rafikov 2004). In this paper we consider that the population of planetesimals is represented by spherical bodies of a single size. It is worth mentioning that although we allow for the three drag regimes according to the above-mentioned criterions, for the ranges of planetesimal sizes considered in this work (100 − 0.1 km) and for the kind of interaction we are mostly interested in here (protoplanet-planetesimal interactions), planetesimals are found to be mainly in the quadratic regime. Therefore, in most of the cases, for determining the solids accretion rate the effect of the gas drag is governed by Eqs. (34), (35).

The planetesimal eccentricities and inclinations are excited by the presence of a protoplanet. Ohtsuki et al. (2002) studied the evolution of the mean square orbital eccentricities and inclinations and introduced semi-analytical formulae to describe the stirring produced by the protoplanet $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{e}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{VS}\mathit{,}\mathit{M}}& \mathrm{=}& \\ \frac{\mathrm{d}{\mathit{i}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{VS}\mathit{,}\mathit{M}}& \mathrm{=}& \end{array}$where b is the full width of the feeding zone of the protoplanet in terms of their Hill radii (here we adopt b ~ 10), and P_VS and Q_VS are given by $\begin{array}{ccc}{\mathit{P}}_{\mathrm{VS}}& \mathrm{=}& \left[\frac{\mathrm{73}\stackrel{}{{{\mathit{e}}^{\mathrm{2}}}_{\mathrm{˜}}}}{\mathrm{10}{\mathrm{\Lambda }}^{\mathrm{2}}}\right]\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathrm{10}{\mathrm{\Lambda }}^{\mathrm{2}}\mathit{/}\stackrel{}{{{\mathit{e}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{)}\\ & & \\ {\mathit{Q}}_{\mathrm{VS}}& \mathrm{=}& \left[\frac{\mathrm{4}\begin{array}{}\mathrm{˜}\\ {\mathit{i}}^{\mathrm{2}}\end{array}\mathrm{+}\mathrm{0.2}\begin{array}{}\mathrm{˜}\\ \mathit{i}\end{array}\stackrel{}{{{\mathit{e}}^{\mathrm{3}}}_{\mathrm{˜}}}}{\mathrm{10}{\mathrm{\Lambda }}^{\mathrm{2}}\begin{array}{}\mathrm{˜}\\ \mathit{e}\end{array}}\right]\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathrm{10}{\mathrm{\Lambda }}^{\mathrm{2}}\stackrel{}{{{\mathit{e}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{)}\\ & & \mathrm{+}\left[\frac{\mathrm{72}{\mathit{I}}_{\mathrm{QVS}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}{\mathit{\pi }\begin{array}{}\mathrm{˜}\\ \mathit{e}\end{array}\begin{array}{}\mathrm{˜}\\ \mathit{i}\end{array}}\right]\ln \mathrm{(}\mathrm{1}\mathrm{+}{\mathrm{\Lambda }}^{\mathrm{2}}\mathrm{)}\mathit{,}\end{array}$with $\mathrm{\Lambda }\mathrm{=}\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}\mathrm{(}\stackrel{}{{{\mathit{e}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{+}\begin{array}{}{}^{\mathrm{˜}}\\ {\mathit{i}}^{\mathrm{2}}\end{array}\mathrm{)}\mathit{/}\mathrm{12}$. The functions I_PVS(β) and I_QVS(β) can be approximated for 0 < β ≤ 1 by $\begin{array}{ccc}{\mathit{I}}_{\mathrm{PVS}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}& \mathrm{\simeq }& \\ {\mathit{I}}_{\mathrm{QVS}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}& \mathrm{\simeq }& \end{array}$(Chambers 2006). The excitation that the protoplanet produces on the planetesimals weakens with the increase in the distance between the protoplanet and the planetesimals, i.e. farther away planetesimals are less excited. Here we follow the approach of Guilera et al. (2010) and consider that the effective stirring is given by $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{e}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{VS}\mathit{,}\mathit{M}}^{\mathrm{eff}}& \mathrm{=}& \\ \frac{\mathrm{d}{\mathit{i}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{VS}\mathit{,}\mathit{M}}^{\mathrm{eff}}& \mathrm{=}& \end{array}$where f(Δ) ensures that the perturbation of the protoplanet is confined to its neighbourhood, $\mathit{f}\mathrm{\left(}\mathrm{\Delta }\mathrm{\right)}\mathrm{=}{\left[\mathrm{1}\mathrm{+}{\left(\frac{\mathrm{\Delta }}{\mathit{n}{\mathit{R}}_{\mathrm{H}}}\right)}^{\mathrm{5}}\right]}^{-1}\mathit{,}$(50)with Δ =  | a_M − a_m | , where a_M is the semi-major axis of the protoplanet and a_m is the semi-major axis of the planetesimal. Although the functional form is arbitrary, the scale on which the stirring acts is similar to the one found in N-body calculations (excluding the effects of resonances). For this work we have chosen n = 5 to limit the perturbation of the planet to its feeding zone. In the future, with the aid of N-body calculations, we plan to obtain a better semi-analytical function to characterise the extent of the planetary perturbation.

We also consider that the planetesimals’ eccentricities and inclinations are stirred by their mutual interactions. For a population of planetesimals of equal mass m, the evolution of their eccentricities and inclinations are well described by $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{e}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{VS}\mathit{,}\mathit{m}}& \mathrm{=}& \\ \frac{\mathrm{d}{\mathit{i}}^{\mathrm{2}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{VS}\mathit{,}\mathit{m}}& \mathrm{=}& \end{array}$(Ohtsuki et al. 2002), with ${\mathit{h}}_{\mathit{m}}\mathrm{=}{\left(\frac{\mathrm{2}\mathit{m}}{\mathrm{3}{\mathit{M}}_{\mathit{\star }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(53)In this case P_VS and Q_VS are evaluated with the reduced eccentricity and inclination relative to the planetesimal mass (i.e. $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathrm{=}\mathrm{2}\mathit{e}\mathit{/}{\mathit{h}}_{\mathit{m}}$, $\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}\mathrm{=}\mathrm{2}\mathit{i}\mathit{/}{\mathit{h}}_{\mathit{m}}$). There is no dynamical friction term in Eqs. (51), (52), as it vanishes when a single mass population of planetesimals is considered (Ohtsuki et al. 2002). Although strictly speaking we have two populations of planetesimals (rocky and icy bodies depending on whether they are inside or beyond the ice line) that have the same size but not the same mass (because of the difference in their density), the region where the two types of planetesimals are present at the same time is very narrow.

In this work we neglect changes in the mass of planetesimals due to fragmentation and changes in the surface density owing to planetesimal drifting.

3.1. Comparison with previous accretion rate of solids

In previous works (A05, M09), the accretion rate of solids has been treated in a very simple way which led to an underestimation of the formation timescale of planets. In those works, the prescription for the planetesimals’ eccentricities and inclinations was the same as in Pollack et al. (1996), who assumed that the planetesimals’ inclinations depend only on planetesimal-planetesimal interactions. Under this assumption, the reduced value of the planetesimals’ inclination, $\begin{array}{}{\mathrm{˜}}^{}\\ \mathit{i}\end{array}$, was prescribed as $\begin{array}{}\mathrm{˜}\\ \mathit{i}\end{array}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{E}}}{\sqrt{\mathrm{3}}\mathrm{\Omega }{\mathit{R}}_{\mathrm{H}}}\mathit{,}$(54)where v_E is the escape velocity from the surface of a planetesimal. This means that the planetesimal’s inclination $\mathit{i}\mathrm{=}\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}{\mathit{R}}_{\mathrm{H}}\mathit{/}\mathit{a}$ is constant independently of the mass of the planet. On the other hand, eccentricities are assumed to be controlled by both planetesimals and protoplanet stirring, its value given by $\begin{array}{}\mathrm{˜}\\ \mathit{e}\end{array}\mathrm{=}\max \mathrm{\left(}\mathrm{2}\begin{array}{}\mathrm{˜}\\ \mathit{i}\end{array}\mathit{,}\mathrm{2}\mathrm{\right)}\mathit{.}$(55)Therefore, if $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathrm{=}\mathrm{2}\begin{array}{}{}^{\mathrm{˜}}\\ \mathit{i}\end{array}$, the protoplanet is growing according to the runaway regime because e and i would be independent of the mass of the protoplanet. If $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{e}\end{array}\mathrm{=}\mathrm{2}$, the eccentricity of planetesimals would be affected by the presence of the protoplanet, so to a certain extent the stirring of the embryo is taken into account. However, this condition corresponds to the protoplanet-planetesimal scattering in the shear-dominated regime. Ida & Makino (1993) showed that the shear-dominated regime lasts for only a few thousand years, after which planetesimals are strongly stirred by the protoplanet. During the shear-dominated period, eccentricities and inclinations of planetesimals in the vicinity of the protoplanet remain low. This leads to an accretion scenario that is much faster than that corresponding to the oligarchic regime (the oligarchic regime usually occurs in the dispersion-dominated regime). To clearly show the difference between this quasi-runaway accretion of solids and the oligarchic regime, we performed two simulations that are identical in all parameters except for the prescription of e and i. The planet is assumed to form in situ at 6 AU. Accreted planetesimals are 100 km in radius. No disc evolution is considered, so simulations are stopped when the planet reaches one Jupiter mass. For the quasi-runaway regime we use Eqs. (55) and (54) to calculate e and i, while for the oligarchic growth we solve the differential equations presented in the previous section. Figure 2 (left panel) shows the ratio of oligarchic and quasi-runaway eccentricities and inclinations as a function of the mass of the planet. In the case of the eccentricity, the one corresponding to the oligarchic regime is  ~4 times higher than in the quasi-runaway regime. The oligarchic inclination is several tens of times higher than the quasi-runaway regime. As a consequence, the accretion rate of solids is much lower in the oligarchic than in the quasi-runaway regime because in the oligarchic regime planetesimals are more excited and more difficult to accrete. A comparison between the two accretion rates of solids is shown in Fig. 2 (right panel). Owing to the lower accretion rate in the oligarchic regime, while in the quasi runaway regime it takes less than 1 Myr to form the planet, in the oligarchic regime formation is much longer, taking 3.25 × 10⁷ years.

Fig. 2 Left panel: the ratio between eccentricities (red) and inclinations (blue) as a function of the planet’s mass considering both the oligarchic and the quasi-runaway growth. Clearly, the planetesimals’ excitation is several times higher in the oligarchic growth than in the quasi-runaway growth. The right panel shows the impact of the planetesimals’ excitation in the accretion rate of solids: for the quasi-runaway growth (yellow) the accretion rate is much higher than in the oligarchic growth (green).

4. Results

In the previous section we introduced the main characteristics of the planet formation model, with special focus on the differences in the physical and numerical model with respect to A05. In the following sections we will concentrate on the impact that the accretion rate of solids has on the formation of giant planets. Here we will consider the formation of isolated planets, i.e. only one planet per disc. The computations of planetary system formation will be presented in other papers (Alibert et al. 2012; Carron et al. 2012).

As we described in Sect. 3, the treatment of the evolution of eccentricities and inclinations of planetesimals intends to minimise the assumptions on their values (keeping in mind that it is not an N-body calculation, but the adopted formulas reproduce N-body results of the planetesimals’ accretion rates and excitation). To consider a realistic accretion rate that is not too computationally expensive, Thommes et al. (2003) considered that the planetesimals’ eccentricities and inclinations can be estimated assuming that the stirring produced by the protoplanet is instantaneously balanced by the gas drag. The approximation to the equilibrium values of e (e_eq) can be derived by equating the stirring timescale and the damping timescale, resulting in ${\mathit{e}}_{\mathrm{eq}}\mathrm{=}\mathrm{1.7}\frac{{\mathit{m}}^{\mathrm{1}\mathit{/}\mathrm{15}}{\mathit{M}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\rho }}_{\mathit{m}}^{\mathrm{2}\mathit{/}\mathrm{15}}}{{\mathit{b}}^{\mathrm{1}\mathit{/}\mathrm{5}}{\mathit{C}}_{\mathrm{D}}^{\mathrm{1}\mathit{/}\mathrm{5}}{\mathit{\rho }}_{\mathrm{gas}}^{\mathrm{1}\mathit{/}\mathrm{5}}{\mathit{M}}_{\mathit{\star }}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{a}}^{\mathrm{1}\mathit{/}\mathrm{5}}}\mathrm{·}$(56)The equilibrium value for i (i_eq) is assumed to be half the value of e_eq, as this relationship has been shown to be a good approximation in the high-velocity cases (Ohtsuki et al. 2002): ${\mathit{i}}_{\mathrm{eq}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{e}}_{\mathrm{eq}}\mathit{.}$(57)However, it is not clear whether planetesimals are always in equilibrium, especially if we consider that depending on their mass, planetesimals are differently affected by gas drag, and that during its formation, a planet migrates and the protoplanetary disc evolves. On the other hand, if equilibrium is attained, it is interesting to compare the equilibrium values obtained by explicitly solving Eqs. (31), (32) with the approximations given by Eqs. (56), (57) based on timescale considerations.

Fortier at al. (2007, 2009) and Benvenuto et al. (2009) assumed, for simplicity, the equilibrium approximation of Thommes et al. (2003) in their in situ, giant planet formation models. However, Chambers (2006) found that significant deviations from equilibrium occur at the very beginning of the growth of the embryo, but eventually equilibrium is attained. In the cases he discussed, these deviations do not seem to have a noticeable effect on the final mass of the planet as long as no time restriction for the lifetime of the disc is assumed. On the other hand, Guilera et al. (2010, 2011) did not use any approximations and explicitly calculated e and i by solving the corresponding time evolution differential equations in their giant planet formation model. No study comparing the equilibrium approximation and the explicit calculation of e and i has been made using a self-consistent giant planet formation model. Moreover, out-of-equilibrium effects can be important, not only at the beginning of the formation of a planet. When planets migrate, they can enter regions where planetesimals are in principle cold (low values of e and i), or already excited if there is another planet growing in the neighbourhood. Depending on the ratio between the stirring and the migration timescale, one can expect some cases where the equilibrium approximation may not be accurate.

We study the formation of giant planets considering a self-consistent model for the interplay between the disc evolution, accretion by the growing planets and gas-driven migration. In this paper, we focus our study on the formation of single planets, looking in detail at the dependence of planetary growth upon the planetesimal size and the differences in the final results between the equilibrium approximation and the explicit calculation of e and i, all together with planetary migration. We perform this analysis in four steps. First we consider the formation of a 1 M_⊕ planet, neglecting the presence of an envelope, of planetary migration and of disc evolution (Sect. 4.1). Second, we compute the full formation of a planet (except migration), which is assumed to be over when the gaseous component of the disc disappears. We keep the in situ formation hypothesis to facilitate the analysis (Sect. 4.2). Third, we allow for gas driven migration during the formation of the planet (Sect. 4.3). Fourth, we generalise the examples presented in the third step by considering a wide range of plausible protoplanetary discs and initial locations for the embryo to have an overview of all possible outcomes (Sect. 4.4).

Fig. 3 In situ formation of a 1 M⊕ planet at 6 AU for different radii of accreted planetesimals (red, 100 km; blue, 10 km; green, 1 km; orange, 0.1 km). The solid line is for the equilibrium approximation of eccentricities and inclinations (see Eqs. (56), (57)). The dashed and dotted lines are for the explicit calculation of e and i solving the differential equations Eqs. (31), (32). Two different initial values, e0, i0, were considered: the dashed line represents the case where e0 and i0 are given by Eqs. (59), (58), while the dotted line assumes e0, i0 to be the initial equilibrium values of the equilibrium case. The top-left panel depicts the mass growth of the protoplanet as a function of time. The top-right panel shows the time evolution of the eccentricity. The bottom-left panel shows the eccentricity as a function of the mass of the embryo. Clearly, in all cases, equilibrium is attained after an out of equilibrium state. The bottom-right panel plots the accretion rate of solids as a function of the embryo mass. Note that, as a function of the planet mass, the difference in the accretion rate is only a consequence of the difference in e and i.

4.1. Formation of a 1 M_⊕ planet

We analyse first the formation of a 1 M_⊕ planet (we stop the calculation when this mass is reached), neglecting the presence of an envelope, planetary migration, and disc evolution. The point of this case is to focus on the initial stages of the accretion and analyse the importance of the size of the accreted planetesimals and the implications for the growth of the embryo when considering the equilibrium approximation. This is because we aim to clearly show the consequences of these assumptions that we neglect other physical processes that act simultaneously. This means, for example, that for a fixed mass of the embryo and independently of the elapsed time, the state of the protoplanetary disc is the same in terms of surface density (solids and gas). The same applies for the capture radius of the planet. In this example, the embryo is assumed to be located at 6 AU, where the initial solids surface density is Σ_m = 10   g   cm^-2 and the density of the nebular gas is ρ_gas = 2.4    × 10^-9   g   cm^-3. For this disc the snow line is at 3.5 AU. The initial surface density profile of the disc is given by Eq. (8) with γ = 0.9 and a_C = 127 AU. The initial mass of the embryo is 0.01 M_⊕. For the equilibrium approximation we adopt Eqs. (56), (57) to calculate the values of e and i. For the explicit calculation of the eccentricity and inclination of planetesimals we solve Eqs. (31), (32). To solve Eqs. (31), (32), initial conditions for e and i must be given. We consider two possibilities for the initial conditions that we think bracket the parameter space. On one hand, we consider that the planetesimal disc is initially cold, and the planetesimals’ eccentricities and inclinations are given by the equilibrium value between their mutual stirring and the gas drag. These values can be derived by equating the stirring timescale and the damping timescale, which results in $\begin{array}{ccc}{\mathit{e}}_{\mathrm{eq}}^{\mathit{m}\mathrm{-}\mathit{m}}& \mathrm{=}& \\ {\mathit{i}}_{\mathrm{eq}}^{\mathit{m}\mathrm{-}\mathit{m}}& \mathrm{=}& \end{array}$This means that we are assuming that the embryo instantaneously appears in an unperturbed planetesimal disc. The other extreme situation is to assume a hot disc, where planetesimals are already excited by the embryo and their initial eccentricities and inclinations are those corresponding to the value of equilibrium between the stirring of the embryo (0.01 M_⊕) and the gas drag, approximated by Eqs. (56), (57). The initial values of e for these cases are given in Table 2, where ${\mathit{e}}_{\mathrm{eq}\mathit{,}\mathrm{0}}^{\mathit{m}\mathrm{-}\mathit{m}}$ corresponds to Eq. (58) and e_eq,0 to Eq. (56). The initial values of i are assumed to be e/2 (i₀ = e₀/2). We perform calculations for the equilibrium approximation and the two sets of initial conditions for the explicit calculation for four radii of the accreted planetesimals: 100, 10, 1, and 0.1 km.

Figure 3 shows the results of the simulations mentioned above. The top-left panel depicts the mass growth of the planets as a function of time. What is first evident from this plot is the timescale difference in the formation of the 1 M_⊕ planet depending on the size of the accreted planetesimals: while for accreted planetesimals of 100 km it takes  ~10⁷ yr to grow from 0.01 M_⊕ mass to 1 M_⊕, it takes  ~10⁵ yr in the case of 0.1 km planetesimals. This difference in the growth timescale is entirely due to the fact that large planetesimals are less damped by gas drag than the smaller ones. Hence, while they are stirred up by the massive embryos to about the same value, they keep higher e and i values (bottom-left panel), making the accretion process much slower. We recall that the disc does not evolve and planets do not migrate, therefore differences depend only on the planetesimal size. Note that for a fixed mass of the embryo the accretion rate of solids differs by two orders of magnitude between the two extreme cases (bottom-right panel).

If we now turn our attention to the differences in growth rate for a fixed planetesimal size, but for different approaches in the calculation of the eccentricities and inclinations, we see that adopting the equilibrium approximation may lead to significant differences in the mass of the planet, more evident when the accreted planetesimals are small (top-left panel, compare the solid line with the dashed or the dotted lines of the same colour). For 100 km and 10 km planetesimals we do not see deviations from equilibrium when the initial conditions for e and i are given by Eqs. (56), (57) (red and blue dotted lines, top-right panel). If the initial values for e and i are given by Eqs. (58), (59), the equilibrium values given by Eqs. (56), (57) are reached in an almost negligible fraction of the formation timescale. We conclude that the equilibrium approximation for e and i gives results that nicely agree with more complex calculations, as far as the planetesimal size is relatively big or the evolution time is sufficiently long. However, as shown in the same picture, the growth of a small planet is very long when considering such massive planetesimals, and the formation of a gas giant under these conditions is highly compromised.

The situation becomes more critical when we consider smaller planetesimals. Deviations from equilibrium are evident regardless of the initial values adopted for e and i, especially for r_m = 0.1 km. Nevertheless, equilibrium is always attained, although the equilibrium values for e and i are lower than those given by Eqs. (56), (57), especially for the inclination. As we can see from Fig. 4, the approximation β_eq = i/e = 1/2 is poor for smaller planetesimals. Note that the real equilibrium value of β depends upon the planetesimal size. While β_eq = 1/2 is a good approximation for planetesimals larger than  ~1 km, this is not the case for smaller planetesimals, which tend to have inclination values lower than 1/2   e. We find that for small planetesimals, the velocity regime is at the limit between the high- and medium-velocity regime ($\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{e}\end{array}$ , $\begin{array}{}{\mathrm{˜}}^{}\\ \mathit{i}\end{array}$  ≲  2), therefore eccentricities are more effectively excited than inclinations (see Ohtsuki et al. 2002). The low values of e, i, and β increase the accretion rate, as can be seen in the bottom-right panel of Fig. 3, speeding up the formation of the embryo (relative to the equilibrium approximation case).

Fig. 4 Inclination-eccentricity ratio for the out-of-equilibrium cases of Fig. 3. Red is for 100 km planetesimals, blue for 10 km, green for 1 km, and orange for 0.1 km. The dotted line represents the case of an initially hot disc, while the dashed line represents an initially cold one. The grey line is the standard equilibrium value, βeq = i/e = 1/2.

4.2. Planet formation without migration

Fig. 5 In situ planetary formation at 6 AU for different radii of the accreted planetesimals (red, 100 km; blue, 10 km; green, 1 km; orange, 0.1 km). The same line code as in Fig. 3 is adopted. The formation of the planet ends when the gas component of the disc is dispersed.

In this section we analyse the complete formation of a planet for the same cases as before. Now, the planet accretes mass (solids and gas) as long as the disc is still present. This means that we allow for disc evolution (Sect. 2.1), therefore planet formation is considered to be over when the disc disperses. The presence of the gaseous envelope is taken into account for calculating the capture radius. However, we still neglect the migration of the embryo to facilitate the analysis of the dependence on the planetesimals’ e and i.

As has already been shown in other works, the presence of the envelope increases the capture radius, speeding up the formation of the planet. The enhancement in the capture cross section makes the accretion of solids more effective. This effect is more noticeable for small planetesimals. Compared to the cases of Sect. 4.1, the formation timescale of a 1 M_⊕ can be reduced to  ~35% (for the smallest planetesimals, r_m = 0.1 km), due only to the enhancement in the capture cross section. The presence of an atmosphere, even if its mass is negligible compared to the total mass of the planet, has to be considered for embryos as small as 0.1 M_⊕ because it plays an important role for the accretion of solids.

Figure 5 is the analog of Fig. 3 for the complete formation of the planets. As the disc disperses after 6 Myr, in the cases where the accretion of planetesimals is slow (r_m = 100,10 km), the final masses of the planets are lower than 1 M_⊕. The differences in growth observed in the previous section as a function of the planetesimal size have dramatical consequences for the final mass of the planet, which can be  ~0.1 M_⊕ if the accreted planetesimals have a radius of 100 km,  ~0.8 M_⊕ for planetesimals of 10 km,  ~1200 M_⊕ (3.7 Jupiter masses, M_J) for 1 km planetesimals, and 7100 M_⊕ (22 M_J) for 0.1 km planetesimals. These numbers show how differences in the accretion rate of solids (here regulated by the size of the accreted planetesimals, bottom-right panel) impact on the final mass of a planet and the non linear aspect of planet formation in the core accretion model: once the critical mass is attained, the high accretion of gas rapidly leads to massive planets. The fact that bigger planets form when the accreted planetesimals are small is a consequence of the gas drag, which operates in two ways that combine positively: nebular gas drag is more effective in damping the planetesimal’s eccentricities and inclinations when planetesimals are small (bottom-left panel) and atmospheric gas drag is able to deflect the more distant planetesimals’ trajectories, therefore enlarging the capture radius of the planet. In fact, for the cases of 1 km and 0.1 km accreted planetesimals, embryos grow to become big giants. When a massive solid embryo is formed, it triggers the accretion of gas, leading to the formation of a gas giant planet (green and orange lines).

Indeed, planets can end up being very massive if they enter the runaway phase of gas. As explained before, during the attached phase, the accretion of gas is a result of solving the equations presented in Sect. 2.3. The planet will remain attached to the disc until its accretion rate exceeds the maximum amount of gas that the disc can deliver. When this condition is not longer satisfied, it goes into the detached phase. During the detached phase, the maximum accretion rate is given by Eq. (20). As in M09 and M12, here we follow the results of Kley & Dirksen (2006) for the accretion of gas during the detached phase. Planet-disc interactions can lead to eccentric instabilities, which means that the planet can enter regions that are outside its gap and have full access to the gas present in the disc, with no limitation for accretion except for the ability of the disc to supply it. Although it is not clear whether all planets can suffer from an eccentric instability, for the sake of simplicity we assume this is the general situation in our simulations. On the other hand, other works (e.g. Lissauer at al. 2009) include a limitation for the accretion of gas when the planet opens a gap in the disc, because they consider that planets are in circular orbits. Comparatively, this assumption leads to less massive planets.

The equilibrium approximation and the explicit calculation of e and i also bring differences in the final mass of the planet. While for large planetesimals (100 and 10 km) the equilibrium approximation works fine, as shown in the previous section, for smaller planetesimals this might not be the case. For r_m = 1 km, the final mass of the planet in the equilibrium approximation is 537 M_⊕, while with the explicit calculation considering an initially hot disc it is 950 M_⊕ and with an initially cold disc, 1200 M_⊕. For r_m = 0.1 km, the final mass of the planet in the equilibrium approximation is 4745 M_⊕, while with the explicit calculation considering an initially hot disc it is 6745 M_⊕ and with an initially cold disc, 7100 M_⊕. Differences in mass are the consequence of a timing effect. When the planet grows faster, the damping produced by the nebular gas is more effective because the gas density in a younger disc is higher than in an older one. Then, the faster an embryo grows, the more it profits from the nebular gas drag. Moreover, the cross-over mass (the mass of the core for which the mass of the envelope equals the mass of the core), which in this case is  ~30 M_⊕, is achieved earlier when the accretion rate of solids is higher. When a protoplanet reaches the cross-over mass its growth is dominated by the accretion of gas that is already in the runaway regime. If the protoplanet starts the accretion of gas when the disc is younger and more massive, it provides a larger reservoir of gas which translates, in the end, into a bigger planet.

It is important to remark that Eqs. (56), (57) are obtained assuming that the stirring timescale is equal to the damping timescale (Thommes et al. 2003), and the equilibrium approximation (i = 0.5   e). On the other hand, Eqs. (31), (32) explicitly solve the coupled evolution of e and i given a set of initial conditions. This means that out-of-equilibrium states are allowed (e.g. as seen at the beginning of the calculations) and equilibrium states are reached naturally and without a fixed ratio between e and i. Clearly, the explicit resolution of the differential equations is a better physical approach. We here included calculations with the equilibrium approximation just for comparison. Results show that the equilibrium approximation works fine for larger planetesimals, but overestimates the excitation of smaller ones, making their accretion less effective. This delays the whole process of planet formation. Because the planet is embedded in a disc with a finite (and short) lifetime, this delay impacts on the final mass of the planet. When explicitly solving the equations of e and i, initial conditions for these quantities are needed. This is a problem because we cannot be certain about the state of the disc at the beginning of our calculations. So assumptions for the initial values cannot be avoided. The initially cold disc favours the growth of a planet allowing for high accretion rates in the first thousands of years (Fig. 5, bottom-right panel). In the initially hot disc this effect is not present because departures from equilibrium are smaller and therefore equilibrium is attained faster (we recall that the equilibrium attained when solving the differential equations can be different from the equilibrium approximation for small planetesimals). However, final results do not strongly depend upon the initial choices. The differences in the final mass are 20% and 5% for r_m = 1 km and r_m = 0.1 km, respectively. Given all the uncertainties of the model, these differences are acceptable.

Fig. 6 Full planetary formation, now allowing for the migration of the embryos, for different radii of the accreted planetesimals (red, 100 km; blue, 10 km; green, 1 km; orange, 0.1 km). The same line code as in Fig. 3 is adopted. Embryos are initially located at 6 AU. The left panel shows the cumulative mass of the embryo (solids plus gas) as a function of time (in logarithmic scale). The right panel depicts the migration path of the protoplanet (in linear scale). For the accretion of 1 km planetesimals, the embryos are lost in the Sun (indicated with crosses in the plot). The same happens for the planet that grows by accretion of 0.1 km planetesimals in the equilibrium approximation (orange solid line). When planets are lost, calculations are stopped.

4.3. Planet formation with migration

We now complete our analysis by including the migration of the protoplanets. The migration model is the one presented in Dittkrist et al. (in prep.) and Mordasini et al. (2010). Figure 6 shows the total mass of the planet and its semi-major axis as a function of time for the same cases as in the previous sections. Clearly, the situation is very different from the in situ hypothesis. To start with, in all cases the final location of the planet is far from its initial location. Indeed, the three calculations for 1 km planetesimals result in lost planets (planets cross the inner edge of the disc and are considered to fall into the central star). The same is the case for the equilibrium approximation of 0.1 km planetesimals.

Fig. 7 Mass versus semi-major axis for the cases shown in Fig. 6. Note that the accretion rate of solids plays a major role not only in the growth of the planet but also in its migration path.

Fig. 8 Protoplanet’s growth and migration timescales (τgrowth = M/ | Ṁ |  and τmig = a/ | ȧ |  respectively) as a function of the protoplanet’s mass for the four different sizes of the accreted planetesimals of Fig. 6 (corresponding to the dashed line in Fig. 6). The peaks in the migration timescale correspond to the changes in the sense of migration of the protoplanet. From the two bottom panels, it can be seen that when the protoplanet’s mass is  ≳10 M⊕, both timescales become comparable and eventually the migration timescale becomes shorter than the accretion timescale, leading to a fast migration of the protoplanet with very little accretion. For 1 km planetesimals the protoplanet growth is too slow to gain enough mass to change the type of migration before they are lost in the star (when its mass is  ~30 M⊕). For 0.1 km planetesimals, the protoplanet can grow more massive because its accretion rate is higher than in the previous case, then it is able to open a gap in the disc and migration switches to type II, preventing it from falling in the star.

For the cases where planets are not lost in the star, we find that the final masses for accreted planetesimals of r_m =  100 km and 10 km are independent of the way e and i were computed (~0.2 M_⊕ and 0.7 M_⊕, respectively). The final locations are also similar for the different considerations for e and i (~1.7 AU for both r_m = 100 km and 10 km). For a fixed planetesimal size, the migration paths for the different assumptions on the eccentricities and inclinations are somehow shifted in time but are similar if we consider them as a function of mass (Fig. 7).

In the case of r_m = 1 km, protoplanets are lost in the central star when their mass is  ~20 M_⊕, most of which is in the solid core. At  ~8 M_⊕ the protoplanet undergoes inward migration in the adiabatic saturated regime⁵ (see Paardekooper et al. 2010, 2011; Mordasini et al. 2010; Dittkrist et al.). This timescale turns out to be shorter than the accretion timescale. The protoplanet covers Δa ≃  5 AU in  ~7 × 10⁵ yr. In this time it doubles its mass, which, however, is not big enough for the planet to enter type II migration. Planet migration in general slows down when the protoplanet is able to open a gap in the disc (type II migration), which, as a rule of thumb, happens when the planet mass is  ~100 M_⊕. This situation is not reached in this case, where accretion is too slow compared to migration, resulting in the loss of the forming planet. For r_m =  0.1 km the differences between the equilibrium approximation and the explicit calculation of e and i are more dramatic: adopting the equilibrium approximation leads to the loss of the planet in the central star (for the same reason as in the previous case, the growth rate is very slow) while for the explicit calculation of e and i, although the planet ends in an orbit very close to the central star (~0.2 AU), the previous growth of the embryo is fast enough to enable a large accretion of gas. In this case the planet reaches a mass that enables it to switch to type II migration. The planet decelerates its migration speed until it stops. The final masses are 13 M_J for an initially hot disc and 15 M_J for an initially cold disc. The fate of an embryo (becoming a high- or a low-mass planet, surviving or being lost in the central star), as it is shown here, depends upon the size of the accreted planetesimals and on the assumptions we adopt to describe their dynamics, because these strongly impact on the accretion rate of solids and therefore on the whole formation process through the regulation of the growth timescale. Here we have shown the interplay between the evolution of the protoplanetary disc, the growth of the protoplanet, and the operation of migration. In all the cases we are considering the same disc, that globally evolves with time in the same way. Differences in the local evolution of the disc arise due to accretion (solids and gas accreted by the planet are removed from the disc) and ejection of planetesimals (when the planet is massive enough). That in Fig. 7, for the same mass of a protoplanet, the location in the disc can be different is a consequence of this interplay: protoplanets reach a certain location earlier or later in the evolution of the disc, depending on their growth rate. The state of the disc and the mass of the protoplanet at that moment determines its migration rate. Accordingly, independently of what regulates the growth of the planets, the examples here show that planet’s growth and migration rate are tightly coupled.

Figure 8 shows a comparison between the migration timescale (τ_mig = a/ | ȧ | ) and the protoplanet’s growth timescale (τ_growth = M/ | Ṁ | ) for the four sizes of the accreted planetesimals we have considered, solving explicitly the equations for e and i and adopting an initially cold planetesimal disc. For accreted planetesimals of 100 km and 10 km, the jump in the growth rate of the planets at M ≳ 0.1 M_⊕ corresponds to the crossing of the ice line. Planetesimals in the inner region (inside the ice line) are denser and the gas drag is less effective on them (see Eqs. (34), (35) and (37)), therefore their random velocities are higher and accretion rates are lower. Also, the solids surface density is lower, which contributes to decrease the accretion rate. The peaks in the migration timescale correspond to changes in the sense of migration⁶.

When the accreted planetesimals have a radius of 1 km, the migration timescale becomes lower than the growth timescale when the mass of the protoplanet is  ~10   M_⊕. The planetesimals relative velocities are very high in the neighbourhood of the protoplanet and accretion becomes more difficult as the protoplanet increases its mass. The protoplanet grows slowly and migrates fast. This situation is never reverted and the protoplanet is lost in the central star. Because the migration rate is high and the accretion rate is not enough to counteract it, the protoplanet migrates without limit, almost at a constant mass, until it reaches the central star.

For accreted planetesimals of 0.1 km, the protoplanet’s migration timescale is shorter than the protoplanet’s accretion timescale (just as in the former case) for the critical mass interval of a few tens to about one hundred of Earth masses. However, in this case accretion proceeds fast enough for the protoplanet to start runaway accretion of gas before reaching the inner edge of the disc. Therefore the planet is able to accrete mass very fast (gas accretion dominates the growth of the planet), which means that it becomes massive enough to enter type II migration. The migration rate (in type II with gap opening) decreases as the protoplanet grows in mass. Therefore, being in the runaway phase is what saves the planet from falling into the star.

The situation can be summarised as follows: when the growth of the planet is dominated by the accretion of solids in the oligarchic regime (before gas runaway accretion), the growth timescale is proportional to M^1/3 (Ida & Makino 1993). However, for planets massive enough to trigger rapid gas accretion, the growth timescale is much shorter. On the other hand, typical migration time scales with M^-1 in type I and is independent of the planet’s mass in type II, being generally much longer than in type I. As a consequence, if a planet migrates in type I and is dominated by the oligarchic growth, it is very likely to be lost in the star. It is only if the planet succeeds to become massive enough to start the runaway gas accretion, and quasi-simultaneously enter type II migration, that it can brake before reaching the central star. Therefore, there is a critical mass range between  ~10−100 M_⊕: a planet in this mass range is likely to be undergoing inward type I migration and a decreasing growth rate, which at this stage is dominated by the accretion of solids. Indeed, because it is massive, the protoplanet excites the random velocities of planetesimals, making its accretion difficult. Therefore the growth of the planet is very slow at this stage, because planetesimals cannot be effectively accreted. Although this is the pathway to the runaway accretion of gas, the transition between being dominated by the accretion of solids and the accretion of gas, even though it is fast, is not immediate. On the other hand, inward type I migration in this mass range is still very fast. Therefore, if planets enter the runaway gas phase, they would probably grow massive, change to type II migration, and avoid a fatal end in the central star. If they remain very small, type I migration is not very harmful. Finally, if they grow to about Neptune mass while the disc is still young, but do not manage to grow fast enough, their migration accelerates and they end in the central star.

4.4. Exploring the parameter space of initial conditions

One may wonder if the examples presented in the previous section are representative of a general case. To answer this question we applied our planet formation model to a variety of protoplanetary discs (including different lifetimes, masses, metallicities, etc.) and initial locations for the embryo.

In this framework, we performed sets of more than 10 000 simulations. Each set considers a different size for the accreted planetesimals (r_m =  100, 10, 1 and 0.1 km). For these simulations we adopted the disc models described in Sect. 2.1.2. The initial mass of the embryo is, in all the cases, 0.01 M_⊕. The initial location of the embryo is varied between 0.2 and 30 AU. For an embryo to start at a certain location, we checked that the initial mass of solids at its location is greater than the mass of the embryo. We subtracted the mass of the embryo from the initial mass of solids of its feeding zone, which means we are assuming instantaneous, in situ formation for it. Planets grow by accreting solids and gas, and can migrate in the disc, according to the model we presented in the previous sections. The evolution of the eccentricity and inclination of planetesimals are calculated solving the differential equations (Eqs. (31), (32)) all throughout the disc for every time step. As initial conditions for e and i we assume that they are given by Eqs. (58), (59), what we called before an initially cold disc.

The parameter space explored for the initial mass and lifetime of protoplanetary discs is schematically shown in Fig. 9. The rectangle surrounded by a black box shows the whole set of initial conditions studied. Note however that the distribution of these parameters is likely not uniform, and not all combinations have the same probability of occurrence: long-lived discs as well as very low and very massive discs are unlikely. We found that giant planets⁷ do not form for the initial conditions corresponding to the grey region in Fig. 9. In fact, for a population of 100 km and 10 km planetesimals, giant planets do not form under any of the initial conditions considered in this work. For the 100 km population of planetesimals this is because the accretion time to form a massive solid core is longer than the discs’ lifetimes. For 10 km planetesimals, some giant planets would form if migration were not at work (see Sect. 5). However, according to disc-planet interaction models, planets do migrate and the migration timescale turns out to be much shorter than the accretion timescale. Planets are lost in the central star before they are able to accrete enough solids to trigger the runaway accretion of gas (as seen in the previous section, this would allow them to grow faster and switch from type I to type II migration, therefore preventing their loss in the central star).

However, the situation reverts when the accreted planetesimals are smaller. Figure 9 shows that massive discs favour the formation of giant planets. The coloured regions depict the characteristics of discs where, for a particular initial position, the embryo succeeded in growing to a giant planet. This does not mean that a giant planet will form in any location of the disc, but that formation is possible at certain locations. The yellow region shows the parameter space of giant planets that succeeded to survive in the disc accreting 1 km planetesimals. Clearly, there is a dependence on the discs lifetime: less massive but longer-lived discs favour the formation of giant planets. In red is plotted the situation of r_m = 0.1 km accreted planetesimals. Because smaller planetesimals are accreted faster and more efficiently, the disc parameter space that leads to the formation of giant planets is bigger. Interestingly, in this case there is no dependence on the lifetime of the disc. When the amount of solids present is large enough, accretion to form a core with a critical mass proceeds so fast that it is always shorter than the lifetime of the discs considered here. If the planetesimal disc is massive, accretion can be fast enough to enable protoplanets to reach the cross-over mass. Massive protoplanets switch from fast type I migration to a much slower type II, therefore decelerating and eventually braking before reaching the central star, preventing this way the loss of planets.

Fig. 9 Lifetime vs. initial solid mass of the protoplanetary discs considered in our simulations. The total mass of the disc relates to the solid mass of the disc through the gas-dust ratio. The rectangle surrounded by a black box schematically shows the parameter space explored. The grey rectangle depicts the region where giant planets do not form at all. The region of the parameter space where protoplanets reached the cross-over mass are shown in yellow for protoplanets accreting 1 km planetesimals, and in red for protoplanets accreting 0.1 km planetesimals. For bigger accreted planetesimal (100 and 10 km) no giant planets form.

Very small mass planets are the most abundant outcome in all simulations (regardless of the size of the accreted planetesimals) and their final location could be anywhere in the disc. Planets more massive than 10 M_⊕, in general, start their formation beyond the ice line (where solids are more abundant and feeding zones are larger), and due to migration they reach the inner regions of the disc. The final location of these planets extends between the inner edge of the disc and 10 AU. Planets in the range of 10² to 10³ M_⊕ are less abundant, and their final locations are very different from their initial emplacements. These planets undergo much migration, and those that remain are the ones that were able to accrete gas fast enough to enter the runaway accretion of gas, which prevented their loss in the central star. Planets in the mass range 10−10² M_⊕ are those that undergo the largest net displacement from their original location. Most of the surviving planets in this mass range have masses lower than 20 M_⊕. These planets, which were not able to reach their cross-over masses while they were migrating towards the central star, are preserved because the disc dissipated before they could fall into the star. In the case of r_m = 0.1 km,  ~18% of the simulations lead to lost planets, 85% of which were not able to reach their cross-over mass. This confirms our analysis of the previous section. Planets that are massive enough to undergo rapid inward type I migration (≳10 M_⊕) but whose growth rate is still dominated by the accretion of solids are likely to be lost in the star.

In the previous sections we have noticed that results of giant planet formation depend upon the dynamical model adopted to describe planetesimals’ dynamics. The equilibrium approach has the advantage of being numerically not very time-consuming. However, it can lead to very different results when compared to the case of explicitly solving the differential equations for e and i. When solving the equations, we also have the problem of setting the initial conditions, which are unknown. To test the importance of these assumptions, we performed 10   000 simulations under three different conditions. Each calculation differs from the other only in the treatment of the planetesimals e and i. This means that for a given set of initial conditions, we calculate the formation of the planet three times: once assuming the full equilibrium situation for planetesimals, another solving the explicit differential equations using as initial conditions for e and i the equilibrium value when the stirring timescale of the embryo equals the nebular gas drag timescale (we call this scenario hot disc) and finally also solving the differential equations but using as initial conditions the equilibrium value between mutual planetesimal stirring and gas drag damping (we call this scenario cold disc).

Figure 10 shows the fraction of surviving planets with respect to the total amount of planets formed (surviving planets plus planets lost in the central star). Accreted planetesimals are 0.1 km size. To make this plot we classified the planets according to their final mass into five mass bins. Clearly, in all the cases, the most abundant planets are the low-mass planets (<1 M_⊕). The equilibrium scenario being the one corresponding to the slowest accretion rate, is the one with more planets in the lowest mass bin (<1 M_⊕). At the other extreme, in the case of the cold initial disc, accretion of solids at the beginning is more efficient, embryos grow bigger in a shorter time which, in turn, gives them more chances to continue accreting. That is why in the case of an initial cold disc planets are more massive. This is evident in the histogram: for a cold initial planetesimal disc, the fraction of planets in each mass bin is the highest (except, of course, in lowest mass bin). The amount of planets in the interval of 10 to 10² M_⊕ represents around 2% of the surviving planets. Most of them are in the mass range of 10 to 20 M_⊕, and just a few in the range of 20 to 50 M_⊕. There are no planets in the range of 50 to 10² M_⊕, which shows the dramatic effect of type I migration on these intermediate mass planets. In the mass interval 10² − 10³ M_⊕ a profound decrease in the number of planets can be noticed, independently of the approach used for the planetesimal dynamics. This can be understood as follows: protoplanets with masses greater than 10² M_⊕ are usually in the runaway phase of gas accretion. This means that hundreds of Earth masses can be accreted in a very short time. Therefore planets in the runaway phase easily grow more massive than Jupiter. Only if the disc dissipates during this process of accretion, final planetary masses can be between Saturn and a few Jupiter masses. To a lesser extent, planet migration has also some consequences for depleting this bin: although the mass bin in which planet migration is most effective in eliminating planets is the 10 − 10² M_⊕ range, it is still very important in this mass regime, eliminating around 25% of these planets (especially those closer to the 10² M_⊕ edge of the mass interval).

Fig. 10 Fraction of surviving planets with final mass in a certain mass interval. Ntotal is the sum of all surviving planets (nsurv) plus the planets lost in the central star. Final masses of the formed planets are separated into five mass bins (M < 1   M⊕, 1 < M [M⊕]  < 10, 10 < M [M⊕]  < 100, 100 < M [M⊕]  < 1000, M > 1000   M⊕). The solid component of the planets grows by the accretion of 0.1 km planetesimals. Three different approaches in the treatment of the evolution of the eccentricities and inclinations of planetesimals are considered: equilibrium, and explicit calculation of the differential equations for a “hot” and “cold” initial conditions.

Fig. 11 In situ formation. Fraction of planets with final mass in a certain mass interval. Ntotal represents the total amount of planets formed, independently of their mass. Planets are separated into five mass bins (M < 1   M⊕, 1 < M [M⊕]  < 10, 10 < M [M⊕]  < 100, 100 < M [M⊕]  < 1000, M > 1000   M⊕). Accreted planetesimals have a radius of 100 km (red), 10 km (green), 1 km (blue) and 0.1 km (magenta). Most of the final masses are in the first mass bin. The inset caption represents the magnified histogram for the other mass bins.

5. Discussion

The calculations presented in this paper focus on the formation of planets and on the high degree to which this process is affected by the accretion rate of solids. The accretion rate of solids introduced intends to be realistic while computationally tractable. Some other aspects of the model are therefore rather simplified here; for example, we considered the formation of a single planet. The application of these calculations to the formation of planetary systems will be presented elsewhere (Alibert et al. 2012). We introduced a semi-analytical description for the eccentricities and inclinations of planetesimals. In this work we explicitly calculated the planetesimals’ eccentricities and inclinations, taking into account the stirring of the growing planets, the gas drag from the nebula, and the mutual stirring of the planetesimals themselves. The stirring produced by the growing planet excites planetesimals and makes their accretion more difficult as it grows. On the other hand, gas drag counteracts the stirring effect, an effect more important for smaller planetesimals. We considered three approaches to determine the rms eccentricity and inclination of planetesimals:

• analytical equilibrium calculation, as described byEqs. (56), (57);

• out-of-equilibrium, by solving the time evolution of the rms eccentricity and inclination of planetesimals (Eqs. (31), (32)); starting from a cold planetesimal disc (their excitation state being the results of planetesimal-planetesimal interaction and gas drag only); and

• out-of-equilibrium, starting from a hot disc, where planetesimals are already excited by the planetary embryo.

We have shown that these three approaches lead to different accretion rates, the difference depending on the planetesimal size, and being more important as a result of the migration and disc evolution feedbacks. For large planetesimals, the three approaches lead to similar results, but the accretion rate of solids under this hypothesis is very low, preventing the formation of massive planets. On the other hand, for low-mass planetesimals, the excitation state of planetesimals as derived from the analytical equilibrium calculations, and the excitation state of planetesimals following the second or third approach, even after a long time (when the equilibrium solution of Eqs. (31), (32) is reached), are different, resulting in different accretion rates of solids. In particular, the ratio of the eccentricity to the inclination can deviate substantially from the 1/2 ratio assumed in the first approach.

As a consequence of the size dependence of gas drag, small planetesimals are easier to accrete, leading to a faster formation of planets. We have shown that in the framework of the model hypothesis outlined in Sect. 2, the formation of planets ranging from a fraction of an Earth mass to several Jupiter masses can only be accomplished under the assumption of a population of small planetesimals or, more generally, if planetesimals are, by any process, maintained in a cold dynamical state. Similar conclusions were identified in other works (e.g. Fortier et al. 2007, 2009; Benvenuto et al. 2009), where it was shown that to achieve the formation of the giant planets of the solar system in less than 10 Myr most of the accreted planetesimals have to be small. Those models, however, assumed in situ formation, and did not consider a consistent calculation of the planet’s final mass (the computations were stopped when the masses of the giant planets of the solar system were reached). Although the approach to calculating the planet formation (regarding the accretion of gas and solids) is similar, our work accounts for the migration of the planets and their final masses are determined by the coupled evolution of the disc and the planet.

To compare our results with these previous works, we performed simulations where migration is switched off. We calculated the in situ formation of 10   000 planets that accrete planetesimals of 100, 10, 1 and 0.1 km in radius. Indeed, as can be seen from Fig.  11, small planets are the most abundant ones. Larger planets form when there are small planetesimals to accrete (smaller than r_m = 10 km). Planets in the mass range 10 − 10²   M_⊕ (whose growth is not yet dominated by the accretion of gas) are produced in about the same fraction as giant planets in the mass range  > 10³   M_⊕, whose growth is dominated by the a runaway accretion of gas. The runaway gas accretion leads to the difficulty of forming planets in the mass range of 10² − 10³   M_⊕: only a fine-tuned timing effect, the gas disc disappearing when the planet is in this mass range, can lead to planets in the Saturn-Jupiter domain.

If we focus on the mass bin 10 to 10²   M_⊕, one difference between the migration and in situ scenario is that, while in the former there are no planets in the mass range between 50 and 100 M_⊕ (because they are all engulfed by the central star), in the in situ case planets are distributed in the whole mass spectrum of the bin (however, most of them are in the lowest mass range [10 to 20 M_⊕]). It is also interesting to compare the mass percentage of solids for planets between 5 − 10²   M_⊕ for these two scenarios. As can be seen in Table 3, planets that migrate have a more massive core than planets that formed in situ. In general, when planets migrate, they have access to regions of the disc that have not been depleted of planetesimals. Therefore, the solids surface density is higher and so is the accretion rate, which results in the formation of a more massive core.

We considered that the opacity of the envelope corresponds to the full interstellar medium opacity. However, several works (e.g. Pollack et al. 1996; Podolak 2003; Hubickyj et al. 2005; Movshovitz & Podolak 2008) suggest that opacity in the planet’s envelope should be much lower, leading to a faster formation. As a complementary result, we computed one in situ population and one population including migration (10 000 planets each) reducing the opacity to 2% of that of the interstellar medium. We found that with a reduced opacity the loss of planets decreases, increasing the number of planets that survive in the protoplanetary disc in every mass bin. Nevertheless, the shape of the mass distribution is very similar to that shown in Fig. 10: the fraction of surviving planets decreases with the mass of the planet, with a minimum in the interval 10² − 10³M_⊕ (although with the reduced opacity the fraction of planets in this bin is four times higher than with the full opacity), and rising again in the last mass interval. In Table 4 we show the solids mass fraction for planets with masses between 5 to 100 M_⊕, considering both the in situ and a migration scenario. In both cases, the fraction of solids is smaller than in the corresponding case calculated with full opacity.

Using the full model (including migration and disc evolution), giant planet formation by accretion of 100 km planetesimals is quite unlikely, if not impossible. In our simulations, to actually form giant planets we had to reduce the planetesimal size to the 0.1 − 1 km radius range. Such a conclusion raises the question of the most likely planetesimal size during the epoch of planet formation. Recent studies on planetesimal formation give different conclusions, however. On one hand, models that explain the formation of planetesimals by direct collapse in vortices in turbulent regions (e.g. Johansen et al. 2007) predict a fast formation of very big planetesimals (r_m > 100 km). We note, however, that this formation process may not be totally efficient and only a small fraction of solids initially present in protoplanetary discs is likely to end up in such big planetesimals. The conclusions of Johansen et al. (2007) are consistent with the results of Morbidelli et al. (2009) on the initial function of planetesimals. On the other hand, a recent study of Windmark et al. (2012) shows that direct growth of planetesimals via dust collisions can lead to the growth of 0.1 km planetesimals. In addition, initially small planetesimals show better matches to the observed size distribution of objects in the asteroid belt and among the trans-Neptunian objects (TNOs): Weidenschilling (2011) showed that the size distribution currently observed in the asteroid belt in the range of 10 to 100 km can be better explained by an initial population of 0.1 km planetesimals. Kenyon & Bromley (2012) concluded, by combining observations of the hot and cold populations of TNOs with time constraints on their formation process, that TNOs form from a massive disc mainly composed of 1 km planetesimals. More investigations on the formation of planetesimals, and planetary embryos, are definitely required to test the viability of planetary core formation by accretion of low-mass planetesimals. We note finally that what is important in the work we have presented now is not the initial mass function of planetesimals, but their mass function at the time of planet formation. The two quantities are likely to differ due to planetesimal-planetesimal collisions and the resulting mass growth and/or fragmentation.

6. Conclusions

We presented calculations of planetary formation considering the formation of a single planet at each time, and starting with embryos of 0.01   M_⊕. Our simulations consist of calculating the formation of a planet, including its growth in mass by accretion of solids and gas, its migration in the disc, and the evolution of the disc until the gas component of the disc is dispersed. When the nebular gas is gone, simulations were stopped. Therefore, the subsequent growth of the planets, by accretion of residual planetesimals or collisions among embryos if we were considering planetary systems, was not considered. During their formation, the growth of the planets was calculated self-consistently coupling the accretion of solids and gas. The accretion of solids was computed assuming the particle-in-a-box approximation, and computing the excitation state of planetesimals, which in turn regulates to a great extent the accretion rate of the planet. The accretion of gas was computed by solving the differential equations that govern the evolution of the structure of the planet. Finally, protoplanets grow in an evolving protoplanetary disc, whose density, temperature and pressure was calculated at every time step.

The combination of oligarchic growth (for the solid component of the planet) with the migration of the planet has severe consequences for protoplanets that are able to grow up to a few tens Earth masses: these planets tend to collide with the central star (or at least to migrate to the innermost location of the protoplanetary disc). Indeed, planets that are between 10 and 100 Earth masses are usually undergoing very rapid inward type I migration, but are not massive enough to switch to a slower, type II migration. In our simulations, the only surviving planets in the range of 10 to 20 Earth masses correspond to cases where the gas component of the disc dissipates during their growth, preventing them from falling into the star. In contrast, if the solid core grows fast enough, it enables the accretion of large amounts of gas when the critical mass is reached. At this point, the runaway of gas ensures an extremely quick growth in the mass of the planet, and the planet migration rate decreases.

In the model we presented in this paper, we assumed a uniform population of small planetesimals whose size remains unchanged during the whole formation of the planet. Most probably the initial population of planetesimals in protoplanetary discs is not uniform in size, but follows a size distribution. We have shown, however, that without small planetesimals giant planet formation is difficult to explain, at least in the way we understand it now. However, even with an initial population of small planetesimals, the collisions among themselves are likely to be disruptive as soon as their random velocities start to be excited by a planetary embryo.Therefore, it is also unlikely that an initial population of only small planetesimals can be used to explain the formation of giant planets. Moreover, even assuming this to be true, only a few planets in the range of several tens of Earth masses to a few Jupiter masses can be formed.

In addition, small-mass planetesimals are subject to strong radial drift as a result of gas drag. Planetesimal drift can have positive or negative consequences in the formation of planetary systems, as has been shown by Guilera et al. (2010, 2011). Similarly, fragmentation and coagulation can hasten or delay planet formation as a whole: fragments of smaller mass are easier to accrete but they can also leave a planet’s feeding zone very quickly as a result of gas drag. Finally, in a planetary system, fragments that are not accreted by the embryo that generated them can be accreted by another embryo located in an innermost region. It is not clear what the possible outcomes of adding planetesimals drifting, fragmentation, and many embryos forming in an evolving disc could be. This represents a very important step in the understanding of the first stage of planet formation, however.

Because most of the accretion of solids should, at some point, be dominated by small planetesimals or fragments, our calculations can be understood as a description of that stage. An interesting scenario to analyse, in particular if planetesimals are born massive (Johansen et al. 2007; Morbidelli et al. 2009), would be the following. An initial population dominated by  ~100 km planetesimals would prevent a fast growth of the embryos at the beginning, a time during which the planetary embryos would suffer only of a little migration and the protoplanetary disc would evolve, reducing its gas surface density. Fragments (smaller planetesimals), resulting from collisions between big planetesimals, would start to be generated later (the timescale of fragmentation of 100 km planetesimals affected by the stirring of Moon to Mars mass embryos is of the order 10⁶ years), therefore accelerating the formation of the embryo in a later stage of the disc evolution. The collisional cascade would probably still produce small fragments fast enough to help the growth of an embryo, even if they leave the feeding zone very fast. Therefore, protoplanets could grow by the accretion of fragments, not necessarily generated by themselves, but generated by another distant embryo. Adding these different processes will give us a better insight into the formation process and would help us to constrain, from planetary formation models, possible initial size distributions for planetesimals.

It is also important to mention that we have not considered the possibility of planetesimal-driven migration. Although in general planetesimal-driven migration acts on a longer timescale than type I migration, recently Ormel et al. (2012) found that planetesimal-driven migration can have a mild effect on mid-sized planets in massive planetesimal discs, competing with type I migration.

The main conclusions of our work are that formation of giant planets in the framework of the sequential accretion model needs the presence of unexcited planetesimals. One obvious way to de-excite planetesimals is through gas drag, but this requires this latter to be efficient, which in turn translates into low-mass planetesimals. These planetesimals can be primordial or fragments of originally bigger ones. But, at some point, small boulders are needed to build protoplanetary massive cores before the dissipation of the disc. The combination of migration and oligarchic growth, on the other hand, prevents the formation of

intermediate mass planets. However, this result can change when considering the formation of planets in planetary systems where their gravitational interactions are taken into account. Captures in resonances can prevent planets from colliding with the central star, preserving them in planetary systems. Exploring these effects will be the subject of future works.

Acknowledgments

This work was supported by the European Research Council under grant 239605, the Swiss National Science Foundation. We thank C. Mordasini for the fruitful discussions. We also wish to thank the referee, J. Lissauer, for constructive criticisms which enabled us to improve the original version of this article.

References

All Tables

All Figures

	Fig. 1 Disc model. Gas surface density as a function of radius at different times. For the initial model (red line) the parameters are those of the first row of Table 1, the disc lifetime being 3 Myr. Each line corresponds to a different time (from top to bottom 0 to 3 Myr, every 0.1 Myr).
In the text

Fig. 2 Left panel: the ratio between eccentricities (red) and inclinations (blue) as a function of the planet’s mass considering both the oligarchic and the quasi-runaway growth. Clearly, the planetesimals’ excitation is several times higher in the oligarchic growth than in the quasi-runaway growth. The right panel shows the impact of the planetesimals’ excitation in the accretion rate of solids: for the quasi-runaway growth (yellow) the accretion rate is much higher than in the oligarchic growth (green).

In the text

Fig. 3 In situ formation of a 1 M⊕ planet at 6 AU for different radii of accreted planetesimals (red, 100 km; blue, 10 km; green, 1 km; orange, 0.1 km). The solid line is for the equilibrium approximation of eccentricities and inclinations (see Eqs. (56), (57)). The dashed and dotted lines are for the explicit calculation of e and i solving the differential equations Eqs. (31), (32). Two different initial values, e0, i0, were considered: the dashed line represents the case where e0 and i0 are given by Eqs. (59), (58), while the dotted line assumes e0, i0 to be the initial equilibrium values of the equilibrium case. The top-left panel depicts the mass growth of the protoplanet as a function of time. The top-right panel shows the time evolution of the eccentricity. The bottom-left panel shows the eccentricity as a function of the mass of the embryo. Clearly, in all cases, equilibrium is attained after an out of equilibrium state. The bottom-right panel plots the accretion rate of solids as a function of the embryo mass. Note that, as a function of the planet mass, the difference in the accretion rate is only a consequence of the difference in e and i.

In the text

	Fig. 4 Inclination-eccentricity ratio for the out-of-equilibrium cases of Fig. 3. Red is for 100 km planetesimals, blue for 10 km, green for 1 km, and orange for 0.1 km. The dotted line represents the case of an initially hot disc, while the dashed line represents an initially cold one. The grey line is the standard equilibrium value, βeq = i/e = 1/2.
In the text

	Fig. 5 In situ planetary formation at 6 AU for different radii of the accreted planetesimals (red, 100 km; blue, 10 km; green, 1 km; orange, 0.1 km). The same line code as in Fig. 3 is adopted. The formation of the planet ends when the gas component of the disc is dispersed.
In the text

Fig. 6 Full planetary formation, now allowing for the migration of the embryos, for different radii of the accreted planetesimals (red, 100 km; blue, 10 km; green, 1 km; orange, 0.1 km). The same line code as in Fig. 3 is adopted. Embryos are initially located at 6 AU. The left panel shows the cumulative mass of the embryo (solids plus gas) as a function of time (in logarithmic scale). The right panel depicts the migration path of the protoplanet (in linear scale). For the accretion of 1 km planetesimals, the embryos are lost in the Sun (indicated with crosses in the plot). The same happens for the planet that grows by accretion of 0.1 km planetesimals in the equilibrium approximation (orange solid line). When planets are lost, calculations are stopped.

In the text

	Fig. 7 Mass versus semi-major axis for the cases shown in Fig. 6. Note that the accretion rate of solids plays a major role not only in the growth of the planet but also in its migration path.
In the text

Fig. 8 Protoplanet’s growth and migration timescales (τgrowth = M/ | Ṁ |  and τmig = a/ | ȧ |  respectively) as a function of the protoplanet’s mass for the four different sizes of the accreted planetesimals of Fig. 6 (corresponding to the dashed line in Fig. 6). The peaks in the migration timescale correspond to the changes in the sense of migration of the protoplanet. From the two bottom panels, it can be seen that when the protoplanet’s mass is  ≳10 M⊕, both timescales become comparable and eventually the migration timescale becomes shorter than the accretion timescale, leading to a fast migration of the protoplanet with very little accretion. For 1 km planetesimals the protoplanet growth is too slow to gain enough mass to change the type of migration before they are lost in the star (when its mass is  ~30 M⊕). For 0.1 km planetesimals, the protoplanet can grow more massive because its accretion rate is higher than in the previous case, then it is able to open a gap in the disc and migration switches to type II, preventing it from falling in the star.

In the text

Fig. 9 Lifetime vs. initial solid mass of the protoplanetary discs considered in our simulations. The total mass of the disc relates to the solid mass of the disc through the gas-dust ratio. The rectangle surrounded by a black box schematically shows the parameter space explored. The grey rectangle depicts the region where giant planets do not form at all. The region of the parameter space where protoplanets reached the cross-over mass are shown in yellow for protoplanets accreting 1 km planetesimals, and in red for protoplanets accreting 0.1 km planetesimals. For bigger accreted planetesimal (100 and 10 km) no giant planets form.

In the text

Fig. 10 Fraction of surviving planets with final mass in a certain mass interval. Ntotal is the sum of all surviving planets (nsurv) plus the planets lost in the central star. Final masses of the formed planets are separated into five mass bins (M < 1   M⊕, 1 < M [M⊕]  < 10, 10 < M [M⊕]  < 100, 100 < M [M⊕]  < 1000, M > 1000   M⊕). The solid component of the planets grows by the accretion of 0.1 km planetesimals. Three different approaches in the treatment of the evolution of the eccentricities and inclinations of planetesimals are considered: equilibrium, and explicit calculation of the differential equations for a “hot” and “cold” initial conditions.

In the text

Fig. 11 In situ formation. Fraction of planets with final mass in a certain mass interval. Ntotal represents the total amount of planets formed, independently of their mass. Planets are separated into five mass bins (M < 1   M⊕, 1 < M [M⊕]  < 10, 10 < M [M⊕]  < 100, 100 < M [M⊕]  < 1000, M > 1000   M⊕). Accreted planetesimals have a radius of 100 km (red), 10 km (green), 1 km (blue) and 0.1 km (magenta). Most of the final masses are in the first mass bin. The inset caption represents the magnified histogram for the other mass bins.

In the text