© ESO, 2015

1. Introduction

Since the pioneering discovery of the first exoplanet orbiting a main-sequence star by Mayor & Queloz (1995), the number of discovered exoplanets has increased to over 1000 planets with a plethora of characteristics challenging planet formation models. The latest generation of core accretion models such as the planet formation model in Alibert et al. (see Alibert et al. 2005; Mordasini et al. 2009a,b; Fortier et al. 2013, Alibert et al. 2013, henceforth referred to as A05, M09, M09b, F13, and A13, respectively) or the semianalytical method of Ida and Lin (see Ida & Lin 2004, 2010; Ida et al. 2013) permits the study of planet formation not only for isolated planets, but for planets forming concurrently in the same system. This allows for the study of the effects of mutual perturbation between growing protoplanets. Thus, planetary characteristics, such as eccentricity, and architectural features, such as period ratios, which both are influenced by mutual interactions, can be compared to the distributions for observed exoplanets.

While the distributions of semimajor axes and planetary masses are generally well reproduced by planet formation models (see, e.g., M09 and M09b), the eccentricities of synthetic giant planets are statistically smaller than those for observed exoplanets (see, e.g., A13). This has been attributed to the damping effect of the gas disk during the formation phase. Studies by, e.g., Jurić & Tremaine (2008) showed that observed eccentricities can be reproduced by following the gravitational post-formation evolution of systems, but for their initial conditions the simulations relied on ad-hoc architectures for the initial system setup. It is not clear whether these types of architectures can actually result from the formation process. It is therefore uncertain whether dynamical instabilities suppressed while the gas disk is still present can account for the observed exoplanet eccentricities or if other effects such as e.g. perturbations from stellar fly-bys (see, e.g., Malmberg & Davies 2009) are necessary as a source of eccentric planets. Thommes et al. (2008), using a code coupling a 1D gas disk evolution code with SyMBA, found that for some disk parameters, giant planets with substantial eccentricities can be formed by planet-planet and planet-disk interactions. These planets are also obtained by A13, but their occurrence rates in the planet populations are relatively low and thus insufficient to match observed exoplanet distributions.

N-body simulations of the orbital evolution of (proto-)planets in gas disks (see, e.g., McNeil et al. 2005; Lee & Thommes 2009; Matsumura et al. 2010) found that a large number of planets ends up locked in low-order mean-motion resonances due to resonant migration. The large sample of exoplanets observed by the Kepler spacecraft, however, does not exhibit a similar enrichment in commensurable period ratios. Therefore, in formation models including both N-body dynamics and planet-disk interactions, a mechanism to break out of resonance is necessary to reduce the fraction of commensurable orbits.

In the present paper, we study the impact of post-formation dynamical evolution on planetary systems formed with the model of A13. We investigate whether the post-formation evolution alone can account for the observed statistical increase in eccentricities. Previous studies presented in M09b and A13 followed the formation of planets until the dispersal of the protoplanetary gas disk and typically lasted for up to T_disk ~ 10 Myr. We improve on these results by following the subsequent gravitational evolution of synthetic planetary systems for an additional 100 Myr after the dispersal of the protoplanetary disk, partially bridging the time gap between the formation era and the actual observation epoch. Between the time of disk dispersal and the long-term dynamically stable configurations of mature planetary systems, dynamical instabilities hitherto suppressed by the gas disk through eccentricity and inclination damping can arise, resulting in scatterings, collisions, and ejections of planets from their respective systems.

In Sect. 2, we briefly review our planet formation model and discuss the procedure we employed to compute the post-formation evolution and the planetary populations studied in Sect. 3. We discuss the effect on the evolved eccentricity distributions of a number of initial condition parameters as well as small stochastic perturbations in Sect. 4, while the post-formation evolution effects in general are presented in Sect. 5. Finally, we provide a brief discussion and conclusion.

2. Planet formation model

In the present work, we study the post-formation evolution decoupled from the planet formation process, numerically integrating the orbits of planets generated with a planet formation model. We use the planet formation model described in F13 and A13, based on earlier works (see A05 and M09), to generate populations of planetary systems in a Monte Carlo approach from initial parameters calibrated to observations (given in A13). The formation model self-consistently simulates the evolution of the gas and planetesimal disks, computes the solid and gas accretion together with the internal structure of planets, and accounts for both the gravitational interaction between planets and the interaction of planets with the surrounding gas disk.

The disk structure and evolution is modeled using a 1+1D-model computing the radial and vertical structure of the disk. The time evolution of the disk is obtained from the diffusion equation for the gas disk surface density. A detailed description of the disk structure and evolution model is given in Papaloizou & Terquem (1999) and A05. The interactions between growing protoplanets and the gas disk are computed using the prescriptions of Fogg & Nelson (2007) for Type I migration in the case of an isothermal disk and by setting the eccentricity damping proportional to the semimajor axis damping for planets in Type II migration. For nonisothermal disks, we use the migration and damping prescriptions of Paardekooper et al. (2010, 2011) and Dittkrist et al. (2014). A more detailed description of planet-disk interactions is given in A13.

2.1. Evolution procedure

To follow the long-term evolution of the systems generated with our population synthesis code, we use the symplectic integrator SyMBA (Duncan et al. 1998). The post-formation evolution of multiplanetary population syntheses is then obtained in the following way:

