© ESO, 2012

1. Introduction

The number of known transiting extrasolar planets or planet candidates has recently increased exponentially, thanks both to ground-based observations (e.g., Gillon et al. 2011) and space missions like the COnvection ROtation and planetary Transits satellite (CoRoT; e.g., Léger et al. 2009) and Kepler (Borucki et al. 2011). By combining the radius measurements with radial velocity measurements, which yield the mass of the planet, one gets the planetary mass-radius (M-R) diagram, an observational result of similar importance to the semimajor axis-mass (a-M) diagram.

The latter relation is available owing to the success of ongoing radial velocity surveys (e.g., Mayor et al. 2011). It is a goal of population synthesis models to understand the structure of the a-M distribution because of the multitude of clues it contains for planet formation theory (e.g., Ida & Lin 2004; Mordasini et al. 2009a). A recent comparison of numerous theoretical and observational results obtained mostly by the radial velocity technique can be found in Alibert et al. (2011) and Mordasini et al. (2012a).

The M-R diagram which is now available for a statistically significant number of planets, plays an important role because it allows the mean density of the planet to be derived. This constrains the internal planetary structure, which is of central importance to understand the nature (Leconte et al. 2011) and the formation of the planet. The formation and evolution of the planetary M-R relationship is studied in the companion paper Mordasini et al. (2012b, hereafter Paper II; see also Mordasini et al. 2011c).

Besides transiting planets, the number of planets detected by direct imaging has also increased significantly in the past few years, even though in absolute numbers, far fewer such planets have been found to date. But these discoveries, like the planetary system around HR 8799 (Marois et al. 2008, 2010), have triggered numerous theoretical studies regarding the formation of these objects (e.g., Dodson-Robinson et al. 2009; Kratter et al. 2010). Two points about these planets are particularly interesting: their large semimajor axis a and the fact that we measure the luminosity L of young giant planets at some time t. Both quantities are important to understand the formation mechanism (e.g., Marley et al. 2007; Janson et al. 2011). In the near future, more capable instruments like the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE) at the VLT (Beuzit et al. 2007) or the Gemini Planet Imager (GPI) at Gemini South (McBride et al. 2011) and later the Exoplanet Imaging Camera and Spectrograph (EPICS) at the E-ELT (Kasper et al. 2008) will become operational. We can therefore expect the number of points that we can put in the t-L, the a-L, and (for cases with an independent, dynamical mass determination) the M-L diagram to increase from a handful at the moment to hundreds in a few years, similar to what has happened for the M-R diagram in the past few years.

This shift from an era of discovery to one of a first physical characterization of extrasolar planets by their radii and luminosities has profound implications for planet formation theory. Until recently, planet formation models studied mostly giant planets detected by radial velocity measurement; here only the minimal mass and some orbital elements were known. Now we are confronted with a multitude of different sets of observational data and constraints, each regarding primarily planets of different types: transiting, close-in planets, most of them with a small radius (as found by Kepler, see also Paper II); directly imaged, self-luminous, massive young giant planets at large distances from their host star; and a wide range in masses going from super-Earth to Jovian mass planets at distances varying from close to the star out to several AU found by radial velocity measurements (e.g., Mayor et al. 2011).

The final goal of planet formation theory to develop one single theory that explains the formation (and evolution) of all these very different planets is a challenging one; it will require a great deal of effort in the coming years. Nevertheless, the approach to bring together all these different observational data in a coherent way (which in itself is not a trivial task, cf. Wolfgang & Laughlin 2012) to constrain planet formation and evolution theories seems to be a promising route. In the end, one is not interested in a theory that can explain certain types of planets, but fails for others.

In this and the companion paper, we take a first step in this direction. We present multiple upgrades of our formation model (introduced first in Alibert et al. 2004). The most important addition is that we now couple formation in a self-consistent way with the subsequent evolution at constant mass once the protoplanetary disk is gone. With this approach, we can calculate directly from one single model not only the planet’s mass and semimajor axis, but also the planet’s main physical characteristics like the radius, luminosity, surface gravity, effective temperature, and composition in terms of iron, silicates, ices, and H₂/He.

These basic characteristics are calculated for a planet at any time during its “life”, starting as a tiny sub-Earth mass seed in the protoplanetary disk to a mature, billion-of-years-old planet. A direct coupling of the planet’s formation and its evolution is necessary, as it is well known that the formation has important, direct consequences for the evolution, in particular for the luminosity of giant planets at young ages (“cold” vs. “hot” start models, cf. Marley et al. 2007; Spiegel & Burrows 2012). With this model development, we can now compare our simulations directly with observational constraints coming from radial velocity (and microlensing) transits as well as direct imaging.

Other upgrades are the following: a detailed description for the rate at which gas is accreted by the planet in the disk-limited gas runaway accretion phase, an internal structure model for the solid (iron/silicate and possibly ice) part of the planet, the inclusion of radioactive decay for the luminosity of the core, a new initial profile for the gaseous disk, a new prescription for the photoevaporation of the disk (including external and internal photoevaporation), and finally a realistically low grain opacity for the gaseous envelope (presented in a dedicated work, Mordasini et al. 2012c).

We thus deal in this paper and in Paper II mostly with improvements of the physical description of one planet. In other papers, we addressed upgrades regarding the disk model (Fouchet et al. 2011), the migration of low-mass planets (Dittkrist et al., in prep.), and the effect of the concurrent formation of several embryos in one disk (Alibert et al., in prep.).

1.1. Organization of the paper

The organization of the paper is as follows: in Sect. 2 we give a short overview of the model. In Sect. 3, we describe the modifications to the computational module that describes the gaseous envelope structure of the planet, i.e., extending it to calculate the structure not only during the pre-runaway formation phase as in our previous models, but also during the gas runaway accretion/collapse phase and the subsequent evolutionary phase after the disk is gone. In Sect. 4 we address a related subject: how to calculate the gas accretion rate in the disk-limited regime, i.e., once gas runaway accretion of forming giant planets has started.

In Paper II, we present further upgrades to the planet module, namely, a realistic model for the density of the solid core of the planets and the inclusion of radiogenic heating in it. We also describe briefly some modifications to the protoplanetary disk model. Finally, we use the upgraded model in population synthesis calculations to study the formation and evolution of the planetary M-R diagram, the distribution of planetary radii, and the comparison with observational data.

In the remainder of this first paper, we show specific results obtained with the upgraded model. In Sect. 6 we study the coupled formation and evolution of a Jovian mass planet at 5.2 AU from the Sun. Many of the effects seen during this particular simulation are characteristic for the effects encountered during the formation and evolution of the planets in a general population synthesis calculation (Paper II). In Sects. 7 and 8, we discuss our results concerning the radii and luminosities of giant planets, putting special weight on a comparison with other, more complex models. The luminosity of young Jupiters is further addressed in a dedicated paper (Mordasini et al., in prep.).

2. Combined model of planet formation and evolution

The formation model used here relies on the core accretion paradigm, coupled to standard models of protoplanetary disk evolution and tidal migration of protoplanets.

The basic concept of core accretion is to follow the concurrent growth of an initially small solid core and its surrounding gaseous envelope, embedded in a disk of planetesimals and gas (Perri & Cameron 1974; Mizuno et al. 1978; Bodenheimer & Pollack 1986). Within the core accretion paradigm, giant planet formation happens as a two-step process: first a solid core with a critical mass (of order 10 M_⊕) must form, then the rapid accretion of a massive gaseous envelope sets in. This process is thought to take typically several million years. Two other processes in the protoplanetary disk happen on a similar timescale: the evolution of the protoplanetary disk itself and the orbital migration of the protoplanets due to angular momentum exchange with the gaseous disk. As described in Alibert et al. (2005a), we therefore coupled a classical core accretion model (very similar to Pollack et al. 1996) to models describing these two processes, using for the disk a standard α model and prescriptions for the type I and II migration of the protoplanets. While we use simple, 1D models for the individual processes, we consider the full coupling between them, which alone leads to complex formation scenarios (cf. Mordasini et al. 2009a).

2.1. Limitations: no accretion after disk dissipation, primordial H₂/He envelopes only

From its conception, our model mainly deals with the formation and evolution of relatively massive Neptunian and Jovian planets, even though the vast majority of planets that actually form in a population synthesis are low-mass, failed core (proto-terrestrial) planets (Mordasini et al. 2009b). We currently do not include the giant impact stage (after the gaseous disk has disappeared), during which terrestrial planets acquire their final mass. Therefore, our model is incomplete in the description of the formation of planets less massive than ~10 M_⊕ (see Mordasini et al. 2009a, for a discussion).

Regarding the long-term evolution, we only consider primordial H₂/He envelopes for which we assume that the mass remains constant after the protoplanetary disk is gone. Thus we neglect a possible evaporation of the primordial gaseous envelope as well as the outgassing of a secondary atmosphere. Such atmospheric mass loss for planets close to the host star will be considered in future work. Note that the minimal allowed semimajor axis for model planets is 0.1 AU at the moment, so that no simulated planets are exposed to the intense radiation field occurring on very tight orbits as, for example, CoRoT-7b (a = 0.017 AU), where atmospheric mass loss may be very important (Valencia et al. 2010).

An exact lower boundary in mass for the application of the model is difficult to specify, because we find (in agreement with Rogers et al. 2011) that planets of just a few Earth masses can accrete non-negligible amounts of H₂/He during the presence of the nebula, provided that the grain opacity in the envelope is low, which is probably the case (Movshovitz et al. 2010). In Paper II we study the formation and evolution of a close-in super-Earth (~4 M_⊕) planet with a tenuous 1% primordial H₂/He envelope. We find good agreement of the radius with the study of Rogers et al. (2011), showing that the evolution of this kind of low-mass planet can also be modeled.

We describe in this paper only improvements (or the inclusion of new physical effects) relative to the original model described in Alibert et al. (2005a), to which the reader is referred to for the remaining description.

3. Gaseous envelope calculation during the attached, detached, and evolutionary phase

We now turn to the most important improvement of the model, which is the ability to calculate planetary envelope structures during the entire formation and evolution of the planets. In previous models, we calculated the gaseous envelope only during the phase when the planet had a subcritical core mass less than the one needed to trigger gas runaway accretion (about 10 M_⊕). This is sufficient if one is interested in the mass of the planets only. In order to characterize the planetary structure (and thus to have R and L, too), we now also include the gas runaway phase (which is attained by planets becoming supercritical during the lifetime of the nebula) and the evolution at constant mass after the dissipation of the nebula. This phase is eventually attained by all planets.

The calculation of the gaseous envelope in all phases makes it possible to obtain the core mass more accurately because we can now calculate the capture radius for planetesimals R_capt in the gas runaway accretion phase. This capture radius is necessary to calculate the core accretion rate Ṁ_Z. In previous calculations, we had to extrapolate the capture radius as found in the pre-runaway phase into the runaway phase. The behavior of R_capt is discussed in Sect. 6.2.2.

3.1. Structure equations

The gas accretion rate of the planet (in the early phases before runaway, cf. Sect. 6.2.1) is obtained by solving the one-dimensional, hydrostatic planetary structure equations (e.g., Bodenheimer & Pollack 1986). These equations are similar to those for stellar interiors. The difference is that the energy release by stellar nuclear fusion is replaced by the heating due to impacting planetesimals. These impacts yield the dominant energy source during the early formation stage. The other equations are the standard equations (except for the assumption of a constant luminosity within the envelope, as later discussed) of mass conservation, hydrostatic equilibrium, and energy transfer (e.g., Broeg 2009): $\begin{array}{ccc}& \frac{\mathrm{d}\mathit{m}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathit{\rho } & \frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathit{Gm}}{{\mathit{r}}^{\mathrm{2}}}\mathit{\rho }\\ & \frac{\mathrm{d}\mathit{l}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{0} & \frac{\mathrm{d}\mathit{T}}{\mathrm{d}\mathit{r}}\mathrm{=}\frac{\mathit{T}}{\mathit{P}}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{r}}\mathrm{\nabla }\mathrm{\left(}\mathit{T,P}\mathrm{\right)}\mathit{.}\end{array}$In these equations, r is the radius as measured from the planetary center, m the mass inside r (including the core mass M_Z), l the luminosity at r, and ρ,P, and T are the gas density, pressure, and temperature. The quantity ∇(T,P) is given by $\mathrm{\nabla }\mathrm{\left(}\mathit{T,P}\mathrm{\right)}\mathrm{=}\frac{\mathrm{d}\hspace{0.17em}\ln \hspace{0.17em}\mathit{T}}{\mathrm{d}\hspace{0.17em}\ln \hspace{0.17em}\mathit{P}}\mathrm{=}\min \mathrm{\left(}{\mathrm{\nabla }}_{\mathrm{ad}}\mathit{,}{\mathrm{\nabla }}_{\mathrm{rad}}\mathrm{\right)}\mathit{,}$(3)where the radiative gradient (in radiative zones) in the diffusion approximation is given as ${\mathrm{\nabla }}_{\mathrm{rad}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{64}\mathit{\pi \sigma G}}\frac{\mathit{\kappa lP}}{{\mathit{T}}^{\mathrm{4}}\mathit{m}}\mathit{,}$(4)where κ denotes the opacity (given by Bell & Lin 1994; and Freedman et al. 2008) and σ is the Stefan-Boltzmann constant. The adiabatic gradient (in convective zones) ∇_ad is directly given by the equation of state SCvH (Saumon et al. 1995). We thus assume zero entropy gradient convection and use the Schwarzschild criterion to determine whether a layer is convectively unstable.

