Semi-major axis–mass diagram

Figure 1 shows the planetary semi-major axis–mass (a–M) diagram. It is a classical observational constraint for population synthesis and is still one of the most important observational results. Explaining the structures seen in this plots is one of the goals of planetary population synthesis. The extreme diversity, but also the existence of certain structures in the a–M diagram is visible. For extrasolar planets, the mass–distance diagram has become a representation of similar importance as the Hertzsprung–Russell diagram for stellar astrophysics (Ida & Lin Reference Ida and Lin2004a).

Fig. 1. Semi-major axis–mass diagram of extrasolar planets. The colours show the observational technique that was used for the detection. The figure is not corrected for the various observational biases that favour for the radial velocity and the transit technique the detection of close-in, giant planets. The planets of the Solar System are also shown for comparison. Data from www.exoplanet.eu (Schneider et al. Reference Schneider, Dedieu, Le Sidaner, Savalle and Zolotukhin2011).

In the plot, one can distinguish several groups of planets. There are, for example, massive close-in planets without an equivalent in the Solar System. Such hot Jupiters are found around approximately 1% of solar-like stars (Marcy et al. Reference Marcy, Butler and Fischer2005; Howard et al. Reference Howard, Marcy and Johnson2010; Mayor et al. Reference Mayor, Marmier and Lovis2011). A class of extrasolar planets that has only been detected in the last few years thanks to the progress in the observational precision are low-mass planets with masses between 1 and 30 M _⊕ (Earth masses). These super-Earths and mini-Neptunes seem to be very abundant, since every second FGK star is found to have such a companion with a period of up to 100 days (Mayor et al. Reference Mayor, Marmier and Lovis2011). This result is at least in qualitative agreement with the analysis derived from the KEPLER mission, which also detects an extremely numerous population of planets with small radii ≲4R _⊕ (e.g. Howard et al. Reference Howard, Marcy and Bryson2012; Fressin et al. Reference Fressin, Torres and Charbonneau2013). Since hot Jupiters are much more easily detected by both the radial velocity and the transit method compared to low-mass (respectively small) planets, the number of such low-mass planet is underestimated in Fig. 1 which is not corrected for observational biases.

If we further inspect Fig. 1, we may ask whether it points to a statistically significant deficit of planets with a mass of approximately 40 Earth masses (so-called ‘planetary desert’, Ida & Lin Reference Ida and Lin2004b). This is, among many others, a very interesting question from a theoretical point of view that will be discussed later on (subsection ‘Radial velocity: the planetary initial mass function’, see also Mordasini et al. Reference Mordasini, Mayor and Udry2011).

The semi-major axis distribution and the planetary mass function are two fundamental one-dimensional (1D) statistical distributions that are encoded into the a–M diagram. These distributions are (besides of the radius distribution) of prime interest for statistical studies (e.g. Ida & Lin Reference Ida and Lin2004b; Thommes et al. Reference Thommes, Matsumura and Rasio2008; Mordasini et al. Reference Mordasini, Alibert and Benz2009a) and can be compared to theoretical results with Kolmogorov–Smirnov tests (Mordasini et al. Reference Mordasini, Alibert, Benz and Naef2009b). The planetary initial mass distribution is further addressed in subsection ‘Radial velocity: the planetary initial mass function’.

In the figure, there are also planets that were discovered by direct imaging. In this technique, one measures of course not the mass, but the luminosity of a planet. The conversion of luminosity into mass is model dependent and uncertain as shown by Marley et al. (Reference Marley, Fortney, Hubickyj, Bodenheimer and Lissauer2007) or Spiegel & Burrows (Reference Spiegel and Burrows2012). In Mordasini (Reference Mordasini2013) a new aspect was pointed out which is important for the conversion of luminosity into mass. It is found, perhaps surprisingly at first sight, which the post-formation luminosity of giant planets formed by core accretion depends significantly on the mass of the solid core. We address this finding in subsection ‘Direct imaging: luminosity at young ages’.

Mass–radius diagram

Figure 2 shows the observed mass–radius diagram of the extrasolar planets and compares it with three theoretical mass–radius relationships for planets with different bulk compositions (from Mordasini et al. Reference Mordasini, Alibert and Georgy2012b). The combination of measurements of the radius of a transiting planet (first by Henry et al. Reference Henry, Marcy, Butler and Vogt2000 and Charbonneau et al. Reference Charbonneau, Brown, Latham and Mayor2000) and its mass (via radial velocity) makes it possible to derive the mean density of the planet. In the past years, such combined measurements were made for many exoplanets, so that the planetary mass–radius diagram became knownFootnote ². It is an observational result of similar importance as the semi-major axis–mass diagram. One notes that as in the a–M diagram there is a large diversity, but that there are also clear trends leading, e.g. to regions in the plot that are not populated.

Fig. 2. The observed mass–radius relationship of the extrasolar planets (red points) together with theoretical mass–radius relationships for planets with an Earth-like composition, with an interior consisting purely of water ice, and for planets with a bulk composition roughly like Jupiter. The planets of the Solar System are also shown. Data from www.exoplanets.org Wright et al. (Reference Wright, Fakhouri and Marcy2011).

A surprising observational result was the discovery of numerous ‘inflated’ planets with radii much larger than Jupiter which was not predicted from standard planet evolution theory. It is now clear (Demory & Seager Reference Demory and Seager2011; Miller & Fortney Reference Miller and Fortney2011) that these bloated radii are related to the proximity of the planets to the host star (most currently known transiting planets have small orbital distances of a few 0.1 AU or less due to the decreasing geometrical transit probability with distance). The exact mechanism that leads to the large radii is still not completely understood. Several possible explanations have been put forward in the past, including tidal heating (Bodenheimer et al. Reference Bodenheimer, Laughlin and Lin2003), dissipation of stellar irradiation deep in the atmosphere (Guillot & Showman Reference Guillot and Showman2002), double diffusive convection (Chabrier & Baraffe Reference Chabrier and Baraffe2007), enhanced atmospheric opacities (Burrows et al. Reference Burrows, Hubeny, Budaj and Hubbard2007) and ohmic dissipation (Batygin et al. Reference Batygin, Stevenson and Bodenheimer2011).

The importance of the M–R diagram stems from its information content about the inner bulk composition of planets which is the first very basic geophysical characterization of a planet. This first characterization is found by the comparison of the observed mass and radius with theoretical models of the internal structure (e.g. Fortney et al. Reference Fortney, Marley and Barnes2007; Seager et al. Reference Seager, Kuchner, Hier-Majumder and Militzer2007). In the Solar System, there are three fundamental types of planets, namely, terrestrial, ice giant and gas giant planets. The imprint of the bulk composition on the radius is indicated by three theoretical lines in Fig. 2. Two lines show the theoretical mass–radius relationship for solid planets made of silicates and iron in a 2 : 1 ratio as for the Earth and for pure water planets, whereas the third line shows the M–R relationship for giant planets consisting of hydrogen/helium (H/He) and a solid core of about 10% in mass. Being able to understand and to reproduce this second fundamental figure besides the a–M diagram is another goal of statistical studies of planet formation and evolution. The reason for the importance of the M–R diagram for formation theory is the additional constraints on the formation process that cannot be derived from the mass–distance diagram alone.

An example are the observational constraints coming from the M–R diagram on the radial extent of orbital migration. Efficient inward migration brings ice-dominated, low-density planets from the outer parts of the disc close to the star. These planets can in principle be distinguished from denser planets consisting only of silicates and iron that have presumably formed more or less in situ in the inner hotter parts of the disc. Complications arise from the fact that the mass–radius relationship is degenerate in some parameter space (different bulk compositions can lead to an identical mass–radius relationship, see Rogers & Seager Reference Rogers and Seager2010). Therefore, spectroscopic measurements might be necessary to actually distinguish the two types of planets. Other key questions are: (1) What are the heavy element masses contained in giant planets? This is related to the question about the fundamental formation mechanism of giant planets, core accretion or gravitational instability. (2) Which planets can accrete and also keep primordial H/He envelopes (envelope evaporation)? Figure 2 shows that there are low-mass planets which likely contain important amounts of H/He. They therefore form a class of low-mass, low-density planets without counterpart in the Solar System. (3) Are there correlations between the planetary bulk composition and stellar properties (metallicity) as has been found by Guillot et al. (Reference Guillot, Santos and Pont2006) or Burrows et al. (Reference Burrows, Hubeny, Budaj and Hubbard2007)?

Direct imaging

The first technique besides transits that has recently yielded important new results for planet characterization is the direct imaging technique. The method is technically challenging due to the small angular separation of a very faint source (the planet) from a much brighter one (the host star). The number of planets detected by direct imaging is currently still low. But already these discoveries, like the planets around HR 8799 (Marois et al. Reference Marois, Macintosh and Barman2008) or β Pictoris (Lagrange et al. Reference Lagrange, Bonnefoy and Chauvin2010) have triggered numerous theoretical studies regarding their formation (e.g. Dodson-Robinson et al. Reference Dodson-Robinson, Veras, Ford and Beichman2009; Kratter et al. Reference Kratter, Murray-Clay and Youdin2010). Two points about these planets are interesting: their large semi-major axis and the fact that we directly measure the intrinsic luminosity at young ages in several infrared (IR) bands. Both quantities are important to understand the formation mechanism (core accretion or gravitational instability) and in particular the physics of the accretion shock occurring when the accreting gas hits the planet's surface during formation (e.g. Commerçon et al. Reference Commerçon, Audit, Chabrier and Chièze2011). If the gravitational potential energy of the accreting gas is radiated away, low entropy gas is incorporated into the planet, leading to a faint luminosity and small radius (so-called ‘cold start’, Marley et al. Reference Marley, Fortney, Hubickyj, Bodenheimer and Lissauer2007) while the accretion of high entropy material leads to a ‘hot start’ with a high luminosity and large radius (e.g. Burrows et al. Reference Burrows, Marley and Hubbard1997; Baraffe et al. Reference Baraffe, Chabrier, Barman, Allard and Hauschildt2003). Recently, Spiegel & Burrows (Reference Spiegel and Burrows2012) have shown that the different scenarios result in observable difference in the magnitudes of the young planets. The global model mainly discussed in this work (see Fig. 3) calculates the luminosity during both the formation and evolution phase with a self-consistent coupling. For young giant planets, this is a significant difference compared to purely evolutionary models, since this coupling is necessary to know the entropy in the envelope directly after formation and to correctly predict the luminosity at young ages. Since multi-band photometry can be used to estimate the metal enrichment of a planet and because new direct imaging instruments are currently becoming operational (Sphere and Gpi), it is important that future global models will include better descriptions of the gas accretion shock and better atmospheric models (cf. subsection ‘Atmosphere of the planet’).

Fig. 3. Schematic representation of the workflow in the population synthesis method (Ida & Lin Reference Ida and Lin2004a; Mordasini et al. Reference Mordasini, Alibert and Benz2009a). The 11 computational sub-models of the combined global planet formation and evolution model are based on the core accretion paradigm (see Alibert et al. Reference Alibert, Mordasini, Benz and Winisdoerffer2005a, Reference Alibert, Carron and Fortier2013, Mordasini et al. Reference Mordasini, Alibert and Georgy2012b, Reference Mordasini, Alibert, Klahr and Henningc).

Spectroscopy

Second, one of the most important aspects of the recent observational progress towards characterization are the spectra of a number of exoplanets transiting bright stars (e.g. Richardson et al. Reference Richardson, Deming, Horning, Seager and Harrington2007). The atmosphere represents a window into the composition of a planet and contains a multitude of clues to its formation history. The atmospheric composition depends on the composition of the host star, the nebula properties such as the temperature where the planet formed, the composition of the accreted gas and planetesimals, the size of the planetesimals and their (material) strength, the evolution of the distribution of the chemical compounds inside the planet and so on. Each migration and accretion history will result in a different atmospheric composition, as well as core and total heavy element mass. Jupiter, for example, is enriched in carbon by about a factor three relative to the Sun, while Uranus and Neptune are enriched by a factor ~30 (e.g. Guillot Reference Guillot1999). An enrichment relative to the sun is a natural prediction of the core accretion formation model but not of the competing gravitational instability model. For example, planet formation simulations based on the core accretion paradigm that reproduce (some of) the observed chemical composition such as the enrichment in carbon were presented in Gautier et al. (Reference Gautier, Hersant, Mousis and Lunine2001) and Alibert et al. (Reference Alibert, Mousis, Mordasini and Benz2005b). One should however note that other aspects of the composition of Jupiter are not straightforward to understand neither in the context of the core accretion nor gravitational instability model. In particular, the measurement that the enrichment of Jupiter both in highly volatile argon on the one hand and the more refractory sulphur on the other is similar (again by a factor of 2–4 relative to solar, Owen et al. Reference Owen, Mahaffy and Niemann1999) is difficult to explain in the context of conventional models of trapping of highly volatile gases in amorphous ice and a growth of Jupiter in a region similar to its current position. Several explanations to this conundrum have been put forward like a formation of the planetesimals that enriched Jupiter's envelope at lower temperatures and thus probably at significantly larger orbital distances (Owen et al. Reference Owen, Mahaffy and Niemann1999) or the incorporation of noble gases in the form of clathrate hydrates in crystalline ice (Gautier et al. Reference Gautier, Hersant, Mousis and Lunine2001).

An exoplanet that lately got a lot of attention as an example how spectroscopy constrains planet formation and evolution theory is GJ1214b (Charbonneau et al. Reference Charbonneau, Berta and Irwin2009). It has a radius of about 2.7 Earth radii and a mass of 6.5 Earth masses, so that it is a planet without counterpart in the Solar System between terrestrial planets and ice giants. The measured mass and radius is compatible with different internal compositions, like a rocky core with an H/He atmosphere or a planet dominated by water with a water vapour atmosphere (Rogers & Seager Reference Rogers and Seager2010). Such different compositions and the associated different mean molecular weight lead to different transmission spectra (e.g. Bean et al. Reference Bean, Désert and Kabath2011). Initially observed spectra showed that GJ1214b must either have an atmosphere with high-altitude clouds or hazes, or contain at least 70% water by mass (Bean et al. Reference Bean, Désert and Kabath2011; Berta et al. Reference Berta, Charbonneau and Désert2012). Recent precise observations could rule out cloud-free atmospheres even of high mean molecular mass (Kreidberg et al. Reference Kreidberg, Bean and Désert2014). A clear solar-composition atmosphere is excluded with very high confidence. Combined formation and evolution simulations that keep track of where planetesimals deposit their mass during impacts in the gaseous envelope (Mordasini et al. Reference Mordasini, Alibert and Benz2006) indicate that highly enriched atmospheres are a typical outcome for a planet such as GJ1214b with clear consequences for spectroscopy (Fortney et al. Reference Fortney, Mordasini and Nettelmann2013). But in general, despite the fact that some processes have been identified (Madhusudhan et al. Reference Madhusudhan, Mousis, Johnson and Lunine2011; Öberg et al. Reference Öberg, Murray-Clay and Bergin2011; Thiabaud et al. Reference Thiabaud, Marboeuf and Alibert2014), no self-consistently linked calculations have been made to date that keep track of the chemical composition of both the accreted gas and planetesimals during formation and directly predict the atmospheric composition and spectrum.

Global planet formation and evolution model

The most important element of population synthesis is a global planet formation and evolution model that establishes the link between disc and planetary properties. The sub-models of the global models typically used in population synthesis are described below in the section ‘Global planet formation and evolution models’ and are the main subject of this paper. For the different sub-models, already relatively well-established standard physical descriptions are employed if possible, which are the result of specialized models. Important simplifications are, however, often necessary for computational time restrictions. It is clear that the current global models (and often also the specialized ones) only provide a first, very rough approximation of the complex processes that actually govern planet formation. In this sense, it is likely that the global models will in future undergo important modifications, tracing in this way the developments in the field of planet formation theory. In order to still test the global models as far as possible, dedicated simulations are made for relatively well-known individual planetary systems, in particular the Solar System. But also some extrasolar systems can be studied individually, like, for example, the planetary system around HD 69830 with three Neptunian planets (Alibert et al. Reference Alibert, Baraffe and Benz2006; Lovis et al. Reference Lovis, Mayor and Pepe2006).

A global formation models should output as many observable quantities as possible, since in this case, one can use combined constraints from many techniques. The typical outputs are: the planetary mass, orbital distance and eccentricity (for comparison with radial velocity and microlensing), the radii for comparisons with transit observations, the intrinsic luminosity for comparison with discoveries made with direct imaging, and the atmospheric structure and composition for comparison with spectroscopy. The ability to compare concurrently and self-consistently with many different observational techniques (and therefore different sub-populations of planets) is crucial because global models suffer from the fact that they necessarily rely on a relatively high number of ill-constrained parameters (e.g. the viscosity parameter α, the sizes of the planetesimals, the initial radial slope of the solids etc.). The more comparisons are possible, the better these quantities can be determined, and the less likely it becomes that agreement of the model with observations is only obtained because of a sufficiently high number of unconstrained quantities.