• 1.

  A population of planetary systems is generated using the planet formation model from A13. In this model, planets are followed dynamically until the gas disk mass decreases to M_disk ≤ 10^-5 M_⊙.

• 2.

  The system is then evolved over 100 Myr using the symplectic integrator SyMBA.

The chosen timescale for evolution of 100 Myr is a compromise between the time required to explicitly integrate a system in time and the number of systems we need to integrate to obtain a statistically meaningful result. The chosen integration time of 100 Myr is sufficiently large to allow for dynamical effects (cf. Jurić & Tremaine 2008, where they find a rapid drop in dynamical activity for T> 100 Myr in integrations running up to ~1 Gyr) acting on timescales longer than the disk lifetime such as, e.g., for secular resonances, orbit crossing, and collisions, where the timescales are typically a few Myr.

To reduce computational load and to focus on potentially observable systems, we set a minimum mass of the most massive planet (m_start) below which the system is not integrated, i.e., if no planet with m_p ≥ m_start is present in the system, the system is omitted from long-term integration. For the simulations studied in this paper, we set m_start to 1 M_⊕. This choice is also motivated by results from A13, where a convergence in the number of planets with a given mass is shown for planets with M_p ≥ 5 M_⊕, as a function of the initial number of embryos and initial conditions. For planets with M_p ≤ 5 M_⊕, significant effects from initial conditions on their occurrence rates and properties cannot be excluded. These uncertainties increase with further decreasing mass.

Collisions during the post-formation evolution are treated as perfect mergers, occurring in SyMBA if the separation d_ij between two objects is less than the sum of their physical radii. Note that we do not distinguish between core and envelope radius in our treatment of collisions.

3. Population synthesis

To test the effects of different parameters in the formation model, we generated several populations of planetary systems with the planet formation code (see Sect. 2), using the model described in A13.

A summary of the populations we studied is listed in Table 1. Below, each population is briefly described. We use a criterion based on the radial velocity semiamplitude K to select the planets in each population which we compare to a similarly selected sample of observed exoplanets, with K ≳ 10 m s^-1. We do not, however, apply a cutoff based on orbital periods as the fraction of planets with orbital periods larger than ~20 yr¹ is less than ~3% except for the population ND (see below), where the fraction is of ~13%.

3.1. The reference population

The reference population against which we compare the results of changes to the initial conditions of the planet formation process (populations ND, RD, E, see below) and stochastic perturbations (populations RA and RI, see below) was generated using 10 initial embryo seeds and disk masses, disk lifetimes, and dust-to-gas ratios drawn randomly from the parameter distributions described in A13. The disk masses vary from M_d = 0.001 M_⊙ to M_d = 0.3 M_⊙, with a dust-to-gas ratio scaled to the metallicity of the star. The gas disk masses are derived from Andrews et al. (2010), where a dust-to-gas ratio of f_d/g = 0.01 is assumed. The solid disk mass is computed by scaling the mass of the solid disk according to the randomly drawn stellar metallicity from a list of ~1000 CORALIE targets (see also A13) with a new f_d/g. In our simulations, most disks are low mass, with ~62% of all disks having M_d ≤ 0.03 M_⊙. Disk lifetimes T_disk are chosen from an exponentially decaying cumulative distribution function with characteristic time T = 2.5 Myr (see Mamajek 2009). The photoevaporation rate is then adjusted to yield M_d = 10^-5 M_⊙ at T = T_disk. Finally, embryos are seeded throughout the disk. All embryos have a starting mass of M = 10^-2 M_⊕ and an initial semimajor axis randomly drawn from a log-uniform distribution between 0.1 AU and 20 AU, with an added constraint that embryos are located at least 10 times their mutual Hill radius from each other. After all systems in the population have reached M_disk ≤ 10^-5 M_⊙, the systems containing no planet with M_p ≥ m_start were filtered out, while the remaining systems were integrated for an additional 100 Myr to follow the post-formation evolution under the mutual gravitational interactions of the planets alone (i.e., any residual gas or planetesimal disk and tidal interactions were not included). This resulted in a population of 2834 planets in 344 systems being evolved beyond the formation phase, with 2519 planets remaining after 100 Myr.

3.2. Modified damping timescales – populations ND and RD

As there are some uncertainties about the timescales of inclination and eccentricity damping of planets due to the interactions with protoplanetary disks (see, e.g., Bitsch & Kley 2010), we performed population syntheses with the same initial conditions as in the reference case, but modifying the damping timescales. The resulting populations without any eccentricity or inclination damping (population ND) and with the damping timescales increased by a factor of 10 with respect to the damping timescale prescription given in A13 (population RD) resulted in a population of 1853 planets in 264 systems, with 1492 planets remaining after 100 Myr (population ND) and in a population of 2480 planets in 343 systems, of which 2155 planets remained after 100 Myr (population RD).

3.3. Different initial number of embryo seeds

The number of embryos is a free parameter of the model, and Alibert et al. (2013) have already discussed the effects of changing the number of initial embryo seeds on the distributions of planetary masses and semimajor axes as well as on period ratios of neighboring planet pairs for the formation phase until the disk has dispersed. However, as the number of embryos increases, we also expect an increase in dynamical interactions. We therefore also computed the post-formation evolution of a population with the same initial conditions as the reference case, but with the number of initial embryo seeds increased to 20 (population E20), as well as a population with the same initial conditions as the reference population, but with only five initial embryo seeds (population E5). To maintain a similar total number of planets, a lower number of systems was simulated, resulting in a population of 3241 planets in 228 systems, of which 2238 planets remain after 100 Myr, whereas a slightly larger number of systems was simulated for the population E5, resulting in 2171 planets in 483 systems, of which 2129 planets remain after 100 Myr.

