Global models of planet formation and evolution

Abstract

Despite the strong increase in observational data on extrasolar planets, the processes that led to the formation of these planets are still not well understood. However, thanks to the high number of extrasolar planets that have been discovered, it is now possible to look at the planets as a population that puts statistical constraints on theoretical formation models. A method that uses these constraints is planetary population synthesis where synthetic planetary populations are generated and compared to the actual population. The key element of the population synthesis method is a global model of planet formation and evolution. These models directly predict observable planetary properties based on properties of the natal protoplanetary disc, linking two important classes of astrophysical objects. To do so, global models build on the simplified results of many specialized models that address one specific physical mechanism. We thoroughly review the physics of the sub-models included in global formation models. The sub-models can be classified as models describing the protoplanetary disc (of gas and solids), those that describe one (proto)planet (its solid core, gaseous envelope and atmosphere), and finally those that describe the interactions (orbital migration and N-body interaction). We compare the approaches taken in different global models, discuss the links between specialized and global models, and identify physical processes that require improved descriptions in future work. We then shortly address important results of planetary population synthesis like the planetary mass function or the mass–radius relationship. With these statistical results, the global effects of physical mechanisms occurring during planet formation and evolution become apparent, and specialized models describing them can be put to the observational test. Owing to their nature as meta models, global models depend on the results of specialized models, and therefore on the development of the field of planet formation theory as a whole. Because there are important uncertainties in this theory, it is likely that the global models will in future undergo significant modifications. Despite these limitations, global models can already now yield many testable predictions. With future global models addressing the geophysical characteristics of the synthetic planets, it should eventually become possible to make predictions about the habitability of planets based on their formation and evolution.

Introduction

Thanks to the progress of observational techniques in the past decades, we are the first generation of human beings that has had the technological capabilities to answer the question about the existence of planets around other stars (Mayor & Queloz 1995). Since then, there was an enormous increase in observational data on extrasolar planets. The latest observational results from different detection techniques (Borucki et al. 2011; Mayor et al. 2011; Cassan et al. 2012) even indicate that the presence of planets is the rule rather than the exception, at least around solar-like stars. However, this increase in observational data on extrasolar planets does not mean that we can now fully explain how these planets came into existence from a theoretical point of view. On the contrary, many observational findings on extrasolar planets were not predicted from planet formation and evolution theory, or were even in contrast to it, showing that this field is still in its infancy.

A young method to improve the theoretical understanding of planet formation is planetary population synthesis. It is a statistical method that makes it possible to improve the theoretical understanding of the physics governing planet formation and evolution using statistical comparisons to observational constraints provided by the population of extrasolar planets. With this approach the global effects of many key physical processes occurring during planet formation and evolution can be put to the observational test, something which is notoriously difficult in astronomy, since the objects that are studied are far away and only accessible via the radiation they (or their host star) emit.

In this paper, we review global models of planet formation and evolution that are used in such planetary population synthesis calculations (Ida & Lin 2004a, b, 2005, 2008a, b, 2010; Thommes et al. 2008; Miguel & Brunini 2009; Mordasini et al. 2009a, b, 2012a, b, c; Alibert et al. 2011, 2013; Hellary & Nelson 2012; Benz et al. 2013; Forgan & Rice 2013; Galvagni & Mayer 2013; Hasegawa & Pudritz 2013; Ida et al. 2013). In this context, ‘global’ means that these models can directly predict planetary properties based on properties of the protoplanetary disc in which the planets form. To do so, they unite in one model the results of many specialized models that address a specific important mechanism occurring during planet formation such as accretion or migration. Here we concentrate on the physical description of these processes as included in the global models because this is the key ingredient of the entire population synthesis approach.

Observational motivation

Currently, there are approximately 1500 confirmed exoplanets known¹ that were mostly detected with the spectroscopic radial velocity technique or photometric transit observations. Additionally, there are about four thousand candidates from the Kepler satellite (Borucki et al. 2011) that were detected with extremely precise transit photometry. These detections have revealed an exiting diversity in the properties of planetary companions that was not expected from the structure of our own planetary system, the Solar System. The detections have, however, not only revealed a surprising diversity, but also a number of interesting correlations and structures in the properties of the planets.

These insights were in particular possible thanks to the large number of planets now known. For the first time, this allows us to look at the extrasolar planets no more solely as single objects. Instead, it is possible to look at them as a population that is characterized by a number of statistical properties. Important examples are the distributions of masses, semi-major axes, radii, eccentricities and the relations between theses quantities. Understanding these statistical properties from the point of view of planet formation and evolution theory is one of the fundamental goals of the models presented here. The understanding that can be gained in this way also feeds back into the way we understand our own Solar System. Also the Solar System itself provides a large body of precise observational constraints against which planet formation models must be compared. But developing a theory that is focused on one planetary system can be misleading; the discovery of exoplanets that are very different from any Solar System planet has shown that.