This means that our model builds on two interlinked traditional assumptions about the interiors of (giant) planets. While these assumptions are mainly made for simplicity of the models, they are not necessarily correct (Stevenson 1985; Leconte & Chabrier 2012). First, it is assumed, as just stated, that the interiors are adiabatic due to efficient large-scale convection. In this case, for Jupiter and Saturn at present time, the degree of super-adiabaticity (which is the fractional degree by which the actual temperature gradient is larger than the adiabatic gradient) is very small (10^-8−10^-9, Leconte & Chabrier 2012). There are, however, several situations where the assumption of an adiabatic interior may break down (Saumon & Guillot 2004). Additionally, at early stages of evolution, the planets could be going through regions in the surface gravity-effective temperature parameter space where the formation of H₂ leads to a nonlinear response of the atmospheric structure, which influences the degree of super-adiabaticity (Baraffe et al. 2002). We are currently working on including the mixing length theory (e.g., Kippenhahn & Weigert 1990) to better quantify this issue (see also Rafikov 2007).

The second traditional assumption is that the planet consists of a few well-separated regions, which are chemically homogeneous. In our model, we assume that the envelope consists purely of hydrogen and helium. In reality (Stevenson 1985; Chabrier & Baraffe 2007), it is possible that compositional gradients exist, originating, for example, from the dissolution of the core (Wilson & Militzer 2012) or the destruction of planetesimals in deeper layers (Mordasini et al. 2005). If a stabilizing compositional gradient is present, it is likely that the large-scale convection breaks up into many thin chemically homogeneous convective layers, which are separated by narrow, diffusive interfaces with large compositional gradients (Stevenson 1979). Due to the large number of diffusive interfaces, the efficiency of heat transport is strongly reduced, resulting in a super-adiabatic temperature gradient. As recently shown by Leconte & Chabrier (2012), such a semiconvective model is able to reproduce observational constraints coming from the gravitational moments and atmospheric composition of Jupiter and Saturn. This shows that an adiabatic interior cannot be taken for granted.

3.2. A new, simple method for the calculation of total luminosity and evolutionary sequences

In a planetary population synthesis model, the evolution of thousands of different planets is calculated, covering an extremely wide parameter range of planetary core, envelope masses, accretion rates, and background nebula conditions. As a result, we need a stable and rapid method for the numerical solution of these equations. We thus replaced the ordinary equation for dl/dm = −T∂S/∂t (where S is the entropy) by the assumption that l is constant within the envelope and that we can derive total luminosity L (including solid and gas accretion as well as contraction and release of internal heat) and its temporal evolution by total energy conservation arguments, an approach somewhat similar to Papaloizou & Nelson (2005) and Hartmann et al. (1997). We recall that −dE_tot/dt = L and that, in the hydrostatic case, the total energy (neglecting rotation) is given as $\begin{array}{ccc}{\mathit{E}}_{\mathrm{tot}}& & \mathrm{=}{\mathit{E}}_{\mathrm{grav}}\mathrm{+}{\mathit{E}}_{\mathrm{int}}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{M}}\frac{\mathit{Gm}}{\mathit{r}}\mathrm{d}\mathit{m}\mathrm{+}{\mathrm{\int }}_{{\mathit{M}}_{\mathit{z}}}^{\mathit{M}}\mathit{u}\mathrm{d}\mathit{m}\\ & & \mathrm{=}\mathrm{-}\mathit{\xi }\frac{\mathit{G}{\mathit{M}}^{\mathrm{2}}}{\mathrm{2}\mathit{R}}\mathit{,}\end{array}$where u is the specific internal energy, M the total mass, and R the total radius of the planet. The integration of the gravitational energy includes the core, for which we assume for simplicity a density that is constant as a function of r. This is strictly speaking not self-consistent, since we assume that the core is differentiated, see Paper II. Then, the binding energy of the core is given as $\mathrm{-}\mathrm{(}\mathrm{3}\mathit{/}\mathrm{5}\mathrm{)}\mathit{G}{\mathit{M}}_{\mathit{Z}}^{\mathrm{2}}\mathit{/}{\mathit{R}}_{\mathrm{core}}$. By differentiating this energy with respect to time, assuming a constant density, the resulting luminosity from the growth of the core is equal to the accretion luminosity of planetesimals falling onto the core with a velocity equal to 0 at infinity. On the other hand, the core does not contribute to the internal energy because we do not consider the thermal evolution of the core. In Eq. (6) we introduced a parameter ξ, which represents the distribution of mass within the planet and its internal energy content. This last formula is nothing else than a definition of ξ. For example, a fully convective, nearly isentropic star that can be approximated by a n = 3/2 polytrope would have ξ = 6/7 (Hartmann et al. 1997).

The quantity ξ can be found for any given structure at time t with the equations above. Then one can write $\begin{array}{ccc}\mathrm{-}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{\mathit{E}}_{\mathrm{tot}}& \mathrm{=}\mathit{L}\mathrm{=}{\mathit{L}}_{\mathit{M}}\mathrm{+}{\mathit{L}}_{\mathit{R}}\mathrm{+}{\mathit{L}}_{\mathit{\xi }}& \\ & \mathrm{=}\frac{\mathit{\xi GM}}{\mathit{R}}\mathit{Ṁ}\mathrm{-}\frac{\mathit{\xi G}{\mathit{M}}^{\mathrm{2}}}{\mathrm{2}{\mathit{R}}^{\mathrm{2}}}\mathit{Ṙ}\mathrm{+}\frac{\mathit{G}{\mathit{M}}^{\mathrm{2}}}{\mathrm{2}\mathit{R}}\mathit{\xi ̇}\mathit{,}& \end{array}$where Ṁ = Ṁ_Z + Ṁ_XY is the total accretion rate of solids and gas, and Ṙ is the rate of change of the total radius. This equation corresponds to a generalization of Eq. (6) in Hartmann et al. (1997). All quantities except  can readily be calculated at time t. We now set $\mathit{L}\mathrm{\simeq }\mathit{C}\mathrm{(}{\mathit{L}}_{\mathit{M}}\mathrm{+}{\mathit{L}}_{\mathit{R}}\mathrm{)}\mathit{.}$(9)The factor C corrects approximately for neglecting the L_ξ term, and is obtained by calculating a posteriori, the total energy in the new structure at t + dt, which gives the exact luminosity as L_ex = − [E_tot(t + dt) − E_tot(t)] /dt. By setting C = L_ex/L, one obtains the correction factor, so that exact energy conservation would have occurred. As an approximation, we then use this C for the next time step. Using this method, the estimated luminosity L and the actual luminosity L_ex generally agree very well, provided that dt is small enough¹. With this prescription, we can thus always calculate the total luminosity at t + dt, which is one of the necessary boundary conditions for the envelope calculations, allowing us in the end to construct evolutionary sequences. We will show below that this method leads to results in terms of luminosity and radius evolution that are in very good agreement with traditional calculations based on entropy, such as found in Burrows et al. (1997) and Baraffe et al. (2003).

This method automatically includes the (core) luminosity due to the accretion of planetesimal during the formation phase as well as the release of energy due to the contraction of the core at constant mass (its dominant contribution during the evolution phase for giant planets, Baraffe et al. 2008). To this luminosity we finally add radiogenic luminosity L_radio (see Paper II) to get the total intrinsic luminosity. By planetary luminosity, we always mean the luminosity emitted from the interior, without the contribution from absorbed stellar radiation, which makes up about 40% of the total flux for Jupiter today (Guillot & Gautier 2009).

We plan on further improving our scheme to estimate the luminosity at the next timestep, which is possible by considering the contribution from the core and the envelope separately.

3.3. Hydrostatic approximation

Equation (5) assumes that the planet is always in hydrostatic equilibrium, particularly during the rapid contraction phase at the moment the planet detaches from the nebula (see the left bottom panel in Fig. 2), which we shall call the collapse phase independent of its true nature. To check if this assumption is self-consistent, we performed an order of magnitude estimate. Using the temporal change of the total radius v(R) = dR/dt (which is of order −10 m/s during the fastest contraction), we assumed that the contraction of the layers in the interior is homologous. This allows us to estimate the total kinetic energy as ${\mathit{E}}_{\mathrm{kin}}\mathrm{=}{\mathrm{\int }}_{{\mathit{M}}_{\mathit{Z}}}^{\mathit{M}}\frac{\mathrm{1}}{\mathrm{2}}{\left(\mathit{v}\mathrm{\left(}\mathit{R}\mathrm{\right)}\frac{\mathit{r}}{\mathit{R}}\right)}^{\mathrm{2}}\mathrm{d}\mathit{m}\mathit{.}$(10)Comparing E_kin with the total energy, we find that E_kin is always many orders of magnitude smaller than the total energy during the entire formation and evolution of a giant planet. During the collapse, it is smaller by about six to eight orders of magnitude (and during the subsequent slow long-term contraction even by about 20 to 30 orders). This very large difference makes one think that the hydrostatic approximation is justified. However, our check is only looking at global quantities and, thus, would not catch if in certain layers of the planet hydrodynamic effects become locally important. In addition, one cannot exclude a priori that the usage of the hydrostatic approximation itself forces the evolution of the planet to behave hydrostatically.

The large security margin of six to eight orders of magnitudes combined with the fact that Bodenheimer & Pollack (1986) reached the same conclusion seems, however, to indicate that the evolution will not become a dynamical collapse but remain just a rapid gravitational contraction, validating our approach. Nevertheless we must note that opposite views do exist (Wuchterl 1991).

3.4. Simplification using dl/dr = 0

The method introduced above only yields the total luminosity at the surface and not l(r), which is necessary to solve the structure equations. We use the simplest setting, namely, that dl/dr = 0. This seems at first sight not to be a reasonable approximation as it (formally) means that the complete luminosity originates in the core, which is, of course, not the case except for the early phases. However, the solution of the structure equations is usually very insensitive to the specific form that l takes in the interior. This result can be understood because l enters into the remaining three structure equations only in radiative zones of the planet (see Eqs. (3) and (4)). In convective zones, l does not enter at all. This is only true for the assumption that convection is strictly adiabatic, which might not be the case, as discussed at the end of Sect. 3.1.

However, significant radiative zones (in terms of contained mass) only exist during the very early phases of formation of the planet (Bodenheimer & Pollack 1986). But then the dominant part of the luminosity is generated by the core’s accretion of planetesimals, so that L ≈ L_core, or in other words dl/dr ≈ 0. Consequently, our approximation is close to the exact solution. The core luminosity due to the accretion of planetesimals is given as ${\mathit{L}}_{\mathrm{core}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathit{Z}}}{{\mathit{R}}_{\mathrm{core}}}\mathit{Ṁ}\mathit{Z}\mathit{,}$(11)where M_Z is the mass of the core, R_core its radius, and Ṁ_Z is the accretion rate of planetesimals. This means that we use the sinking approximation (Pollack et al. 1996).

During the collapse phase and the long-term evolution when no planetesimals are accreted any more, l varies strongly inside the envelope, from a very small value (L_radio) at the core to a typically much larger total luminosity at the surface for giant planets due to the cooling and contraction of the envelope. However, in this phase the planet is nearly fully convective (e.g., Bodenheimer & Pollack 1986; Guillot 2005), with a radiative gradient dominant over the adiabatic one by orders of magnitude (Alibert et al. 2005b). As a result, the radial variation of l is irrelevant for the structure. A radiative layer still exists, but only near the very surface of the planet (e.g., Burrows et al. 1997), where it contains negligible amounts of matter (≲ 1%, Guillot & Gautier 2009). Therefore, the layer does not contribute significantly to the contractional luminosity, which comes from the deeper interior (Kippenhahn & Weigert 1990), so that dl/dr = 0 is again a good approximation in those parts where l matters. This also means that the Schwarzschild criterion, which uses l to decide if a layer is convective or not, remains valid.

3.4.1. Tests of the simplification

We performed a number of tests regarding the assumption that dl/dr = 0. We adopted different prescriptions for the form of l as a function of radius (or mass). First, we (arbitrarily) assumed a linear increase with m from l(R_core) = L_core at the core to the total luminosity at the surface.

For the second set of tests, we proceeded as follows. Knowing the planetary internal structure at two timesteps (t and t − dt), we computed the Lagrangian change in the specific entropy S(m,t). From this, we computed l(m,t) = −T(S(m,t) − S(m,t − dt))/dt, where T is the mean temperature for mass shell m between t − dt and t. Then, we computed the planetary internal structure at time t + dt by scaling this l(m,t) onto the new core and total mass, and the new L_core and total luminosity. This rescaled, variable l(m,t + dt) was then used to calculate the structure.

In all cases, there are only very small differences between these calculations and the simulations with a constant l. For example, the time needed to go to runaway gas accretion in a Pollack et al. (1996) J1 like calculation (see Appendix A) is about 2−5% longer in the dl/dr = 0 simulation compared to the variable l case. This is a much smaller difference than those introduced by, for example, the uncertainties regarding the grain opacity or the planetesimal size. For the long-term evolution phase, the differences are virtually vanishing, because the planets are almost fully convective. This large independence of the results on the form of l has also been noticed by others (Fortier 2011, pers. comm.).

Finally, we tested the impact of the simplified luminosity calculation during the attached and collapse phase by comparing it with simulations based on the solution of the full set of structure equations using the traditional entropy method, obtained by a completely independent model (Broeg & Benz 2012). As described in their paper, these authors solve the structure equations in the quasi-hydrostatic approximation on a self-adaptive 1-dimensional grid using the classical Henyey method (Kippenhahn & Weigert 1990). For a simulation similar to the J1 case discussed in Appendix A, we compared the accretion rates of gas and solids, the luminosity, and the radius of the forming planet. Also in this case, we found very good agreement.