3.4. Additional stochastic perturbations

In the Bernese planet formation model, the migration prescriptions assume that the disk is azimuthally symmetric and homogeneous. Moreover, changes due to, e.g., turbulence (see, e.g., Johnson et al. 2006) which can affect the migration of low- to intermediate-mass planets are not included in the model. To simulate the effect of stochastic perturbations on the planets, we computed the post-formation evolution of the reference population, but added a random perturbation within ±1% to the semimajor axis of each planet (population RA) before computing the post-formation evolution. This is not intended as a self-consistent scheme as the perturbation bears no effect on the formation process, but rather as cursory investigation into the effect of semimajor axis perturbation on the subsequent evolution of planetary systems. However, we also computed a second population with stochastic perturbations (the population RI), this time in a self-consistent way, by adding a perturbation to planetary inclinations during the formation phase. The inclinations are perturbed by drawing a random inclination from a uniform distribution in [0°,1°] for planets with i ≤ 1°.

4. Results: effect of initial conditions on the eccentricity distribution

We used the post-formation evolution of a sample of 344 systems with ten initial embryos per system as the baseline against which we compared the effect of altering different aspects of the planet formation model, namely the number of embryos and the efficiency of the eccentricity (and inclination) damping due to interactions with the surrounding gas disk.

4.1. The initial number of embryos

In A13, the impact of the initial number of embryo seeds on the distributions of mass and semimajor axis was studied. The distributions for 5, 10 and 20 initial embryo seeds show a converging trend both with respect to the distribution of planetary masses (for planets with M_p ≳ 10 M_⊕) and the distribution of semimajor axes. In the present work, we compared the distribution of eccentricities between the different populations, both after the initial formation phase which typically lasts a few Myr until the gas disk is dissipated, and after a subsequent gravitational evolution lasting for 100 Myr after the time of disk dissipation. We find that unlike for the mass and semimajor axes, the eccentricities and their evolution differ between the populations with 5 (Fig. 1), 10 (Fig. 2) and 20 (Fig. 3) initial embryo seeds. The median eccentricity at the time of disk dissipation increases with increasing number of embryos from e = 0.019 (E5) to e = 0.030 (R) to e = 0.039 (E20) for planets with K ≥ 10 m s^-1. The full populations including the lower-mass planets follows a similar trend, albeit at much smaller eccentricities, with e = 4.1 × 10^-4, e = 5.8 × 10^-4, and e = 8 × 10^-4. Post-formation evolution affects the different populations in different ways. While all three populations show an increase in eccentricity for the full populations (median eccentricities of e = 4.6 × 10^-4, e = 1.2 × 10^-3, and e = 5.8 × 10^-3, respectively), the populations of planets with K ≥ 10 m s^-1 show a slight decrease in median eccentricity for the population with 5 (e = 0.017) and 10 embryos (e = 0.029), whereas the population with 20 embryos shows a slight increase in median eccentricity (e = 0.042), see Figs. 4 and 5.

We note that in all three populations the median eccentricities (and maximum eccentricities) are larger for the populations restricted to planets with K ≳ 10 m s^-1. On one hand, the eccentricity damping in our model is proportional to the damping in the semimajor axis and thus most efficient for planets with m_p ≲ 1−10 M_⊕, the mass around which the planets enter runaway gas accretion and ultimately transition into the slower Type II migration. Consequently, the more massive planets with K ≳ 10 m s^-1 migrating in Type II are less strongly damped than the lower-mass planets. On the other hand, the eccentricities of our planet populations are due to dynamical interactions, i.e., resonant excitations and secular perturbations as well as close encounters. Previous numerical and analytical studies (see e.g. Ida & Lin 2010; Ida et al. 2013, and references therein) have found that planetary eccentricities from two-body encounters correspond to a Rayleigh distribution multiplied with the mass-weighted escape eccentricity of the encounter, the latter being defined as ${\mathit{e}}_{\mathrm{esc}\mathit{,ij}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{esc}\mathit{,}\mathrm{ij}}}{{\mathit{v}}_{\mathit{K}}}\mathrm{=}\sqrt{\frac{\left({\mathit{m}}_{\mathit{i}}\mathrm{+}{\mathit{m}}_{\mathit{j}}\right)\mathit{a}}{{\mathit{M}}_{\mathit{\star }}\left({\mathit{R}}_{\mathit{i}}\mathrm{+}{\mathit{R}}_{\mathit{j}}\right)}}$, where $\mathit{a}\mathrm{=}\sqrt{{\mathit{a}}_{\mathit{i}}{\mathit{a}}_{\mathit{j}}}$ is the geometric mean semimajor axis of the planets i and j, m_i, m_j, and M_⋆ denote the mass of both planets and the star, and where R_i and R_j are the planetary radii. The resulting maximum post-encounter eccentricity is then given by ${\mathit{e}}_{\max \mathit{,j}}\mathrm{=}\frac{{\mathit{m}}_{\mathit{i}}}{{\mathit{m}}_{\mathit{i}}\mathrm{+}{\mathit{m}}_{\mathit{j}}}\mathrm{ℛ}{\mathit{e}}_{\mathrm{esc}\mathit{,ij}}$, where ℛ denotes the random variate from the Rayleigh distribution. For an encounter between low-mass planets, the escape eccentricity typically remains small, such that encounters between low-mass planets only lead to small increases in eccentricity. In contrast, encounters between a low-mass planet (j) and a massive planet (i) can rapidly result in e_max,j> 1 and thus in the ejection of the low-mass planet without substantially increasing the eccentricity of the massive planet. Encounters between two massive planets, while still potentially resulting in ejections, increase the eccentricities of both planets and thus produce a population of massive planets with some eccentricity. This bias against eccentric, low-mass planets is further compounded in our simulations by a large number of low-mass planets in systems without any massive companion, such that low-mass planets predominately interact with other low-mass planets. The correlation between planetary masses and eccentricities, however, is also found in observed exoplanets (see, e.g., Ida et al. 2013).