The special interest in a statistical population wide approach also comes from the fact that the knowledge about a single extrasolar planet is often limited. For the large majority of the extrasolar planets, still only a few orbital elements (semi-major axis, eccentricity, etc.) and a minimum mass are known (or a radius, but no mass in the case of most Kepler candidates). In order to benefit also from the large number, but individually limited data sets, statistical methods are necessary. Having this ability is important since several future surveys such as the Gaia space mission or the Sphere and Gpi direct imaging surveys (Beuzit et al. 2008; McBride et al. 2011) will yield additional statistical data sets.

For a handful planets this has partially changed in the last few years, and a first rough geophysical characterization of some extrasolar planets has become possible by multi-band photometry or spectroscopy (e.g. Richardson et al. 2007; Konopacky et al. 2013). These are typically planets around bright stars for which both the mass and the radius are known (or alternatively the intrinsic luminosity in the case of directly imaged planets). In this case, more observational constraints can be derived like the mean density or the atmospheric structure and composition. These planets are investigated in detail, and will likely have a special role in the next decade of extrasolar planet study. This observational progress was the motivation for some recent work on the theoretical side discussed in this paper, which is the extension of an existing planet formation model (Alibert et al. 2005a) into a self-consistently coupled planet formation and evolution model (Mordasini et al. 2012a). With such a combined model one can predict all major observable characteristics of a planet. We will discuss this topic and how it can be incorporated into statistical studies in the section ‘From detection to characterization’.

The fundamental assumption behind the planetary population synthesis method is that the observed statistical properties of exoplanets (like the aforementioned distributions) can be explained by the action of the always identical governing physical processes during the formation phase of the planets, but under different initial conditions. The initial conditions for the planet formation process are the properties of the protoplanetary disc which are found to surround most newly born stars (e.g. Haisch et al. 2001; Fedele et al. 2010). Observations of such circumstellar discs show that they also have, as planets, a wide variety of properties in terms of their most important characteristics such as their mass, lifetime or radius (e.g. Andrews et al. 2010). For an individual planetary system, the properties of the protoplanetary disc from which it formed are mostly unknown, except maybe for the dust-to-gas ratio in the disc that is likely correlated with the stellar metallicity that can be measured spectroscopically today. This means that the initial conditions are only known in a statistical sense, which again makes a statistical approach appropriate. It is clear that this fundamental assumptions neglects that stars, and therefore planets, usually do not form in isolation, but in stellar clusters. This environment can influence the formation process. Two examples are the impact of close-by massive stars on the lifetime of a protoplanetary disc (Adams et al. 2004), or the gravitational perturbation by a passing star or temporary binary companion (Malmberg et al. 2011).

Observational constraints

In comparison to other studies in the domain of planet formation theory, population synthesis is directly at the interface between theory and observation, since it is a goal to connect these two domains. It is therefore interesting to briefly discuss two central observational results that are important for the theoretical studies, namely, the semi-major axis–mass and the mass–radius diagram.

It is clear that besides these two results, there is very large number of additional statistical constraints that can be deduced from the extrasolar planets such as the distributions of radii, eccentricities and luminosities, as well as the relations between all these quantities. Further important observational constraints are the frequency and properties of planets as a function of host star properties such as mass and metallicity, the mean spacing between planets in multiple systems, the alignment of the planets among themselves and with the stellar equator, the frequency of planets in mean motion resonances and so on. A recent review on the observed properties of extrasolar planets can be found in Fischer et al. (2014).

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 2004a).

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. 2011).

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. 2005; Howard et al. 2010; Mayor et al. 2011). 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. 2011). 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. 2012; Fressin et al. 2013). 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 2004b). 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. 2011).

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 2004b; Thommes et al. 2008; Mordasini et al. 2009a) and can be compared to theoretical results with Kolmogorov–Smirnov tests (Mordasini et al. 2009b). 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. (2007) or Spiegel & Burrows (2012). In Mordasini (2013) 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. 2012b). The combination of measurements of the radius of a transiting planet (first by Henry et al. 2000 and Charbonneau et al. 2000) 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 known². 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. (2011).

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 2011; Miller & Fortney 2011) 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. 2003), dissipation of stellar irradiation deep in the atmosphere (Guillot & Showman 2002), double diffusive convection (Chabrier & Baraffe 2007), enhanced atmospheric opacities (Burrows et al. 2007) and ohmic dissipation (Batygin et al. 2011).

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. 2007; Seager et al. 2007). 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 2010). 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. (2006) or Burrows et al. (2007)?