Probability distributions for the initial conditions

The second central ingredients for population synthesis are sets of initial conditions. These sets of initial conditions are drawn in a Monte Carlo way from probability distributions. These probability distributions represent the different properties of protoplanetary discs and are derived as closely as possible from results of disc observations. At least three different fundamental disc properties have been considered in past population synthesis studies: The total disc (gas) mass (Beckwith & Sargent Reference Beckwith and Sargent1996; Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2010), the dust-to-gas ratio (assumed to be correlated with the stellar [Fe/H], Santos et al. Reference Santos, Israelian and Mayor2004; Fischer & Valenti Reference Fischer and Valenti2005) and the lifetime of the disc (Haisch et al. Reference Haisch, Lada and Lada2001; Mamajek Reference Mamajek2009; Fedele et al. Reference Fedele, van den Ancker, Henning, Jayawardhana and Oliveira2010). Additionally properties can be in the outer radius of the disc (controlled by the angular momentum of the collapsing cloud) or the initial radial slope of the solid surface density (Kornet et al. Reference Kornet, Stepinski and Rozyczka2001; Miguel et al. Reference Miguel, Guilera and Brunini2011). It is important to note that the derivation of such distributions for the initial conditions is not straightforward, introducing uncertainties in the final populations. For example, the total disc mass is typically found from sub-mm observations of cold dust at large semi-major axes (e.g. Andrews & Williams Reference Andrews and Williams2007; Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2010) to which gas is added at a ratio that is typical for the interstellar medium (usually a factor 100). The disc masses that are derived with this method can be up to one order of magnitude smaller compared to disc mass estimates based on the stellar accretion rate (Hartmann et al. Reference Hartmann, Calvet, Gullbring and D'Alessio1998). This discrepancy can result from a substantial growth of the dust to sizes to which submm observations are no more sensitive (e.g. Andrews & Williams Reference Andrews and Williams2007). Furthermore, the concept of an ‘initial’ protoplanetary disc is in any case questionable since discs form dynamically during the infall of the protostellar core. A more realistic model for the initial conditions would therefore include the early formation phase of the discs based on a simple infall model (e.g. Hueso & Guillot Reference Hueso and Guillot2005).

Synthetic detection bias

For a given set of initial condition, the formation model is used to calculate the final outcome, that is, the planetary system. This step is repeated many times leading to a population of synthetic planets (typically ~10 000 planets). Many of these synthetic planets could not be detected by current observational techniques because their mass (or radius) is too small. In order to make quantitative comparisons with the observations, one must therefore apply a synthetic observational detection bias. This leads to the sub-population of detectable synthetic planets. This group is then compared with a comparison sample of actual exoplanets. Depending on the observational technique, different detection biases are used, meaning that typically, different sub-populations are probed. It is clear that for quantitative comparisons, the selection bias of a given observational survey should be known as accurately as possible for this step. This makes that large, well characterized surveys like, e.g. the Harps high-precision survey (Mayor et al. Reference Mayor, Marmier and Lovis2011), the Kepler satellite (Borucki et al. Reference Borucki, Koch and Basri2011), and also microlensing surveys (e.g. Gould et al. Reference Gould, Dong and Gaudi2010; Cassan et al. Reference Cassan, Kubas and Beaulieu2012) are of particular interest.

Comparisons/statistical tests

For the comparison of the detectable synthetic planets and the actual planets various statistical methods can be used, like, for example, 2D Kolmogorov–Smirnov tests in the a–M or M–R plane. Other quantities that are tested are the detection frequency, and all different 1D distributions. It can further be studied if correlations exist between the initial conditions and the planet properties (Ida & Lin Reference Ida and Lin2004b; Mordasini et al. Reference Mordasini, Alibert, Benz, Klahr and Henning2012a), and if similar correlation exist in reality. The most important observed correlation is the one between the stellar metallicity and the frequency of giant planets. Giant planets are much more frequent around high metallicity stars (Gonzalez Reference Gonzalez1997; Santos et al. Reference Santos, Israelian and Mayor2004; Fischer & Valenti Reference Fischer and Valenti2005), a correlation that can be reproduced with formation models based on the core-accretion theory (Ida & Lin Reference Ida and Lin2004b;Mordasini et al. Reference Mordasini, Alibert, Benz and Naef2009b).

Depending on the results of this procedure, one can judge if the theoretical model is able to reproduce certain observed properties. In the ideal case, one single population should reproduce all observational constraints coming from many different techniques in a statistically significant fashion. In reality, there will always be differences between the model and the observations. This is, however, not an issue but instead the modus operandi of population synthesis, because the reasons for these differences are then analysed, so that various physical descriptions of the processes occurring during formation and evolution can be tested. This can have the consequence that an existing sub-model must be modified or even abandoned for being inconsistent with observations, or that new physical mechanisms must be added to the theoretical model. This is the fundamental process by which statistical studies improve the understanding of planet formation and evolution. Clearly, we currently still stand at the beginning of this process, even if the global formation models discussed below already concentrate a non-negligible amount of physics in one framework.

Predictions

In case of a satisfactory agreement between theory and observation (at least for a given aspect), one can return to the full underlying synthetic population and make predictions about planets or planetary properties that currently cannot be observed yet, like very low-mass planets and on the longer term, their habitability. The capability of population synthesis to produce output for direct falsifiability with future observations, that is, its predictive power is the strength of the method. Besides that, such predictions are also useful to estimate the yield of future instruments and surveys.

Vertical structure of the protoplanetary disc

A first sub-model calculates the structure and evolution of the gaseous protoplanetary disc (Lynden-Bell & Pringle Reference Lynden-Bell and Pringle1974; Papaloizou & Terquem Reference Papaloizou and Terquem1999; Alibert et al. Reference Alibert, Mordasini, Benz and Winisdoerffer2005a). The gaseous disc model yields the ambient properties in which the planets form. The ambient pressure and temperature serve as outer boundary conditions for the calculation of the structure of the gaseous envelope of the planets. The structure of the disc is also very important for the orbital migration of the protoplanets, since the direction and migration rate depend on the radial slopes of the temperature and the gas surface density (subsection ‘Orbital migration’). A good compromise for the numerical description of the discs that are in reality very complex, 3D structures driven by (magneto-) hydrodynamical processes (e.g. Flock et al. Reference Flock, Dzyurkevich, Klahr, Turner and Henning2011) is provided by α-viscosity models (Shakura & Sunyaev Reference Shakura and Sunyaev1973) that describe the disc as a rotating viscous fluid. In the model of Alibert et al. (Reference Alibert, Mordasini, Benz and Winisdoerffer2005a), the disc is described in the usual 1+1D approximation, that is, the disc has a vertical and a radial structure, but is assumed to be axisymmetric.

The vertical structure as a function of height z above the disc midplane is obtained by solving the coupled equations of hydrostatic equilibrium, energy conservation and energy transfer in the diffusion approximation for the radiative flux described by the equations (e.g. Papaloizou & Terquem Reference Papaloizou and Terquem1999)

(1)$$\displaystyle{1 \over {\rm \rho}} \displaystyle{{\partial P} \over {\partial z}} = - {\rm \Omega} ^2 z,\quad \displaystyle{{\partial F} \over {\partial z}} = \displaystyle{9 \over 4}{\rm \rho \nu \Omega} ^2, \quad F = \displaystyle{{ - 16{\rm \sigma} T^3} \over {3{\rm \kappa \rho}}} \displaystyle{{\partial T} \over {\partial z}}.$$

The first equations is the hydrostatic equilibrium in a thin disc at a pressure P and density ρ for a Keplerian frequency ${\rm \Omega} = \sqrt {GM_* /r^3} $ where G is the gravitational constant, M _* is the mass of the star and r is the orbital distance. The second and third equations state that the energy liberated by viscous dissipation (characterized by a turbulent viscosity ν) causes a flux of energy F that is transported via radiative diffusion to be radiated away at the surface of the disc. The other quantities in the equation are the temperature T and the opacity κ, whereas σ is the Stefan–Boltzmann constant.

The turbulent viscosity is calculated using the classical α-parametrization of Shakura & Sunyaev (Reference Shakura and Sunyaev1973) as (c _s is the sound speed) ${\rm \nu} = {\rm \alpha} c_{\rm s}^2 {\rm /\Omega} $ in which all the complex physics about the angular transport processes in protoplanetary discs is hidden. The current understanding about these processes is actually still quite poor; for a recent review, see Turner et al. (Reference Turner, Fromang and Gammie2014). In the context of the vertical disc structure, we note that the α-parametrization of viscosity leads to viscous heating that is concentrated around the disc midplane. This is a consequence of the vertically constant α. Direct radiation-magnetohydrodynamic simulation where the turbulence is due to the magneto-rotational instability (MRI, Balbus & Hawley Reference Balbus and Hawley1991) lead in contrast to a less concentrated heating, and thus to a different vertical temperature profile than in an α model (Flock et al. Reference Flock, Fromang, González and Commerçon2013). Such a difference is particularly strong when the turbulence is primarily occurring above the midplane as it is the case in an MRI active layer above a midplane deadzone (e.g. Dzyurkevich et al. Reference Dzyurkevich, Turner, Henning and Kley2013). This (and other differences, see Turner et al. Reference Turner, Fromang and Gammie2014) show that the α-parametrization of the turbulence might not be sufficient for a realistic description of protoplanetary discs.

Viscous dissipation is the dominant heating mechanism in the inner parts of the disc (e.g. Chambers Reference Chambers2009). At larger orbital distances, the irradiation of the host star becomes dominant. This heating source can be incorporated into the surface boundary condition of the temperature structure T _s as (Barrière-Fouchet et al. Reference Barrière-Fouchet, Alibert, Mordasini and Benz2012)

(2)$$T_{\rm s}^4 = T_{{\rm s},{\rm vis}}^4 + T_{{\rm s},{\rm irr}}^4, $$

where T _s,vis would be the temperature due to viscous heating only. The temperature due to the stellar irradiation is approximated as (Hueso & Guillot Reference Hueso and Guillot2005)

(3)$$T_{{\rm s},{\rm irr}} = T_* \left[ {\displaystyle{2 \over {3{\rm \pi}}} \left( {\displaystyle{{R_*} \over r}} \right)^3 + \displaystyle{1 \over 2}\left( {\displaystyle{{R_*} \over r}} \right)^2 \left( {\displaystyle{H \over r}} \right)\left( {\displaystyle{{{\rm d}\,{\rm ln\;} H} \over {{\rm d}\,{\rm ln}\; r}} - 1} \right)} \right]^{1/4}, $$

where T _*, R _* and H are the stellar temperature, radius and the vertical pressure scale height, respectively. Barrière-Fouchet et al. (Reference Barrière-Fouchet, Alibert, Mordasini and Benz2012) set the flaring angle d ln H/d ln r for simplicity to the equilibrium value of 9/7 (Chiang & Goldreich Reference Chiang and Goldreich1997) which, however, means that the flaring angle and the possible effects of shadowing are not described in a self-consistent way. The solution of the vertical structure equations yield the disc midplane temperature and pressure, the vertical scale height H and the vertically averaged viscosity.

Radial structure of the protoplanetary (gas) disc

The evolution of gas surface density Σ as a function of distance r and time t is described by the classical viscous evolution equation of Lynden-Bell & Pringle (Reference Lynden-Bell and Pringle1974) (first term on the right-hand side) supplemented by the effects of mass loss by photoevaporation ${\rm \dot \Sigma} _w (r)$ and accretion onto the planet ${\rm \dot \Sigma} _{{\rm planet}} (r)$

(4)$$\displaystyle{{\partial \Sigma} \over {\partial t}} = \displaystyle{1 \over r}\displaystyle{\partial \over {\partial r}}\left[ {3r^{1/2} \displaystyle{\partial \over {\partial r}}\left( {r^{1/2} {\rm \nu \Sigma}} \right)} \right] - {\rm \dot \Sigma} _{\rm w} (r) - {\rm \dot \Sigma} _{{\rm planet}} (r).$$

The mass loss due to photoevaporation (for a recent review on disc dispersal mechanisms, see Alexander et al. Reference Alexander, Pascucci, Andrews, Armitage and Cieza2013) has two origins (Mordasini et al. Reference Mordasini, Alibert and Georgy2012b): first external photoevaporation due to far-ultraviolet (F _UV) radiation coming from massive stars in the vicinity of the host star (e.g. Matsuyama et al. Reference Matsuyama, Johnstone and Hartmann2003). This drives a wind approximately outside of a gravitational radius R _g,I given as

(5)$$R_{{\rm g},{\rm I}} = \displaystyle{{GM_*} \over {c_{{\rm s,I}}^2}}, $$

where c _s,I is the speed of sound in the heated layer (T _I≈1000 K) of dissociated neutral hydrogen. For a solar-like star, R _g,I is about 140 AU. The reduction of the surface density for a disc with total radius r _max is given as

(6)$${\rm \dot \Sigma} _{w,{\rm ext}} = \left\{ {\matrix{ 0 & {{\rm for}\; r \lt {\rm \beta} _{\rm I} R_{g,{\rm I}},} \cr {\displaystyle{{\dot M_{{\rm wind},{\rm ext}}} \over {{\rm \pi} (r_{{\rm max}}^2 - {\rm \beta} ^2 R_{g,I}^2 )}}} & {{\rm otherwise},} \cr}} \right.$$

that is, it occurs outside of an effective gravitational radius β_IR _g,I. The total rate $\dot M_{{\rm wind},{\rm ext}} $ is an input parameter, and one of the Monte Carlo random variables in a population synthesis calculation. Its distribution function is chosen in a way that the resulting distribution of lifetimes of the synthetic discs is in agreement with the observed distribution of lifetimes of actual protoplanetary discs (Haisch et al. Reference Haisch, Lada and Lada2001). Physically, it depends on the number of and distance to massive stars in the vicinity of the host star (Adams et al. Reference Adams, Hollenbach, Laughlin and Gorti2004).

The second contribution to ${\rm \dot \Sigma} _{\rm w} (r)$ is due to the internal photoevaporation driven by the extreme ultraviolet (EUV) radiation of the host star itself. This ionizing radiation heats the surface layers to a temperature of ~10⁴ K which launches a wind with a velocity c _s,II outside around β_IIR _g,II which corresponds to approximately 7 AU for a solar-like star. The mass loss rate is

(7)$${\rm \dot \Sigma} _{{\rm w},{\rm int}} = \left\{ {\matrix{ 0 & {{\rm for}\; r \lt {\rm \beta} _{{\rm II}} R_{{\rm g},{\rm II}},} \cr {2c_{{\rm s},{\rm II}} n_0 (r)m_{\rm H}} & {{\rm otherwise},} \cr}} \right.$$

where m _H is the mass of the ionized hydrogen. The density of ions at the base of the wind n ₀(r) as a function of distance is approximately (Hollenbach et al. Reference Hollenbach, Johnstone, Lizano and Shu1994)

(8)$$n_0 (r) = n_0 (R_{14} )\left( {\displaystyle{r \over {{\rm \beta} _{{\rm II}} R_{{\rm g},{\rm II}}}}} \right)^{ - 5/2}, $$

which means that most of the wind is originating close to the effective gravitational radius. The density n ₀ at a normalization radius R ₁₄ is found from radiation-hydrodynamic simulations (Hollenbach et al. Reference Hollenbach, Johnstone, Lizano and Shu1994). It increases with the square root of the ionizing photon luminosity of the central star.

The last term in the master equation for the evolution of Σ is calculated by assuming that the gas accreted by a planet at a rate $\dot M_{XY} $ is removed from an annulus around it with a width equal to the planet's Hill sphere so that

(9)$${\rm \dot \Sigma} _{{\rm planet}} = \displaystyle{{\dot M_{XY}} \over {2{\rm \pi} aR_{\rm H}}}. $$

The Hill sphere radius R _H of a planet of mass M and semi-major axis a around a star of mass M _* is given as (M/3M _*)^1/3a. One finds that this term is important for the global evolution of the disc only for massive planets undergoing gas runaway accretion.

A surface density at time t=0 must be specified as initial conditionFootnote ³. Early population synthesis calculations (Ida & Lin Reference Ida and Lin2004a; Mordasini et al. Reference Mordasini, Alibert, Benz and Naef2009b) typically used an initial radial profile inspired by the minimum mass Solar nebula (MMSN) (Weidenschilling Reference Weidenschilling1977; Hayashi Reference Hayashi1981) which is a power-law for Σ falling as r^{− 3/2}. Recent sub-millimetre observation of protoplanetary discs (e.g. Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2010) rather indicate profiles like

(10)$${\rm \Sigma} (r,t = 0) = {\rm \Sigma} _0 \left( {\displaystyle{r \over {R_0}}} \right)^{ - {\rm \gamma}} \exp {\rm \;} \left[ { - \left( {\displaystyle{r \over {R_{\rm c}}}} \right)^{2 - {\rm \gamma}}} \right],$$

where the mean values of the parameters γ and R _c are approximately 0.9 and 30 AU, respectively (Andrews et al. Reference Andrews, Wilner, Hughes, Qi and Dullemond2010). Such profiles were used as initial conditions in more recent syntheses like Mordasini et al. (Reference Mordasini, Alibert and Georgy2012b) and Alibert et al. (Reference Alibert, Carron and Fortier2013) even though the parameters were strictly speaking determined at distances outside of ~20 AU, that is outside of the main planet formation region. The initial surface density Σ₀ at a normalization distance R ₀ determines the initial total mass and is another Monte Carlo random variable in population syntheses.

Figure 4 shows a typical example of the evolution of the gas surface density under the action of viscosity and photoevaporation (but without an embedded planet). The inner disc edge is fixed to 0.1 AU, whereas the outer radius is free to spread or shrink. The initial mass of the disc is similar to the MMSN.

Fig. 4. Example of the evolution of the gas surface density Σ as a function of time and radius a. The uppermost line shows a state close to the beginning of the simulation and is similar to the initial condition given by equation (10). Subsequent profiles are plotted 2×10⁴ years apart (figure from Mordasini et al. Reference Mordasini, Alibert and Georgy2012b).

The global models of Ida & Lin (Reference Ida and Lin2004a) and Hellary & Nelson (Reference Hellary and Nelson2012) use descriptions of the protoplanetary gas disc that are more parameterized. Typically, one assumes a power-law dependency of the gas surface density as a function of radius all the time (i.e. not only as initial condition) and an exponential decrease on a uniform timescale τ_disc. For a power-law exponent p _Σ, the gas surface density is then given as

(11)$${\rm \Sigma} (r,t) = {\rm \Sigma} _0 \left( {\displaystyle{r \over {R_0}}} \right)^{ - p_{\rm \Sigma}} {\rm e}^{ - t/{\rm \tau} _{{\rm disc}}}. $$

The radial temperature profile of the disc is described in an analogous way, which is as a power law with a constant exponent p _T. These exponents must be carefully chosen since the migration of low-mass planets particularly in non-isothermal parts of the disc depends on a sensitive way on the local values of p _Σ and p _T (Lyra et al. Reference Lyra, Paardekooper and Mac Low2010; Kretke & Lin Reference Kretke and Lin2012; Dittkrist et al. Reference Dittkrist, Mordasini, Klahr, Alibert and Henning2014, subsection ‘Orbital migration’). In 1+1D α models, the exponents are functions of distance and time (they depend via the vertical structure on opacity transitions), which leads in particular to convergence zones that act as traps for migrating planets (e.g. Paardekooper et al. Reference Paardekooper, Baruteau, Crida and Kley2010, see also subsection ‘Orbital migration’).

This does, however, not mean that 1+1D α models are already sufficient to catch all important effects found in 3D magnetohydrodynamical simulations. Regarding the general evolution and structure of the disc itself, an important assumption in the models presented here is the one of a radially constant α. The results obtained from the radial disc structures (such as the surface density or temperature) depend significantly on the magnitude and assumed radial profile of α (e.g. Bell et al. Reference Bell, Cassen, Klahr and Henning1997; Kretke & Lin Reference Kretke and Lin2010). It is clear that in actual protoplanetary discs, various processes occurring at different radii like the strength of the MRI, the radial and vertical extent of a possible low-turbulence deadzone (e.g. Dzyurkevich et al. Reference Dzyurkevich, Turner, Henning and Kley2013), or the occurrence of layered accretion (Gammie Reference Gammie1996) produce radial variations in the effective α. The use of one constant α must therefore be seen as a strong simplification.

Regarding planet migration, changes in the effective α can significantly affect planetary migration tracks by slowing or stopping the migration in planet traps (e.g. Matsumura et al. Reference Matsumura, Pudritz and Thommes2007). In the context of global planet formation models, this was studied in Ida & Lin (Reference Ida and Lin2008b) and Hasegawa & Pudritz (Reference Hasegawa and Pudritz2013). Furthermore, Uribe et al. (Reference Uribe, Klahr, Flock and Henning2011) find additional torques in 3D magnetohydrodynamical simulations where the turbulence is given by the MRI. They are directly related to the presence of the magnetic fields so that they cannot be described in a 1+1D α model.

The output of the disc structure model is used in several other sub-models. The ‘atmosphere’ model, for example, needs the midplane pressure and temperature in the nebula to calculate the boundary conditions for the planet's interior. The ‘migration’ models need, as mentioned, the surface density and temperature gradients, the surface density itself, the vertical scale height H and the viscosity to calculate the planetary migration rate.

Disc of solids (planetesimals, fragments)

A third sub-model describes the structure and evolution of the disc of small solid bodies (planetesimals) from which the bigger protoplanets grow. The division in ‘small’ and ‘big’ bodies is made since runaway growth leads to a bimodal size distribution (Weidenschilling et al. Reference Weidenschilling, Spaute, Davis, Marzari and Ohtsuki1997). The model of the disc of solids yields as a function of time and orbital distance the typical size s _solid (or size distribution) of the bodies, their dynamical state (random velocities or, equivalently, inclinations and eccentricities) and their surface density Σ_solid. These quantities control the core accretion rate of a forming protoplanet.

In a protoplanetary disc, the solids are initially in the form of tiny dust grains that grow in time to form kilometre-sized planetesimals (by a process that is currently debated, e.g. Cuzzi et al. Reference Cuzzi, Hogan and Shariff2008). Destructive collisions can again reduce the size of the bodies. Additionally, solid bodies drift radially through the disc at a rate that depends among others things on the size of the bodies. The dynamical state of the small bodies is influenced by several non-linear processes such as dynamical friction, viscous stirring or damping by gas drag (e.g. Ida & Makino Reference Ida and Makino1993; Ormel et al. Reference Ormel, Dullemond and Spaans2010). The interaction of all this processes make that the evolution of the disc of solids is in reality a very complex problem. Some aspects (in particular the viscous stirring by the protoplanets) were recently included in global models by Fortier et al. (Reference Fortier, Alibert, Carron, Benz and Dittkrist2013). In earlier population syntheses, Mordasini et al. (Reference Mordasini, Alibert, Benz and Naef2009b) used a very simple model for the solids where the disc of planetesimals only changes due to the accretion and ejection of planetesimals by the protoplanet. The size of the planetesimals is uniform in time and space (usually 100 km). It is furthermore assumed that the accretion and ejection (for massive protoplanets) homogeneously reduces the surface density of planetesimals in the feeding zone of the planet, which is taken to have a radial width equal to W _feed=B _LR _H where B _L≈4–5 (Lissauer Reference Lissauer1993). This ignores the potential formation of a gap in the disc of planetesimals (Shiraishi & Ida Reference Shiraishi and Ida2008). The decrease of the planetesimal surface density is then simply (Thommes et al. Reference Thommes, Duncan and Levison2003)

(12)$$\displaystyle{{{\rm d}\Sigma _{{\rm solid}}} \over {{\rm dt}}} = - \displaystyle{{(3M_{\rm *} )^{1/3}} \over {6{\rm \pi} a^2 B_{\rm L} M^{1/3}}} \displaystyle{{{\rm d}M_{Z,{\rm tot}}} \over {{\rm dt}}},$$

where dM _Z,tot/dt denotes the total change of surface density of planetesimals due to the presence of the protoplanet, that is, the sum of the rate at which the planet accretes and ejects planetesimals. For the initial condition it is assumed that the surface density of planetesimals is proportional to the initial surface density of the gas times the dust-to-gas ratio f _D/G as

(13)$${\rm \Sigma} _{{\rm solid}} (r,t = 0) = f_{{\rm D}/{\rm G}} f_{{\rm R}/{\rm I}} {\rm \Sigma} (r,t = 0).$$

The dust-to-gas ratio is another Monte Carlo random variable in population synthesis calculation, and set under the simplifying assumption that it is directly given by the stellar [Fe/H] as f _D/G=f _D/G,⊙×10^[Fe/H] where f _D/G,⊙ is the solar metal fraction. The factor f _R/I represents the increase of the surface density due to the condensation of water ice outside of the iceline, and set to (Hayashi Reference Hayashi1981)

(14)$$f_{{\rm R}/{\rm I}} = \left\{ {\matrix{ 1 & {{\rm for}\,\; T_{{\rm mid}} (t = 0) \lt 170K\; \; (a \gt a_{{\rm ice}} ),} \cr {1/4} & {{\rm otherwise}\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \,(a \lt a_{{\rm ice}} ),} \cr}} \right.$$

which means that the initial midplane temperature profile determines the location of the iceline a _ice. More recent works on the composition of the solar nebula (like Lodders Reference Lodders2003) indicate that the change of the surface density at the iceline is probably rather about a factor of 2 than 4 as assumed by Hayashi (Reference Hayashi1981). Figure 5 shows the initial surface density and its evolution due to a growing protoplanet. The initial profile simply traces the initial profile of the gas (equation 10) except for the sudden increase at the iceline. The planet initially forms at about 15 AU, and then migrates inwards depleting the planetesimal disc down to a distance of about 5 AU. Due to migration, the expansion of the feeding zone occurs mainly towards the interior part of the disc. It occurs on a shorter timescale than for in situ formation, which can reduce the overall formation timescale of a giant planet (Alibert et al. Reference Alibert, Mordasini and Benz2004).

Fig. 5. Surface density of planetesimals. The black line shows the initial planetesimal surface density in the entire disc. The increase at the iceline is visible. The red lines show the reduction of Σ_solid due to the presence of an accreting and migrating planet.

Other global planet formation models (Ida & Lin Reference Ida and Lin2004a; Hellary & Nelson Reference Hellary and Nelson2012) use similar prescriptions for the initial surface density distribution. The model of Ida & Lin (Reference Ida and Lin2008b) is more complex because it takes into account that a magnetohydrodynamically inactive dead zone in the protoplanetary disc can lead to a local density maximum of planetesimals in the vicinity of the iceline. But it is clear that all global models treat to date the processes occurring during the early condensation of the solids and the subsequent growth phase from dust to planetesimals in a simplistic way (or not at all).

A complete discussion of the key physical processes associated with the formation and growth of planetesimals is beyond the scope of this work (see Johansen et al. Reference Johansen, Blum and Tanaka2014 for a recent review). There is currently no generally accepted complete theory for the growth from μm-sized dust to sizes where gravitational focusing sets in, allowing efficient runaway growth (of order of 1–1000 km, Ormel & Okuzumi Reference Ormel and Okuzumi2013). An incomplete list of important issues that must be addressed includes: (i) the basic mode of growth which is two-body coagulation or self-gravitational instability (Johansen et al. Reference Johansen, Oishi and Low2007; Birnstiel et al. Reference Birnstiel, Dullemond and Brauer2010); these are not mutually exclusive processes (Johansen et al. Reference Johansen, Blum and Tanaka2014), (ii) the mode of concentration of small bodies in the turbulent disc on large (Johansen et al. Reference Johansen, Klahr and Henning2006) or small scales (Cuzzi et al. Reference Cuzzi, Hogan and Shariff2008), (iii) the strong dependency of planetesimal formation rates on the dust-to-gas ratio (e.g. Johansen et al. Reference Johansen, Youdin and Mac Low2009; Cuzzi et al. Reference Cuzzi, Hogan and Bottke2010) so that planetesimals may form preferentially at certain disc locations rather than uniformly (as assumed in the global models), (iv) the material properties of dust grains and planetesimals and their behaviour during mutual collisions which are very poorly constrained (Güttler et al. Reference Güttler, Blum, Zsom, Ormel and Dullemond2010; Leinhardt & Stewart Reference Leinhardt and Stewart2012), (v) the possibility that planetesimal formation may be inefficient or happen over an extended period of time (several Myrs) based on observations of chondrite parent bodies (Bizzarro et al. Reference Bizzarro, Baker and Haack2004), in contrast to the assumptions in the global models, (vi) the possibility that much of the solid mass existed in mm-to-m size particles that drift radially, significantly modifying the surface density profile, the solid-to-gas and rock-to-ice ratios over time (e.g. Cuzzi & Zahnle Reference Cuzzi and Zahnle2004; Ciesla & Cuzzi Reference Ciesla and Cuzzi2006), also in contrast to the simplifications in the global models, (vii) the mechanism that planetesimal–planetesimal collisions may rapidly convert most of the solid mass into smaller objects (e.g. Inaba et al. Reference Inaba, Wetherill and Ikoma2003; Chambers Reference Chambers2008) so that at least three types of bodies must be included in the solid disc model (embryos, planetesimals and fragments, see Ormel & Kobayashi Reference Ormel and Kobayashi2012) and finally that the timescale and even existence of runaway growth to form protoplanets depends sensitively on the turbulence level in the disc due to eccentricity excitation by turbulent density fluctuations (e.g. Nelson Reference Nelson2005; Ida et al. Reference Ida, Guillot and Morbidelli2008; Ormel & Okuzumi Reference Ormel and Okuzumi2013).

Once a clearer and more unified picture of planetesimal formation arises, simplified version of it will again be included in the global models to understand the observational consequences. Some steps towards this can be found in the aforementioned works and, e.g. Kornet et al. (Reference Kornet, Wolf and Rozyczka2007); Carter-Bond et al. (Reference Carter-Bond, O'Brien and Lauretta2010); Birnstiel et al. (Reference Birnstiel, Klahr and Ercolano2012) or Chambers (Reference Chambers2014).

Planetary core accretion rate

The solid core of a protoplanet grows by the accretion of background planetesimals and by the collision with other protoplanets in the case that several planets form concurrently. The simplest way to describe the collisional growth of protoplanet due to the accretion of planetesimals is a Safronov-type rate equation (Safronov Reference Safronov1969) for the accretion rate of the core of mass M _Z

(15)$$\displaystyle{{{\rm d}M_Z} \over {{\rm dt}}} = {\rm \Omega \Sigma} _{{\rm solid}} R_{{\rm capture}}^2 F_{\rm G} (e,i).$$

The capture radius R _capture for planetesimals is larger than the core radius due to the presence of the gaseous envelope. It is calculated in the ‘infall’ model, while the Σ_solid is an output of the model discussed in the final section. The gravitational focussing factor F _G(e,i), that takes into account three-body effects (planet, planetesimal and star), is given, for example, by Greenzweig & Lissauer (Reference Greenzweig and Lissauer1992). Its magnitude depends on the dynamical state of the planetesimal swarm. This is the key quantity determining dM _Z/dt since it gives raise to different growth regimes such as runaway, oligarchic or orderly growth; see, e.g. Rafikov (Reference Rafikov2003). In the population syntheses of Mordasini et al. (Reference Mordasini, Alibert, Benz and Naef2009b) and (Reference Mordasini, Alibert and Georgy2012b), the original expressions of Pollack et al. (Reference Pollack, Hubickyj and Bodenheimer1996) are used to estimate the eccentricities and inclinations. These equations predict low random velocities, allowing rapid growth, which is probably difficult to achieve for 100 km planetesimals in reality (e.g. Thommes et al. Reference Thommes, Duncan and Levison2003). An updated description with a more detailed description of the interactions between protoplanet and planetesimals can be found in Fortier et al. (Reference Fortier, Alibert, Carron, Benz and Dittkrist2013), but it is likely that the prescriptions for the core accretion rate will be further significantly modified in view of new result of specialized models: as understood recently (e.g. Ormel & Klahr Reference Ormel and Klahr2010; Lambrechts & Johansen Reference Lambrechts and Johansen2012; Morbidelli & Nesvorny Reference Morbidelli and Nesvorny2012; Chambers Reference Chambers2014) the accretion rate of decimetre-sized pebbles instead of km-sized planetesimals can be very high, especially also at larger semi-major axes. Such high rates are necessary to form a ~10 M _⊕ core during the presence of the nebula in a MMSN disc which is otherwise far from simple to achieve; see, e.g. Ormel & Okuzumi (Reference Ormel and Okuzumi2013) or Chambers (Reference Chambers2014). The high accretion rates are due to a strongly increased capture cross-section for small particles due to drag with the disc gas during the encounter with the protoplanetary core. Understanding the global predictions of pebble accretion when coupled to global planet formation models and population syntheses is an important task for future studies.

The formation model of Ida & Lin (Reference Ida and Lin2004a) also uses a Safronov-type rate equation but describes the dynamical state of the swarm based on more recent results for the oligarchic growth regime (Ida & Makino Reference Ida and Makino1993; Ormel et al. Reference Ormel, Dullemond and Spaans2010). As they assume smaller planetesimals (km-sized instead of 100 km), the resulting core growth timescales are nevertheless relatively similar (Mordasini et al. Reference Mordasini, Alibert and Benz2009a). The global model of Hellary & Nelson (Reference Hellary and Nelson2012) is much more detailed because the solid accretion rate is found by direct N-body integration of a few ten protoplanets plus a several thousand planetesimals. The downside of this is that the associated long computation time makes the calculation of statistically sufficiently large population of synthetic planets difficult, at least for the moment.

Internal structure of the planetary gas envelope

The sub-models discussed up to this point describe the protoplanetary disc and the growth of the core. The ‘envelope’ model (together with the internal structure model of the core) deals with the protoplanet itself, calculating the internal 1D (spherically symmetric) radial structure of the gaseous envelope (H/He) of the planet. During the formation phase, the model in particular yields the gas accretion rate $\dot M_{XY} $. Low-mass cores can only bind tenuous atmospheres, while cores more massive than roughly ten Earth masses can trigger rapid runaway gas accretion, so that a giant planet forms. After the formation phase the long-term evolution, that is, the contraction and cooling at constant mass are calculated, which yield the radius and intrinsic luminosity of a planet (recent reviews on planetary internal structures and the temporal evolution can be found in Baraffe et al. (Reference Baraffe, Chabrier, Fortney and Sotin2014) and Chabrier et al. Reference Chabrier, Johansen, Janson and Rafikov2014). The internal structure is found by integrating the standard equations for planetary interiors (e.g. Bodenheimer & Pollack Reference Bodenheimer and Pollack1986) which are the equations of mass conservation, hydrostatic equilibrium, energy conservation and energy transfer:

(16)$$\displaystyle{{\partial m} \over {\partial r}} = 4{\rm \pi} r^2 {\rm \rho,} \quad \displaystyle{{\partial P} \over {\partial r}} = - \displaystyle{{Gm} \over {r^2}} {\rm \rho,} $$

(17)$$\displaystyle{{\partial l} \over {\partial r}} = 4{\rm \pi} r^2 {\rm \rho} \left( {{\rm \varepsilon} - P\displaystyle{{\partial V} \over {\partial t}} - \displaystyle{{\partial u} \over {\partial t}}} \right),\quad \displaystyle{{\partial T} \over {\partial r}} = \displaystyle{T \over P}\displaystyle{{\partial P} \over {\partial r}}\nabla (T,P).$$

In these equations, r is the radius as measured from the planetary centre, m the gas mass inside r, ρ, P, T, u the gas density, pressure, temperature and specific internal energy, V the specific volume 1/ρ and t the time. The gradient ∇=d ln T/d ln P is either the radiative or the adiabatic gradient, whichever is shallower. The adiabatic gradient is directly given by the equation of state, whereas the radiative gradient is

(18)$$\nabla _{{\rm rad}} = \displaystyle{3 \over {64{\rm \pi \sigma} G}}\displaystyle{{{\rm \kappa} lP} \over {T^4 m}}.$$

In this equation, σ is the Stefan–Boltzmann constant, whereas κ is the Rosseland mean opacity. In the outer layers of the envelope, it is dominated by the opacity due to small grains suspended in the gas (Podolak Reference Podolak2003; Movshovitz et al. Reference Movshovitz, Bodenheimer, Podolak and Lissauer2010; Cuzzi et al. Reference Cuzzi, Estrada and Davis2014, subsection ‘Impact of grain opacity on the planetary radius distribution’). These are the same equations as for stellar interiors with the difference that the luminosity is not due to nuclear burning of hydrogen but due to the contraction and cooling of the gas, and the energy deposition by impacting planetesimals that is found with the ‘infall’ sub-model. The burning of deuterium in sufficiently massive objects can also be included in the energy source term ε (Mollière & Mordasini Reference Mollière and Mordasini2012; Bodenheimer et al. Reference Bodenheimer, D'Angelo, Lissauer, Fortney and Saumon2013). Note that the models of Mordasini et al. (Reference Mordasini, Alibert, Klahr and Henning2012c) and Alibert et al. (Reference Alibert, Carron and Fortier2013) simplify the energy equation by assuming a luminosity that is radially constant.

Figure 6 shows an example of the radial envelope structure inside a growing giant planet that forms in situ at 5.2 AU, with initial conditions similar to the classical J1 simulations of Pollack et al. (Reference Pollack, Hubickyj and Bodenheimer1996). The left end of the lines corresponds to the core–envelope interface at R _core, while the outer radius approximately corresponds to the Hill sphere radius. The radial structure is shown at the ‘crossover’ point which is the moment when the core and envelope mass are equal. One sees that the pressure and density increase by many orders of magnitude across the envelope. At the outer edge, the density is of order of 10⁻¹¹ g cm⁻³, which is typical for the outer parts of the solar nebula, while at the core–envelope interface, the density approaches an almost fluid-like value of 0.17 g cm⁻³. Also the temperature at this position is already quite high with T _c≈18 000 K. The plot also indicates the position of the protoplanet's capture radius for 100 km planetesimals. It is about ten times larger than the core radius, leading to a core accretion rate that is about a factor 100 larger compared to the case of an envelope-free core (equation 15).

Fig. 6. Radial envelope structure inside a growing giant planet at ‘crossover’. The four panels show the pressure, temperature, density and mass as function of radius. The values at the core–envelope boundary are also given. The black dot on the lines indicates the position of the capture radius for 100 km planetesimals.

At a small core mass, the gas accretion rate found by solving the internal structure equations $\dot M_{{\rm KH}} $ is small, because the Kelvin–Helmholtz timescale is long (Ikoma et al. Reference Ikoma, Nakazawa and Emori2000). Once the crossover mass is reached which corresponds to the critical core mass in the early strictly static calculations (like Perri & Cameron Reference Perri and Cameron1974; Mizuno et al. Reference Mizuno, Nakazawa and Hayashi1978), $\dot M_{{\rm KH}} $ starts to increase rapidly, and becomes at some moment (typically when the total mass of the planet is of order of 100 M _⊕) larger than the maximal rate $\dot M_{{\rm XY},{\rm max}} $ at which the protoplanetary disc can supply gas to the planet. This means that the envelope of the planet now contracts so quickly that (at least formally) an empty shell between the planet's envelope and the background nebula develops. The planet therefore detaches from the nebula and contracts rapidly, becoming much more centrally condensed. The gas now freely falls from the Hill sphere onto the planet, where it is accreted through an accretion shockFootnote ⁴. The planetary gas accretion rate is thus

(19)$$\dot M_{XY} = \min {\rm \; (}\dot M_{{\rm KH}}, \dot M_{XY,{\rm max}} {\rm )}.$$

The disc-limited gas accretion rate $\dot M_{XY,{\rm max}} $ is controlled by a number of processes, namely, the local availability of gas, the gas flux through the disc that decreases itself in time and the effect of gap formation (Lubow et al. Reference Lubow, Seibert and Artymowicz1999; D'Angelo et al. Reference D'Angelo, Durisen, Lissauer and Seager2011). Initially, the planet accretes rapidly the gas in its vicinity, so that the gas accretion rate can be approximate by a Bondi-type rate (D'Angelo & Lubow Reference D'Angelo and Lubow2008)

(20)$$\dot M_{XY,{\rm max},{\rm R}} \approx \displaystyle{\Sigma \over H}{\rm \Omega} R_{{\rm gc}}^3, $$

where R _gc is the capture radius for gas which is approximately the smaller of (a fraction of) the Hill and the Bondi radius ${GM}{\rm /}{c}_{\rm s}^2 $ where M is the planet's mass and c _s the sound speed in the background nebula. After the local reservoir has been exhausted, the global flux of gas in the disc towards the planet starts to limit the gas accretion. This rate is now given by (e.g. Mordasini et al. Reference Mordasini, Alibert, Klahr and Henning2012c)

(21)$$\dot M_{XY,{\rm max},{\rm F}} = k_{{\rm Lub}} \left( {3{\rm \pi \nu \Sigma} + 6{\rm \pi} r\displaystyle{{\partial {\rm \nu \Sigma}} \over {\partial r}}} \right),$$

where k _Lub is a factor that describes the fraction of the gas streaming viscously through the disc towards the planet that is eventually accreted onto it (Lubow et al. Reference Lubow, Seibert and Artymowicz1999; Lubow & D'Angelo Reference Lubow and D'Angelo2006). Note that there is considerable uncertainty about the efficiency of gas accretion in the disc-limited regime (Benz et al. Reference Benz, Ida, Alibert, Lin and Mordasini2013). This directly influences the predictions for the upper end of the planetary mass function, but also for type II migration and the resulting formation of Hot Jupiters.

Other global models of planet formation describe the accretion of gas in the phase when it is limited by the envelope's contraction in a parameterized way without solving the internal planet structure (Ida & Lin Reference Ida and Lin2004a; Thommes et al. Reference Thommes, Matsumura and Rasio2008; Hellary & Nelson Reference Hellary and Nelson2012). Instead, a parametrization of the Kelvin–Helmholtz timescale τ_KH is used, and the gas accretion rate is written in an equation of the form

(22)$$\dot M_{XY} = \displaystyle{M \over {{\rm \tau} _{{\rm KH}}}}, $$

where the Kelvin–Helmholtz timescale is itself a function of the planet mass M and the opacity κ (Ikoma et al. Reference Ikoma, Nakazawa and Emori2000), typically approximated as τ_KH=kM^−pκ^−q. A potential limitation of such an approach is the dependency on the parameters k, p and q that are not well constrained (p≈1.9–3.5, Miguel & Brunini Reference Miguel and Brunini2008). This can have directly visible consequences in population syntheses, for example, in the predicted planetary mass function (Miguel & Brunini Reference Miguel and Brunini2008). When no internal structures are calculated, no direct predictions for transit or direct imaging searches can be made, especially at early times when the formation still directly influences these quantities. On the other hand, the computational cost is much reduced. It should be noted that also the solution of the 1D structure equations is (currently) not completely free from prespecified parameters. This is due to the fact that in global models, the opacity κ in equation (18) is typically not found ab initio as in specialized models, but approximated as the interstellar medium (ISM) opacity multiplied with some prespecified reduction factor (cf. subsection ‘Impact of grain opacity on the planetary radius distribution’). These specialized models (Podolak Reference Podolak2003; Movshovitz & Podolak Reference Movshovitz and Podolak2008; Movshovitz et al. Reference Movshovitz, Bodenheimer, Podolak and Lissauer2010) find the grain opacity self-consistently by solving the Smoluchowski equation in each atmospheric layer including the effects of grain growth, settling and vaporization. This is computationally time consuming, therefore Mordasini et al. (Reference Mordasini, Klahr, Alibert, Miller and Henning2014) have recently tried to parametrize these results for use in population synthesis models by deriving the reduction factor of the ISM opacity that leads to gas accretion timescales that agree with the results of the specialized models. The problem with this approach is that one global reduction factor cannot capture the dependency of the mechanisms governing the grain dynamics on planetary properties such as the core and envelope mass.

The ‘envelope’, ‘infall’ and ‘accretion rate’ sub-models of the global formation model shown in Fig. 3 are very similar to classical 1D giant planet formation codes like Bodenheimer & Pollack (Reference Bodenheimer and Pollack1986); Pollack et al. (Reference Pollack, Hubickyj and Bodenheimer1996) or Lissauer et al. (Reference Lissauer, Hubickyj, D'Angelo and Bodenheimer2009) which are the conceptual origin of this model (see Helled et al. Reference Helled, Bodenheimer and Podolak2013 for a recent review on giant planet formation). Despite this origin, most planets that form in a population synthesis are in fact low-mass planets that do not trigger gas runaway accretion. These planets can, however, still be modelled with the restriction that primordial H/He envelopes are the only type of envelope that can currently be considered.

The extension into an evolutionary model (Mordasini et al. Reference Mordasini, Alibert, Klahr and Henning2012c) can bring a global model close to classical models of (giant) planet long-term evolution (cooling and contraction) like Burrows et al. (Reference Burrows, Marley and Hubbard1997) or Baraffe et al. (Reference Baraffe, Chabrier, Barman, Allard and Hauschildt2003). The set of equations for the internal structure during the evolution after the dissipation of the protoplanetary nebula remains the same as during formation, but they are solved with different outer boundary conditions as described in the ‘atmosphere’ section below.

Figure 7 shows the long-term evolution of the internal structure of a close-in hot Jupiter planet represented by its pressure–temperature profile. The figure shows a temporal series of profiles, including both the atmosphere and the complete interior. The upper end of the lines thus corresponds to the outer radius of the planet, whereas the lower is the core–envelope boundary. The uppermost line shows the state shortly after the end of formation when the final mass has been reached, whereas the lowest line corresponds to the state after 5 Gyr. The gradual cooling of the interior is obvious. The atmospheric model is the semi-grey solution of Guillot (Reference Guillot2010) discussed below. The formation of a deep radiative zone that is characteristic for strongly irradiated giant planets is visible (Guillot & Showman Reference Guillot and Showman2002). In contrast to the interior, the temperature at the outer radius is nearly constant since it is dominated by stellar irradiation, which was assumed to be constant (the temporal evolution of the host star and its luminosity is neglected).

Fig. 7. Evolution of the interior and atmosphere over gigayears for a close-in giant planet. The uppermost line shows the pressure–temperature profile shortly after the end of formation. The lowermost line is the state after 5 Gyr. The blue part of the lines corresponds to convective regions, while red indicates radiative zones. The mass of the planet is approximately one Jovian mass while the semi-major axis is 0.04 AU.

Atmosphere of the planet

This sub-model provides the structure of the outer part of the (proto)planet and therefore also the boundary conditions necessary to solve the internal structure equations. These boundary conditions depend on the stage of formation or evolution a planet is in. Three major phases can be distinguished (e.g. Bodenheimer et al. Reference Bodenheimer, Hubickyj and Lissauer2000).

During the first attached (or nebular) phase which applies to low-mass, sub-critical planets embedded in the protoplanetary nebulae, the envelope bound to the protoplanet smoothly transitions into the unbound, background conditions in the nebula. The structure extends out to a radius R that is approximately the smaller of the Hill or Bondi radius (Lissauer et al. Reference Lissauer, Hubickyj, D'Angelo and Bodenheimer2009). No atmosphere in the classical sense existsFootnote ⁵, and the boundary conditions are

(23)$$T^4 = T_{{\rm neb}}^4 + T_{{\rm int}}^4, \quad P = P_{{\rm neb}}, $$

(24)$${\rm \tau} = {\rm max}\left( {{\rm \rho} _{{\rm neb}} {\rm \kappa} _{{\rm neb}} R,2/3} \right),\quad T_{{\rm int}}^4 = \displaystyle{{3{\rm \tau} L_{{\rm int}}} \over {8{\rm \pi \sigma} R^2}}. $$

The pressure is equal to the background midplane pressure in the nebula (provided by the disc model) P _neb. The surface temperature contains the contributions from the nebula midplane temperature T _neb and from the intrinsic luminosity of the planet L _int. The effect of the optical depth τ of the background nebula is taken into account with the approximation of Papaloizou & Nelson (Reference Papaloizou and Nelson2005).

The second detached (or transitional) phase applies to massive cores M _core≳10M _⊕, total mass M≳100M _⊕) that have triggered runaway gas accretion, that is, where the gas accretion rate due to the contraction of the envelope $\dot M_{{\rm KH}} $ would be larger than the maximal rate at which the nebula can supply gas to the planet $\dot M_{XY,{\rm max}} $, as described above. The planet has now a free outer radius that collapses rapidly from initially the Hill sphere radius to about 2–3 R (for cold accretion, i.e. low entropies, Marley et al. Reference Marley, Fortney, Hubickyj, Bodenheimer and Lissauer2007). The gas falls approximately at free fall velocity v _ff from the Hill sphere onto the planet where it shocks. The boundary conditions are then (Bodenheimer et al. Reference Bodenheimer, Hubickyj and Lissauer2000)

(25)$$\dot M_{XY} = \dot M_{XY,{\rm max}}, \quad v_{{\rm ff}}^2 = 2GM\left( {\displaystyle{1 \over R} - \displaystyle{1 \over {R_{\rm H}}}} \right),$$

(26)$$P = P_{{\rm neb}} + \displaystyle{{\dot M_{XY}} \over {4{\rm \pi} R^2}} v_{{\rm ff}} + \displaystyle{{2g} \over {3{\rm \kappa}}}, \quad {\rm \tau} = {\rm max}\left( {{\rm \rho} _{{\rm neb}} {\rm \kappa} _{{\rm neb}} R,2/3} \right),$$

(27)$$T_{{\rm int}}^4 = \displaystyle{{3{\rm \tau} L_{{\rm int}}} \over {8{\rm \pi \sigma} R^2}}, \quad T^4 = (1 - A)T_{{\rm neb}}^4 + T_{{\rm int}}^4. $$

The outer pressure therefore contains the contributions from the background nebula, the accretion shock and from the standard Eddington expression for the photospheric pressure due to the material residing above the τ=2/3 surface. The temperature contains the contributions from the nebula and the intrinsic luminosity, where A is the albedo of the planet. With these boundary conditions one describes the accreting planet as a scaled-down version of a stellar hydrostatic core undergoing spherical accretion as in the classical work of Stahler et al. (Reference Stahler, Shu and Taam1980). Hydrodynamic simulations indicate that in reality, the actual infall geometry and thus the boundary conditions are more complicated than in this 1D picture, making 3D radiation-hydrodynamic calculations necessary (e.g. Klahr & Kley Reference Klahr and Kley2006; Ayliffe & Bate Reference Ayliffe and Bate2009).

The third evolutionary (or isolated) phase starts after the protoplanetary disc has disappeared. The planet now evolves at constant mass (if we neglect further accretion or mass loss through envelope evaporation for close-in planets, see subsection ‘Atmospheric escape’). Mordasini et al. (Reference Mordasini, Alibert and Georgy2012b) model this phase with a simple grey atmosphere so that (e.g. Guillot Reference Guillot2005)

(28)$$P = \displaystyle{{2g} \over {3{\rm \kappa}}}, \quad T_{{\rm int}}^4 = \displaystyle{{L_{{\rm int}}} \over {4{\rm \pi \sigma} R^2}}, $$

(29)$$T_{{\rm equi}} = 280\,{\rm K}\left( {\displaystyle{a \over {1\,{\rm AU}}}} \right)^{ - 1/2} \left( {\displaystyle{{M_*} \over {M_ \odot}}} \right),\quad T^4 = (1 - A)T_{{\rm equi}}^4 + T_{{\rm int}}^4. $$

The equilibrium temperature due to stellar irradiation is calculated assuming that the star with mass M _* is on the main sequence where the luminosity approximately scales as ${M}_*^4 $. Other evolutionary models like Burrows et al. (Reference Burrows, Marley and Hubbard1997); Baraffe et al. (Reference Baraffe, Chabrier, Barman, Allard and Hauschildt2003, Reference Baraffe, Chabrier and Barman2008) in contrast couple the interior calculation to full non-grey atmospheres.

For a giant planet at a few AU where the irradiation flux from the host star is rather low, the grey atmosphere and the full non-grey atmospheres lead to similar cooling curves (Bodenheimer et al. Reference Bodenheimer, Hubickyj and Lissauer2000; Mordasini et al. Reference Mordasini, Alibert and Georgy2012b). For Hot Jupiter planets and in general all strongly irradiated planets, grey atmospheres lead, however, to too low temperatures deep in the atmosphere, so that the cooling timescale is underestimated (Guillot & Showman Reference Guillot and Showman2002). A better atmospheric model than the grey atmosphere is the semi-grey approximation of Guillot (Reference Guillot2010). This model provides the temperature as a function of optical depth in an atmosphere that transports both an intrinsic heat flux and receives an outer irradiation flux. The model is parametrized by two mean opacities, one in the visual and one in the IR. The mean temperature as a function of IR optical depth τ is then

(30)$$\eqalign{T^4 = & \displaystyle{{3T_{{\rm int}}^4} \over 4}\left( {\displaystyle{2 \over 3} + {\rm \tau}} \right) + \displaystyle{{3T_{{\rm equi}}^4} \over 4}\left( {\displaystyle{2 \over 3} + \displaystyle{2 \over {3{\rm \gamma}}} \left[ {1 + \left( {\displaystyle{{{\rm \gamma \tau}} \over 2} - 1} \right)e^{ - {\rm \gamma \tau}}} \right]} \right. \cr & + \left. {\displaystyle{{2{\rm \gamma}} \over 3}\left( {1 - \displaystyle{{{\rm \tau} ^2} \over 2}} \right)E_2 ({\rm \gamma \tau} )} \right),} $$

where γ denotes the ratio of the mean opacity in the visual κ_v to the mean opacity in the thermal infrared κ_th, whereas E ₂ is an exponential integral. This expression provides a fair approximation of detailed irradiated atmosphere models (e.g. Showman et al. Reference Showman, Fortney and Lian2009) and therefore also gives more realistic cooling curves. This is important for the comparison of the radii of synthetic and actual transiting exoplanets. A better atmospheric structure also makes it possible to calculate the transit radius that is bigger than the τ≈1 radius to the grazing observational geometry more accurately (e.g. Hansen Reference Hansen2008).

Figure 8 shows a typical example of the atmospheric structure of a Hot Jupiter found with the semi-grey model. The nominal pressure–temperature profile is calculated with the Rosseland mean opacity in the IR for a solar-composition gas from Freedman et al. (Reference Freedman, Marley and Lodders2008) and the parameter γ is set to 0.4. For the black line, the same ratio is used, but the opacity in the thermal domain is fixed to 0.01 cm² g⁻¹. A strong reduction of the opacity in the visual leads to a deeper penetration of the irradiation into the planet (blue line), while an enhanced optical opacity leads to a hotter outer atmosphere with a temperature inversion, while the deep atmosphere is cooler (green line).

Fig. 8. Atmospheric pressure–temperature profiles for a Hot Jupiter found with the semi-gray approximation of Guillot (Reference Guillot2010). The dotted line is the nominal model, while other lines are calculated with different opacities in the visual and thermal domain.

Besides the semi-grey atmospheres for strongly irradiated planets that are studied with transit observations, the grey atmosphere should also be replaced with a full coupling of the internal structure calculations with the non-grey atmospheres of, e.g. Allard et al. (Reference Allard, Homeier and Freytag2011) for non (or weakly)-irradiated giant planets, following the procedure explained in Chabrier & Baraffe (Reference Chabrier and Baraffe1997). This will not only provide more accurate cooling curves for giant planets observed with direct imaging, but also their magnitudes in the different observational bands instead of the total luminosity only. This will make it possible to directly compare the predictions from population syntheses with the discoveries of new direct imaging instruments such as Sphere or Gpi.

Infall of planetesimals into the protoplanet's envelope

This sub-model calculates the interaction of planetesimals and the gaseous envelope of the protoplanet during the formation phase (Mordasini et al. Reference Mordasini, Alibert and Benz2006). A similar sub-model was described by Podolak et al. (Reference Podolak, Pollack and Reynolds1988) for the classical giant planet models of Pollack et al. (Reference Pollack, Hubickyj and Bodenheimer1996). This sub-model links the accretion of solid (5.4) and the envelope structure (5.5). Two quantities are the main output of the sub-model: the protoplanet's capture radius R _capture for planetesimals that enters the solid accretion rate (equation 15), and the radial energy and mass deposition profiles that enter into the calculation of the envelope structure in particular via the equation for the luminosity. It also yields how much high-Z material is deposited in the envelope to enrich it, and how much of it directly reaches the core. This is important in the context of the (atmospheric) composition of planets (Fortney et al. Reference Fortney, Mordasini and Nettelmann2013). The sub-model calculating the interaction (see Mordasini et al. Reference Mordasini, Alibert and Benz2006, for a short overview) includes gravity and gas drag, thermal ablation as for shooting stars and aerodynamical disruption inspired by the destruction of comet Shoemaker-Levy 9 in Jupiter's atmosphere (e.g. Zahnle & Mac Low Reference Zahnle and Mac Low1994; Crawford et al. Reference Crawford, Boslough, Trucano and Robinson1995) or the Tunguska event (Chyba et al. Reference Chyba, Thomas and Zahnle1993).

To model the interaction, a set of coupled ordinary differential equations are integrated numerically, giving the position in three dimensions r, velocity ${\rm \dot r}$, mass M _pl and radius R _pl of the impacting planetesimal as a function of time t. Infalling planetesimals are accelerated by the gravity of the protoplanet and slowed down by gas drag. The equation of motion for the planetesimal is therefore given by (planetocentric reference frame, 2 body approximation)

(31)$$M_{{\rm pl}} {\bf \ddot r} = - \displaystyle{{Gm(r)M_{{\rm pl}}} \over {r^2}} \cdot \displaystyle{{\rm r} \over r} - \displaystyle{1 \over 2}C_D {\rm \rho} \,\dot r^2 \displaystyle{{{\rm \dot r}} \over {\dot r}}S,$$

where C _D is the drag coefficient that can be written as a function of the Reynolds and the Mach number, while ρ is the density of the gas through which the planetesimal is plowing. It is obtained from the calculations of the internal structure of the planet as is the mass m(r) of the protoplanet inside of the position of the planetesimal. S is the cross-section of the planetesimal, equal to ${\rm \pi} {R}_{{\rm pl}}^2 $ for a spherical planetesimal. Note that it is only initially assumed that the planetesimal is spherical, later on it can get distorted due to aerodynamic forces.

As the planetesimal flies through the envelope, pressure and temperature are increasing. Eventually, two effects can lead to its destruction: thermal ablation and mechanical mass loss. These effects determine how deep the planetesimal is able to penetrate, thus determining where in the planet's envelope the energy and the debris of the planetesimals are deposited. The mass loss rate due to thermal ablation can be written in its simplest form as (Opik Reference Opik1958)

(32)$$\displaystyle{{{\rm d}M_{{\rm pl}}} \over {{\rm dt}}} = - \displaystyle{1 \over 2}C_H {\rm \rho} \,\dot r^3 S\displaystyle{1 \over {Q_{abl}}}. $$

Q _abl is the amount of energy needed to heat body material from its initial temperature to the point where melting or vaporization occurs plus the specific heat needed for this phase change. The heat transfer coefficient C _H is an a priori unknown function that can vary over many orders of magnitude (Svetsov et al. Reference Svetsov, Nemtchinov and Teterev1995). It depends on the velocity, envelope conditions, flow regime, shape of the body, etc. and gives the fraction of the incoming kinetic energy flux of the gas that is available for ablation. Note that in order to compute this fraction, the hydro- and thermodynamic state of the flow surrounding the impacting planetesimals needs to be known. The most important flow regime for massive impactors is a hypersonic highly turbulent continuum flow. This means that a strong shock wave forms. The radiation field generated by the shock is the dominant heating source leading to thermal ablation, as temperatures in excess of 30 000 K are reached (e.g. Zahnle Reference Zahnle1992). Therefore, in these conditions, the bow shock temperature must be computed by solving the shock jump conditions for a non-ideal gas (cf. Chevalier & Sarazin Reference Chevalier and Sarazin1994).

Besides thermal ablation, mechanical effects can result in a rapid destruction of the planetesimal by fragmenting it into a large number of small pieces that eventually get thermally ablated. Mechanical destruction is of prime importance for massive (km-sized) impactors as was understood in the hydrodynamical simulations for SL-9 (e.g. Svetsov et al. Reference Svetsov, Nemtchinov and Teterev1995). The main effect to consider is the large pressure difference between the front of the planetesimal (stagnation point) and the back where the pressure almost vanishes. This pressure difference (if large enough) leads to a lateral spreading of the body (the so-called ‘pancake’ model, introduced by Zahnle Reference Zahnle1992) and eventually its disruption by fragmentation. The rate of lateral spreading is given as

(33)$$\displaystyle{{{\rm d}^2 R_{{\rm pl}}} \over {{\rm dt}^2}} = C_{\rm S} \displaystyle{{\rm \rho} \over {{\rm \rho} _{\rm b}}} \displaystyle{{\dot r^2} \over {R_{{\rm pl}}}}, $$

where C _S is a coefficient of order unity (Chyba et al. Reference Chyba, Thomas and Zahnle1993), whereas ρ_b is the mean density of the fluidized impactor.

The impactor acts approximately as a fluid for ram pressure exceeding the internal tensile strength by a large factor. The disruption into many fragments is then due to Rayleigh–Taylor (RT) instabilities that grow due to the deceleration of the front side of the body (a denser fluid) by the shocked gas (a less dense fluid). Such instabilities are seen to develop in hydrodynamical simulations (e.g. Mac Low & Zahnle Reference Mac Low and Zahnle1994; Korycansky et al. Reference Korycansky, Zahnle and Low2000). Mordasini et al. (Reference Mordasini, Alibert and Benz2006) describe this process as a fragmentation cascade due to growing RT fingers. In the non-linear stage, the height until which fingers have grown into the body h _inst can be estimated as (Sharp Reference Sharp1984; Youngs Reference Youngs1989)

(34)$$\displaystyle{{{\rm d}^2 h_{{\rm inst}}} \over {{\rm dt}^2}} = 2{\rm \alpha} \left( {\displaystyle{{{\rm \rho} _{\rm b} - {\rm \rho} _{{\rm stagn}}} \over {{\rm \rho} _{\rm b} + {\rm \rho} _{{\rm stagn}}}}} \right)\ddot r,$$

where the expression in brackets is the Atwood number, α a parameter ≈0.1 (Sharp Reference Sharp1984; Youngs Reference Youngs1989), while ρ_stagn is the gas density at the stagnation point. The combined actions of flattening and growth of RT fingers leads to a fast destruction of the impactor in a terminal explosion, similar as in the Tunguska event (Chyba et al. Reference Chyba, Thomas and Zahnle1993).

Figure 9 from Mordasini et al. (Reference Mordasini, Alibert and Benz2006) illustrates the output of the ‘infall’ sub-model. It shows the fate of stony planetesimals of various initial radii in protoplanetary envelopes with masses between 0.001 and 100 M _⊕. The basic result is that the larger the impactor and the lower the envelope mass, the more likely it is that the planetesimal can penetrate to the core, as one expects. The plot shows that the detailed structure is, however, more complex. There are intermediate sized planetesimals (radii ~100–1000 m) that penetrate through surprisingly massive envelopes. They are too massive for efficient purely thermal ablation, but too small to undergo intense fragmentation. Big bodies (≳100 km) on the other hand are protected by their self-gravity from intense fragmentation. It is clear that these results depend on the material properties of the impactors, which are typically not well constrained. Icy planetesimal are more prone to destruction in the envelope (Podolak et al. Reference Podolak, Pollack and Reynolds1988), for large impactors mainly due to their lower tensile strength (by about a factor 100, Chyba et al. Reference Chyba, Thomas and Zahnle1993). If planetesimals have properties similar to the parent body of the Shoemaker-Levy 9 impactors which was a nearly strength-less rubble pile hold together only by self-gravity (Asphaug & Benz Reference Asphaug and Benz1994), the ability of planetesimals to reach the solid core would be even more reduced.

Fig. 9. Fate of planetesimals in the envelope of a growing giant planet. The plot shows whether they can penetrate to the solid core or if they get completely destroyed in the envelope as a function of the envelope mass of the protoplanet and the initial radius of the impacting planetesimal (thick solid line). Planetesimals inside the roughly squared region limited by the thin line in the upper right corner undergo aerodynamic fragmentation (figure from Mordasini et al. Reference Mordasini, Alibert and Benz2006).

The population synthesis framework shown in Fig. 3 is currently the only global planet formation and evolution model that explicitly addresses the interaction of planetesimals and the protoplanetary atmosphere. The coupling with the rest of the model is, however, at the moment not yet self-consistent. Despite the fact that the mass deposition profiles are calculated, it is assumed in the ‘envelope’ sub-model that all accreted solids immediately sink to the core (‘sinking approximation’ Pollack et al. Reference Pollack, Hubickyj and Bodenheimer1996). The output of the ‘infall’ sub-model is thus not used to calculate the (radially changing) composition of the gas, and linked to that, its opacity. Both are known to be important for the formation process (e.g. Movshovitz et al. Reference Movshovitz, Bodenheimer, Podolak and Lissauer2010; Hori & Ikoma Reference Hori and Ikoma2011), and compositional gradients even have the potential to lead to a semi-convective interior instead of the usually assumed fully convective state in giant planets, with important consequences for the cooling and inferred bulk composition (Leconte & Chabrier Reference Leconte and Chabrier2012). Including these effects in future work along the lines demonstrated by, e.g. Hori & Ikoma (Reference Hori and Ikoma2011) and Iaroslavitz & Podolak (Reference Iaroslavitz and Podolak2007) is therefore important, especially for a more accurate theoretical characterizations of planets in terms of their bulk and atmospheric composition, a quantity that can potentially be measured by spectroscopy (e.g. Bean et al. Reference Bean, Désert and Kabath2011; Bonnefoy et al. Reference Bonnefoy, Boccaletti and Lagrange2013; Fortney et al. Reference Fortney, Mordasini and Nettelmann2013; Konopacky et al. Reference Konopacky, Barman, Macintosh and Marois2013).

Internal structure of the solid core

This sub-model calculates the radius of the solid core as a function of its mass, bulk composition (iron, silicate and ice mass fraction) and external pressure due to the surrounding gaseous envelope which is important for giant planets. It is necessary for the correct prediction of the total radius of solid planets without a sizable H/He envelope. But also for giant planets, it is necessary to calculate the core radius correctly to obtain accurate total radii during the evolutionary phase (e.g. Mordasini et al. Reference Mordasini, Alibert, Klahr and Henning2012c).

The sub-model for the solid core used in the global model of Fig. 3 was originally developed in Figueira et al. (Reference Figueira, Pont and Mordasini2009) and assumes a differentiated planet consisting of concentric shells of iron, silicates, and if the planet accreted outside of the iceline, ices. As described in Fortney et al. (Reference Fortney, Marley and Barnes2007), the radius is found by solving the 1D internal structure equations that are in principle the same as for the gaseous envelope. The situation is, however, considerably simplified by assuming that the density of the solid material ρ is (in contrast to gas) approximately independent of temperature, so that the density is function of pressure P only. The system of equations to solve is then just

(35)$$\displaystyle{{{\rm d}m} \over {{\rm d}r}} = 4{\rm \pi} r^2 {\rm \rho,} \quad \displaystyle{{{\rm d}P} \over {{\rm d}r}} = - \displaystyle{{Gm} \over {r^2}} {\rm \rho,} $$

which are the equations of mass conservation and hydrostatic equilibrium (m denotes the core mass inside a radius r, and G is the gravitational constant).

While Figueira et al. (Reference Figueira, Pont and Mordasini2009) used the more accurate tabulated equations of state (EoS) of Fortney et al. (Reference Fortney, Marley and Barnes2007), the global model of Mordasini et al. (Reference Mordasini, Alibert and Georgy2012b) employs the simpler, but more widely applicable EoS of Seager et al. (Reference Seager, Kuchner, Hier-Majumder and Militzer2007) that is a modified polytropic EoS with the material parameters ρ₀, c and n

(36)$${\rm \rho} (P) = {\rm \rho} _0 + cP^n. $$

For silicates, the parameters appropriate for perovskite (MgSiO₃) are used. An advantage of this EoS is that it approaches at sufficiently high pressures (giant planet cores) the correct asymptotic limit which is the EoS of a completely degenerate, non-relativistic electron gas (Zapolsky & Salpeter Reference Zapolsky and Salpeter1969). The mean density of a core of a massive giant planet can in principle reach very high values exceeding 100 g cm⁻³, with potentially observable consequences (Charpinet et al. Reference Charpinet, Fontaine and Brassard2011).

Figure 10 illustrates the output of this sub-model by showing the radius as a function of mass for low-mass solid planets (no external pressure). Three different bulk compositions are shown: Earth-like with a 2 : 1 silicate-to-iron ratio, pure water ice and a mixture of both types of materials. The plot shows that the radius follows to good approximation a power-law, as noted by Valencia et al. (Reference Valencia, O'Connell and Sasselov2006).

Fig. 10. Radius of low-mass solid planets as a function of core mass and composition. The lines show planets having an Earth-like composition (red solid line), pure water ice planets (blue dashed) and a mixed composition with 50% water ice, 33% silicates and 17% iron (green dotted line).

A second output of this sub-model is the time-dependent radiogenic luminosity of the solid core that is due to the decay of long-lived radionuclides in the mantle. As described in Mordasini et al. (Reference Mordasini, Alibert and Georgy2012b) the radiogenic luminosity can be modelled according to the law of radioactive decay assuming that the mantle has a chondritic composition (Urey Reference Urey1955; Lowrie Reference Lowrie2007). The radiogenic core luminosity is very small compared to the luminosity due to the cooling of the gaseous envelope of a giant planet. For core dominated planets it can, however, become the dominant internal heat source during the evolution over gigayears, and therefore affect the contraction of the planet (e.g. Nettelmann et al. Reference Nettelmann, Fortney, Kramm and Redmer2011). Much can be improved in the basic description of the interior of solid planets presented here, like the inclusion of a more accurate EoS, the calculation of the radial temperature structure or the thermal cooling of the interior (e.g. Valencia et al. Reference Valencia, O'Connell and Sasselov2006; Lopez & Fortney Reference Lopez and Fortney2013b). At early times, the effect of cooling of the core could be particularly important for rocky planets with cores retaining residual heat from the accretion process, especially if the core contains heat from a violent, short timescale build-upFootnote ⁶. Also the process of core formation (differentiation of rocky planets into a metallic core and silicate mantle) provides an important heating source, the timescale of which is however not yet entirely clear (Rubie et al. Reference Rubie, Nimmo, Melosh, Schubert and Stevenson2007). Further improvements could be to include the dissolution of the solid core in giant planets (e.g. Guillot et al. Reference Guillot, Stevenson, Hubbard and Saumon2004; Wilson & Militzer Reference Wilson and Militzer2012) and the outgassing of secondary atmospheres given the core's composition acquired during formation (e.g. Elkins-Tanton & Seager Reference Elkins-Tanton and Seager2008). The latter point is particularly relevant as the James Webb Space Telescope should make it possible to characterize the atmospheres of a few low-mass planets (e.g. Belu et al. Reference Belu, Selsis and Morales2011).

Atmospheric escape

The Kepler satellite has discovered a very large population of close-in low-mass planets, in agreement with result from high-precision radial velocity searches. Owing to their proximity to the host star and low surface gravity, these low-mass planets are sensitive to envelope mass loss due to atmospheric escape of the primordial H/He envelope (e.g. Lammer et al. Reference Lammer, Selsis and Ribas2003; Baraffe et al. Reference Baraffe, Selsis and Chabrier2004; Erkaev et al. Reference Erkaev, Kulikov and Lammer2007; Murray-Clay et al. Reference Murray-Clay, Chiang and Murray2009; Owen & Jackson Reference Owen and Jackson2012). This means that for such planets, atmospheric escape is important on a population level, shaping the statistical properties like the radius distribution (Lopez et al. Reference Lopez, Fortney and Miller2012; Lopez & Fortney Reference Lopez and Fortney2013a; Owen & Wu Reference Owen and Wu2013). This is due to the fact that the presence of even a very tenuous H/He envelope (low-mass fraction) has a large impact on the total radius (Adams et al. Reference Adams, Seager and Elkins-Tanton2008; Rogers et al. Reference Rogers, Bodenheimer, Lissauer and Seager2011; Mordasini et al. Reference Mordasini, Alibert, Klahr and Henning2012c).

Therefore, if one wants to connect predictions of a formation model (which yields the H/He mass after formation) with observations by the Kepler satellite, it is necessary to include envelope evaporation during the evolutionary phase. In principle, to find the escape rate, it is necessary to solve the radiation-hydrodynamic equations describing the flow of the upper atmosphere under the effects of heating by UV and X-ray irradiation on a (potentially multidimensional) grid (Murray-Clay et al. Reference Murray-Clay, Chiang and Murray2009; Owen & Jackson Reference Owen and Jackson2012).

A simplified description of the mass loss rate of close-in planets suitable for global models consists of the following elements: first a description for the incoming stellar EUV (and X-ray) flux as a function of time t and orbital distance of the planet a is needed. The EUV flux can be approximated as (Ribas et al. Reference Ribas, Guinan, Güdel and Audard2005)

(37)$$F_{{\rm UV}} \approx F_{{\rm UV},0} \left( {\displaystyle{t \over {1\; \,{\rm Gyr}}}} \right)^{ - 1} \left( {\displaystyle{a \over {1\,\; {\rm AU}}}} \right)^{ - 2}, $$

where F _UV,0 depends on the host star type (Lecavelier Des Etangs Reference Lecavelier Des Etangs2007), while the power-law index for the temporal decrease depends somewhat on the wavelength interval that is considered. This equation is valid for stars older than ~100 Myr, while at younger ages, the flux saturates at a maximum value.

Envelope evaporation can either be driven by X-rays or EUV, depending on the relative position of the ionization front and the X-ray sonic point, which allows us to identify which mechanism is dominant (Owen & Jackson Reference Owen and Jackson2012).

In the EUV regime, two sub-regimes exist (Murray-Clay et al. Reference Murray-Clay, Chiang and Murray2009). At lower EUV fluxes, most of the incoming energy flux goes into pdV work that lifts gas out of the planet's potential well, while radiative losses and internal energy changes are small. Therefore one can write the evaporation rate with an equation of the form (so-called energy-limited approximation, Watson et al. Reference Watson, Donahue and Walker1981)

(38)$$\displaystyle{{{\rm d}M_{{\rm UV},{\rm e} - {\rm lim}}} \over {{\rm dt}}} = \displaystyle{{{\rm \epsilon} _{{\rm UV}} {\rm \pi} F_{{\rm UV}} R_{{\rm UV}}^3} \over {GMK_{{\rm tide}}}}. $$

In this equation, M is the planet's mass, ϵ_UV is an efficiency factor (that hides the complex physics), R _UV the planetary radius where EUV radiation is typically absorbed (estimated as in Murray-Clay et al. Reference Murray-Clay, Chiang and Murray2009), and K _tide is a factor to take into account that mass only needs to reach the Hill sphere to escape (Erkaev et al. Reference Erkaev, Kulikov and Lammer2007).

At high EUV fluxes, most of the UV heating is lost via radiative cooling, so that an equilibrium between photoionization and radiative recombination is established. The evaporation rate in this radiation-recombination limited regime can be approximated as (Murray-Clay et al. Reference Murray-Clay, Chiang and Murray2009)

(39)$$\displaystyle{{{\rm d}M_{{\rm UV},{\rm rr} - {\rm lim}}} \over {{\rm dt}}} = 4{\rm \pi \rho} _{\rm s} c_{\rm s} r_{\rm s}^2, $$

where ρ_s and c_s are the density and speed of sound at the sonic point at a radius r _s. These quantities can be estimated as described in Murray-Clay et al. (Reference Murray-Clay, Chiang and Murray2009).

The mass loss rate in the X-ray-driven regime can (roughly) be estimated with an analogous equation as in the energy-limited UV regime. These equations are then solved together with the equations for the internal structure and the atmosphere, yielding the evolution of the planet's envelope mass and radius as a function of time. When coupled with population synthesis calculations, this allows us to study the global effects of envelope evaporation, for example, on the planetary radius distribution (Jin et al. Reference Jin, Mordasini and Parmentier2014).

It should, however, be noted that the efficiency factors (both for EUV and X-ray-driven evaporation) are in reality not constants, but depend on the planet mass, radius and magnitude of the heating flux. This means that the mass loss rates found with constant factors should be seen as rough estimates, in particular for the mass loss history of an individual object (Owen & Wu Reference Owen and Wu2013). On the other hand, calculations of the evolution of an entire population of planets under the influence of atmospheric escape indicate that the statistical consequences do not vary drastically when the efficiency factors are varied within reasonable limits (Jin et al. Reference Jin, Mordasini and Parmentier2014).

Orbital migration

The discovery of a Jovian planet at an orbital distance of only 0.05 AU from its star by Mayor & Queloz (Reference Mayor and Queloz1995) was a surprise because theoretical planet formation models had rather predicted (e.g. Boss Reference Boss1995) that giant planets should be found several AU away, as it is the case in the Solar System. The mechanism that was underestimated was orbital migration, that is, the radial displacement of planets. Several mechanisms can cause orbital migration including classical disc migration (Goldreich & Tremaine Reference Goldreich and Tremaine1980), migration due to the planetesimal disc (e.g. Levison et al. Reference Levison, Thommes and Duncan2010; Ormel et al. Reference Ormel, Ida and Tanaka2012) or Kozai migration with tidal circularization (Kozai Reference Kozai1962; Fabrycky & Tremaine Reference Fabrycky and Tremaine2007). Another possible mechanism producing close-in planets is planet–planet scattering followed by tidal circularization (e.g. Ivanov & Papaloizou Reference Ivanov and Papaloizou2004; Papaloizou & Terquem Reference Papaloizou and Terquem2006; Beaugé & Nesvorný Reference Beaugé and Nesvorný2012). Support for such a dynamical origin of Hot Jupiters arises from the distribution of the orbital distances of these planets relative to the Roche limit (Ford & Rasio Reference Ford and Rasio2006; Valsecchi & Rasio Reference Valsecchi and Rasio2014, but see Rice et al. Reference Rice, Veljanoski and Collier Cameron2012 for a contrasting view) or the observation of highly inclined planetary orbits relative to the stellar equator via measurements of the Rossiter–McLaughlin effect (see, e.g. Winn et al. Reference Winn, Johnson and Peek2007; Triaud et al. Reference Triaud, Collier Cameron and Queloz2010). Such misaligned orbits are not expected from classical disc migration (but see also Batygin Reference Batygin2012). The relative importance of the various mechanisms that lead to Hot Jupiters is currently debated. A recent comparison of the distribution of the spin–orbit angles indicates (Crida & Batygin Reference Crida and Batygin2014) that both disc migration and dynamical mechanisms contribute.

Here we concentrate on classical disc migration which is the consequence of the gravitational interaction of the gaseous protoplanetary disc and embedded planets. This mechanism was the first to be included in most global planet formation models. The dynamical effects in multi-body systems which can potentially also lead to Hot Jupiters were in contrast only addressed recently in population syntheses (Alibert et al. Reference Alibert, Carron and Fortier2013; Ida et al. Reference Ida, Lin and Nagasawa2013). The main result of the large body of studies addressing disc migration (see Kley & Nelson Reference Kley and Nelson2012, for a recent review) is that angular momentum is being transferred between disc and planet via torques that lead in most cases to a loss of angular momentum for the planet (inward migration). The angular momentum J of a planet of mass M in orbit around a star of mass M _* at a semi-major axis a, and the migration rate da/dt under the action of a total torque Γ_tot are given as

(40)$$J = M\sqrt {GM_* a}, \quad \displaystyle{{{\rm d}a} \over {{\rm dt}}} = 2a\displaystyle{{\Gamma _{{\rm tot}}} \over J}.$$

Two types of migration are distinguished depending upon the mass of the planet, respectively its impact on the disc structure. Low-mass planets (M≲10–100 M _⊕) lead to angular momentum fluxes that are much smaller than the background viscous angular momentum transport and therefore do not affect the surface density of the disc in a significant way. They migrate in type I migration where the total torque is found by summing up the contributions from the inner and outer Lindblad torques plus the corotation torque (e.g. Ward Reference Ward1986; Tanaka et al. Reference Tanaka, Takeuchi and Ward2002). The total torque Γ_tot can be expressed in a form like (Paardekooper et al. Reference Paardekooper, Baruteau, Crida and Kley2010)

(41)$${\rm \Gamma} _{{\rm tot}} = \displaystyle{1 \over {\rm \gamma}} (C_0 + C_1 p_{\rm \Sigma} + C_2 p_{\rm T} ){\rm \Gamma} _0, \quad {\rm \Gamma} _0 = \left( {\displaystyle{q \over h}} \right)^2 {\rm \Sigma} a^4 {\rm \Omega} ^2 $$

where γ is the adiabatic index of the gas, q=M/M _*, h the disc aspect ratio, Σ the gas surface density at the planet's position, and Ω the Keplerian frequency. The constants C _i depend on the local thermodynamical regime in the disc. The quantities p _Σ and p _T are the local radial power-law slopes of the gas surface density and temperature profile in the protoplanetary disc. The work of Tanaka et al. (Reference Tanaka, Takeuchi and Ward2002) assumed a locally isothermal disc which results in fast inward migration for typical disc conditions. These equations were used in several older population synthesis calculations (e.g. Ida & Lin Reference Ida and Lin2008a; Thommes et al. Reference Thommes, Matsumura and Rasio2008; Mordasini et al. Reference Mordasini, Alibert, Benz and Naef2009b; Alibert et al. Reference Alibert, Mordasini and Benz2011) where it was found – not surprisingly – that the migration rates needed to be reduced by large factors to bring the synthetic results in better agreement with observations. In the meantime, it was understood (e.g. Baruteau & Masset Reference Baruteau and Masset2008; Casoli & Masset Reference Casoli and Masset2009; Kley et al. Reference Kley, Bitsch and Klahr2009; Paardekooper et al. Reference Paardekooper, Baruteau, Crida and Kley2010) that in more realistic non-isothermal discs, there are several sub-types of type I migration, some of which lead to outward migration. The different sub-types can be identified by the comparison of four characteristic timescales, namely, the U-turn, the viscous, the libration and the cooling timescale (Dittkrist et al. Reference Dittkrist, Mordasini, Klahr, Alibert and Henning2014). These more realistic type I descriptions were used in several recent population synthesis simulations (Hellary & Nelson Reference Hellary and Nelson2012; Dittkrist et al. Reference Dittkrist, Mordasini, Klahr, Alibert and Henning2014).

In each type I sub-regime, the migration rate and direction depend besides the planet's mass on the local radial slopes of the surface density and temperature that are given by the disc model (subsection ‘Radial structure of the protoplanetary (gas) disc’). These slopes change due to opacity transitions which has the consequence in the adiabatic sub-regime that there are special places in the disc with zero net torque and a negative dΓ_tot/da. This means that inwards of these positions, migration is directed outwards, and outwards of them, it is directed inwards, so that planets in this convergence zone migrate from both sides towards the migration traps (Lyra et al. Reference Lyra, Paardekooper and Mac Low2010; Sándor et al. Reference Sándor, Lyra and Dullemond2011; Kretke & Lin Reference Kretke and Lin2012; Dittkrist et al. Reference Dittkrist, Mordasini, Klahr, Alibert and Henning2014).

Figure 11 shows the semi-major axis of six growing planets undergoing non-isothermal migration. All planets migrate in the same protoplanetary disc, but do not interact mutually. The planetary embryos start at different locations ranging from 3 to 8 AU. One sees that at the beginning, the planets starting inside of 6 AU migrate quickly outwards, while those starting further out migrate inwards, so that all planets reach the convergence point where the total torque vanishes. If the disc itself would not evolve, migration would stop at this point. Owing to disc evolution the point of vanishing torque moves itself inwards. This inward migration is much slower than isothermal type I migration as it happens on a viscous timescale (Paardekooper et al. Reference Paardekooper, Baruteau, Crida and Kley2010). During this time, the planets grow by accreting planetesimals and gas.

Fig. 11. Semi-major axis as a function of time for six migrating and growing (but mutually non-interacting) protoplanets in the same viscously evolving disc. The plot shows the rapid migration towards the convergence point and the slow inward migration while the planets are captured into it.

After a few 10⁵ years, the planets leave the convergence point because the corotation torque saturates, so that they are back at faster inward migration. Shortly afterwards, the planets are sufficiently massive that they transition into type II migration (see below) which is again slower. The planets eventually stop as the local disc mass becomes smaller than the planet mass, so that the inertia of the planet prevents further rapid migration. Simulations that include the gravitational interaction and collisions between the protoplanets find that the convergence zone can be the place of rapid growth as it concentrates a lot of matter in one region. It has therefore the potential to be a preferred place for rapid core formation for (giant) planets (e.g. Horn et al. Reference Horn, Lyra, Mac Low and Sándor2012).

The second main type of disc driven migration is type II migration of sufficiently massive planets. The gravitational interaction of such massive planets with the protoplanetary discs repels the gas in an annulus around the planet, so that a gap forms (Lin & Papaloizou Reference Lin and Papaloizou1986). According to the so-called viscous criterion for gap formation, the torque due to the presence of the planet (that repels the gas) must be larger than the background torque in the disc due to turbulent viscosity (that tries to fill up the gap). The additional thermal criterion demands that the disc vertical scale height should be smaller than the planetary Hill sphere radius, so that no sharp, unstable density gradients arise. Crida et al. (Reference Crida, Morbidelli and Masset2006) derived a combined criterion that was used in Mordasini et al. (Reference Mordasini, Alibert and Georgy2012b) while earlier population syntheses partially only used the thermal criterion.

Type II migration itself comes in two sub-regimes (Armitage Reference Armitage2007). Disc-dominated type II migration occurs if the mass of the planet M is much smaller than the local disc mass (~Σa ²). The planet then acts as a relay that communicates the viscous torques in the disc across the gap by tidal torques. The planet's migration is then locked to the evolution of the disc itself. The migration rate is thus equal to the radial velocity of the gas v _r,gas

(42)$$\displaystyle{{{\rm d}a} \over {{\rm dt}}} = v_{{\rm r},{\rm gas}} = - \displaystyle{3 \over {\Sigma \sqrt r}} \displaystyle{\partial \over {\partial r}}\left( {{\rm \Sigma \nu} \sqrt r} \right)$$

(43)$$ = - \displaystyle{{3{\rm \nu}} \over {2r}} - \displaystyle{3 \over \Sigma} \displaystyle{\partial \over {\partial r}}({\rm \Sigma \nu} ),$$

but note that hydrodynamic simulations find a more complex behaviour depending on the planet mass and gas surface density (Edgar Reference Edgar2007). This means that inside of the radius of maximum viscous couple (or velocity reversal, Lynden-Bell & Pringle Reference Lynden-Bell and Pringle1974), the planet migrates inwards (which is the normal case) while in the outer parts where the disc is spreading, it moves outwards.

Massive planets in the inner disc (or at late times into the disc evolution) are more massive than the local disc mass, so that planet-dominated migration occurs. The migration rate is then given as

(44)$$\displaystyle{{{\rm d}a} \over {{\rm dt}}} = \left( {\displaystyle{{2\Sigma a^2} \over M}} \right)^{k_{\rm p}} v_{{\rm r,gas}}, $$

where k _p is equal to 1 and 1/2 in the fully and partially suppressed case (Alexander & Armitage Reference Alexander and Armitage2009). The resulting slow-down of planets is important for the final semi-major axis distribution of giant planets (Mordasini et al. Reference Mordasini, Alibert, Benz and Naef2009b).

The rough approximation to estimate the planetary type II migration rate based simply on the radial velocity of the gas can be replaced by the direct summation of the torques according to the original impulse approximation (Lin & Papaloizou Reference Lin and Papaloizou1986; Alexander & Armitage Reference Alexander and Armitage2009). On a longer timescale, it might be desirable for global models to transition to hydrodynamic 2D disc models, because they allow us to capture phenomena that are difficult to model in 1D. An example is the outward migration of two giant planets locked into mean motion resonances (Masset & Snellgrove Reference Masset and Snellgrove2001).

The global formation model of Hellary & Nelson (Reference Hellary and Nelson2012) uses a similar description of non-isothermal type I migration and type II migration as presented here. Also Thommes et al. (Reference Thommes, Matsumura and Rasio2008) use the isothermal type I migration rates of Tanaka et al. (Reference Tanaka, Takeuchi and Ward2002), but obtain the type II migration rate based on the more accurate impulse approximation (Lin & Papaloizou Reference Lin and Papaloizou1986). Their disc model is a 1D viscously evolving model, similar as presented in subsection ‘Radial structure of the protoplanetary (gas) disc’. Ida & Lin (Reference Ida and Lin2008a) study the global effects of type I migration by population synthesis calculations. They use the migration description of Tanaka et al. (Reference Tanaka, Takeuchi and Ward2002) but reduce the rate by an arbitrary factor ≤1, an approach very similar to the one of Mordasini et al. (Reference Mordasini, Alibert, Benz and Naef2009b). This is an example how specialized models can be compared with observations thanks to population synthesis calculations.

Interaction between several (proto-)planets

From the oligarchic growth regime (e.g. Ida & Makino Reference Ida and Makino1993) it is expected that throughout the disc massive bodies (‘oligarchs’) form with radial spacings equal to a few mutual Hill spheres. For such a concurrent formation of many protoplanets, multiple effects arise due to the interaction of the planets: regarding accretion, the planets compete for the accretion of gas and solids, excite the random velocities of the planetesimals with consequences for the solid accretion rate of neighbouring protoplanets, and increase (for insufficient damping by the disc of gas or planetesimals) the eccentricities of the massive bodies which leads to collisions among the protoplanets and/or their ejection. Also the orbital migration is modified, since migrating planets can get caught into mean-motion resonances, which can, for example cause outward migration of giant planets (Masset & Snellgrove Reference Masset and Snellgrove2001). The modification of the surface density of the gas disc due to an already formed giant planets affects the migration of subsequently forming cores (e.g. Masset et al. Reference Masset, Morbidelli, Crida and Ferreira2006). The population synthesis calculations of Mordasini et al. (Reference Mordasini, Alibert, Benz and Naef2009b and Reference Mordasini, Alibert and Georgy2012b) in contrast assumed that only one planet per disc forms. This one-embryo-per-disc approximation is probably the most severe limitation of the first generation of global models, especially for low-mass planets since the observations show (e.g. Mayor et al. Reference Mayor, Marmier and Lovis2011) that such planets are mostly found in multiple planet systems. Therefore it is likely that the formation of each planet was influenced by the presence of other bodies. For single giant planets, the impact is probably less severe.

In the context of global planet formation models, this limitation was recently addressed in Ida & Lin (Reference Ida and Lin2010), Ida et al. (Reference Ida, Lin and Nagasawa2013) and Alibert et al. (Reference Alibert, Carron and Fortier2013). In the latter work, an explicit N-body integrator was added ‘on top’ of the existing model that calculates the gravitational interaction and the collisions of concurrently forming protoplanets. Typically ten protoplanets per disc are considered in order to keep the computational time on a level that still allows us to conduct planetary population synthesis calculations. The effects of the gravitational interaction with the gaseous disc (eccentricity damping) are modelled as in Fogg & Nelson (Reference Fogg and Nelson2007), and it is assumed that the planetesimal surface density is uniform in overlapping feeding zones. The planetesimals are still described via a surface density and a mean dynamical state, and not as individual (super-) particles as it is the case in the model of Hellary & Nelson (Reference Hellary and Nelson2012). In this model, both the protoplanets and the planetesimals (represented by super particles) interact via an explicit N-body integration. The population synthesis models of Ida & Lin (Reference Ida and Lin2004a and Reference Ida and Lin2008b) also initially used the one-embryo-per-disc approximation but allowed for several generations of subsequently forming planets. In Ida & Lin (Reference Ida and Lin2010) and Ida et al. (Reference Ida, Lin and Nagasawa2013) a new semi-analytical approach was presented to describe the gravitational interaction of several protoplanets in a statistical way based on orbit crossing timescales. The advantage of such an approach is the small computational cost that is several orders of magnitude lower than direct N-body integrations. Here, one can clearly see the different origins of various global formation models that are used for statistical studies, like giant planet formation as in Pollack et al. (Reference Pollack, Hubickyj and Bodenheimer1996) for the Alibert, Mordasini and Benz models, N-body calculations for the Hellary & Nelson (Reference Hellary and Nelson2012) models, or a dedicated code for a rapid calculations of statistical results as in the case of the Ida and Lin models. An overview of the numerous effects that arise from the dynamical interactions of (proto) planets and external perturbers can be found in Davies et al. (Reference Davies, Adams and Armitage2013).

Illustrative output: formation tracks

An illustrative output of the population synthesis framework represented in Fig. 3 is shown in Fig. 12. It shows formation tracks in the mass–distance plane for the one-embryo-per-disc approximation. The non-isothermal type I migration model is used. Planetary embryos are inserted at a given starting distance into protoplanetary disc of varied properties with an initial mass of 0.6 Earth masses. They then grow by accreting planetesimals and gas, and concurrently migrate due to the interaction with the gas disc. The distribution of the final positions of the planets (at the moment the protoplanetary disc goes away) can be compared with the observed semi-major axis–mass distribution.

Fig. 12. Theoretical planetary formation tracks that show how planetary seeds (initial mass 0.6 Earth masses) concurrently grow and migrate. The colours indicate the different types of orbital migration (type I: brown: locally isothermal; red: adiabatic, unsaturated corotation torque; blue: adiabatic, saturated coronation torque; green: type II). The position of the planets at the moment in time that is shown (4.9 Myr) is indicated by black symbols. Some planets have reached the inner border of the computational disc at 0.1 AU.

One can see that the outcome of the formation process is of a high diversity, despite the fact that always exactly the same formation model is used. This is a basic outcome similar to the observational result. In the figure, one can, for example, find tracks that lead to the formation of hot Jupiters. Most embryos, however, remain at a low mass, since they cannot accrete a sufficient amount of planetesimals to start rapid gas accretion and become a giant planet. At an orbital distance of 0.2–1 AU, an overdensity of low-mass planets (M≲5M _⊕) can be seen. These are planets that are captured in the inner convergence zone (cf. Fig. 11). One also notes that almost all giant planets are inside of 1 AU, which is not in agreement with observations. These points to too rapid inward orbital migration in the model, meaning that the theoretical description of this process must be further improved. It is a typical result that the synthetic mass distribution (discussed in the next section) is in better agreement with the observational data than the synthetic semi-major axis distribution (e.g. Mordasini et al. Reference Mordasini, Alibert, Benz and Naef2009b).

Radial velocity: the planetary initial mass function

Among the many outputs that can be compared with observations, one of the most fundamental results of population synthesis is a prediction for the distribution of planetary masses. It is obvious that the planetary mass function has many important implications, including the question about the frequency of habitable extrasolar planets. In the left panel of Fig. 13, the planetary mass function is shown as derived from the high-precision radial velocity Harps survey of FGK dwarfs (Mayor et al. Reference Mayor, Marmier and Lovis2011). It makes clear that below a mass of approximately 30 Earth masses, there is a strong increase in the frequency. In particular, low-mass planets of with masses less than ~10M _⊕ are very frequent. The right panels shows the predicted mass function from an early population synthesis calculations of planets around a 1 M _⊙ star (Mordasini et al. Reference Mordasini, Alibert, Benz and Naef2009b, see also Ida & Lin Reference Ida and Lin2004a). The synthetic model predicted a large population of low-mass planets. A feature that is very interesting from a theoretical point of view is that also in the theoretical curve, there is a strong change in the frequency at a similar mass as in the observations. This is explained by the fact that in this mass domain (beyond the critical core mass), planets start to accrete nebular gas in a rapid, runaway process (Mizuno et al. Reference Mizuno, Nakazawa and Hayashi1978). They then quickly grow to masses M≳100M _⊕. It is unlikely that the protoplanetary disc disappears exactly during the short time during which the planet transforms from a Neptunian into a Jovian planet. This makes that planets with intermediate masses ≳30M _⊕ are less frequent (‘planetary desert’, first found in Ida & Lin Reference Ida and Lin2004a). The ‘dryness’ of the desert depends directly on the rate at which planets can accrete gas during the runaway phase (see Mordasini et al. Reference Mordasini, Mayor and Udry2011), while the mass where the frequency drops represents the mass where runaway gas accretion starts, that is the critical core mass which is of order of 15M _⊕, so that the total mass is about 30M _⊕ (core and envelope mass are approximately equal when the rapid accretion of gas sets in, Pollack et al. Reference Pollack, Hubickyj and Bodenheimer1996). We thus see how the comparison of synthetic and actual mass function constrains the core accretion model. We further note that while qualitatively, model and observation agree in the basic result that low-mass planets are very frequent, quantitatively the number of detectable planet at 1 and 0.1 m s⁻¹ is clearly underestimated in the model. This could partially be related to the one-embryo-per-disc approximation that is used in this early simulation (see Benz et al. Reference Benz, Ida, Alibert, Lin and Mordasini2013 for an updated synthetic mass function).

Fig. 13. Comparison of the observed and the synthetic planetary mass distribution. The left panel shows the distribution of planetary masses as found with high-precision radial velocity observations (Mayor et al. Reference Mayor, Marmier and Lovis2011). The black line gives the raw count, whereas the red line corrects for the observational bias against the detection of low-mass planets. The right panel shows the planetary mass function as found in an early population synthesis calculation (figure adapted from Mordasini et al. Reference Mordasini, Alibert, Benz and Naef2009b). The black line gives the full underlying population, while the blue, red and green lines are the detectable synthetic planets at a low (10 m s⁻¹), high (1 m s⁻¹) and very high (0.1 m s⁻¹) radial velocity precision.

Astrometry and microlensing: exploring different sub-populations

The observational constraints from astrometric observations are in principle similar to those from the radial velocity technique, with the difference that the actual mass is measured, and that the detection sensitivity increases with semi-major axis. To date, detecting extrasolar planets with this technique has proven difficult to achieve (e.g. Sahlmann et al. Reference Sahlmann, Segransan and Queloz2011), but is the ultimate goal of a number of ground based (e.g. Prima, Launhardt et al. Reference Launhardt, Queloz and Henning2008) and space-based missions (like the proposed Neat satellite, Malbet et al. Reference Malbet, Léger and Shao2011). In any case, the Gaia satellite is predicted to discover a very large number of extrasolar giant planets (several thousands at intermediate orbital distances of ~1–4 AU, Casertano et al. Reference Casertano, Lattanzi and Sozzetti2008). Since these discoveries will result from an unbiased, magnitude-limited survey with a well-defined detection bias (similar to the Kepler satellite), they will be extremely useful for statistical studies of giant planet formation.

Also results from the microlensing technique are important for statistical studies, since they probe the sub-population of low-mass planets at a few AU of orbital distance that is not accessible to other techniques. Microlensing is in some sense an extreme statistics-only method, because in its simplest form it only yields the distance of the planet in units of the star's Einstein radius and the ratio of the planet's mass to the (unknown) host star mass, meaning that very few physical information about an individual planet is revealed (this can change if a number of effects like the microlensing parallax, the orbital motion of the planet, or finite source effects can be measured; see, e.g. Gaudi Reference Gaudi2012).

For statistical studies, this is not necessarily a problem, provided that the number of detections is sufficiently high and that the observational detection bias is well known. A number of studies (e.g. Gould et al. Reference Gould, Dong and Gaudi2010; Cassan et al. Reference Cassan, Kubas and Beaulieu2012) have already used the microlensing discoveries to derive the frequency and a power-law exponent for the planetary mass function. From a theoretical point of view it seems rather unlikely that two parameters (normalization and power-law slope) are sufficient to describe the planetary mass function over three orders of magnitude in mass (from a few 10⁰ to a few 10³M _⊕) as it is currently made in the observational studies due to the low number of microlensing planets. From the core accretion paradigm one expects that at least four parameters are needed (this will still be a rough approximation only), because there are two different fundamental types of planets (solid planets and gas giants) for which different physical mechanisms determine the mass. The two types should therefore come with separate slopes and normalizations as one can already deduce from Fig. 13. As the number of microlensing planets is expected to increase thanks to ground and space-based observations with satellites like Euclid (Penny et al. Reference Penny, Kerins and Rattenbury2013) or Wfirst (Goullioud et al. Reference Goullioud, Content and Kuan2012), it will become possible to test this prediction observationally.

Transits

In the past few years, there was a very rapid increase of observational data coming from photometric observations, in particular thanks to the Kepler satellite. This new observational data provides important new impulses to planet formation and evolution theory, since it adds constraints that go beyond the position of a planet in the mass–distance diagram. Extended comparisons of theoretical and observational results derived from transit observations can be found, for example, in Howard et al. (Reference Howard, Marcy and Bryson2012), Mordasini et al. (Reference Mordasini, Alibert and Georgy2012b), Lopez & Fortney (Reference Lopez and Fortney2013b) or Marcy et al. (Reference Marcy, Isaacson and Howard2014) to name just a few. Here we concentrate on two important results.

Synthetic mass–radius diagram

The first is the mass–radius diagram. The observed mass–radius relationship was shown in subsection ‘Mass–radius diagram’. Here we discuss its synthetic counterpart. Figure 14 shows a comparison of the mass–radius relationship of actual and synthetic planets as found in a population synthesis that combines planet formation and evolution (modified from Mordasini et al. Reference Mordasini, Alibert and Georgy2012b). The global shape of the planetary mass–radius relationship can be understood from the core accretion paradigm and the basic properties of matter as expressed in the EoS: Low-mass planets can only accrete tenuous H/He envelopes since their Kelvin–Helmholtz timescale of envelope contraction during the formation phase is long compared to the typical lifetime of a protoplanetary disc (e.g. Ikoma et al. Reference Ikoma, Nakazawa and Emori2000). Therefore, the top left corner in the M–R plane remains empty, as no low-mass, gas-dominated planets come into existence.

Fig. 14. Mass–radius diagram of synthetic planets with a primordial H/He envelope at an age of 5 Gyr and a semi-major axis between 0.1 and 5 AU together with actual planets in- and outside of the Solar System with a well-known mass and radius, and a semi-major axis of at least 0.1 AU. The colours indicate the mass fraction of heavy elements Z in the synthetic planets. The black symbols, for example, correspond to solid-dominated low-mass planets which contain at most 1% of H/He, while the most massive planets (dark yellow) consist of at least 99% H/He. The other colours are: red: 1<Z≤5%, green: 5 <Z≤20%, blue: 20<Z≤40%, cyan: 40<Z≤60%, magenta: 60<Z≤80%, yellow: 80<Z≤95% and brown: 95<Z≤99% (figure modified from Mordasini et al. Reference Mordasini, Alibert and Georgy2012b).

Also the bottom right and middle part remains empty. For typical protoplanetary discs, this is simply a consequence of the fact that such discs do not contain enough mass in metals to form a solid Jovian-mass planet. For the MMSN, for example, one can roughly estimate a total amount of heavy elements of M _disc,MMSN×Z _⊙≈0.013M _⊙×0.015 (Hayashi Reference Hayashi1981; Lodders Reference Lodders2003) which corresponds to about 64 M _⊕. However, for a very massive metal-rich disc, we can in contrast have a disc gas mass that is still self-gravitationally stable of about 10% of the stellar mass and a metallicity [Fe/H]=0.5. Then one estimates for a 1 M _ε star a total mass of heavy elements of order of 1M _⊙×0.1×Z _⊙×10^0.5, which gives about 1600 M _⊕ (for a similar discussion, see Baraffe et al. Reference Baraffe, Chabrier and Barman2008). Obviously, it is not clear how much of this mass can be incorporated into one planet and at what point of the formation process this should happen.

In this context, it useful to note that a number of transiting giant planets seem to have extreme amounts of metals in their interior as estimated from their mass–radius relation: for Hat-P-20b Leconte et al. (Reference Leconte, Chabrier, Baraffe and Levrard2011) estimate of order of 340 M _⊕ of metals (water) in the interior; about 300–1000 M _⊕ of solids seem to reside in Corot-20b (Deleuil et al. Reference Deleuil, Bonomo and Ferraz-Mello2012, using the internal structure model of Guillot & Morel Reference Guillot and Morel1995). This indicates that large amounts of solids are indeed available in some discs and that some planets even manage to accrete a significant fraction of the totally available metal inventory. From these considerations it appears that just from the availability of metals in protoplanetary discs, it is not excluded that in principle a, say, 100 M _⊕ planet made entirely of ices could form. However, no such planet has been observed to date: all known planets in this mass range (M≳100M _⊕) have a mass–radius relation that shows that they contain significant amounts of H/He (see Fig. 2).

This is expected from the core accretion model since massive cores necessarily cause runaway gas accretion (e.g. Papaloizou & Terquem Reference Papaloizou and Terquem1999), so that the final composition of the forming planet contains significant amounts of envelope gas, at least if the cores form during the presence of the gaseous nebula. Thus, no massive purely solid planets come into existence that would populate the bottom right and middle parts of the mass–radius diagram. One further notes that the synthetic and most actual planets (both in- and outside of the Solar System) populate similar loci in the mass–radius plane.

From the position of a planet in the mass–radius relationship it is possible to deduce (within limits due to the degeneracies, cf. Rogers & Seager Reference Rogers and Seager2010) the bulk composition of a planet. The plot shows that depending on the mass range, there can be many different associated radii for one mass, reflecting a large diversity in interior compositions (in this case, fraction of heavy elements versus H/He). These different compositions are in turn due to the different formation histories. It is, for example, found that for the low planetesimal random velocities assumed here, planets at large distances typically contain a higher fraction of solid elements, since the mass of planetesimals available to accrete (the isolation mass) increases for typical disc models (radial slope of the planetesimal surface density) with distance. This correlation is preserved even under the action of disc migration. It might not be preserved if scattering is responsible for the formation of close-in planets. This could open a new possibility to distinguish different modes of formation for Hot Jupiters.

Impact of grain opacity on the planetary radius distribution

The second example of a synthetic result that can be compared with transit observations is the planetary radius distribution. It can serve as an example of how planet population synthesis can be used to study the global consequences of a specific physical mechanism (see Hasegawa & Pudritz Reference Hasegawa and Pudritz2014; Mordasini et al. Reference Mordasini, Klahr, Alibert, Miller and Henning2014).

Figure 15 compares the observed distribution of radii of planets inside of 0.27 AU as found by the Kepler satellite (Howard et al. Reference Howard, Marcy and Bryson2012) with the radius distribution in three different population syntheses. The three calculations are identical except for the value of f _κ. This parameter describes the reduction factor of the opacity due to grains suspended in the protoplanetary atmosphere during the formation phase relative to the interstellar value. A value of f _k=1 means that the full interstellar opacity is used (Bell & Lin Reference Bell and Lin1994), while f _κ=0 means that we are dealing with a grain-free gas where only molecular opacities contribute (Freedman et al. Reference Freedman, Marley and Lodders2008). These opacities are used when calculating the internal structure of the planets (subsection ‘Internal structure of the planetary gas envelope’) where the magnitude of the opacity is important for the rate at which planets can accrete primordial H/He envelopes (e.g. Ikoma et al. Reference Ikoma, Nakazawa and Emori2000). At low opacities, the liberated gravitational potential energy of the accreted gas can be more promptly radiated away, allowing the envelope to contract faster, so that new gas can be accreted. Specialized microphysical models of grain evolution predict that the opacity should be strongly reduced in protoplanetary atmospheres relative to the ISM because grains grow rapidly in the denser atmosphere and then settle into the deeper parts of the envelope where they are vaporized (Podolak Reference Podolak2003; Movshovitz & Podolak Reference Movshovitz and Podolak2008; Movshovitz et al. Reference Movshovitz, Bodenheimer, Podolak and Lissauer2010).

Fig. 15. Impact of the opacity due to grains in the protoplanetary gas envelope during the formation phase on the planetary radius distribution after 5 Gyr of evolution. In all panels, the blue line with error bars shows the bias corrected distribution of radii found by the Kepler satellite (data set from Borucki et al. Reference Borucki, Koch and Basri2011) according to the analysis of Howard et al. (Reference Howard, Marcy and Bryson2012). The black dotted lines are the preliminary analysis by the same authors of the updated Kepler data set of Batalha et al. (Reference Batalha, Rowe and Bryson2012). The red line shows the synthetic distribution for an opacity due to grains in the protoplanetary envelope equal to the ISM grain opacity, reduced relative to the ISM by a factor 0.003 and for vanishing grain opacity (left to right).

Considering the radius bin in Fig. 15 that contains the giant planets at about 1 Jovian radius (≈11 Earth radii), we see that with full-grain opacities, there are too few synthetic planets relative to the observations. With vanishing grain opacities, the efficiency of giant planet formation is on the contrary too high in the model relative to the data. A relatively good agreement with the observations is found with an f _κ=0.003 (middle panel), which is the value derived from fitting gas accretion timescales found with detailed grain growth models (Movshovitz et al. Reference Movshovitz, Bodenheimer, Podolak and Lissauer2010). Even if this comparison is preliminary (it is unclear whether the fitting value found for a specific simulation can generalized to the entire population, see Mordasini et al. Reference Mordasini, Klahr, Alibert, Miller and Henning2014), it shows that the grain opacity has statistically visible population-wide consequences. This opens the possibility to observationally test microphysical grain models.

Direct imaging: luminosity at young ages

The direct imaging technique measures the intrinsic luminosity of young Jupiters. This is interesting for planet formation and evolution theory because it is an observable quantity that is determined by the entropy of the gas in the interior of the planet. The entropy is in turn determined by the structure of the accretion shock of gas during the gas runaway accretion phase. Even more fundamentally, the entropy state could be related to the basic giant planet formation mechanism (core accretion versus gravitational instability, see, e.g. Galvagni et al. Reference Galvagni, Hayfield and Boley2012; Spiegel & Burrows Reference Spiegel and Burrows2012).

The number of planets detected by this method has increased significantly in the last few years (see Marleau & Cumming Reference Marleau and Cumming2014, for an overview). In absolute numbers it is still low compared to the radial velocity and transit method, but this should change in the coming years. It is nevertheless already now possible to use statistical methods to study the discoveries made by this method as shown in Fig. 16.

Fig. 16. Comparison of theoretical results and combined observational constraints on the nature and formation of β Pictoris b from radial velocity observations and direct imaging (adapted from Bonnefoy et al. Reference Bonnefoy, Boccaletti and Lagrange2013). The plot shows the mass–distance diagram of synthetic planets predicted for the specific properties of the host star. The blue points represent all synthetic planets while the large green circles and red open symbols indicate planets that agree at an age of 12 Myr within two respectively one σ with the constraints derived from direct imaging (T _eff and L). The grey-shaded region corresponds to the mass domain that agrees with the limits derived from radial velocity observations.

The figure show the mass–distance diagram for a synthetic population that is customized for the star β Pictoris which is a young (~12 Myr) A5 V star visible to the naked eye. Customized means that the quantities that are known for this specific star (its mass and metallicity) are fixed at the observed values, in contrast to normal population syntheses where these quantities are Monte Carlo variables. Additionally, the properties of the synthetic planets are studied at the actual age of the star.

β Pictoris is orbited by a directly imaged companion at a semi-major axis of about 8–10 AU (Lagrange et al. Reference Lagrange, Gratadour and Chauvin2009, Reference Lagrange, Bonnefoy and Chauvin2010). The recent analysis of the near-IR spectral energy distribution by Bonnefoy et al. (Reference Bonnefoy, Boccaletti and Lagrange2013) shows that the companion has a luminosity of log(L/L _ε)=−3.87±0.08, an effective temperature T _eff=1700±100 K and a surface gravity log(g)=4.0±0.5 (cgs units). According to ‘hot start’ models (e.g. Burrows et al. Reference Burrows, Marley and Hubbard1997; Baraffe et al. Reference Baraffe, Chabrier, Barman, Allard and Hauschildt2003) that assume an arbitrary (high) value of the entropy after formation without considering the actual formation phase, these values correspond to a mass of the companion of ${\rm 9}_{ {- \rm 2}} ^{{ + {\rm 3}}} $ Jovian masses. Radial velocity observations (Lagrange et al. Reference Lagrange, De Bondt and Meunier2012) show that the mass of the companion must be less than 10 and 25 Jovian masses for semi-major axes of 8 and 12 AU, respectively. The allowed mass domain is indicated in Fig. 16.

The mass–distance diagram as a function of stellar mass was studied in the past with models that predict only these quantities (Ida & Lin Reference Ida and Lin2005; Alibert et al. Reference Alibert, Mordasini and Benz2011). The new result shown in Fig. 16 is that the theoretical model predicts besides mass and semi-major axis also the luminosity and effective temperature. Thus it becomes possible to combine the constraints from both direct imaging and radial velocity. The plot shows that there are indeed a number of planets that agree in all four fundamental properties (a, M, L and T _eff) with the observed values, indicating that β Pictoris b could have formed by core accretion and that it has a mass in the planetary mass domain of about ten Jovian masses. Note that this result directly depends on the assumption in the formation model that the planetesimal random velocities are low as in Pollack et al. (Reference Pollack, Hubickyj and Bodenheimer1996) and that gap formation does not lead to a reduction of the gas accretion rate (Kley & Dirksen Reference Kley and Dirksen2006).

The simulations are conducted assuming that during gas runaway accretion, the gas accretion shock radiates all potential energy liberated by the infalling gas (‘cold accretion’). Still there are planets that have a luminosity that agrees with the observed value of about log(L/L _ε)=−3.9. This might seems surprising at first, because under the same assumption, Marley et al. (Reference Marley, Fortney, Hubickyj, Bodenheimer and Lissauer2007) had in contrast found that the post-formation luminosities in the relevant mass domain are always less than log(L/L _⊙)=−5. While investigating the reason for the discrepancy, it was found that the difference between the simulations stems from different core masses M _core: in the Marley et al. (Reference Marley, Fortney, Hubickyj, Bodenheimer and Lissauer2007) simulations, the core masses are less than 19 M _⊕, whereas the synthetic planets that agree with the observations in Fig. 16 have much more massive cores exceeding ~100 M _⊕. For identical core masses, the two models agree well. This insight from the population synthesis led to a systematic dedicated study of the dependency of the post-formation entropy and luminosity on the core mass (Mordasini Reference Mordasini2013). This is therefore another example how global models feed back into specialized ones. It was found that the post-formation luminosity of massive giant planets is very sensitive to the core mass due to a self-amplifying mechanism (see Mordasini Reference Mordasini2013 and Bodenheimer et al. Reference Bodenheimer, D'Angelo, Lissauer, Fortney and Saumon2013 for details).

It is currently unclear whether such massive cores can actually form, especially at large orbital distances. Several factors play a role, like for example the settling of dissolved planetesimal material towards the centre of the planet (Iaroslavitz & Podolak Reference Iaroslavitz and Podolak2007) and the timescale of core accretion which might get speeded up if mainly small bodies are accreted (see, e.g. Dodson-Robinson et al. Reference Dodson-Robinson, Veras, Ford and Beichman2009; Ormel & Klahr Reference Ormel and Klahr2010; Kratter & Murray-Clay Reference Kratter and Murray-Clay2011; Lambrechts & Johansen Reference Lambrechts and Johansen2012). At least for some close-in transiting giant planets, very large core masses (≳100M _⊕) have been inferred as mentioned from the observed mass–radius relationship (Miller & Fortney Reference Miller and Fortney2011; Deleuil et al. Reference Deleuil, Bonomo and Ferraz-Mello2012).