4.2. Eccentricity damping

A13 presented results for population syntheses with modified eccentricity damping timescales and found that including no eccentricity and inclination damping (ND) or increased damping timescales (RD) naturally led to populations with significantly larger eccentricities. We computed the post-formation evolution of the same populations to study the effect of reduced eccentricity damping during formation on the final architecture of planetary systems (Figs. 6 and 7). We found that in both cases with modified damping, the median eccentricities of the populations of planets with K ≳ 10 m s^-1 decrease over time, from e = 0.148 to e = 0.129 (ND) and from e = 0.102 to e = 0.095. For the full populations, the median eccentricity decreases in the case of no damping, dropping to e = 0.066 from e = 0.122 at the time of disk dispersal, whereas the median eccentricity increases from e = 1.0 × 10^-3 at the time of disk dispersal to e = 2.2 × 10^-3 after 100 Myr (Fig. 8).

Contrary to our original expectation that eccentricities would grow throughout the post-formation evolution, we find that in both populations, the median eccentricity for planets with K ≳ 10 m s^-1 decreases over 100 Myr. This is explained by the fact that starting with larger post-formation eccentricities and inclinations, the systems are more dynamically active, which leads to more collisions and ejections. This typically only affects planets with eccentricities close to or larger than unity, as we consider a planet to be ejected from the system if its heliocentric distance r_p ≥ 1000 AU and only very few planets are scattered onto wide orbits with moderate eccentricities. This eventually leads to a decrease in eccentricity not fully compensated by eventual excitation of the remaining body in a two-body encounter (see Sect. 4.1). In addition to the loss of eccentric planets through ejections, collisions also reduce the eccentricities. This can be more readily understood by evaluating the angular momentum deficit of a planet, introduced by Laskar (1997) and defined as ${\mathrm{𝒞}}_{\mathit{i}}\mathrm{=}{\mathit{m}}_{\mathit{i}}\sqrt{\mathit{G}{\mathit{M}}_{\mathit{\star }}{\mathit{a}}_{\mathit{i}}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathit{i}}^{\mathrm{2}}}\cos {\mathit{i}}_{\mathit{i}}\right)\mathit{,}$(1)with semimajor axis a_i, eccentricity e_i and inclination i_i, and with m_i and M_⋆ denoting the planet and host star mass, respectively. The AMD characterizes the degree of nonlinearity in a given orbit. For purely secular systems, the total angular momentum deficit, given by the sum over all individual , is conserved. Laskar (2000) notes that in a fully elastic collision, the total angular momentum deficit of the colliding pair decreases. Consequently, the eccentricity of the resulting planet tends to be smaller than the eccentricities of the parent bodies. As the planets in our simulations are generally on well-separated orbits, only eccentric planets undergo collisions, thus again removing the more eccentric planets from the population.

The relatively large fraction of collisions is due to large initial populations of planets with Safronov numbers Θ < 1, for which the preferred outcome of close encounters are collisions rather than ejections. These large fractions (>94% for both ND and RD) are due to a large number of both close-in and low-mass planets for which the local orbital escape velocity is larger than the surface escape velocity, such that a scattering by the planet is not likely to eject the scattered planet from the system.

4.3. Stochastic perturbations

In previous studies (see A13), we found that the fraction of planet pairs with period ratios close to the nominal locations of low-order mean motion resonances are significantly larger than for observed exoplanets. As resonances can stabilize otherwise unstable systems and thus have an important effect on the eccentricitiy evolution of planetary systems, stochastic perturbations to orbital elements from, e.g., density fluctuations in the disk can perturb planets out of resonance and thereby increase the likelihood of strong dynamical interaction (scatterings and collisions). To evaluate the effect of small stochastic perturbations, we computed the post-formation evolution for two modified versions of the reference population (R). In the first case, we added a post facto perturbation δa ∈ [−0.01,0.01]a to the semimajor axis of each planet at the end of their formation, but before starting the post-formation evolution (population RA). In the second case, we ran the formation code and added a random, uniformly distributed inclination i ∈ [0°,1°] for each planet reaching an inclination of i ≤ 1° during the formation phase while the gas disk was present (population RI). Contrary to the semimajor axis perturbation, the population with the perturbed inclination is self-consistent, i.e., the perturbation effects were applied to the full formation process rather than just perturbing the final orbital elements.