From detection to characterization

The above three questions are an observational motivation to extend a global planet formation model into a coupled formation and evolution model as in Mordasini et al. (2012c). It then becomes possible to calculate radii based on the bulk composition obtained during formation, which makes new observational constraints usable for population synthesis. Such a model can then be used to study the population-wide mass–radius relationship and to compare with the radius distribution found by the Kepler satellite (Mordasini et al. 2012b). Compared to other, well-established planet evolution models (e.g. Burrows et al. 1997; Baraffe et al. 2003), the evolutionary model of Mordasini et al. (2012c) is, however, still significantly simplified in several aspects as discussed in subsection ‘Atmosphere of the planet’ In view of future observations yielding very precise radii (e.g. by the photometric CHaracterizing ExOPlanet Satellite Cheops, Broeg et al. 2013) it will probably be necessary to find more accurate physical descriptions also in global models.

The mass–radius diagram represents, in a prototypical way, the transition of the focus from pure exoplanet detection to beginning exoplanet characterization in the past few years. Besides the M–R relationship, there was recent observational progress towards characterization in two other domains.

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. 2008) or β Pictoris (Lagrange et al. 2010) have triggered numerous theoretical studies regarding their formation (e.g. Dodson-Robinson et al. 2009; Kratter et al. 2010). 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. 2011). 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. 2007) while the accretion of high entropy material leads to a ‘hot start’ with a high luminosity and large radius (e.g. Burrows et al. 1997; Baraffe et al. 2003). Recently, Spiegel & Burrows (2012) 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 2004a; Mordasini et al. 2009a). 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. 2005a, 2013, Mordasini et al. 2012b, c).

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. 2007). 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 1999). 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. (2001) and Alibert et al. (2005b). 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. 1999) 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. 1999) or the incorporation of noble gases in the form of clathrate hydrates in crystalline ice (Gautier et al. 2001).

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. 2009). 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 2010). Such different compositions and the associated different mean molecular weight lead to different transmission spectra (e.g. Bean et al. 2011). 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. 2011; Berta et al. 2012). Recent precise observations could rule out cloud-free atmospheres even of high mean molecular mass (Kreidberg et al. 2014). 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. 2006) indicate that highly enriched atmospheres are a typical outcome for a planet such as GJ1214b with clear consequences for spectroscopy (Fortney et al. 2013). But in general, despite the fact that some processes have been identified (Madhusudhan et al. 2011; Öberg et al. 2011; Thiabaud et al. 2014), 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.

Elements of the population synthesis method

Planetary population synthesis as a suitable method for statistical studies of planet formation and evolution was introduced in the pioneering work of Ida & Lin (2004a). A similar framework, but intended for more quantitative comparisons was established in Mordasini et al. (2009a). The basic idea is to run a global planet formation model for varying initial conditions. With this framework the population-wide, statistical consequences of a theoretical description of a specific physical mechanism can be studied and compared with the population of actual extrasolar planets. Examples are specialized models of type I migration or of grain growth in protoplanetary atmospheres, see subsection ‘Transits’. The possibility to test specialized models is an important aspect of population synthesis. These specialized models are typically more complex and contain more subtleties than their simplified counterpart embedded in a global model. But if the simplified counterpart is still able to capture the essence of the original specialized model, then population synthesis is often the only possibility to test them observationally. A framework for population synthesis typically consists of the following elements that are shown in the flowchart of the method in Fig. 3.

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. 2006; Lovis et al. 2006).

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 1996; Andrews et al. 2010), the dust-to-gas ratio (assumed to be correlated with the stellar [Fe/H], Santos et al. 2004; Fischer & Valenti 2005) and the lifetime of the disc (Haisch et al. 2001; Mamajek 2009; Fedele et al. 2010). 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. 2001; Miguel et al. 2011). 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 2007; Andrews et al. 2010) 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. 1998). 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 2007). 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 2005).

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. 2011), the Kepler satellite (Borucki et al. 2011), and also microlensing surveys (e.g. Gould et al. 2010; Cassan et al. 2012) 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 2004b; Mordasini et al. 2012a), 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 1997; Santos et al. 2004; Fischer & Valenti 2005), a correlation that can be reproduced with formation models based on the core-accretion theory (Ida & Lin 2004b;Mordasini et al. 2009b).

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.

Global planet formation and evolution models