The tests we have presented may seem complicated, and one may wonder why we do not simply solve the whole set of equations, replacing Eq. (2) by its standard form dl/dm = −T∂S/∂t. We do not do this because population synthesis requires the computation of thousands of formation models, starting from different initial condition. Such a computation is done automatically, using programs generating the initial conditions, then computing the formation model, and finally compiling the important results (e.g., mass, semimajor axis, radius, and luminosity of the formed planets). In this situation, it is important that only very few simulations not converge, since nonconverging simulations could introduce a bias in the results (e.g., if all simulations leading to massive planets were not to converge, one would predict only low-mass planets). We have therefore chosen not to use the standard Newton-Raphson to compute the planetary internal structure, but instead to use a shooting method, which converges extremely well (less than 0.3% of our simulations in the synthesis in Paper II has problems converging using this method). However, the shooting method is possible in practice only if one knows the luminosity profile a priori, which is the case assuming a constant luminosity or in the two sets of tests presented above.

3.5. Cold versus hot start conditions

Finally, one must take into account that part of the released gravitational energy of newly accreted material is not incorporated into the planetary structure but lost as L_acc in the accretion shock on the planet’s surface (Bodenheimer et al. 2000) or in the surrounding circumplanetary disk (Papaloizou & Nelson 2005). For matter falling in free fall onto the planet, the released luminosity at the surface of the planet is approximately given by ${\mathit{L}}_{\mathrm{acc}}\mathrm{=}\frac{\mathit{GM}}{\mathit{R}}\mathit{Ṁ}\mathit{XY}\mathit{,}$(12)where Ṁ_XY is the gas accretion rate. The same amount of energy is also lost if the matter spirals onto the planet through a circumplanetary disk, with the difference that some fraction (half of it for a Keplerian accretion disk) is already lost in the disk.

We can now consider two limiting scenarios (see Spiegel & Burrows 2012, for intermediate “warm start” simulations). First, as in Papaloizou & Nelson (2005), we assume that all energy liberated in the shock is radiated away from the planet and does not contribute to the luminosity inside the planet’s structure. This is the physically likely scenario (Stahler et al. 1980; Chabrier 2010, pers. comm.). The opposite extreme is to assume that no energy at all is radiated away at the shock. While physically unlikely, it is of interest to consider this case. In absence of radiative losses at the shock, our simulations become similar to the so-called “hot start” calculations. This is because there is no mechanism for these “hot start” models that can act as a sink of entropy, which is the central role of the radiative losses at the accretion shock. As Fortney et al. (2005) and Marley et al. (2007) demonstrated, the smaller initial entropy of models that assume radiative losses at the shock has very important consequences for the luminosity of young Jupiters (“cold start” models, see Sect. 8.2.2).