The population RA has the same median eccentricities as the reference population for planets with K ≳ 10 m s^-1, with medians e = 0.030 at the time of disk dispersal and e = 0.029 after 100 Myr. For the full population, the median eccentricity after 100 Myr is slightly higher than in the reference population, with e ≃ 2 × 10^-3 versus e = 1.2 × 10^-3. Similarly, the eccentricities for the population with perturbed inclinations (RI) are comparable to the eccentricites in the reference population. For planets with K ≳ 10 m s^-1, the median eccentricity after 100 Myr is even slightly smaller than in the reference population, with e = 0.028 versus e = 0.029. We conclude that the stochastic perturbations we applied to planetary inclinations or semimajor axes have negligible influence on the final eccentricities of planet populations.

4.4. Comparison with observations

As already stated in A13, we made no attempt to have the synthetic systems reproduce the observed systems. Hence, we did not expect the evolved synthetic population’s characteristics to match the observed one either. Nevertheless, we compared the eccentricity distributions of the simulated planet populations to the eccentricity distribution of observed exoplanets² (see Fig. 9), limiting all samples to planets with K ≥ 10 m s^-1. We find that the eccentricities in the populations with full eccentricity and inclination damping are significantly smaller than for observed exoplanets. The population where the eccentricity damping timescale was increased by a factor of 10 compared to the nominal value used provides the best agreement³ with observed exoplanet eccentricities. Nevertheless, the median eccentricity of e = 0.095 after 100 Myr is still smaller than the median eccentricity of observed exoplanets (e = 0.15). A possible source of this difference is the inclusion of planets around stars with no further detected planetary mass companions, but possible distant stellar companions, in the sample of exoplanets. Excluding all exoplanets in systems with multiple stellar components (list taken from Roell et al. 2012), however, still yields a median eccentricity of e = 0.14 and does not appreciably change the eccentricity distribution except for the removal of a few highly eccentric planets such as, e.g., HD 20782 b (see Jones et al. 2006). A possible source of additional eccentricity not yet included in our model are perturbations by stellar fly-bys, which can increase planetary eccentricities both by direct perturbations and by destabilizing planetary systems, thus increasing the incidence of strong dynamical interactions, see, e.g., Malmberg et al. (2011).

5. Results: post-formation evolution effects

In addition to the change in planetary eccentricities, the post-formation also potentially affects the distributions of other quantities. We compared the mass, semimajor axis, and period ratio distributions as well as the chemical composition of planets at the time of disk dispersal and after 100 Myr.

The fraction of planets removed by the post-formation evolution strongly depends on both the number of embryos and the efficiency of eccentricity and inclination damping. In the reference population, ~11% of planets present at the end of the formation phase are removed by collisions and ejections. For the population E5, this amount drops to <2% while for E20, it increases to ~30% mostly because of a significant increase in the number of collisions compared to the reference population. Reducing the efficiency of eccentricity and inclination damping by a factor of 10 (population RD) results in a slight increase in the fraction of removed planets (~13%) through an increase in the fraction of ejected planets (~3% for RD versus <1% for R), while completely removing the eccentricity and inclination damping results in ~23% of all planets removed through collisions and ejections, with marked increase in both types of removal. The situation is similar for the fraction of systems experiencing collisions and/or ejections in each population (see Table 2), with a sharp increase for increasing number of embryos and an increase for decreasing damping. Any change in the distribution of mass (Figs. 10, 11, 12), semimajor axis, or period ratios is therefore likely most pronounced in the populations E20 and ND. The addition of stochastic perturbations results in similar fractions of lost planets and active systems in the case of inclination perturbations (RI), but in substantially different fractions for perturbations of the semimajor axes (RA). For these systems, the fraction of lost planets increases by almost a third to f_lost = 16%, while the fraction of active systems, i.e., systems with either collisions or ejections, increases from f_sys,active = 36% in the reference population to f_sys,active = 57% in the perturbed population. The large increase in activity despite the rather small perturbation suggests that many systems in our reference population are close to instability. Most likely, such instabilities are stabilized by the presence of mean-motion resonances, which are then broken by the semimajor axis perturbations. The lack of a significant increase in eccentricity for the perturbed population, however, implies that the increased activity predominately affects the lower-mass planets.

5.1. The planetary mass function