We now come to the description of the physics of global models, which is the central subject of this paper. On the observational side, usually only the initial conditions (the protoplanetary discs) and the final outcomes of the planet formation process (the planets) are accessible to observations, even if in future direct imaging and Alma might allow us to observe planet formation as it happens for a limited number of cases (e.g. Quanz et al. 2013). With global models it is possible to bridge this gap at least on the theoretical side. Many elements of modern planet formation theory that hold to this day were first developed by Safronov (1969). Other works that laid the foundations for the theoretical descriptions used here include Shakura & Sunyaev (1973); Lynden-Bell & Pringle (1974); Perri & Cameron (1974); Mizuno et al. (1978); Goldreich & Tremaine (1979); Hayashi (1981); Bodenheimer & Pollack (1986); Lin & Papaloizou (1986); and Lissauer (1993).

Global numerical models of planet formation (like Ida & Lin 2004a, 2008a, b; Alibert et al. 2005a; Thommes et al. 2008; Hellary & Nelson 2012; Mordasini et al. 2012c) try to cover all major mechanisms that govern the formation and evolution process, and try to follow the formation process from its beginning to its end. The various mechanisms like accretion or migration must be treated in an interlinked way because they happen on similar timescales and feed back into each other. Ideally, the models would start with a protoplanetary disc at a very early stage when the solids are in the form of micrometre-sized dust grains (or in principle even earlier with the collapse of the cloud), and yield as an output full-blown planetary systems with fully characterized planets at an age of several billions of years. In reality this is impossible, but it means that also the evolution of the planets over long timescales should be modelled, since in most cases, we observe planets a long time after they have formed. This applies both to the long-term evolution of the internal structure of a planet (its cooling and contraction) and to the secular evolution of the orbits due to gravitational interactions and tides.

The fact that the models must follow the planetary systems during the entire formation process is the reason why global models are typically 1D (or at most 2D), for example, in the description of the protoplanetary disc (assumed to be axisymmetric) or the internal structure of the planets (assumed to be spherically symmetric). More realistic 2D or 3D hydrodynamical simulation can at least currently not be used to simulate thousands of different initial conditions (protoplanetary discs) over their entire lifetime.

Most global models used to date in population synthesis calculations are based on the core accretion paradigm (Perri & Cameron 1974; Mizuno et al. 1978). Core accretion states that first, solid cores form, some of which later accrete massive gaseous envelops to become giant planets (bottom-up process). The remaining cores collide to form both ice giants and terrestrial planets (for an overview of this sequential picture of planet formation, see, e.g. Papaloizou & Terquem 2006; Mordasini et al. 2010). First statistical considerations based on the competing gravitational instability model where giant planets form directly from a gravitational instability in the protoplanetary gas disc (Cameron 1978; Boss 1997) were recently made in Janson et al. (2011, 2012); Forgan & Rice (2013) and Galvagni & Mayer (2013).

Global models address the different physical processes in a number of interlinked sub-models. Figure 3 lists the 11 sub-models of the combined formation (Alibert et al. 2005a, 2013) and evolution model (Mordasini et al. 2012b, c), on which we concentrate in the following. However, we also compare with other global models and review important ongoing and future work. Each sub-model is relatively simple, but the interaction of them leads to a considerable complexity. We next discuss the physics included in the different sub-models. They can be split in three classes: models for the protoplanetary disc, for one (proto) planet and for the interactions (migration and N-body interaction).

Vertical structure of the protoplanetary disc

A first sub-model calculates the structure and evolution of the gaseous protoplanetary disc (Lynden-Bell & Pringle 1974; Papaloizou & Terquem 1999; Alibert et al. 2005a). 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. 2011) is provided by α-viscosity models (Shakura & Sunyaev 1973) that describe the disc as a rotating viscous fluid. In the model of Alibert et al. (2005a), 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 1999)

(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 (1973) 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. (2014). 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 1991) lead in contrast to a less concentrated heating, and thus to a different vertical temperature profile than in an α model (Flock et al. 2013). 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. 2013). This (and other differences, see Turner et al. 2014) 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 2009). 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. 2012)

(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 2005)

(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. (2012) set the flaring angle d ln H/d ln r for simplicity to the equilibrium value of 9/7 (Chiang & Goldreich 1997) 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 (1974) (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. 2013) has two origins (Mordasini et al. 2012b): 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. 2003). 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. 2001). Physically, it depends on the number of and distance to massive stars in the vicinity of the host star (Adams et al. 2004).

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. 1994)