The luminosity within the planetary structure L_int for the two cases is given as ${\mathit{L}}_{\mathrm{int}}\mathrm{=}\{\begin{array}{c}\\ \mathit{L}\mathrm{-}{\mathit{L}}_{\mathrm{acc}}& \mathrm{“Coldstart”}\\ \mathit{L}& \mathrm{“Hotstart}\mathrm{\left(}\mathrm{accreting}\mathrm{\right)}\mathrm{”,}\end{array}$(13)while an observer would see the total luminosity L coming from both the intrinsic and accretional luminosity². We have added the suffix “accreting” to our “hot start” calculations to distinguish them from traditional “hot start” calculations (Burrows et al. 1997; Baraffe et al. 2003), where one starts with a fully formed planet (i.e., at the final mass) with some high, arbitrary entropy. It is clear that our treatment is a strong simplification of the exact physics of the accretion shock (e.g., Stahler et al. 1980). As already mentioned by Marley et al. (2007), three-dimensional, radiation-hydrodynamic calculations resolving the shock physics would be necessary to get more accurate results. We will show in a dedicated paper (Mordasini et al., in prep.) how our results for the luminosities depend on a number of parameters that should be indicative of general trends for the behavior of the luminosity of newly formed giant planets.

3.6. Possible issues for “hot” accretion and dl/dr = 0

The following concern occurs when “hot” accretion is simulated with the dl/dr = 0 simplification in the envelope. When the internal energy of the infalling gas (or a part of it) is incorporated into the growing object rather than radiated away, important modifications of the internal structure can occur, as illustrated by several papers studying the effects of accretion onto low-mass stars (Hartmann et al. 1997; Mercer-Smith et al. 1984; Prialnik & Livio 1985). These studies show that under some circumstances the heating of the surface layers by the rapid accretion of hot, high-entropy gas can reduce the temperature gradient to a degree that convection stops. This means that a radiative layer develops in the interior, even at considerable depth.

Prialnik & Livio (1985) studied in detail the reaction of an initially fully convective 0.2 M_⊙ star to the accretion of “cold” and “hot” gas at a wide range of accretion rates. For the case of “cold” accretion, they found that, independent of the accretion rate, the objects remain fully convective and contract as they grow, indicating that the dl/dr = 0 simplification is appropriate in this situation. Such a behavior of the radius is also found in our simulations of giant planet formation in the “cold” assumption presented in Sect. 8.2.2. If the gas keeps a fraction α_h of L_acc (Eq. (12)), different scenarios occur, depending on α_h, and the gas accretion rate. In this paper, we study the two limiting cases corresponding to α_h = 1 for completely “hot” and α_h = 0 for completely “cold” accretion (see Eq. (13)).

Prialnik & Livio (1985) find that if the gas accretion rate is very small (Ṁ_XY ≲ 10^-10 M_⊙/yr ≈ 3 × 10^-5 M_⊕/yr), the star also remains fully convective for “hot” accretion. If both the accretion rate and α_h are high, a significant part of the entire star becomes radiative and the entire star expands by a large amount. For intermediate values of α_h and Ṁ_XY (including planetary accretion rates expected in a protoplanetary disk ≲ 10^-2 M_⊕/yr), ~10−30% of the mass of the star becomes radiative. At the same time, the outer radius, including the newly accreted matter, expands, while the matter originally contained in the star contracts. An expansion of the star undergoing “hot” accretion is also seen in the simulations of Mercer-Smith et al. (1984).

As discussed by Prialnik & Livio (1985; see also Hartmann et al. 1997), these results can be understood by comparing the accretion timescale of “hot” material with the thermal relaxation timescale of the newly accreted material. If the accretion timescale is short compared to the relaxation timescale, thermal equilibrium cannot be established, which would be needed for the star to remain convective. In the “hot start” simulations (Sect. 8.2.1), we also see that the radii of the planets reinflate during the gas runaway accretion phase, after the initial collapse. These results will be discussed in detail in Mordasini et al. (in prep.). It is clear that the presence of a deep radiative zone, which we cannot catch with the dl/dr = 0 simplification, could cause strong departures from the results obtained with this assumption. This would also affect the luminosity and radius of the planets after formation.

The results of Prialnik & Livio (1985) were obtained for spherically symmetric accretion. In reality, most of the matter would probably be accreted via a circumplanetary disk (e.g., Lubow et al. 1999), so that only a fraction of the planetary surface would be undergoing direct accretion, while the rest would be free to radiate. As for stars, it seems likely that this could be important for an accurate criterion, when convection is inhibited by “hot” accretion (Hartmann et al. 1997). Three-dimensional radiation hydrodynamic simulations seem necessary to clarify the situation. We will investigate this issue in simulations with a variable luminosity in future work. For the moment, this concern means that our results for the luminosity of young giant planets formed with a “hot start” (and especially at high gas accretion rates) must be regarded with caution.

3.7. Inflation mechanism

At the moment, we do not include any mechanism possibly responsible for the so-called “radius anomaly” (Leconte et al. 2011) which is observed for many transiting Hot Jupiters. These planets have radii clearly larger than expected from standard internal structure modeling as performed here. There are several possible physical mechanisms leading to this bloating, such as tidal heating (Bodenheimer et al. 2001), dissipation of stellar irradiation deep in the interior (Showman & Guillot 2002), ohmic dissipation (Batygin & Stevenson 2010), and double diffusive convection (Chabrier & Baraffe 2007). Not all these mechanisms are fully understood, and we leave their study for the moment to dedicated works (e.g., Leconte et al. 2011). But we note that population synthesis calculations, including the effects of various such mechanisms, could be a fruitful way to understand which one reproduces best the observed radius distribution (see Paper II for the predicted synthetic radius distribution for a ≥ 0.1 AU).

3.8. Boundary conditions

In order to solve the structure equations, boundary conditions must be specified, which should also provide a continuous transition between different phases. For the formation and evolution of a giant planet, three different fundamental phases must be distinguished: attached, detached, and evolutionary. Lower mass planets stay in the first phase until the nebula disappears and then directly pass into the evolutionary phase. Five boundary conditions are given, namely, core radius, luminosity, total planetary radius (in the attached phase) or total planetary mass (in the detached and evolutionary phase), and surface temperature and pressure. The four differential equations with five boundary conditions only have a solution for a given total planetary mass (attached phase) or radius (detached and evolutionary phase).

3.8.1. Attached (or nebular) phase

At low masses, the envelope of the protoplanet is attached continuously to the background nebula, and the conditions at the surface of the planet are the pressure P_neb and approximately the temperature T_neb in the surrounding disk. Total radius R is given in this regime by about the minimum of the Hill radius R_H and the Bondi (or accretion) radius R_A. The gas accretion rate is found by the solution of the structure equations and is given by the ability of the envelope to radiate away energy so that it can contract (i.e., its Kelvin-Helmholtz timescale), enabling new gas to stream in Pollack et al. (1996). Numerically, for the shooting method, one iterates in this phase on the mass of the envelope. We use $\begin{array}{ccc}\mathit{R}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{A}}}{\mathrm{1}\mathrm{+}{\mathit{R}}_{\mathrm{A}}\mathit{/}\mathrm{\left(}{\mathit{k}}_{\mathrm{liss}}{\mathit{R}}_{\mathrm{H}}\mathrm{\right)}}\mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{neb}}& & \\ \mathit{\tau }\mathrm{=}\max \left({\mathit{\rho }}_{\mathrm{neb}}{\mathit{\kappa }}_{\mathrm{neb}}\mathit{R,}\mathrm{2}\mathit{/}\mathrm{3}\right){\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathrm{=}\frac{\mathrm{3}\mathit{\tau }{\mathit{L}}_{\mathrm{int}}}{\mathrm{8}\mathit{\pi \sigma }{\mathit{R}}^{\mathrm{2}}}& & \\ {\mathit{T}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}_{\mathrm{neb}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathit{l}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{L}}_{\mathrm{int}}\mathit{.}& & \end{array}$The fitting formula for the radius of the planet in Eq. (14) from Bodenheimer et al. (2000) reduces to R_A for R_A ≪ k_lissR_H and to k_lissR_H for k_lissR_H ≪ R_A, where $\begin{array}{ccc}{\mathit{R}}_{\mathrm{A}}\mathrm{=}\frac{\mathit{GM}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}& {\mathit{R}}_{\mathrm{H}}\mathrm{=}{\mathrm{(}{\frac{\mathit{M}}{\mathrm{3} {\mathit{M}}_{\mathrm{\ast }}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{a}\mathit{.}& \end{array}$(17)The sound speed c_s in the equation for the Bondi radius is calculated with the background nebula conditions. As found by hydrodynamic simulations by Lissauer et al. (2009), only the inner part of the material inside the planet’s Roche lobe is permanently bound to the envelope of the planet, while the upper part participates in the general gas flow in the protoplanetary disk. We therefore follow these authors in setting k_liss = 1/3.

Regarding the temperature, Eq. (14) reflects the fact that the planet has to be hotter in order to radiate into the nebula (Papaloizou & Terquem 1999). Equation (15) means that the planet emits energy from the interior and is concurrently embedded in the protoplanetary nebula with a certain temperature, where T_neb is the midplane background nebula temperature as given by our disk model.

3.8.2. Detached (or transition) phase

The solution of the structure equations together with the specified boundary conditions leads to the well-known result that the gas accretion rate starts to increase exponentially once the core has grown to a mass of the order of 10 M_⊕ (see Sect. 6). This means that at some moment the gas accretion rate obtained in this way (i.e., through the planet’s Kelvin-Helmholtz timescale) must become higher than some external maximal possible gas supply. At this point, the planet enters the second phase and contracts to a radius that is much smaller than the Hill sphere radius. This is the “detached” regime of high-mass, runaway gas accretion planets. The planet adjusts its now free radius to the boundary conditions that are given by an accretion shock for matter falling at free fall velocity v_ff from about the Hill sphere radius to the planet’s surface. Probably more realistic boundary conditions would be those appropriate for the interface to a circumplanetary disk (Papaloizou & Nelson 2005). The maximal possible gas accretion rate Ṁ_XY is given by how much gas is supplied by the disk and can flow through the gap onto the planet Ṁ_max (e.g., Lubow et al. 1999). The calculation of this rate is specified in Sect. 4. The boundary conditions are now (cf. Bodenheimer et al. 2000; Papaloizou & Nelson 2005) $\begin{array}{ccc}\mathit{Ṁ}\mathit{XY}\mathrm{=}\mathit{Ṁ}\mathit{XY,}\max {\mathit{v}}_{\mathrm{ff}}^{\mathrm{2}}\mathrm{=}\mathrm{2}\mathit{GM}\left(\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{R}}_{\mathrm{H}}}\right)& & \\ \mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{neb}}\mathrm{+}\frac{\mathit{Ṁ}\mathit{XY}}{\mathrm{4}\mathit{\pi }{\mathit{R}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{v}}_{\mathrm{ff}}\mathrm{+}\frac{\mathrm{2}\mathit{g}}{\mathrm{3}\mathit{\kappa }}\mathit{\tau }\mathrm{=}\max \left({\mathit{\rho }}_{\mathrm{neb}}{\mathit{\kappa }}_{\mathrm{neb}}\mathit{R,}\mathrm{2}\mathit{/}\mathrm{3}\right)& & \\ {\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathrm{=}\frac{\mathrm{3}\mathit{\tau }{\mathit{L}}_{\mathrm{int}}}{\mathrm{8}\mathit{\pi \sigma }{\mathit{R}}^{\mathrm{2}}}{\mathit{T}}^{\mathrm{4}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{)}{\mathit{T}}_{\mathrm{neb}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathit{.}& & \end{array}$The pressure in Eq. (19) has three components: the ambient pressure (which, however, soon after the start of the collapse becomes very small compared to the other terms), the ram pressure due to the accretion shock ρv², and the standard Eddington expression for the photospheric pressure due to the material residing above the τ = 2/3 surface. In Eq. (20), A is the albedo assumed to have the same value of 0.343 as Jupiter (Guillot 2005), g is the gravitational acceleration GM/R², and for the luminosity we still have l(R) = L_int. In future work, we will use an albedo that depends on planet parameters (Cahoy et al. 2010). Numerically, one now iterates in this (and the next) phase on the radius R.

For the “hot start (accreting)” models, which (artificially) assume no radiative losses at all at the shock, it is not completely obvious whether or not we should include the pressure due to the accretion shock. However, we found that the dominant term for P is the photospheric pressure for all accretion rates we considered, and that the results are similar both with and without the ram pressure. For a gas accretion rate of 10^-2 M_⊕/yr, the photospheric pressure is always a factor 15 or more larger than the ram pressure, while for 10^-1 M_⊕/yr, it is still larger by a factor of at least about two. This dominance of photospheric pressure over ram pressure stems from low opacities at the typically relevant temperatures (~1500−3000 K) and densities at the surface of a giant planet undergoing gas runaway accretion. This temperature range approximately corresponds to the interval where grains have already evaporated, but where atomic and molecular opacities are still low, resulting in low κ ~ 0.01 cm²/g, even for full interstellar grain opacities (Bell & Lin 1994; Freedman et al. 2008). The dominance of the photospheric pressure is increased even more by the low grain opacity reduction factor f_opa = 0.003, which is normally used in the simulations. In simulations that are instead made with a full interstellar grain opacity and Ṁ_XY = 10^-1 M_⊕/yr, the ram pressure becomes dominant by a factor of a few for a short time at the beginning of the detached phase.

3.8.3. Evolutionary (or isolated) phase

The last phase starts when the gaseous disk disappears, so that the planet evolves at constant mass. During this phase, we use simple gray stellar photospheric boundary conditions in the Eddington approximation and write (Chandrasekhar 1939; Guillot 2005) $\begin{array}{ccc}\mathit{P}\mathrm{=}\frac{\mathrm{2}\mathit{g}}{\mathrm{3}\mathit{\kappa }}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathrm{=}\frac{{\mathit{L}}_{\mathrm{int}}}{\mathrm{4}\mathit{\pi \sigma }{\mathit{R}}^{\mathrm{2}}}& & \\ {\mathit{T}}_{\mathrm{equi}}\mathrm{=}\mathrm{280}\hspace{0.17em}\mathrm{K}{\left(\frac{\mathit{a}}{\mathrm{1} \mathrm{AU}}\right)}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\left(\frac{{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{M}}_{\mathrm{\odot }}}\right){\mathit{T}}^{\mathrm{4}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{)}{\mathit{T}}_{\mathrm{equi}}^{\mathrm{4}}\mathrm{+}{\mathit{T}}_{\mathrm{int}}^{\mathrm{4}}\mathit{.}& & \end{array}$Here we still have l(R) = L_int, which now also equals the total luminosity L as L_acc = 0 (Eq. (13)). This procedure means that our outer planetary radius R corresponds to the Rosseland mean τ = 2/3 surface. For Eq. (22), we have assumed that the planet is rotating quickly and redistributing the heat from stellar irradiation over its full surface, and that the stellar luminosity is proportional to ${\mathit{M}}_{\mathrm{\ast }}^{\mathrm{4}}$. This applies for roughly solar-like stars on the main sequence. In future work, we will include boundary conditions that are also accurate for planets close to the parent star undergoing intense irradiation (Guillot 2010; Heng et al. 2012). For the moment, we recall that the minimal allowed distance for planets in the model is 0.1 AU, where irradiation, at least for giant planets, is not yet a very important factor.

It is clear that the boundary conditions for the long-term evolution are described in a much simpler way than in the Burrows et al. (1997) or Baraffe et al. (2003) models, which employ proper non-gray atmospheres (Chabrier & Baraffe 1997). We have, however, found that the general agreement in terms of total luminosity or radius is very good (cf. Sects. 7 and 8) and certainly sufficient for our purpose of population synthesis. This is in agreement with Bodenheimer et al. (2000), who also found that his cooling curves using simple Eddington boundaries and those from Burrows et al. (1997) agree well.

4. Disk-limited gas accretion rate of the planet in the runaway regime

The second improvement of the model presented in this first paper regards the planetary gas accretion rate in the gas runaway accretion phase. It occurs for supercritical cores with a mass larger than about 10 M_⊕. In this phase, the gas accretion rate is no longer limited by the planet’s envelope structure, but rather by the rate at which gas can be provided by the disk. The gravitational interaction of the protoplanet and the surrounding disk determines the accretion rate across the growing gap (e.g., D’Angelo et al. 2010).

The actual gas accretion rate Ṁ_XY must be given by the smaller of the two rates: $\mathit{Ṁ}\mathit{XY}\mathrm{=}\mathrm{Min}\left(\mathit{Ṁ}\mathrm{Tkh}\mathit{,}\mathit{Ṁ}\mathit{XY,}\max \right)\mathit{,}$(23)where Ṁ_Tkh denotes the gas accretion rate as found through the solution of the structure equations in the attached phase, which is determined by the ability of the envelope to radiate away the potential energy of the accreted matter, while Ṁ_XY,max is the disk-limited rate.

In our previous models, we used for the disk-limited rate simply $\mathit{Ṁ}\mathit{XY,}\max \mathrm{=}\mathit{Ṁ}\mathrm{disk}\mathit{,}\mathrm{equi}\mathrm{=}\mathrm{3}\mathit{\pi }\begin{array}{}\mathrm{˜}\\ \mathit{\nu }\end{array}\mathrm{\Sigma }\mathit{,}$(24)which is the viscous accretion rate in a protoplanetary disk in equilibrium (i.e., where the mass flux is constant as a function of distance from the star). This expression, however, neglects three points: (i) the presence of the planet acting as a mass sink locally brings the disk out of equilibrium; (ii) only a fraction of the gas flux through the disk will be accreted onto the planet, while some part will flow past it (Lubow & D’Angelo 2006); (iii) there can be (initially) a local reservoir of gas around the planet that can be accreted faster than on a viscous timescale.

4.1. Non-equilibrium flux

The inner part of a viscous accretion disk evolves in absence of a planet relatively quickly into equilibrium. However, the sucking action of the planet and photoevaporation at late stages disturb the purely viscous evolution and bring the disk out of equilibrium. Then, the mass flux at a distance r from the star is given as $\mathit{Ṁ}\mathrm{disk}\mathrm{=}\mathrm{3}\mathit{\pi }\begin{array}{}\mathrm{˜}\\ \mathit{\nu }\end{array}\mathrm{\Sigma }\mathrm{+}\mathrm{6}\mathit{\pi r}\frac{\mathit{\partial }\begin{array}{}\mathrm{˜}\\ \mathit{\nu }\end{array}\mathrm{\Sigma }}{\mathit{\partial r}}\mathit{,}$(25)where the second non-equilibrium term becomes dominant close to the planet when it is in runaway gas accretion. In order to model the mass flux onto the planet, we have to consider different flow geometries of the gas in the gas-feeding zone around the planet. We look at the outer flux Ṁ_O at R_O = R_p + 1/2R_out and the inner flux Ṁ_I at R_I = R_p−1/2R_out. Here, R_p is the orbital distance of the planet and R_out is its gas capture radius, which is given in Eq. (14).

4.1.1. Flow geometries

Four different flow geometries where we count a mass flux towards the star as positive have to be considered. First, we can have a situation where both Ṁ_O and Ṁ_I are positive. This is the regular situation. It takes place when the planet is inside the radius of maximum viscous couple R_MVC (or radius of velocity reversal). At R_MVC, the flow in the disk changes from inward accretion to outward spreading. Second, we can have a situation where both Ṁ_O and Ṁ_I are negative. This is the case if the planet is outside R_MVC. The third case occurs when the planet is exactly at R_MVC. Then, decretion on both sides happens. This could, in principle, mean that only a very small amount of gas can be delivered to the planet. But since we take into account the local reservoir of gas (Sect. 4.2), this regime has no practical meaning. The fourth flow geometry occurs when gas flows towards the planet from both sides. This regime occurs due to the sucking action of the protoplanet at the beginning of runaway gas accretion, and then again near the end of the disk lifetime, when the global fluxes through the disk are small.

For these four geometries, the flow-limited maximal gas accretion rate Ṁ_XY,max,F is calculated as $\mathit{Ṁ}\mathit{XY,}\max \mathit{,}\mathrm{F}\mathrm{=}\{\begin{array}{c}\\ \mathit{Ṁ}\mathrm{O}& \mathrm{for}{\mathit{R}}_{\mathrm{p}}\mathit{<}{\mathit{R}}_{\mathrm{MVC}}\\ \mathit{Ṁ}\mathrm{I}& \mathrm{for}{\mathit{R}}_{\mathrm{p}}\mathit{>}{\mathit{R}}_{\mathrm{MVC}}\\ \mathrm{Abs}\mathrm{(}\mathit{Ṁ}\mathrm{I}\mathrm{-}\mathit{Ṁ}{\mathrm{O}}^{\mathrm{)}}& \mathrm{for}{\mathit{R}}_{\mathrm{p}}\mathrm{=}{\mathit{R}}_{\mathrm{MVC}}\\ \mathrm{Abs}\mathrm{(}\mathit{Ṁ}{\mathrm{I}}^{\mathrm{)}}\mathrm{+}\mathit{Ṁ}\mathrm{O}& \mathrm{fortwo}\mathrm{-}\mathrm{sidedaccretion}\mathit{.}\end{array}$(26)For the first two cases, we have assumed that the planet migrates more slowly than the gas flows. This is due to the fact that planets that are in gas runaway accretion typically are also in the planet-dominated type II migration regime, where this condition is fulfilled. For the third case, we consider the net flux, and for the last case, the sum of both fluxes is used.

4.2. Local reservoir

A local reservoir of gas, which is already in the vicinity of the planet, can be accreted at a rate larger than Ṁ_XY,max,F, at least for some time. This is relevant at the moment when the planet just passes into runaway gas accretion. Such a higher accretion rate is important for the depth of the so-called “planetary desert” (Ida & Lin 2004): the higher the possible accretion rates, the lower the number of planets of intermediate masses (ca. 10−100 M_⊕, see the discussions in Mordasini et al. 2009a, 2011b).

In order to calculate the accretion rate allowed by the local reservoir, we use an approach directly based on the local gas surface density, and not on the fluxes. Following D’Angelo & Lubow (2008) and Zhou & Lin (2007), we assume a spherical, homogeneous, unimpeded accretion flow for a planet with a cross section σ_cross in a gas of density ρ ≈ Σ/H (taken as the mean over the gas feeding zone) and a relative velocity v_rel, $\mathit{Ṁ}\mathit{XY,}\max \mathit{,}\mathrm{R}\mathrm{=}\mathit{\rho }{\mathit{\sigma }}_{\mathrm{cross}}{\mathit{v}}_{\mathrm{rel}}\mathit{.}$(27)Gas is captured at some radius R_gc, so that ${\mathit{\sigma }}_{\mathrm{cross}}\mathrm{\approx }\mathit{\pi }{\mathit{R}}_{\mathrm{gc}}^{\mathrm{2}}$. For simplicity, we set R_gc equal to the radius defined in Eq. (14). The relative velocity is given approximately as v_rel ≈ R_gcΩ. Therefore, the gas accretion rate in runaway, limited by the local reservoir Ṁ_XY,max,R, is approximately given as $\mathit{Ṁ}\mathit{XY,}\max \mathit{,}\mathrm{R}\mathrm{\approx }\frac{\mathrm{\Sigma }}{\mathit{H}}\mathrm{\Omega }{\mathit{R}}_{\mathrm{gc}}^{\mathrm{3}}\mathit{.}$(28)As shown by D’Angelo & Lubow (2008), such a simple description of the gas accretion rate is in fairly good agreement with results of 3D hydrodynamic simulations. It is, however, strictly speaking only valid in the intermediate phase, where the planet’s Hill sphere R_H is still smaller than the disk scale height H. The accretion flow is no longer spherical and homogeneous once the planet has grown so massive that R_H > H. Instead, it has a more complex geometry due to the tidal interaction of the planet and the disk, leading to gap formation and accretion streamers (Zhou & Lin 2007). However, we find that Ṁ_XY,max,R is larger than Ṁ_XY,max,F and thus the relevant quantity only during an intermediate phase just after the start of the disk-limited gas accretion phase (as illustrated in Fig. 1). In this phase, planets typically have masses still so small that R_H < H, so that the simple prescription is valid.

The final accretion rate in runaway is calculated as $\mathit{Ṁ}\mathit{XY,}\max \mathrm{=}{\mathit{k}}_{\mathrm{lub}}\mathrm{Max}\left(\mathit{Ṁ}\mathit{XY,}\max \mathit{,}\mathrm{R}\mathit{,}\mathit{Ṁ}\mathit{XY,}\max \mathit{,}\mathrm{F}\right)\mathit{,}$(29)where we have introduced a factor k_lub that takes into account that some gas flows past the planet. Following Lubow & D’Angelo (2006), we use values of k_lub of 0.75 to 0.9. Finally, we also make sure that not more gas is accreted than the amount present in the feeding zone M_feed, i.e., that Ṁ_XY,max < M_feed/dt. This condition, however, is automatically fulfilled with the above equations.

4.3. Example of disk-limited gas accretion rate

In Fig. 1, we give an example that illustrates the calculation of the gas accretion rate in runaway described in the last section. The figure shows gas accretion rates of a growing giant planet as well as accretion rates in the disk during the phase when the planet starts gas runaway accretion. The planet is at about 10 AU in a gas- and solid-rich disk. The parameter k_lub is set to 0.9, and the viscosity parameter α in the disk is 0.007.

The thick solid black line shows the accretion rate of the planet Ṁ_XY. The thin solid black line shows the disk-limited accretion rate Ṁ_XY,max. Shortly after 1.33 Myr, the accretion rate found by the solution of the structure equations Ṁ_Tkh increases very rapidly, so that the disk-limited rate is quickly reached. In this phase, Ṁ_XY,max is given as k_lubṀ_XY,max,R. This is because Ṁ_XY,max,R (shown by the magenta dashed-dotted line) is larger than Ṁ_XY,max,F (shown by the red dotted line). We thus find that the local reservoir allows for an accretion rate higher than that given by the viscous rate, namely about 0.03 M_⊕/yr. A typical peak accretion rate of a few times 10^-2 M_⊕/yr has also been found in other studies (e.g., Lissauer et al. 2009). However, later on (starting at about 1.35 Myr), we note that both Ṁ_XY,max,R and Ṁ_XY,max,F converge on approximately the same value. Also, M_feed/dt (cyan long-dashed-dotted line) then takes approximately the same value.

From the curves showing the mass flux inside and outside of the planet, Ṁ_I (blue dashed line) and Ṁ_O (green long-dashed line), we can deduce in which flow geometry we are. Initially, both fluxes are positive, which means that we are in the standard case, where the planet is inside the radius of maximum viscous couple R_MVC and the disk is globally flowing towards the star. At about 1.325 Myr, however, we see that Ṁ_I becomes negative, meaning that the flow geometry changes into the “two-sided accretion” case. This is a consequence of the accretion onto the planet at a rate higher than the viscous one. Later on (at about 1.6 Myr, not shown in the figure), the flow geometry changes back into the standard geometry, and Ṁ_I is then about ten times smaller than Ṁ_O. This corresponds to the 10% of gas that is allowed to flow past the planet. Now, a new equilibrium has established, with the planet acting as an additional sink.

5. Results

We address three different subjects. First, we present a combined formation and evolution calculation for a Jovian mass planet at 5.2 AU. We show the evolution of the most important observable quantities, such as mass, radius, and luminosity in addition to the evolution of the central temperature and pressure. The calculation can be compared with the calculations of, e.g., Pollack et al. (1996), Bodenheimer et al. (2000), or Lissauer et al. (2009). Second, we compare the planetary M-R relation for giant planets (see also Paper II) and the influence of the core mass on the radius as found in our model with several other studies. Third, we use our model to study the luminosities of young Jupiters. The results for “cold starts” can be compared with Marley et al. (2007) or Spiegel & Burrows (2012).

6. Example of a combined in situ formation and evolution simulation

In this section, we use the extension of the planet module to calculate an illustrative example of a combined formation and evolution simulation.

These simulations are strongly simplified in comparison to the calculations performed for the complete population synthesis models (Paper II). In the full model, the moment of the onset of limited gas accretion (and thus detachment from the nebula), the disk-limited runaway gas accretion rates (Sect. 4), and the final mass of the planet result in a self-consistent way from the evolution of the gaseous disk (Alibert et al. 2005a). Here, the evolution of the disk is switched off. Instead, we fix the maximal allowed Ṁ_XY to some prespecified value and artificially switch off accretion on a short timescale when the mass of the planet approaches the desired value. Migration is also completely switched off, otherwise, it is a central component of the formation model as it can lead to much more complex accretion histories (Mordasini et al. 2009a). But this way we allow for direct comparison with numerous previous works that used similar assumptions as, e.g., Pollack et al. (1996), Lissauer et al. (2009), and Movshovitz et al. (2010).

In Table 1, we list the most important settings used for the simulations. Other settings, for example, a planetesimal radius of 100 km, are the same as in Alibert et al. (2005a), except when mentioned.

6.1. Grain opacity reduction factor f_opa

A new parameter in our models is the reduction factor of the grain opacity in the envelope relative to the ISM grain opacity, f_opa. We determined this factor in Mordasini et al. (2012c) in the following way: Movshovitz et al. (2010) coupled for the first time self-consistently a giant planet formation simulation with the concurrent calculation of the grain growth and settling inside the planet’s gaseous envelope. They found very low resulting effective opacities and thus very short durations of the (plateau) phase II (Pollack et al. 1996). Such grain evolution calculations are currently beyond the scope of our modelbecause they are computationally extremely expensive. But to still get an approximative representation of this important effect, we compared the duration of phase II as found in our model with the duration found by Movshovitz et al. (2010) for different values of f_opa and the three planetesimal surface densities considered in Movshovitz et al. (2010). There is a a very clear relationship between f_opa, the planetesimal surface density, and the duration of phase II. We have found that on the mean (for the three planetesimal surface densities) the duration of phase II becomes the same in both models if we reduce the interstellar grain opacities by a factor 0.003.

This strong reduction is clearly even smaller than the previously studied “low opacity case” (Pollack et al. 1996; Hubickyj et al. 2005), which arbitrarily assumed a f_opa = 0.02. Two points should be considered. First, it is clear that the calculations of Movshovitz et al. (2010) are simplified in such aspects as the use of spherical grains in opacity calculations instead of aggregates, or the neglect of icy grains. In addition, a uniform grain opacity reduction factor cannot reproduce the detailed structure found by Movshovitz et al. (2010) for opacity as a function of radius in the envelope. It nevertheless leads to envelope structures that are at least similar, especially in the deeper layers. These two points mean that using a uniform grain opacity reduction factor, calibrated via the simulations of Movshovitz et al. (2010), is a strong simplification of the actual processes in a protoplanetary envelope. However, it represents a better approximation than using arbitrarily reduced values. The reader is referred to Mordasini et al. (2012c) for details.

6.2. Jovian mass planet at 5.2 AU

For this simulation, we use the same initial planetesimal surface density, disk pressure, and temperature as Pollack et al. (1996). This case has been studied frequently in the literature, and we consider it our nominal model. The luminosity is calculated as in the “cold start” assumption, i.e., we assume that the shock luminosity does not contribute at all to the luminosity inside the planetary structure.

The formation of Jupiter is of particular importance as a benchmark for any planet formation model. This is because for no other giant planet does an equally high number of detailed observational constraints exist. In Table A.1 in the appendix, we list the most important quantities characterizing the planet such as the total mass, core mass, or luminosity at particular important moments in time. The data in these tables can be compared with similar data in previous studies mentioned above.

The simulations presented in this section differ from those in Pollack et al. (1996) in a number of points, making them more realistic. The differences are that we include planetesimal ejection, core density is variable (Paper II), and f_opa = 0.003, which leads to a strongly reduced formation timescale. In Appendix A, we also present simulations where we set these parameters to the same value as in Pollack et al. (1996) for quantitative comparisons with their simulations.

6.2.1. Formation phase: mass

In Fig. 2, we show the most important planetary properties as a function of time during the formation and at the very beginning of the evolutionary phase once the final mass is reached.

The top left panel shows the characteristic three stages labeled with I, II, and III (Pollack et al. 1996) of such classical in situ calculations (all three occur within the attached phase described in Sect. 3.8). First, a rapid build-up of the core occurs until the isolation mass is reached, then a plateau phase follows (characterized by a slow increase of the envelope mass, which allows further core growth), and then the transition to gas runaway accretion is observed.

The crossover point when the core and envelope masses are the same is reached relatively quickly at 0.82 Myr due to the strongly reduced grain opacity. This duration is in good agreement with Movshovitz et al. (2010). The crossover mass is M_cr = 16.3 M_⊕. This is in very good agreement with Pollack et al. (1996) with M_cr = 16.2 M_⊕ or Hubickyj et al. (2005), who obtained M_cr = 16.1 M_⊕ for the same initial planetesimal surface density.

Shortly after the crossover point, the gas accretion rate increases strongly because of beginning runaway accretion, so that it hits the maximal allowed value of 10^-2 M_⊕/yr at t = 0.923 Myr. The total mass of the planet is then 86.7 M_⊕, and τ_KH is as short as ~1640 yrs. At this moment, the second detached phase begins and the collapse of the envelope starts, which is in fact a rapid, but still hydrostatic contraction as discussed in Sect. 3.3 (Bodenheimer & Pollack 1986). This phase is indicated with the label D in the plot.

At a high gas accretion rate of 0.01 M_⊕/yr, a mass of 1 M_X is approached quickly, so that at 0.935 Myr we start to artificially ramp down the accretion rates (see top right panel of Fig. 2). Shortly afterwards, at about 0.961 Myr, the final mass of the planet is almost reached (99.5% of the final mass). This marks the beginning of the evolutionary phase, labeled with an E. Note that the accretion rate of planetesimals also increases strongly during the runaway phase because of the fast expansion of the planet’s feeding zone. It reaches a maximum of about 7 × 10^-4 M_⊕/yr and then starts to decrease before we ramp it down artificially³. This is due to the fact that ejection of planetesimals becomes increasingly important as the mass of the planet grows (Alibert et al. 2005a) and because the capture radius of the planet shrinks (bottom left panel). The final core mass that is obtained is therefore very similar to the one that is found without stopping the accretion artificially, as becomes clear when looking at Fig. 6.

The final total mass of the planet is 316.6 M_⊕, while the final mass of the core is 32.9 M_⊕. In this paper, we assume that all accreted planetesimals can reach the core. In reality, once the envelope has grown above a certain mass (about 1−2 M_⊕ for 100 km planetesimals, Mordasini et al. 2005), planetesimals cannot penetrate to the core of the planet but get destroyed in the envelope. Therefore, the 32.9 M_⊕ value should be associated with the total mass of heavy elements contained in the planet rather than the core mass.

The theoretical estimates for the actual core and total heavy element mass in Jupiter vary. Uncertainties in the EOS of hydrogen and helium as well as in the efficiency of convection make an exact determination of the enrichment difficult. Based on the work of Saumon & Guillot (2004), the total heavy element mass was estimated by Guillot & Gautier (2009) to lie between 10 to 42 M_⊕. Besides classical equations of state like the EOS SCvH, which are based on the free energy minimization method, EOS based on ab initio quantum molecular dynamics calculations have become available in the past few years. However, the results for the core and total heavy element mass of Jupiter derived from this method vary: Militzer et al. (2008) find a massive core of M_core = 16 ± 2 M_⊕ and M_Z = 20 ± 4 M_⊕, which is quite different from the result of Nettelmann et al. (2012) of M_core = 0−8 M_⊕ and M_Z = 28−32 M_⊕. Finally, Leconte & Chabrier (2012), who again use the EOS SCvH but allow for semiconvection due to stabilizing compositional gradients, find high M_Z = 41−63 M_⊕ as the interiors are hotter. These results demonstrate that the investigations regarding the composition of Jupiter have not yet settled on a definitive conclusion.

In our model, the accretion of the core occurs in two steps: about 10−15 M_⊕ are accreted in phase I and phase II, while the rest is accreted in the gas runaway accretion phase. The mass of the envelope is already high, and planetesimals should be destroyed in the envelope. Therefore, we can estimate that the actual core mass should be of the order 12 M_⊕ (this is when the envelope mass reaches 2 M_⊕), while the remaining 21 M_⊕ of metals should be dissolved in the envelope, unless they sink to the core-envelope interface (Pollack et al. 1996).

6.2.2. Formation phase: radius

The bottom left panel of Fig. 2 shows the evolution of the total radius, the core radius, and the capture radius. The total radius is initially (during the attached stage) very large, as it is approximately equal to one-third of the Hill sphere radius. When the limiting gas accretion rate is hit, the planet detaches from the nebula and collapses to a radius of initially about 2 to 3 Jovian radii. The shape of the curve showing the total radius is very similar to that found by Lissauer et al. (2009).

The capture radius (which is found by integrating the orbits of planetesimals inside the envelope, see Alibert et al. 2005a) is during phase II much larger than the core radius, namely, by a factor between ~10 and 60. This is very important for the core growth, as the capture radius enters quadratically into the growth rate. After the envelope has started to collapse, R_capt becomes equal or slightly less than the total radius. This is due to the fact that in the collapse the envelope structure becomes much more compact. Before the collapse, it is characterized by very large scale heights of thousands to several ten thousands of kilometers, where planetesimals get captured in a soft catch. Durning the collapse, the planet develops a well-defined surface (i.e., a much smaller scale height near the surface) at which planetesimals are captured.

Looking at the core radius, we note that in the gas runaway accretion phase when the planet grows very quickly in mass, it first increases a bit in size, but then decreases, even though the core grows in this phase by about 20 M_⊕. The reason is that the external pressure on the core by the accreted gas increases so much that the resulting strong increase of the core density over-compensates the mass increase (see Paper II for a dedicated discussion).

The bottom right panel of Fig. 2 shows the core luminosity L_core, the total luminosity L, and the luminosity within the planetary structure L_int. The difference L − L_int is thus the shock/accretion luminosity L_acc. The first peak in the curve is due to the rapid accretion of the core and the second to the combined effects of runaway gas accretion and envelope collapse. During the second peak, for about 2 × 10⁴ years, the accretional luminosity is larger than L_int by up to a factor ~3. For more massive planets discussed below, this ratio can become much larger. The second luminosity maximum occurs at 0.94 Myr, with L = 2.6 × 10⁶ L_X corresponding to log (L/L_⊙) = −2.65. Lissauer et al. (2009) find in their simulations with the same limiting gas accretion rate peak luminosities of about log (L/L_⊙) = −2.35, which is similar to our value.

We note that the classical Pollack et al. (1996) picture with three well-defined phases might not be a realistic formation scenario for Jupiter. Several factors, including migration (Alibert et al. 2004), a slower core growth (Fortier et al. 2007), and the opening of a gap in the planetesimal disk (Shiraishi & Ida 2008), might suppress in particular phase II. We have two reasons for investigating the classical scenario: we want to compare our model with previous studies and we are primarily interested in the collapse and evolutionary phase, which should be similar in modified scenarios.

6.2.3. Evolution phase: radius and luminosity

Figure 3 shows the radius and luminosity of the planet as a function of time, including the long-term evolution over Gigayears when the mass is constant. The evolution now occurs slowly by gradual contraction and cooling. For comparison, the temporal evolution of the radius and the luminosity of a 1 M_X planet as found by Baraffe et al. (2003) and Burrows et al. (1997) is also shown. The temporal zero point in these models has been associated with the time the planet reaches its final mass in our calculations. Baraffe et al. (2003) and Burrows et al. (1997) use non-gray atmospheric boundary conditions and the standard method based on the −T∂S/∂t term to calculate the evolution. As we can seen, there is good agreement among the models. The differences between our model and the other two ones are of the same order as the mutual differences between Baraffe et al. (2003) and Burrows et al. (1997). The different internal compositions (mass of heavy elements) have a certain influence on the exact value of R, so that we cannot expect exactly identical results.

After 4.6 Gyr, the planet has a radius of 0.99 R_X and a luminosity of 1.13 L_X, corresponding to log (L/L_⊙) = −9.01. At the same age, Burrows et al. (1997) found a radius of 1.11 R_X and a luminosity of 1.52 L_X. Larger values are expected, as the models in Burrows et al. (1997) are core-less models.

6.2.4. Central conditions

In Fig. 4, left panel, we show various central conditions during both the formation and evolution phase. With “central” conditions, we mean here the envelope-core boundary and not the actual center of the planet.

The solid red line in the left panel shows the central temperature. It belongs to the scale on the left. Initially, the temperature at the bottom of the envelope is about 3000 K, rising to values between 10 000 and 20 000 K during phase II. When gas runaway accretion starts and the envelope collapses, there is a sharp upturn of T_cent. The maximal value of 69 700 K is reached at 0.98 Myr, afterward it declines again. At an age of 4.6 Gyr, T_cent has fallen to about 17 600 K, which is in good agreement with the estimates of Guillot (1999), which range from 15 000 to 21 000 K.

The second quantity for which the scale is shown on the left is the central pressure, indicated by the black dotted line. It is found to have a value of 4317 GPa at the current age of the Solar System. Guillot (1999) gives a value of 4000 GPa.

The other two lines in the figure belong to the y-axis on the right. The green long-dashed line is the mean density of the core, while the blue dashed line is the gas density at the envelope-core boundary. As described in Paper II, the core consists of iron, followed by a silicate layer and finally an ice layer, with an ice fraction expected for a body forming outside the iceline. During phase II, ρ_core lies between 3 and 4 g/cm³. This is close to the typically assumed constant value of 3.2 g/cm³ in other formation calculations (e.g., Movshovitz et al. 2010). During the collapse, the external pressure exerted onto the core increases dramatically, so that the core density increases to about 10 g/cm³, leading to the aforementioned shrinking of the core size, despite its growth in mass. At 4.6 Gyr, we find ρ_core = 14.3 g/cm³. The central gas density quantitatively follows a similar pattern, starting of course at a smaller value of about 0.1 g/cm³ during phase II. At 4.6 Gyr is has reached a fluid-like density of about 4 g/cm³.

6.2.5. Surface conditions

The right panel of Fig. 4 shows the conditions at the surface of the planet. The solid red line is the surface temperature T_surf calculated as ${\mathit{T}}_{\mathrm{surf}}^{\mathrm{4}}\mathrm{=}{\mathit{T}}^{\mathrm{4}}\mathrm{+}\frac{{\mathit{L}}_{\mathrm{acc}}}{\mathrm{4}\mathit{\pi \sigma }{\mathit{R}}^{\mathrm{2}}}\mathit{,}$(30)where T is defined in Eq. (20). It thus contains the contributions from the intrinsic luminosity, the accretional luminosity (which is, however, important only for a short period; see Fig. 2 bottom right panel), and the absorbed stellar light. When the accretional luminosity is negligible, T_surf = T, which is the temperature used as the outer boundary condition for the structure. When both L_int and L_acc are negligible, T_surf = (1 − A)T_equi. This happens for a Jovian planet at 5.2 AU after ≳ 10 Gyr, when the intrinsic luminosity of the planet becomes small compared to the incident stellar flux. The background temperature is indicated by the black dotted line. During the attached and detached phase, the background temperatures is equal to the (constant) nebula temperature T_neb = 150 K, while in the evolutionary phase, it is given by the equilibrium temperature with stellar irradiation T_equi at 5.2 AU (Eq. (22)). The two temperatures belong to the left scale. In phase I, T_surf can become about a factor 2.5 larger than T_neb, as L is quite large in this phase. In phase II, T_surf ≈ T_neb, in agreement with Papaloizou & Nelson (2005). During the collapse, there is a strong increase of T_surf, which reaches a maximal value of 2050 K. At the moment the final mass is reached, T_surf ≈ 1350 K. At the current age of Jupiter, T_surf = 126 K, while the measured value is about 124 K (Guillot & Gautier 2009). We see that at late times the decrease of T_surf starts to flatten out as it approaches the constant value of (1 − A)T_equi. In the figure, the dashed-dotted red line is (L_int/(4πσR²))^1/4 i.e., the temperature due to the intrinsic flux only. This quantity continues to fall. The nomenclature for the different temperatures is unfortunately not homogeneous in the literature: our quantity T_surf corresponds during the evolutionary phase to T_therm, according to the definition in Baraffe et al. (2003), but to T_eff in Fortney & Nettelmann (2010).

In addition, we show the pressure P at the surface of the planet. This curve belongs to the right axis. During the attached phase, it is simply equal to the (constant) nebula pressure (Eq. (14)). Then, at the moment that the planet detaches, P increases as P_edd = 2/3 g/κ becomes large. The dynamic/ram pressure ρv² also contributes to P in this phase (Eq. (19)). The ram pressure remains, however, always 1−2 orders of magnitude smaller than P_edd, which is related to the low grain opacity (Sect. 3.8.2). At late times, after the accretion is over, we see that P has a peculiar shape, which is given by opacity transitions. At 4.6 Gyr, the pressure is about 0.18 bars.

6.2.6. Atmospheric pressure-temperature profile

In Fig. 5, we show the pressure as a function of temperature for the uppermost part of the planet down to a pressure of 1 kbar at different moments in time. The evolutionary sequence starts at the first structure (at the right) at a time of 1 Myr (i.e., shortly after the planet has reached the final mass), while the last structure (at the left) corresponds to 4.6 Gyr. Structures are shown at irregular intervals between these times. We recall that the surface of the planet and thus the upper end of the lines corresponds to the τ = 2/3 surface. During the evolutionary sequence, the radius of the planet shrinks from 2.36 to 0.99 R_X, which causes an increase of gravitational acceleration g from 444 at the beginning to 2546 cm/s² at 4.6 Gyr. In the figure, the thick black lines indicate radiative zones, while the blue dotted line shows the p − T structure measured by the Galileo probe (Lodders & Fegley 1998). We find a very good agreement.

The figure can be approximately compared with Fig. 2 in Burrows et al. (1997). These authors use non-gray atmospheres, but indicate the location which corresponds to our outer surface. One must also take into account that they keep the gravity at a constant value of 2200 cm/s², so that we cannot expect to see identical results, especially at early times. But a comparison shows that the two p − T profiles nevertheless have a similar general structure. In both works, there exists, for example, a detached radiative zone at temperatures around 1500 K.

In summary, we see from the in situ formation and evolution calculation for Jupiter that our combined formation and evolution model, despite the simplifications, leads to a planet with basic properties in very good agreement with the most important observational constraints coming from Jupiter. This new simple method to calculate evolutionary sequences leads to results that are equivalent to those of the traditional entropy-based method.

7. Radii

We now generalize the results from the last chapter concerning the radius and study the radii of giant planets of different masses and different compositions. The goal is to validate our model by comparing it with existing work (e.g., Fortney et al. 2007; or Baraffe et al. 2008).

We do this by performing the same combined formation and evolution calculations as shown for Jupiter, the only difference being that we vary the initial planetesimal surface density (which will eventually lead to different core masses) as well as the moment when we terminate the gas accretion, so that we get different total masses. Otherwise, the calculations are identical, so that, for example, our calculations apply for planets at a distance of 5.2 AU from a solar-like star.

Such calculations are shown in Fig. 6. The plot shows the total and core mass as a function of time for planets of final masses of 0.15, 1, 2, 5, and 10 M_X. The lines for the 1 Jupiter mass planet are the same as in Sect. 6. The gas accretion rate in runaway is 0.01 M_⊕/yr in all cases. Therefore, more massive planets reach their final mass later. It is interesting to note that all planets have a final core mass of about 33 M_⊕, except for the lowest mass planet (total mass 0.14 M_X = 44.5 M_⊕) which has a core of about 18 Earth masses. The lower core mass in this case is due to the fact that for this planet, gas and solid accretion are (externally) ramped down before it reaches the gas runaway phase. The nearly identical core mass of all other planets is in contrast not externally imposed. It is rather a natural consequence of two effects. First, the capture radius starts to decrease at the moment when the planet collapses. This happens always at the same mass, independent of the final mass. Second, with increasing mass, the planet becomes efficient in ejecting planetesimals rather than accreting them (Sect. 6.2.1, Ida & Lin 2004).

Figure 7 shows the internal pressure-temperature structure of these planets (and of a 20 M_X planet) at an age of four to five Gyr. The left end of the lines corresponds to the surface of the planet, while the right end corresponds to the envelope-core interface. The 0.14 M_X planet has near the surface a significant radiative zone. Otherwise, the planets are nearly fully convective and characterized by a single adiabat. Near a pressure of about 1 Mbar, we see a change in slope, which comes from the molecular to metallic transition of hydrogen, also visible in a similar figure in Guillot & Gautier (2009).

7.1. Mass-radius relation

The M-R relation of (giant) planets has been studied for a long time (e.g., Zapolsky & Salpeter 1969). The interest in this relation lies in its connection to the composition of the planet and the state of matter in its interior. For recent reviews, see Chabrier et al. (2009) and Fortney et al. (2010).

The general result qualitatively found by Zapolsky & Salpeter (1969) is that the M-R relationship in the giant planet regime is characterized by a local maximum in R. This behavior can be understood with polytropic models with a polytropic index that increases with mass. The change is in turn due to the increasing importance of the degeneracy pressure of electrons relative to the classical coulomb contribution of ions with planet mass (e.g., Chabrier et al. 2009). Another general result is that the radius of giant planets decreases with core mass and increases with proximity to the star (e.g., Fortney et al. 2007; Baraffe et al. 2008).

In Fig. 8, we show the mass-radius relation for planets with masses between 0.14 and 20 M_X. These are the same models as in Figs. 6 and 7. Except for the lowest mass planet with a core of about 18 M_⊕, the core mass is approximately 33 M_⊕. The radii are shown at 0.1, 1, 4.6, and 10 Gyr. The planets were calculated using the cold start assumption, which, however, plays only a role at t = 0.1 Gyr and for planets more massive than 5 M_X. In the top right corner of the figure, we show with the magenta dotted line the radius for simulations performed with the “hot start (accreting)” scenario. The radii are as expected larger, by about 0.07 R_X for the most massive planet. For all other cases, the imprint of the formation is lost by ~0.1 Gyr and the radii are virtually identical.

The radii decrease with time, as expected. At all times, there is a maximum of the radius at a mass of about 3−5 M_X, in agreement with, e.g., Chabrier et al. (2009) or Fortney et al. (2010). For comparison, we also show the results from a more detailed model by Fortney et al. (2007). The very light blue dotted lines show the radius as found by them at an age of 4.5 Gyr for planets with a core of 25 M_⊕ (and thus similar as here) at orbital distances of 1 AU (upper line) and 9.5 AU (lower line). We see that these lines bracket our red line at almost all masses, which is exactly what we expect when the two models yield very similar results. Only for the most massive planets do ≳ 10 M_X we see that our models yield slightly larger radii, by a few 0.01 R_X. Also the lowest mass planet has in our simulations a somewhat larger radius than the counterpart in Fortney et al. (2007). This is, however, not a surprise because it has a smaller core by about 6 M_⊕, which is relevant due to the small total mass of about 44 M_⊕ (cf. Fortney et al. 2007).

7.2. Influence of the core mass

The fact that a more massive core for a given total mass leads to a smaller radius has important implications for the study of planets. It allows us, for example, to infer (within certain limits) the bulk composition of transiting extrasolar planets, which in turn constrains formation models (see Paper II). For example it is clear that giant planets formed by the core accretion mechanism must be enriched in heavy elements, as in the case for Jupiter. For the competing formation model, the disk instability model, the situation is in contrast less clear (Boley et al. 2011).

Another important application is to connect the inferred core masses of transiting extrasolar planet with their formation environment. In this context, it was found that the core mass of giant planets and the metallicity of their host star are positively correlated (Guillot et al. 2006; Burrows et al. 2007). This is reproduced by planet population synthesis calculations based on the core accretion paradigm (Mordasini et al. 2009b). Recent findings (Miller & Fortney 2011) even indicate that all transiting giant planets contain a minimum amount of 10−15 M_⊕ of metals, as predicted earlier by the core accretion mechanism (Mordasini et al. 2011c).

It is therefore important to study the dependence of total radius on the core mass and to compare this with previous studies. Figure 9 shows the radius as a function of core mass for planets with a total mass of 1 M_X. The red solid line shows the result from this work for five planets at an age of 4.6 Gyr. We see that for an increase of the core mass from ~10 to ~50 M_⊕ there is a reduction of the radius by about 0.07 R_X. We also show the results from Fortney et al. (2007) and Baraffe et al. (2008). Our result for a = 5.2 AU is located between the curves for 1 and 9.5 AU from Fortney et al. (2007) at 4.5 Gyr, with very similar slopes of the curves. From this we conclude that the two models yield very similar results.

An interesting aspect of the combined formation and evolution calculations is that we can directly relate the radius (respectively the core mass) with the initial surface density of the planetesimals during formation, which is not possible for purely evolutionary models. In the current case, the five different core masses correspond to planetesimal surface densities of 3.5, 5, 7.5, 10 and 15 g/cm². We must, however, take into account that there is no unique mapping between planetesimal surface density and core mass, because other factors like migration or the distance where the planet forms also influence this quantity.

7.3. Comparison and limitations

We conclude from the different comparisons that our upgraded model makes it possible to obtain evolutionary sequences and thus radii with very good accuracy starting with self-consistent initial conditions.

It is, however, also clear that we have mostly addressed relatively “simple” cases of massive, gas-dominated planets at large distances from the star. In the solar systems, there are planets that could not be very accurately modeled with a simple model as used here. It is, for example, well known that for Saturn, He demixing and settling is necessary to explain the measured luminosity (e.g., Hubbard et al. 1999). We must also expect that for low-mass, core-dominated planets, possibly with strongly enriched envelopes, we would get less accurate cooling curves and radii, as we do not model the thermodynamics of the core, assume that all solids are in the core (cf. Baraffe et al. 2008), and also ignore effects of the composition on the opacity. Other limitations apply for strongly irradiated planets. These effects will be added in later modifications of our model.

8. Luminosities

The other fundamental observable quantity apart from the radius that we obtain from our upgraded model is the planetary luminosity. It is obviously central for direct imaging observations. For transit observations, the discovery of giant planets with radii much larger than expected from standard internal structure modeling was a surprise. Regarding the direct imaging technique, there is on theoretical grounds an open question: what is the luminosity of young Jupiters? Currently, two classes of models are discussed in the literature: “hot start” models (Burrows et al. 1997; Baraffe et al. 2003), and “cold start” models” (Fortney et al. 2005; Marley et al. 2007). In the “hot start” model, one assumes as the initial state an already fully formed planet at an (arbitrarily) high entropy state, and thus a large radius and luminosity. In the “cold start” scenario, one assumes that the planet is gradually built up by accretion of gas through an accretion shock on the surface of the planet. The radiative losses of the liberated gravitational energy at the shock lead to a lower entropy, and thus to lower luminosity and radius.

8.1. Relationship of cold/hot start and core accretion/gravitational instability

The importance of the two different initial states stems from the fact that the luminosity of young planets is an observable quantity establishing links to the physical mechanisms occurring during the formation. The currently known directly imaged planets seem to be incompatible (Janson et al. 2011; Marley et al. 2012) with the “cold start” model of Marley et al. (2007), since the observed temperatures and luminosities are higher than predicted by this model. It should be kept in mind that the model of Marley et al. (2007) studies the limiting case of a shock that is 100% efficient in radiating away the accretional energy (and makes a number of assumptions concerning the formation process, cf. below). This limiting case is studied because detailed, possibly multidimensional radiation-hydrodynamic calculations of the shock properties have not yet been conducted. Recently, Spiegel & Burrows (2012) have thus studied a broad range of intermediate “warm-start” models and showed that they lead to significant differences in the brightness of the planets, depending on the mass and spectral band.

The related question about the initial thermodynamic state of newly formed brown dwarfs and stellar objects has been studied in various works (e.g., Prialnik & Livio 1985; Hartmann et al. 1997; Baraffe et al. 2009; Commerçon et al. 2011). These studies indicate that, depending on the fraction of accretion luminosity radiated away from or absorbed by the protostar, the formation through gravitational collapse can lead to bright, extended objects (“hot accretion” leading to a “hot start”) or faint and compact objects (“cold accretion” leading to a “cold start”). The latter is the case if the accretion shock during the formation process is supercritical, i.e., if all the accretion shock energy is radiated away, as found recently by Commerçon et al. (2011) for accretion on the first Larson core.

For giant planets, two formation mechanisms are currently discussed. If a “cold” or “hot” start can preferentially be associated with one of the two formation mechanisms, then the understanding of the luminosity of young Jupiters has important implications for the understanding of giant planet formation as a whole.

8.1.1. Core accretion

The first proposed formation mechanism is core accretion, which is the paradigm underlying this work. A number of observational findings indicate that core accretion represents the dominant formation mechanism, at least for planets inside a few AU (e.g., Mordasini et al. 2012a). Whether core accretion can also explain the directly imaged, massive planets at orbital distances ≫ 10 AU is the subject of an ongoing debate (e.g., Dodson-Robinson et al. 2009; Kratter et al. 2010). The core accretion, “cold accretion” simulations of Marley et al. (2007) lead as mentioned to very low initial luminosities, incompatible with the observed values (Janson et al. 2011; Marley et al. 2012). We will, however, show in Mordasini et al. (in prep.) that this result (found for a specific set of parameters during the formation process) does not mean that planets formed by core accretion in a “cold accretion” scenario must necessarily (for all reasonable parameters) have luminosities as low as in Marley et al. (2007) (see also the recent calculations of Spiegel & Burrows 2012). Furthermore, the results in the next two sections show that the formation of giant planets via core accretion with “hot accretion” leads to high initial luminosities, which very quickly converge on the luminosity of classical “hot start” calculations (Burrows et al. 1997; Baraffe et al. 2003). We thus see that, depending on currently ill-constrained physical mechanisms occurring during the formation phase (in particular the nature of the accretion shock), the luminosity of planets formed by core accretion can probably vary significantly.

8.1.2. Gravitational instability

For the second formation mechanism, the gravitational instability (GI) model, the question about the associated luminosity might also be relatively complex, as indicated by the explorative discussion in this section. More solid results are expected once multidimensional, radiation-hydrodynamic calculations following the entire formation by gravitational instability from the first instability in the disk to a compact planet at the final mass become feasible.

In a qualitative way, we could expect that the post-formation luminosity of planets formed by GI depends on the specific way these objects grow in mass. Hydrodynamic simulations show that the initial mass of the clump that becomes self-gravitationally unstable is typically 3−20 M_X (e.g., Forgan & Rice 2011; Helled & Bodenheimer 2011; Zhu et al. 2012). These masses roughly correspond to the ones expected analytically from the mass contained in the disk patch with a length scale similar to the most unstable mode. The objects that form at this stage are self-gravitating, relatively cold, extended objects. Their radius is comparable to the Hill sphere radius (e.g., Zhu et al. 2012). For example, a 10 M_X clump at 50 AU from a 1 M_⊙ star, which may represent a typical outcome (e.g., Janson et al. 2011), has a Hill sphere radius of about 7.4 AU. These clumps consist of molecular hydrogen and, once formed, contract slowly (e.g., Helled & Bodenheimer 2011). In their nature, they correspond to the first Larson core of star formation.

The mass contained in such an object has not gone through an accretion shock, since it forms the spontaneous, self-gravitational binding of an extended part of the disk, typically within an over-dense spiral arm (e.g., Boley et al. 2010). We therefore expect the specific entropy of the gas in the clump to be high. This is confirmed by numerical simulations. Hayfield (2012, pers. comm.; see also Galvagni et al. 2012) finds that the specific entropy in a clump taken from the radiation-hydrodynamic simulations of Boley et al. (2010) is about 15 k_B/baryon at the moment when the clump should undergo the second collapse (cf. below). If this is representative of the entropy of the final object, it corresponds to a “hottest start” (Spiegel & Burrows 2012). This allows a first supposition to be made: if the final mass of a GI planet is the same (or similar) as the mass of the initially unstable clump (an approximation that is regularly made, e.g., Janson et al. 2011), then we may expect these objects to have a “hot start”, simply because there is no accretion shock involved during their formation.

However, since the clump necessarily forms in a massive disk, it appears more likely that it continues to accrete mass (in principle, gap formation can stop accretion, e.g., Zhu et al. 2012). In this case, we have to distinguish two sub-scenarios. As the first core contracts quasi-statically, its interior heats up, until a central temperature of about 2000 K is reached. Then, molecular hydrogen dissociates, and a dynamical collapse of the entire clump occurs (e.g., Bodenheimer et al. 1980). This corresponds to the formation of the second core in star formation. The timescale until this second collapse occurs is estimated to be roughly 10³ to 10⁴ years, with more massive clumps collapsing earlier (Helled & Bodenheimer 2011).

In the first sub-scenario, the clump grows to its final mass predominately before it collapses. Accretion in this phase is efficient, since the cross section is approximately given by the large Hill sphere radius. This results in high to very high accretion rates of 10^-7 to 10^-4 M_⊙/yr (Boley et al. 2010; Zhu et al. 2012). The object thus may grow quickly out of the planetary mass regime (Kratter et al. 2010). Independently of this issue, for the “hot” vs. “cold accretion” question, it is relevant in which way the gas is incorporated into the clump (in particular if it is accreted via a supercritical shock). The escape speed from the Hill sphere of a 10 M_X object at 50 AU is about 1.5 km s^-1. For comparison, in the calculation of the accretion onto the first (stellar) core in Commerçon et al. (2011), the shock velocity is about 2.4 km s^-1, which is comparable. Because we are in a shearing disk, both Boley et al. (2010) and Zhu et al. (2012), however, find that the relevant velocity with which they can reproduce numerically obtained accretion rates is not the free fall-velocity, but the Keplerian shear velocity across the Hill sphere. The velocity with which the newly accreted matter joins the clump is thus of order ΩR_H, where Ω is the Keplerian frequency. In the example, this corresponds to about 0.6 km s^-1. Whether or not this is sufficient to produce a supercritical shock must be evaluated for realistic velocities, pre-shock densities, and temperatures. If it is not, as we may suspect due to the rather low velocity, then the material accreted in this phase should be added at a high specific entropy. This leads to the second supposition: if a GI planet accretes most of its final mass before it undergoes the second collapse, and the accretion onto the AU sized clump is relatively gentle, we may again expect a “hot start”.

In the second sub-scenario, the clump accretes most of its final mass after it has undergone the second collapse, which leads to an object with a size of a few Jovian radii (Helled & Bodenheimer 2011). Now, the accreted gas hits the protoplanet at a much higher velocity of several 10 km s^-1 which is similar to the runaway gas accretion phase in the core accretion scenario. This leads to the third supposition: the larger the fraction of the final mass in a GI planet accreted after the clump has undergone the second collapse, the “colder” it should be, provided that the accretion shock in the post-collapse phase is supercritical. The effect that the higher the fraction of the final mass processed through the shock, the lower the initial luminosity, is seen in the simulations of Marley et al. (2007), where it causes more massive planets to be less luminous than their lighter counterparts (“luminosity inversion”).

This discussion (see also the related discussion in Spiegel & Burrows 2012) shows that with the currently limited knowledge both about the specific way matter is accreted and the involved thermodynamics in general, it is difficult to firmly associate a “hot/cold start” with a specific formation mechanism. Some of the general trends for the impact of formation on the initial luminosity of core accretion planets will be presented in a dedicated paper (Mordasini et al., in prep.).

8.2. “Hot start” and “cold start” models

In the present paper, we will only compare our results with already published studies, similar as for the radii in the previous section. Figure 10 shows (total) luminosity (including accretional luminosity) as a function of time for planets of 1, 2, 5, and 10 M_X.

The left panel shows the result using “cold start” boundary conditions, i.e., it is assumed that all accretional energy is radiated away at the shock, while no losses are assumed to occur in the “hot start (accreting)” case (Eq. (13)). The cold start, M = 1 M_X simulation is the same one as discussed in detail in Sect. 6, while the other masses correspond to the simulations shown in Fig. 6. We also plot in the figures the L(t) found by Burrows et al. (1997) and Baraffe et al. (2003), for comparison⁴. Both these models are classical “hot start” simulations. We associate the “t = 0” moment of these models with the moment when in our simulations, the planets reach their final mass. For the “cold start” simulation, this moment corresponds to the sharp drop of L particularly well visible for the 10 and 5 M_X cases, when L_acc vanishes.

8.2.1. “Hot start”: comparison with Burrows et al. (1997) and Baraffe et al. (2003)

Focusing first on the right panel with the “hot start (accreting)” models, we see a good agreement between our model and the other two. The differences between our (simpler, gray atmosphere) model and the two more complex ones is of the same order as the differences mutually between Burrows et al. (1997) and Baraffe et al. (2003). This is in agreement with Bodenheimer et al. (2000), who also find very good agreement of their gray atmosphere models and Burrows et al. (1997). It is, however, clear that the precise shape of L(t) (there, is for example, a small bump in our cooling curves when log L/L_⊙ ≈ −6.25) depends directly on the opacities, so that we expect our predictions to be somewhat less accurate during the long-term evolution.

The assumption of no radiative losses at the shock, but a gradual building up of the planet (as in our simulations) leads to results very similar to those of the classical “hot start” scenario, where one starts with a fully formed planet. The physical reason for the similarity is that in both cases, no entropy sink exists, and the ensuing short Kelvin-Helmholtz timescales result in the initial conditions being quickly “forgotten” in a rapid convergent evolution, as discussed at the end of this section. That means that core accretion with a radiatively completely inefficient shock leads to very similar planetary luminosities at young ages, as occurs in classical “hot start” simulations.

This result certainly underlines the importance of a detailed description of the shock structure, as already pointed out by Marley et al. (2007), because we see that the shock structure has an influence that is as important as the formation mechanism. Progress on the shock structure is an important task for future studies (see Commerçont et al. 2011, for a recent study of the stellar case).

The assumption of dl/dr = 0 in the envelope for the “hot start” scenario could lead to significant departures from the actual interior structure and the associated luminosity of the planets during the gas runaway accretion phase. This is because we cannot catch the effects of a possible deep radiative zone, as discussed in Sect. 3.6. A modified structure during the runaway accretion phase would also affect the luminosities at the beginning of the evolutionary phase, i.e., just after the end of formation.

The comparison of the luminosity at this moment as found in our simulations, and as assumed in Burrows et al. (1997) and Baraffe et al. (2003), indicates that we have luminosities that are similar (for the 10 M_X case) or higher (for lower masses). A higher luminosity means that the Kelvin-Helmholtz timescale is very short, so that the models converge on a timescale of only 10⁵−10⁶ yrs. Such a rapid convergence of “hot start” and “hotter start” models has been found by Baraffe et al. (2002, 2009; see also Marley et al. 2007, for a discussion). The effect of a variable luminosity on “hot start” objects, calculated as described in Sect. 3.4.1, will be studied in future work. This should allow a better description of the properties of very young “hot start” objects.

8.2.2. “Cold start”: comparison with Marley et al. (2007)

A comparison of the left and the right panel makes it immediately clear that we have recovered the important findings of Marley et al. (2007) (which were already visible in the simulations of Bodenheimer et al. 2000), namely that the luminosity of young planets formed according to the “cold start” assumption can be substantially lower, at least for high-mass planets.

Focusing on the left panel with the “cold start” simulations, we see that we have also recovered the most important effects found by Marley et al. (2007): a “luminosity inversion”, which means that the post-accretional luminosity (i.e., the luminosity shortly after the final mass of the planet is reached) is the highest for the lowest-mass planet. The picture is somewhat complicated by the fact that the low-mass planets start with the highest L but also cool the most quickly because their KH-timescale GM²/(RL) is the shortest. The reason for the “luminosity inversion” is that the higher the total mass, the higher the fraction of gas in the envelope that has been processed through the entropy-reducing shock. The total mass of the planets at the moment they start to collapse is about 87 M_⊕, with a core mass of 22 M_⊕ and an envelope of 65 M_⊕ (cf. Table A.2). This is equal for all planets (their evolution is identical up to the point where, for the lower mass planets, accretion is shut off). This means that for a 1 M_X planet, about 20% of its final envelope mass is accreted in the attached phase, without going through the entropy-reducing shock. The ratio is proportionally smaller for higher-mass planets, leading to the inversion.

Another identical result is that the higher the mass, the longer it takes until “hot start” and “cold start” luminosities become equal, a result that can be understood from the KH-timescales. For the 10 M_X planet, it takes several 100 Myr, similar to Marley et al. (2007).

Quantitatively, there are also some differences: we find post-accretional luminosities for the massive (5 and 10 M_X) planets that are of order log L/L_⊙ ≈ −4.1 to −4.3, while Marley et al. (2007) find log L/L_⊙ ≈ −5.7, i.e., fainter. The reason for this difference will be investigated in detail in Mordasini et al. (in prep.).

9. Summary

We have extended our formation model to a self-consistently coupled formation and evolution model for planets with a primordial H₂/He envelope. In this paper, we described the modification we made to the computational module that calculates the internal structure in order to make this extension possible. We then compared the results concerning planetary evolution with existing work. We found good agreement with more complex models. In a companion paper we describe further extensions and improvements relevant during the evolution of the planets, like an internal structure model for the solid core assuming a differentiated planet, which also gives realistic radii for planets without significant atmospheres or the radiogenic luminosity of the planet’s mantle. In addition, we discuss extensively our results on planetary radii as obtained in population synthesis calculations.

9.1. Jupiter: coupled in situ formation and evolution

We simulated the combined formation and evolution of Jupiter in the framework of classical core accretion models without migration and disk evolution (Pollack et al. 1996). We showed that the upgraded model with the new, simple, and rapid method to calculate evolutionary sequences leads to results that are in very good agreement with those of several other models (Lissauer et al. 2009; Burrows et al. 1997; Baraffe et al. 2003). We found that the nominal model for Jupiter leads to the formation of a planet that has properties at 4.6 Gyr that are in excellent agreement with observed values in terms of internal composition, radius, luminosity, and surface pressure-temperature profile.

9.2. Planetary radii

We studied the M-R diagram of giant planets and the influence of the core mass on the total radius. We found that the upgraded model yields planetary radii of giant planets that are in very good agreement with more sophisticated models (Fortney et al. 2007; Baraffe et al. 2008). We must assume that our cooling curves are probably less accurate for low-mass, core-dominated planets, possibly with heavily enriched envelopes. Nevertheless, we find in Paper II that the model yields a radius of a close-in, super-Earth planet with a tenuous H₂/He envelope that agrees to better than 10% with the more detailed simulations of Rogers et al. (2011). We thus can calculate planetary radii for a very wide range of planets, which is our goal for the planet population synthesis calculations for which the model is eventually intended (see Paper II).

9.3. Planetary luminosities

We used our upgraded model to study the luminosity of giant planets of different masses as a function of time. We made calculations both for the “cold start” and the “hot start (accreting)” scenario. For the “hot start (accreting)” models, we found cooling curves very similar to those in Burrows et al. (1997) and Baraffe et al. (2003), despite the fact that we use simple gray boundary conditions. In the “cold start” calculations, we recovered the result from Marley et al. (2007) that the luminosities of massive young planets are much smaller, at least for the chosen parameters. In a dedicated work (Mordasini et al., in prep.), we will revisit the luminosity of young Jupiters, using a large suite of combined formation and evolution calculations.

10. Conclusion

In earlier versions of the model, we calculated the internal structure of the gaseous envelope only during the attached, pre-gas runaway accretion phase. This is sufficient if one is only interested in the final mass of the planets, and thus sufficient for comparisons with planets found by the radial velocity method. Now we also calculated the structure during the gas runaway accretion phase (also the phase when the planet’s radius collapses) and during the evolutionary phase at constant mass over Gigayear timescales. With that, we now know all major quantities (mass, semimajor axis, radius, luminosity, composition) that characterize planets during their entire formation and evolution. This allows to compare our population synthesis models directly and consistently with results coming from all major observational techniques used to detect and characterize extrasolar planets (see Paper II for a comparison of the synthetic and actual M-R relationship, the predicted distribution of radii, and the comparison with Kepler).

Extrasolar planet research has entered an era in which massive amounts of observational data regarding very different sub-populations of planets are brought to us from different techniques. All these observations should be explained consistently by planet formation and evolution theory, but it is a difficult task to unite the different elements and observational constraints into one consistent global picture.

We think that a fruitful approach to this situation is to work with theoretical models that are able to make testable predictions in a consistent way for all important observational techniques. With the work presented in this paper, we have taken a step in this direction by testing and improving theoretical formation models using the wealth of data coming from observations.

---

Acknowledgments

We thank Christopher Broeg for dedicated comparison calculations. We thank Jérémy Leconte for answering many questions related to giant planet evolution. We thank Willy Benz, Gilles Chabrier, Isabelle Baraffe, Chris Ormel, Andrea Fortier, and Tristen Hayfield for helpful discussions. We thank an anonymous referee for questions that helped to illuminate the impact of our assumptions. Christoph Mordasini acknowledges the support as an Alexander von Humboldt fellow. This work was supported in part by the Swiss National Science Foundation and the European Research Council under grant 239605. Computations were made on the BATCHELOR cluster at MPIA.

References

Appendix A: Comparison calculation with Pollack et al. (1996)

In this appendix, we present a detailed comparison calculation with the classical simulation of Pollack et al. (1996). In particular, we simulate the cases J1, J1a, J1b, and J1c. Simulation J1 is the nominal case, while for the other cases, solid accretion gets (arbitrarily) shut off at specified times (see Fig. A.1). The consequence of this is a reduced core luminosity, which in turn leads to a faster gas accretion due to the reduced thermal support of the envelope, an effect described by Pollack et al. (1996).

This effect is relevant in our population synthesis calculations in two situations: during planet migration and/or the concurrent formation of several planets. The reason for the first situation is that with the updated (non-isothermal type I) migration model, protoplanets are often seen to first migrate outwards and then back in. On their way back in, they move into the region of the planetesimal disk, which they have cleared from planetesimals while migrating outwards (see Mordasini et al. 2011a). In the second situation, one planet can migrate into a region that has previously been cleared by another growing planet. In both situations, a faster gas accretion will result. From this, we see that migration and accretion can interact.

With this appendix, we want to check whether our simplified luminosity calculation reproduces the effect found by P96, and in general compare quantitatively our results with the published ones.

The J1 initial conditions in particular mean that the initial planetesimal surface density Σ_s,0 is 10 g/cm² as in Sect. 6. The grain opacities are now, however, at 100% of the interstellar value (f_opa = 1). The density of the core is constant at 3.2 g/cm³ as in Pollack et al. (1996), and the nebular boundary conditions are also the same as in this work. For the J1a, J1b, and J1c cases, we switch off planetesimal accretion at the same moments as in Pollack et al. (1996). The parameters used for the simulations are given in Table A.1. It is clear that the two simulations still cannot yield identical results, as they differ in some other aspects, like a different equation of state or a different model for the planetesimal-envelope interaction.

In Fig. A.1, we show the core mass, envelope mass, and total mass for the four models, and in Table A.2, we give the values of the most important quantities at specific moments in time. We first focus on the nominal J1 case. The crossover point, i.e., the moment when the core and the envelope mass are identical and equal to M_cr, is found to occur at t_cr = 7.49 Myr at M_cr = 16.62 M_⊕. This is very close to Pollack et al. (1996), with t_cr = 7.58 Myr and M_cr = 16.17 M_⊕. The very good agreement in t_cr is partially a chance result, because it is well known that changing, for example, just the opacity tables or the EOS can lead to variations of several million years. This is illustrated by the result of Hubickyj et al. (2005) that an updated version of the model used in Pollack et al. (1996) leads for identical initial conditions to a crossover times t_cr of only 6.07 Myr. At crossover, the core and gas accretion rate are 2.37 × 10^-6 and 1.19 × 10^-5 M_⊕/yr, which are both in very good agreement with Pollack et al. (1996).

The limiting gas accretion rate of 0.01 M_⊕/yr is reached at 8.02 Myr, about half a million years after crossover. The core mass is then about 22.90 M_⊕, and the envelope mass 75.8 M_⊕. In Pollack et al. (1996), a similar, but somewhat smaller gas accretion rate is reached 0.42 Myr after crossover, at M_Z = 21.5 M_⊕ and M_XY = 64.4 M_⊕. In contrast to Pollack et al. (1996), we also took the calculation through the detached end evolutionary phase. The second luminosity peak occurs shortly after the limiting gas accretion is reached, at 8.036 Myr, when the planet has already contracted to a radius of 4.49 R_X and grown to a mass of 267 M_⊕. The luminosity is then log (L/L_⊙) = −2.77. Lissauer et al. (2009) find peak luminosities for planets accreting with the same limiting gas accretion rate of log (L/L_⊙) ≈ −2.4 to −2.3.

The case J1a (red lines) is identical to case J1, except that at t = 1.5 Myr, the solid accretion rate is artificially ramped down to zero on a short timescale as in Pollack et al. (1996). This results in a faster growth of the envelope, i.e., an increase of the gas accretion rate. Our model thus reproduces this important effect. In this simulation, crossover occurs already at t_cr = 3.38 Myr and M_cr = 12.98 M_⊕. In Pollack et al. (1996), the values are 3.32 Myr, and 12.24 M_⊕, which shows that there is also quantitatively very good agreement. The results for the other two cases, J1b and J1c, are comparable, even though the differences in timescales are somewhat larger.

In summary, we see that there is also quantitatively a very good agreement between the Pollack et al. (1996) simulations and the ones presented here.

All Tables

All Figures