The effect of the post-formation on mass distributions of the populations studied is negligible for planets with K ≥ 10 m s^-1 (see Fig. 13), with both the distributions at the time of disk dispersal and after 100 Myr (Fig. 15) in solid agreement with the mass distribution of observed exoplanets when taking the observed exoplanet masses to be their true masses. For exoplanets detected by radial velocity measurements, however, planetary masses are in general only known up to a factor sini, i being the angle between the observer and the orbital plane of the planet. Including a random line-of-sight inclination drawn from i ~ sini (see. e.g., Ho & Turner 2011) to compute the msini, we find that our planets have a smaller median msini than observed exoplanets (see Fig. 14). Moreover, using a two-sample Kolmogorov-Smirnov test, all populations studied show a better agreement with the observed mass distribution for the true population rather than for the sini corrected distribution (cf. Table 3). Indeed, the null hypothesis that the msini are drawn from the same distribution as observed exoplanet msini cannot be rejected only for the population E5. Comparing the true masses to observed msini, the null hypothesis can only be rejected for the population ND, whereas the remaining populations are in solid agreement with the observed distribution. As m ≥ msini, the (massive) planets in our populations appear to be somewhat less massive overall than observed exoplanets. In addition to the planets remaining after 100 Myr, we also tracked the masses of planets ejected both during the formation phase and the subsequent post-formation evolution. Microlensing surveys claim a large number of free-floating planets with m ~ M_Jup, with, e.g., Sumi et al. (2011) finding that free-floating planets in the Jovian mass range are about twice as common as main-sequence stars. Even though the planetary mass functions in our simulations are comparable to the planetary mass function for observed exoplanets, our simulations generally produce only few ejections (see above) and an even lower fraction of ejected planets with masses m ≳ 100 M_⊕. Of all planets with m ≥ 100 M_⊕, only ≲7% are ejected from simulations with the full eccentricity and inclination damping (populations R, E5, E20 and RA), and ~19% and ~43, respectively, for the population with reduced (RD) and without (ND) eccentricity damping. With the fraction of planets with m ≥ 100 M_⊕ between 5% and 10% in our simulations, ejections during the formation phase and the post-formation evolution alone cannot explain a large population of Jovian-mass free-floating planets.

5.2. Semimajor axis distribution

Similar to the mass distributions, the distributions of the semimajor axes remain almost unchanged by the post-formation evolution for planets with K ≳ 10 m s^-1. For the full populations, the changes increase with increasing number of embryos and with decreasing eccentricity damping, with the population with 20 initial embryo seeds exhibiting the largest changes to the semimajor axis distribution (Fig. 16). With the number of scattered planets being smaller than the number of planets undergoing collisions, the changes in the distributions are an effect of removing a part of the population rather than of changing the semimajor axes of the planets. This is also evidenced by the lack of change in the distribution of period ratios (see Sect. 5.4).

5.3. The chemical composition of planets

For the same disk parameters and evolutions as in the planet formation code, Marboeuf et al. (2014a,b) and Thiabaud et al. (2014) respectively computed the distribution of volatile and refractory elements from chemical condensation sequences, thus allowing for the computation of the chemical compositions of planets in the formation model. Using the initial compositions computed by Marboeuf et al. (2014a,b) and Thiabaud et al. (2014) in the case of a nonirradiated disk⁴, we computed the chemical compositions after 100 Myr for the reference population and the population with an increased number of embryos by following the collisional evolution of the planets.

Comparing the initial and final distribution of mass fractions of water for the planetary cores f_H2O,core (see Fig. 17), we find that post-formation evolution has a minimal impact on f_H2O,core as the population median remains at f_H2O,core = 0.284. The spread in f_H2O,core decreases somewhat as the lower quartile increases from f_H2O,core = 0.241 to f_H2O,core = 0.255, indicating that in our simulation, post-formation evolution reduced the amount of “dry” planets. Comparing the evolution of the distribution of f_H2O,core at different mass and semimajor axis ranges (Figs. 18 and 19, respectively), we find that planets with m_p ≤ 1 M_⊕ and planets with a_p ≤ 0.5 AU exhibit significant change in f_H2O,core, whereas more massive planets and planets on wider orbits show almost no change in f_H2O,core. However, for close-in, low-mass planets, the fraction of “dry” planets shows a slight increase after 100 Myr.

The lack of change in f_H2O,core for massive planets is primarily due to the very low frequency of collisions (and, to a lesser degree, ejections) involving massive planets. Note that we do not take into account any potential loss of water due to collisions either during the formation phase or during the post-formation evolution.

5.4. Period ratios

The distributions of period ratios of planets with K ≥ 10 m s^-1 remain virtually unchanged by the post-formation evolution for all planet populations we simulated. For the full populations, the post-formation evolution slightly increases the fraction of pairs with large period ratios (T_out/T_in ≳ 10). This effect appears to be an increasing function of the number of initial embryo seeds: for 5 embryos, the effect of the post-formation is negligible (Fig. 21), for 10 embryos, the effect is small (Fig. 20), and for 20 embryos, the effect becomes somewhat important (Fig. 22). Conversely, the fraction of planet pairs close to low-order mean-motion resonances decreases with increasing number of embryos, both for the full populations and for the populations restricted to planets with K ≥ 10 m s^-1. Reducing the efficiency of eccentricity and inclination damping by increasing the damping timescale by a factor of 10 does not significantly change the period ratio distributions compared to the reference case (see Fig. 23). Completely removing the damping, however, produces a starkly different distribution with no discernable population of pairs in low-order mean motion resonances (Fig. 24).

The population with a random perturbation to the semimajor axes (population RA), while showing an overal similar distribution of period ratios for planets with K ≥ 10 m s^-1, has a lower fraction of planets within ± 1% of the nominal locations of low-order mean motion resonances after 100 Myr compared to the reference population. The ratio of the fraction of systems within one percent of the 2:1 mean motion resonance in the reference population at the end and at the beginning of the post-formation evolution, for instance, is f_2:1,fin/f_2:1,ini = 0.17/0.18 ≈ 0.94, whereas the same ratio for the population RA is f_2:1,fin/f_2:1,ini = 0.10/0.14 ≈ 0.71 comparable to the values for the population with more initial embryo seeds (E20), which has f_2:1 ≃ 0.15 initially and f_2:1 ≃ 0.11 after the post-formation evolution. Compared to f_2:1 ≃ 0.05 for observed exoplanet pairs (Figs. 25 and 26), this is still considerably larger and points to our model either lacking strong dynamical interactions that could force the probably resonant planet pairs away from resonance, or having an overly efficient migration producing (too) many pairs in resonance. Both possibilities would also affect the distribution of eccentricities discussed in Sect. 4.4.