(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. 1994). 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 condition³. Early population synthesis calculations (Ida & Lin 2004a; Mordasini et al. 2009b) typically used an initial radial profile inspired by the minimum mass Solar nebula (MMSN) (Weidenschilling 1977; Hayashi 1981) which is a power-law for Σ falling as r^{− 3/2}. Recent sub-millimetre observation of protoplanetary discs (e.g. Andrews et al. 2010) 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. 2010). Such profiles were used as initial conditions in more recent syntheses like Mordasini et al. (2012b) and Alibert et al. (2013) 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. 2012b).

The global models of Ida & Lin (2004a) and Hellary & Nelson (2012) 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. 2010; Kretke & Lin 2012; Dittkrist et al. 2014, 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. 2010, 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. 1997; Kretke & Lin 2010). 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. 2013), or the occurrence of layered accretion (Gammie 1996) 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. 2007). In the context of global planet formation models, this was studied in Ida & Lin (2008b) and Hasegawa & Pudritz (2013). Furthermore, Uribe et al. (2011) 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. 1997). 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. 2008). 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 1993; Ormel et al. 2010). 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. (2013). In earlier population syntheses, Mordasini et al. (2009b) 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 1993). This ignores the potential formation of a gap in the disc of planetesimals (Shiraishi & Ida 2008). The decrease of the planetesimal surface density is then simply (Thommes et al. 2003)

(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 1981)

(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 2003) 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 (1981). 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. 2004).

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 2004a; Hellary & Nelson 2012) use similar prescriptions for the initial surface density distribution. The model of Ida & Lin (2008b) 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. 2014 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 2013). 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. 2007; Birnstiel et al. 2010); these are not mutually exclusive processes (Johansen et al. 2014), (ii) the mode of concentration of small bodies in the turbulent disc on large (Johansen et al. 2006) or small scales (Cuzzi et al. 2008), (iii) the strong dependency of planetesimal formation rates on the dust-to-gas ratio (e.g. Johansen et al. 2009; Cuzzi et al. 2010) 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. 2010; Leinhardt & Stewart 2012), (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. 2004), 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 2004; Ciesla & Cuzzi 2006), 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. 2003; Chambers 2008) so that at least three types of bodies must be included in the solid disc model (embryos, planetesimals and fragments, see Ormel & Kobayashi 2012) 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 2005; Ida et al. 2008; Ormel & Okuzumi 2013).

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. (2007); Carter-Bond et al. (2010); Birnstiel et al. (2012) or Chambers (2014).

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 1969) 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 (1992). 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 (2003). In the population syntheses of Mordasini et al. (2009b) and (2012b), the original expressions of Pollack et al. (1996) 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. 2003). An updated description with a more detailed description of the interactions between protoplanet and planetesimals can be found in Fortier et al. (2013), 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 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012; Chambers 2014) 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 (2013) or Chambers (2014). 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 (2004a) 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 1993; Ormel et al. 2010). As they assume smaller planetesimals (km-sized instead of 100 km), the resulting core growth timescales are nevertheless relatively similar (Mordasini et al. 2009a). The global model of Hellary & Nelson (2012) 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. (2014) and Chabrier et al. 2014). The internal structure is found by integrating the standard equations for planetary interiors (e.g. Bodenheimer & Pollack 1986) 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 2003; Movshovitz et al. 2010; Cuzzi et al. 2014, 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 2012; Bodenheimer et al. 2013). Note that the models of Mordasini et al. (2012c) and Alibert et al. (2013) 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. (1996). 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. 2000). Once the crossover mass is reached which corresponds to the critical core mass in the early strictly static calculations (like Perri & Cameron 1974; Mizuno et al. 1978), $\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 shock⁴. 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. 1999; D'Angelo et al. 2011). 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 2008)

(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. 2012c)

(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. 1999; Lubow & D'Angelo 2006). Note that there is considerable uncertainty about the efficiency of gas accretion in the disc-limited regime (Benz et al. 2013). 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 2004a; Thommes et al. 2008; Hellary & Nelson 2012). 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. 2000), 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 2008). This can have directly visible consequences in population syntheses, for example, in the predicted planetary mass function (Miguel & Brunini 2008). 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 2003; Movshovitz & Podolak 2008; Movshovitz et al. 2010) 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. (2014) 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 (1986); Pollack et al. (1996) or Lissauer et al. (2009) which are the conceptual origin of this model (see Helled et al. 2013 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. 2012c) can bring a global model close to classical models of (giant) planet long-term evolution (cooling and contraction) like Burrows et al. (1997) or Baraffe et al. (2003). 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 (2010) discussed below. The formation of a deep radiative zone that is characteristic for strongly irradiated giant planets is visible (Guillot & Showman 2002). 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. 2000).

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. 2009). No atmosphere in the classical sense exists⁵, 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 (2005).

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. 2007). 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. 2000)

(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. (1980). 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 2006; Ayliffe & Bate 2009).

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. (2012b) model this phase with a simple grey atmosphere so that (e.g. Guillot 2005)

(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. (1997); Baraffe et al. (2003, 2008) 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. 2000; Mordasini et al. 2012b). 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 2002). A better atmospheric model than the grey atmosphere is the semi-grey approximation of Guillot (2010). 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. 2009) 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 2008).

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. (2008) 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 (2010). 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. (2011) for non (or weakly)-irradiated giant planets, following the procedure explained in Chabrier & Baraffe (1997). 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. 2006). A similar sub-model was described by Podolak et al. (1988) for the classical giant planet models of Pollack et al. (1996). 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. 2013). The sub-model calculating the interaction (see Mordasini et al. 2006, 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 1994; Crawford et al. 1995) or the Tunguska event (Chyba et al. 1993).

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 1958)

(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. 1995). 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 1992). 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 1994).

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. 1995). 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 1992) 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. 1993), 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 1994; Korycansky et al. 2000). Mordasini et al. (2006) 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 1984; Youngs 1989)

(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 1984; Youngs 1989), 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. 1993).

Figure 9 from Mordasini et al. (2006) 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. 1988), for large impactors mainly due to their lower tensile strength (by about a factor 100, Chyba et al. 1993). 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 1994), 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. 2006).

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. 1996). 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. 2010; Hori & Ikoma 2011), 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 2012). Including these effects in future work along the lines demonstrated by, e.g. Hori & Ikoma (2011) and Iaroslavitz & Podolak (2007) 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. 2011; Bonnefoy et al. 2013; Fortney et al. 2013; Konopacky et al. 2013).

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. 2012c).

The sub-model for the solid core used in the global model of Fig. 3 was originally developed in Figueira et al. (2009) 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. (2007), 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. (2009) used the more accurate tabulated equations of state (EoS) of Fortney et al. (2007), the global model of Mordasini et al. (2012b) employs the simpler, but more widely applicable EoS of Seager et al. (2007) 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 1969). 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. 2011).

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. (2006).

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. (2012b) the radiogenic luminosity can be modelled according to the law of radioactive decay assuming that the mantle has a chondritic composition (Urey 1955; Lowrie 2007). 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. 2011). 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. 2006; Lopez & Fortney 2013b). 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-up⁶. 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. 2007). Further improvements could be to include the dissolution of the solid core in giant planets (e.g. Guillot et al. 2004; Wilson & Militzer 2012) and the outgassing of secondary atmospheres given the core's composition acquired during formation (e.g. Elkins-Tanton & Seager 2008). 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. 2011).

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. 2003; Baraffe et al. 2004; Erkaev et al. 2007; Murray-Clay et al. 2009; Owen & Jackson 2012). 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. 2012; Lopez & Fortney 2013a; Owen & Wu 2013). 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. 2008; Rogers et al. 2011; Mordasini et al. 2012c).

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. 2009; Owen & Jackson 2012).

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. 2005)

(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 2007), 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 2012).

In the EUV regime, two sub-regimes exist (Murray-Clay et al. 2009). 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. 1981)

(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. 2009), and K _tide is a factor to take into account that mass only needs to reach the Hill sphere to escape (Erkaev et al. 2007).

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. 2009)

(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. (2009).

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. 2014).

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 2013). 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. 2014).

Orbital migration

The discovery of a Jovian planet at an orbital distance of only 0.05 AU from its star by Mayor & Queloz (1995) was a surprise because theoretical planet formation models had rather predicted (e.g. Boss 1995) 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 1980), migration due to the planetesimal disc (e.g. Levison et al. 2010; Ormel et al. 2012) or Kozai migration with tidal circularization (Kozai 1962; Fabrycky & Tremaine 2007). Another possible mechanism producing close-in planets is planet–planet scattering followed by tidal circularization (e.g. Ivanov & Papaloizou 2004; Papaloizou & Terquem 2006; Beaugé & 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 2006; Valsecchi & Rasio 2014, but see Rice et al. 2012 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. 2007; Triaud et al. 2010). Such misaligned orbits are not expected from classical disc migration (but see also Batygin 2012). 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 2014) 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. 2013; Ida et al. 2013). The main result of the large body of studies addressing disc migration (see Kley & Nelson 2012, 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 1986; Tanaka et al. 2002). The total torque Γ_tot can be expressed in a form like (Paardekooper et al. 2010)