6. Discussion

1. Post-formation evolution effects: eccentricity In the simulations we present, evolving various populations obtained with the Bernese planet formation code for an additional 100 Myr after the time of disk dispersal affects the planetary eccentricities. Populations with low initial eccentricities show an overall increase in median eccentricity resulting from increased eccentricities for low-mass planets, whereas massive planets exhibit minimal increase or even decrease in eccentricity. Contrary to the expectation that post-formation evolution leads to increased eccentricities, initially eccentric populations experience an “eccentricity evaporation” effect where the more eccentric planets are rapidly lost through collisions and ejections, resulting in smaller median eccentricities after the post-formation evolution, but still higher than in populations with low initial eccentricity. The eccentricity distributions of all populations we followed for 100 Myr remain a poor match to observations. To further increase the final eccentricities, additional sources of excitation such as stellar fly-bys (see below) would be required in our model.

2. Post-formation evolution effects: mass, semimajor axis and composition While the eccentricity distributions exhibit some change over the course of the post-formation evolution, the distributions of semimajor axes and planetary masses as well as the distribution of period ratios are largely unchanged by post-formation effects. This might be partly because of the relatively low dynamical activity in the populations we studied, as there is a clear correlation between the dynamical activity (gauged by proxy of using the number of ejections and collisions and the fraction of systems exhibiting either) and the amount of change between distributions at the time of disk dispersal and after 100 Myr for both semimajor axis and mass distributions. This low amount of change in distributions, even in the most dynamically active population for planets with K ≳ 10 m s^-1, implies that the observed distributions of exoplanet semimajor axes and masses are a close approximation to their state at the end of their respective disks’ lifetimes, barring external perturbations. Similarly, the chemical compositions of massive planets and planets on wide orbits exhibit minimal changes from the post-formation evolution effect as collisions (and the resulting changes in chemical composition) primarily occur between planets with m_p ≤ M_⊕. Comparisons of masses and semimajor axes between planet formation results and observed exoplanets, such as M09b, therefore are a warranted approach.

3. Number of initial embryos While the initial number of embryo seeds has a minimal influence on the distribution of planetary masses and only small influence on the distribution of semimajor axes after 100 Myr for planets with K ≥ 10 m s^-1 in our simulations, the impact on planetary eccentricities and period ratios is very important. For period ratios, increasing the initial number of embryos decreases the fraction of planet pairs with period ratios close to the nominal location of low-order mean-motion resonances, but this fraction remains relatively unchanged by the post-formation evolution.

   Fig. 26 Comparison of the period ratios after 100 Myr for different planet populations. The black solid line denotes the distribution for Kepler planet candidates (data taken from www.exoplanets.org), the red dashed line for the reference population (R), the blue dot-dashed line for the population with 20 intial embryo seeds (E20), and the cyan dotted line for the population with five initial embryo seeds (E5).

4. Integration time The 100 Myr integration time is sufficiently long for dynamical instabilities to arise and resolve themselves through collisions and scatterings, and studies of dynamically active systems by, e.g., Jurić & Tremaine (2008) find that dynamical activity strongly decreases after 100 Myr. In Fig. 27, we show the distributions of the timing of dynamical events during the post-formation evolution in the populations we studied. The median event times are τ = 1.28 Myr for the reference population, τ = 2.24 Myr for RD, τ = 3.69 Myr for ND, τ = 4.89 Myr for E5, τ = 0.41 Myr for E20 and τ = 0.39 Myr for RA. In all populations except for E5, where the very low number of dynamical events with n_event = 36 makes statistics somewhat unreliable, ≳90% of all events happen within the first 50 Myr of post-formation integration time, leading us to conclude that for most systems, 100 Myr represent a sufficiently long time to transition from the formation phase to long-term stable configurations.

   Fig. 27 Timing of ejections and collisions in the populations R (upper left), RA (upper right), RD (middle left), ND (middle right), E5 (lower left) and E20 (lower right).

5. Free-floating planets Observations of microlensing events have resulted in claims that, in addition to exoplanets bound to their host stars, a large population of free-floating or rogue planets exist, see, e.g., Sumi et al. (2011). While this claim is disputed by, e.g., Veras & Raymond (2012), a sizeable population of planets ejected from their birth systems might still exist. In our populations, planets are ejected both during the formation phase and during the post-formation evolution. Unlike the large fraction of Jovian mass planets claimed by microlensing observations, the ejected planets in our populations have median masses m ≤ 2 M_⊕, and the fraction of ejected planets with m ≳ M_Jup versus the total number of planets with m ≳ M_Jup in a given population are ≲5% for the populations with full eccentricity and inclination damping, ~19% for the population with reduced damping and ~33% for the population without damping. We therefore do not produce the large number of free-floating Jovian planets expected from microlensing surveys. If a large number of free-floating planets would be confirmed, this would imply either a different formation mechanism for such planets, e.g. through gravitational instabilities in the outer disk (see Boss 2011), or a significantly larger amount of dynamical interactions during and after the formation phase.