(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. (2002) 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 2008a; Thommes et al. 2008; Mordasini et al. 2009b; Alibert et al. 2011) 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 2008; Casoli & Masset 2009; Kley et al. 2009; Paardekooper et al. 2010) 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. 2014). These more realistic type I descriptions were used in several recent population synthesis simulations (Hellary & Nelson 2012; Dittkrist et al. 2014).

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. 2010; Sándor et al. 2011; Kretke & Lin 2012; Dittkrist et al. 2014).

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. 2010). 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. 2012).

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 1986). 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. (2006) derived a combined criterion that was used in Mordasini et al. (2012b) while earlier population syntheses partially only used the thermal criterion.

Type II migration itself comes in two sub-regimes (Armitage 2007). 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 2007). This means that inside of the radius of maximum viscous couple (or velocity reversal, Lynden-Bell & Pringle 1974), 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 2009). The resulting slow-down of planets is important for the final semi-major axis distribution of giant planets (Mordasini et al. 2009b).

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 1986; Alexander & Armitage 2009). 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 2001).

The global formation model of Hellary & Nelson (2012) uses a similar description of non-isothermal type I migration and type II migration as presented here. Also Thommes et al. (2008) use the isothermal type I migration rates of Tanaka et al. (2002), but obtain the type II migration rate based on the more accurate impulse approximation (Lin & Papaloizou 1986). Their disc model is a 1D viscously evolving model, similar as presented in subsection ‘Radial structure of the protoplanetary (gas) disc’. Ida & Lin (2008a) study the global effects of type I migration by population synthesis calculations. They use the migration description of Tanaka et al. (2002) but reduce the rate by an arbitrary factor ≤1, an approach very similar to the one of Mordasini et al. (2009b). 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 1993) 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 2001). 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. 2006). The population synthesis calculations of Mordasini et al. (2009b and 2012b) 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. 2011) 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 (2010), Ida et al. (2013) and Alibert et al. (2013). 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 (2007), 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 (2012). 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 (2004a and 2008b) also initially used the one-embryo-per-disc approximation but allowed for several generations of subsequently forming planets. In Ida & Lin (2010) and Ida et al. (2013) 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. (1996) for the Alibert, Mordasini and Benz models, N-body calculations for the Hellary & Nelson (2012) 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. (2013).

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. 2009b).

Comparisons with observation

In this section, we discuss important selected comparisons between theoretical and observed statistical properties. Thanks to the coupling of planet formation and evolution in the global model as shown in Fig. 3, it is possible to compare with all major observational techniques.

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. 2011). 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. 2009b, see also Ida & Lin 2004a). 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. 1978). 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 2004a). The ‘dryness’ of the desert depends directly on the rate at which planets can accrete gas during the runaway phase (see Mordasini et al. 2011), 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. 1996). 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. 2013 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. 2011). 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. 2009b). 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. 2011), but is the ultimate goal of a number of ground based (e.g. Prima, Launhardt et al. 2008) and space-based missions (like the proposed Neat satellite, Malbet et al. 2011). 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. 2008). 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 2012).

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. 2010; Cassan et al. 2012) 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. 2013) or Wfirst (Goullioud et al. 2012), 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. (2012), Mordasini et al. (2012b), Lopez & Fortney (2013b) or Marcy et al. (2014) 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. 2012b). 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. 2000). 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. 2012b).

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 1981; Lodders 2003) 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. 2008). 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. (2011) 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. 2012, using the internal structure model of Guillot & Morel 1995). 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 1999), 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 2010) 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 2014; Mordasini et al. 2014).

Figure 15 compares the observed distribution of radii of planets inside of 0.27 AU as found by the Kepler satellite (Howard et al. 2012) 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 1994), while f _κ=0 means that we are dealing with a grain-free gas where only molecular opacities contribute (Freedman et al. 2008). 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. 2000). 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 2003; Movshovitz & Podolak 2008; Movshovitz et al. 2010).

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. 2011) according to the analysis of Howard et al. (2012). The black dotted lines are the preliminary analysis by the same authors of the updated Kepler data set of Batalha et al. (2012). 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. 2010). 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. 2014), 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. 2012; Spiegel & Burrows 2012).

The number of planets detected by this method has increased significantly in the last few years (see Marleau & Cumming 2014, 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. 2013). 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. 2009, 2010). The recent analysis of the near-IR spectral energy distribution by Bonnefoy et al. (2013) 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. 1997; Baraffe et al. 2003) 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. 2012) 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 2005; Alibert et al. 2011). 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. (1996) and that gap formation does not lead to a reduction of the gas accretion rate (Kley & Dirksen 2006).

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. (2007) 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. (2007) 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 2013). 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 2013 and Bodenheimer et al. 2013 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 2007) and the timescale of core accretion which might get speeded up if mainly small bodies are accreted (see, e.g. Dodson-Robinson et al. 2009; Ormel & Klahr 2010; Kratter & Murray-Clay 2011; Lambrechts & Johansen 2012). 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 2011; Deleuil et al. 2012).