6. Comparison with previous studies The low number of ejected giant planets is at odds with assumptions about the initial conditions in studies reproducing the observed exoplanet eccentricity distributions from giant planet scattering, such as Jurić & Tremaine (2008), Chatterjee et al. (2008) or Malmberg & Davies (2009; referred to as JT08, C08, and MD09, respectively). In all three studies, the initial planetary systems are more massive, with multiple planets with masses between 0.1 M_Jup and 10 M_Jup, and more compactly spaced, with median separations in mutual Hill radii $\mathit{b}\mathrm{=}\frac{{\mathit{a}}_{\mathrm{O}}\mathrm{-}{\mathit{a}}_{\mathrm{I}}}{{\mathit{R}}_{\mathrm{H}\mathit{,}\mathrm{IO}}}$, where ${\mathit{R}}_{\mathrm{H}\mathit{,}\mathrm{IO}}\mathrm{=}\frac{{\mathit{a}}_{\mathrm{I}}\mathrm{+}{\mathit{a}}_{\mathrm{O}}}{\mathrm{2}}{}^{\mathrm{3}}\sqrt{\frac{{\mathit{m}}_{\mathrm{I}}\mathrm{+}{\mathit{m}}_{\mathrm{O}}}{\mathrm{3}{\mathit{M}}_{\mathit{\star }}}}$, ranging from b = 2 in MD09 to b = 4.4 in C08 and active systems in JT08. In contrast, the planets with K ≳ 10 m s^-1 in our reference population have a (slightly) larger median separation of b = 5.2. The timescale for orbit crossing for a given planet pair strongly depends on the orbital separation (see, e.g., Zhou et al. 2007), such that even the moderate increase in separation between planets results in a strongly reduced amount of dynamical activity in our simulations. Moreover, the assumptions about planetary masses, especially on the frequent occurrence of multiple Jovian-mass planets in a single system, are not reproduced by any population in our simulations. Indeed, the fraction of systems with more than one planet with m_p ≥ 30 M_⊕ ≃ 0.1 M_Jup in our reference population is ~22%, the fraction of systems with more than two such planets is only ~3%. If we set the mass limit to m_p ≥ 1 M_Jup, the respective fractions drop to 11% and 1%. Clearly, the assumptions about the initial conditions in C08, JT08, and MD09 are not reproduced by our planet formation model. As in Thommes et al. (2008), we do find an initial population of eccentric giant planets from the formation process. The populations RD and ND (with reduced and no eccentricity damping) in particular produce several eccentric giant planets and show the best initial agreement with observed exoplanet eccentricities. For both populations, however, the final eccentricities are reduced by the efficient removal of eccentric planets by means of ejections and collisions, resulting in poorer agreement with observed exoplanet distributions than for the initial distributions.

7. External perturbation as a source of dynamical activity With the exception of the population without eccentricity and inclination damping, our populations have typical eccentricities that are much smaller than those of observed exoplanets. Despite uncertainties about whether some of the observed highly eccentric exoplanets could be better explained by multiplanet systems with smaller eccentricities, see, e.g., Anglada-Escudé et al. (2010) or Wittenmyer et al. (2013), the post-formation evolution of our populations is insufficient to reproduce observed eccentricities, thus requiring additional sources of eccentricity. In our simulations, we consider the forming systems to be isolated, whereas star formation and, consequently, planet formation typically takes place in clusters. Studies by, e.g., Malmberg et al. (2007, 2011) have shown that in their reference cluster, a solar-mass star on average has four close encounters⁵ and ~75% have at least two close encounters. Roughly half of these close encounters take place during the first 10 Myr, which is comparable to the disk lifetimes (see Mamajek 2009), and ~90% happen within 100 Myr. These stellar fly-bys can affect the dynamical evolution of nascent planetary systems in two ways: They can directly scatter or eject planets in the case of very strong encounters, or by reduce the stability of the system by increasing eccentricities and inclinations, thus potentially forcing planets out of stabilizing mean-motion resonances. Malmberg et al. (2011) found that between ~40% and 60% of their systems underwent ejections or orbit-crossings after a stellar fly-by, implying that the contributions from stellar encounters to the final architecture of planetary systems could be significant. The inclusion of stellar fly-bys in the current planet formation model will be the topic of future studies.

7. Conclusion

Within the first 100 Myr of evolution after the dissipation of the protoplanetary disk, planetary systems can undergo dramatic changes in orbital architecture through planet-planet scattering, ejections, and collisions. These dynamical interactions are necessary to bring the eccentricities of planet populations in the Bernese planet formation model in line with observed exoplanet eccentricities. However, the effect of the post-formation evolution on the eccentricity distributions alone remains insufficient to explain observed eccentricities, requiring further investigation into the assumptions about initial conditions as well as the inclusion of additional sources of dynamical perturbation such as stellar fly-bys in our model.

While the eccentricities clearly show the effects of post-formation evolution, other planetary characteristics such as the mass and semimajor axis distributions remain largely the same after 100 Myr as at the time of disk dispersal, implying that the mass and semimajor axis of a planet are by and large determined by the formation process alone.

---

Acknowledgments

This work has been in part carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF. This work was supported by the European Research Council under grant 239605.

References

All Tables

All Figures