Summary

We have reviewed global planet formation models that are used in population synthesis calculations. Such global models predict the properties of a planetary system based on the properties of a protoplanetary disc. They therefore establish a link between two classes of observable astrophysical objects which are on one hand the protoplanetary discs as the initial conditions for the planet formation process and on the other hand the extrasolar planets that are the final outcome of this process.

Global models have mainly been used in statistical studies. Such studies are a young approach that helps to improve the theory of planet formation by the comparison of theoretical results and statistical observational constraints provided by the entire population of extrasolar planets. With this approach, the global effects of many different physical mechanisms occurring during planet formation and evolution can be assessed, and (often strongly simplified) theoretical descriptions of these processes derived from specialized models can be tested against observations.

The first group of extrasolar planets that was known in sufficient numbers for statistical comparisons were giant planets detected by the radial velocity method. Population synthesis calculations therefore initially concentrated on studying the mass and semi-major axis distribution of this type of planets (Ida & Lin 2004a; Mordasini et al. 2009b). After this first phase of extrasolar planet detection, recent observational progress now also provides a first geophysical characterization of a growing number of planets outside of the Solar System (section ‘From detection to characterization’). This has led to additional observational constraints for global models like the planetary mass–radius relationship (Fig. 2) or the intrinsic luminosities of young giant planets (subsection ‘Direct imaging: luminosity at young ages’).

In order to yield synthetic populations that can be directly compared with all these different techniques, global models unite in one point the essential results of a significant number of specialized sub-models that describe one specific physical mechanism. The global formation and evolution model mainly discussed in this paper (Alibert et al. 2005a, 2013; Mordasini et al. 2012c), for example, consists of 11 sub-models (Fig. 3) addressing the following aspects: (1) the vertical structure of the protoplanetary disc, (2) the radial structure of the protoplanetary disc, (3) the disc of solids (planetesimals), (4) the core accretion rate of the protoplanet, (5) the planetary gas envelope, (6) the atmosphere of the (proto)planet, (7) the infall of planetesimals into the protoplanet's envelope, (8) the internal structure of the solid core, (9) the atmospheric escape (envelope evaporation) during evolution, (10) the orbital migration due to tidal interaction with the protoplanetary gas disc and finally (11) the gravitational interaction between the protoplanets. There are therefore three different classes of sub-models, namely those that describe the protoplanetary disc (1–3), those that describe one (proto) planet (4–9), and those that describe the interactions (10–11). In this paper, we have presented detailed descriptions of the physics included in these different sub-models and addressed their limitations and possible future improvements.

Owing to their nature as meta models, the global models discussed in this review depend directly on the results of many different specialized models, and therefore on the development of the entire field of planet formation theory. There are important uncertainties in this theory regarding even key aspects; therefore it is likely that the global models presented here will in future undergo significant modifications. These aspects include: (1) the formation of planetesimals and the resulting ‘initial’ distribution of solids in the disc; (2) the accretion of the solid core; (3) the opacity in the protoplanetary atmosphere and the associated gas accretion timescale; (4) the efficiency of orbital migration (which is still too rapid even with non-isothermal migration, see Fig. 12) and (5) the magnitude of gas accretion in the runaway phase. The latter two points could be addressed with 2D hydrodynamic simulations which can now be run over long timescales (Zhu et al. 2011) instead of 1D disc models. It is also possible that at some point it becomes necessary to abandon the simplifications that planetary systems form in isolation, because the gravitational interaction with other stars in young stellar clusters could be important (Malmberg et al. 2007).

The description of planetary evolution in the global models should in future include better atmospheric models and address the effects of heavy element settling, core erosion in giant planets and eventually the formation of secondary atmospheres based on the composition acquired during the planets’ formation for low-mass planets. Future global models of planet formation and evolution should therefore also include better descriptions of condensation and disc chemistry, so that the resulting composition of the planets can be predicted in a more detailed and self-consistent way (Thiabaud et al. 2014). On a longer timescale, this should make it possible to predict the habitability of a planet based on its formation.

Despite the current limitations, when used in planetary population synthesis calculations, global models can already now yield many testable predictions for the major observational techniques. This is important in a time where many surveys both from space and ground yield or will soon yield large amounts of additional data on both the global statistics and the detailed physical characteristics of extrasolar planets. Seeking for the theoretical models that best explain these combined data sets will be a promising approach towards a better understanding of planet formation and evolution.