On the formation of compact planetary systems via concurrent core accretion and migration
We present the results of planet formation N-body simulations based on a comprehensive physical model that includes planetary mass growth through mutual embryo collisions and planetesimal/boulder accretion, viscous disc evolution, planetary migration and gas accretion onto planetary cores. The main aim of this study is to determine which set of model parameters leads to the formation of planetary systems that are similar to the compact low mass multi-planet systems that have been discovered by radial velocity surveys and the Kepler mission. We vary the initial disc mass, solids-to-gas ratio and the sizes of the boulders/planetesimals, and for a restricted volume of the parameter space we find that compact systems containing terrestrial planets, super-Earths and Neptune-like bodies arise as natural outcomes of the simulations. Disc models with low values of the solids-to-gas ratio can only form short-period super-Earths and Neptunes when small planetesimals/boulders provide the main source of accretion, since the mobility of these bodies is required to overcome the local isolation masses for growing embryos. The existence of short-period super-Earths around low metallicity stars provides strong evidence that small, mobile bodies (planetesimals, boulders or pebbles) played a central role in the formation of the observed planets.
keywords:planetary systems, planets and satellites: formation, planets-disc interactions, protoplanetary discs.
Both radial velocity (Mayor et al., 2011) and transit surveys have shown conclusively that systems of low mass planets are common around main sequence stars, with the Kepler mission in particular providing some striking examples of short period compact multi-planet systems (Lissauer et al., 2011; Fabrycky et al., 2014). The most recent release of Kepler data contains over 4700 planet candidates, and more than 700 multi-planet systems (Mullally et al., 2015). Approximately 3000 systems show just a single transiting planet candidate, with orbital periods in the range days.
Analysis of the systems properties provides useful insight for understanding how these planets formed and evolved. One noticeable feature of the multi-systems is the paucity of first order mean motion resonances. The period ratio distribution shows features in the vicinity of the 2:1 and 3:2 resonances, suggesting that they have been dynamically important in the past, but relatively few systems are actually in a strict mean motion resonance (Fabrycky et al., 2014). Examples of systems of small planets that are in or very close to resonance, including 3 body resonances or resonant chains, include Kepler 50 (6:5), Kepler 60 (5:4, 4:3) (Steffen et al., 2012), Kepler 221 (displays a 3 body resonance) (Fabrycky et al., 2014). In general the compact multi-planet systems appear to be composed of terrestrial planets, super-Earths and Neptune-like bodies. Mass estimates based on both radial velocity and transit timing variations suggest that there is a strong diversity in the mean densities of these objects, with some being rocky, some appearing to have a mixture of rock and water, and others being of very low density indicating the presence of significant fractions of H/He (Lissauer et al., 2011; Wu & Lithwick, 2013; Marcy et al., 2014; Jontof-Hutter et al., 2015). Kepler 36 provides an example where a pair of neighbouring planets orbiting close to the 7:6 resonance have dramatically different densities, characteristic of a rocky terrestrial inner body and an outer mini-Neptune (Carter et al., 2012). One of the most interesting facts to emerge from the data is the presence of low mass planetary systems around stars with a broad range of metallicities, including stars whose iron contents are factors of smaller than the solar abundance (Buchhave et al., 2014), a result that is supported by radial velocity discoveries of planets around metal-poor M dwarfs, such as Kapteyn’s star (Anglada-Escudé et al., 2014).
A number of ideas have been put forward to explain the formation and early evolution of the compact Kepler and radial velocity systems, which in cases such as Gliese 581 and HD 69830 appear to contain in excess of of solid material within a few tenths of an au (Lovis et al., 2006; Udry et al., 2007). This concentration of solids close to the star led to classical core accretion models combined with disc driven migration being developed using population synthesis codes (Alibert et al., 2006). More recent population synthesis calculations that also include prescriptions for planet-planet interactions have also been presented (Ida & Lin, 2010). N-body simulations, combined with either hydrodynamic simulations or analytic prescriptions for migration and eccentricity/inclination damping of planetary growth, have also been used to examine the origins of such systems (Cresswell & Nelson, 2006, 2008; Terquem & Papaloizou, 2007; McNeil & Nelson, 2009, 2010; Hellary & Nelson, 2012; Cossou et al., 2013; Coleman & Nelson, 2014; Hands et al., 2014). A common outcome of these N-body simulations is the formation of resonant convoys of planets in the presence of convergent migration, an outcome that is not reflected in the Kepler systems. Various ideas have been put forward to explain why the resonances may be unstable, including tidal eccentricity damping followed by separation of the resonance for short period systems (Terquem & Papaloizou, 2007), stochastic migration due to local turbulence (Adams et al., 2008; Rein & Papaloizou, 2009; Rein, 2012) - a process that is likely to only operate close to the star where the disc can be thermally ionised (Umebayashi & Nakano, 1988; Desch & Turner, 2015), resonance breaking due to overstable librations (Goldreich & Schlichting, 2014), orbital repulsion due to nonlinear spiral wave damping in planet coorbital regions (Baruteau & Papaloizou, 2013; Podlewska-Gaca et al., 2012).
The paucity of mean motion resonances in the Kepler data has led to suggestions that the compact systems formed in situ through giant impacts, akin to the final stages of accumulating the terrestrial planets Chambers & Wetherill (1998), after the concentration of small planetesimals in the inner disc followed by their growth into planetary embryos (Hansen & Murray, 2012). Although this model has some success in generating non resonant multiple planet systems with inclinations that are in good agreement with Kepler systems, there are difficulties in explaining how such large amounts of solids become concentrated in the inner disc, and the model fails to reproduce the numbers of single transiting planets detected by Kepler (Hansen & Murray, 2012). An alternative in situ model has been proposed by Chatterjee & Tan (2014) where pebbles/boulders concentrate and form a planet at the pressure maximum generated at the interface between the inner turbulent region of the disc and the dead zone, and exterior planets are spawned in succession by the disc being eroded outwards when the planets reach gap forming masses. While this model may be able to explain some systems, it is not clear that such a model can work for systems such as Kepler 444 and Kepler 186 where the planet masses are likely to be too small to form gaps, or for planetary systems in which the innermost planets orbit further from their stars than the fully active regions are expected to extend.
In this paper we present the results from a suite of N-body simulations using an updated version of the planet formation and protoplanetary disc model presented in Coleman & Nelson (2014), hereafter referred to as CN14. A basic assumption of the model is that the protoplanetary disc contains a population of planetary embryos distributed across a wide range of orbital radii (between 1 - 20 au), which grow through the accretion of boulders or planetesimals, and through mutual collisions, and can accrete gas from the nebula when they reach masses . We refer to this as a distributed core accretion model, in contrast to one where smaller numbers of embryos might form at specific disc locations such as pressure maxima through the trapping and accumulation of solids. We include the most up-to-date prescriptions for migration, self-consistent evolution of the viscous disc, and disc removal by a photoevaporative wind on multi-Myr time scales. The main updates on CN14 include placing the inner boundary of the computational domain close to the star so that we can simulate planets that migrate to regions with orbital periods down to 1 day, addition of an active turbulent region (mimicked as a simple increase in viscosity) where disc temperatures exceed 1000 K, and a magnetospheric cavity close to the star into which planets can migrate. The aim of this work is simply to examine whether or not such a comprehensive model of planet formation is able to produce planetary systems that are similar to those that have been observed, and if so under which set of conditions (disc mass, metallicity, planetesimal/boulder sizes) do these systems form. We emphasise that this is not population synthesis. No attempt is made to select initial conditions from a distribution of possibilities constrained by observations. We are not aiming to reproduce the frequency with which certain types of planetary systems arise, but instead to just examine which conditions allow Kepler-like compact systems to form subject to our model assumptions. We find that the simulations produce a broad range of outcomes that correlate strongly with the amount of solids present initially in the disc, and with the sizes of the boulders/planetesimals that provide the primary feedstock for planetary growth. Short-period compact multi-systems containing both resonant and non resonant planet pairs are one particular outcome of the runs, and these arise within a restricted range of the parameter space that we consider.
2 Physical model and numerical methods
The N-body simulations presented here were performed using the Mercury-6 symplectic integrator (Chambers, 1999), adapted to include the disc models and physical processes described below. We use an updated version of the physical model described in CN14. The main elements of this model are described below, and the implemented updates are outlined in the following subsections. The basic model consists of 52 protoplanets, orbiting within a swarm of thousands of boulders or planetesimals, embedded in a gaseous protoplanetary disc, all orbiting around a solar mass star. For each simulation we adopt a single size for the boulders or planetesimals. We define objects of radius m to be boulders and objects of radius m to be planetesimals. These various sized objects differ from each other and from protoplanets or planetary embryos because they experience gas drag forces that vary with the size.
2.1 Recap on the CN14 model
The basic model from CN14 is comprised of the following elements:
(i) The standard diffusion equation for a 1D viscous -disc model is solved (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974). Temperatures are calculated by balancing blackbody cooling against viscous heating and stellar irradiation. In the presence of a giant planet tidal torques can open a gap in the disc.
(ii) The final stages of disc removal occur through a photoevaporative wind. A standard photoevaporation model is used for most of the disc evolution (Dullemond et al., 2007), but if a large cavity forms in the presence of a gap forming planet, direct photoevaporation of the disc is switched on if the planet sits outside of the innermost radius from where the thermally-driven wind can be launched (Alexander & Armitage, 2009), as outlined in Section 5 of CN14.
(iii) Planetesimals and boulders orbiting in the disc experience aerodynamic drag.
(iv) We use the torque formulae from Paardekooper et al. (2010, 2011) to simulate type I migration due to Lindblad and corotation torques. Corotation torques arise from both entropy and vortensity gradients in the disc, and the possible saturation of these torques is included in the simulations. The influence of eccentricity and inclination on the migration torques, and on eccentricity and inclination damping are included (Fendyke & Nelson, 2014; Cresswell & Nelson, 2008).
(v) Type II migration of planets is included via the impulse approximation of Lin & Papaloizou (1986) if they reach the gap opening mass.
(vi) Gas envelope accretion from the surrounding disc occurs for planets whose masses exceed using fits to detailed 1D models from Movshovitz et al. (2010). Gas accretion occurs at the local viscous supply rate for gap forming planets. Type II migration, and gas accretion rates through the gap, have been calibrated against hydrodynamic simulations as described in CN14.
(vii) The effective capture radius of protoplanets accreting planetesimals is enhanced by atmospheric drag as described in Inaba & Ikoma (2003).
2.2 Model improvements and additions
Active turbulent region
Fully developed magnetohydrodynamic turbulence is expected to arise in regions of the disc where the temperature exceeds 1000 K (Umebayashi & Nakano, 1988; Desch & Turner, 2015). To account for the increased turbulent stress we increase the viscous parameter when the temperature rises above 1000 K using the prescription
where represents the outermost radius with temperature greater than 1000 K, and is the local disc scale height. This transition leads to a maximum in the hottest parts of the disc sitting within from the star at the beginning of the simulations.
Magnetospheric cavity and inner boundary
A rotating star with a strong dipole magnetic field may create an inner disc cavity through magnetic torques repelling the disc, and this can provide an effective mechanism for preventing planets migrating into their host stars (e.g. Lin et al., 1996). We include a cavity in our simulations by assuimg that the outer edge of the cavity is truncated at , corresponding to an orbital period of days, in agreement with the spin periods of numerous T Tauri stars (Herbst & Mundt, 2005). Planets are able to migrate into this region through either type I or type II migration. A planet that has not reached the local gap opening mass halts its migration once it reaches the cavity edge (the assumption here is that strong corotation torques will stop its migration, as shown for migrating circumbinary planets (Pierens & Nelson, 2007), and those migrating in towards a single central star (Benítez-Llambay et al., 2011)). A gap forming planet continues to migrate into the cavity until it reaches the 2:1 orbital resonance with the cavity outer edge, at which point disc torques are switched off. This resonance is located at from the star. It should be noted that a second planet entering the cavity can nudge a planet sitting at the 2:1 resonance location onto a shorter period orbit. The inner boundary of the computational domain is located just inside (corresponding to a day orbit period). Any planets whose semimajor axes are smaller than the boundary radius are removed from the simulation and are assumed to have hit the star. We note that the inner boundary adopted in CN14 corresponded to an orbital period of 20 days.
A summary of the disc and stellar parameters adopted in all simulations is given in Table 1.
|Disc inner boundary||0.02 au|
|Cavity outer boundary||0.05 au|
|Disc outer boundary||40 au|
|Number of cells||1000|
We make a small change to the opacity prescription used in CN14 by assuming that half of the disc solids are in submicron sized dust particles, with the remainder being in planetary embryos and planetesimals/boulders. The opacity used to calculate the thermal diffusion timescale in the disc is thus multiplied by the factor where is the ratio of the disc metallicity to the solar metallicity. for a disc with solar metallicity, for a disc with half the solar metallicity, and 1 for a disc with twice solar metallicity. This modification of the opacity affects both the equilibrium disc temperature and estimates for when the corotation torques acting on planets saturate.
Gas envelope accretion
Once a protoplanet grows to a mass that exceeds 3 M it starts to accrete a gaseous envelope. We have improved on the fits to the 1D giant planet formation models of Movshovitz et al. (2010) used in CN14. In units of Earth masses and Myr, our improved scheme gives a gas accretion rate of:
This scheme allows for the continuation of core growth after a gaseous envelope has been acquired, while allowing the rate of envelope accretion to adapt to the varying core and envelope mass. Figure 1 shows gas accretion onto 3, 10 and 30 cores without the influence of migration or core growth. These are similar to the models in Hellary & Nelson (2012) and CN14, but are in better agreement with the models presented by Movshovitz et al. (2010). Ideally, we would incorporate self-consistent models of gas envelope accretion in the simulations, but unfortunately this is too expensive computationally to run within our current model. While our fits to the Movshovitz et al. (2010) models allows gas accretion to occur at the rates prescribed in that work, these fits do not change according to the local conditions in the disc, or to a time varying planetesimal accretion rate. This is something that we will address in future work.
The gas accretion rate given by eqn. 2 is valid until the planet forms a gap within the disc, after which the gas accretion rate switches to either the value obtained from eqn. 2 or the viscous supply rate given by
whichever is smaller. Here and are the surface density and viscosity of the gas that sits at a distance of 10 Hill radii exterior to the planet’s location. This prescription is chosen because the planet sits in a deep gap and so the supply rate of gas must be evaluated at a location in the disc that sits outside the fully evacuated gap region. The precise value that is quoted here was determined in Section 5.2 of CN14 where different evaluation distances were tested against 2D hydrodynamic simulations, and 10 Hill radii showed the best agreement. We note that our gas accretion routine conserves mass. Gas that is accreted onto the planet is removed from the surrounding disc.
Solid bodies experience aerodynamic drag, reducing semimajor axes whilst simultaneously damping eccentricities and inclinations. Stoke’s drag is applied to planetesimals/boulders (Adachi et al., 1976) when the size of the body is greater than twice the molecular mean free path (). This switches to Epstein drag when the mean free path exceeds roughly half the planetesimal size (Weidenschilling, 1977). Here is given by:
where is the collision cross section, is the gas mean molecular weight, and is the local gas density. When the planetesimal size is greater than we use Stokes’ drag law given as:
Here, a subscript ’’ corresponds to planetesimals, is the internal density of planetesimals, is the planetesimal radius, and is the relative velocity between the gas and planetesimals. is the dimensionless drag coefficient, taken as a function of the Reynolds number () given below
When the planetesimal size is equal to , both drag regimes are equal, thus we transition to the Epstein drag law given as:
When the planetesimal size is smaller than we only use the Epstein drag law.
|Simulation||Disc mass||Metallicity||Planetesimal radius||Formation behaviour|
2.3 Initial conditions
All simulations were run for 10 Myr, allowing the systems of formed planets to continue evolving through scattering and collisions after the dispersal of the protoplanetary discs. A run time of 10 Myr is insufficient for accretion between embryos orbiting at large distances to reach completion, and some of our simulations were halted when systems of planets on longer period orbits were still evolving. This is unavoidable for systems in which large scale migration leads to the formation of short period planets, with longer period planets remaining at larger semimajor axes, since the time steps become prohibitively short for Gyr-run times to be achieved. For this reason, most of our discussion will focus on the short period systems that arise in the simulations as these are dynamically much more mature than the longer period planets.
The runs were all initiated with 52 planetary embryos, of mass , separated by 10 mutual Hill radii, and with semimajor axes between 1 and 20 . These were embedded in a swarm of thousands of planetesimals/boulders, that were distributed with semimajor axes between 0.5 and 20 , and with masses either 10, 20 or 50 times smaller than the embryos, depending on the metallicity of the system. (This varying mass ratio between embryos and planetesimals was implemented to keep the numbers of planetesimals at a number that allowed the simulations to run on reasonable times scales. Between 3000 and 8000 planetesimals/boulders were used and run times for the individual simulations varied between 3 and 9 months.) The effective physical radii of planetesimals were set to either 10 m, 100 m, 1 km and 10 km, such that the primary feedstock of the accreting protoplanets ranged from being boulders to being large planetesimals whose evolution differed principally because of the strengths of the gas drag forces that they experienced. Planetesimals/boulders in our simulations represent a larger group of particles, with realistic masses depending on their physical radii, whose averaged orbits allow them to be approximated as a single massive super-particle with an effective physical radius. Eccentricities and inclinations for protoplanets and planetesimals/boulders were randomized according to a Rayleigh distribution, with scale parameters and , respectively.
Collisions between protoplanets and other protoplanets or planetesimals resulted in perfect sticking. We neglect planetesimal-planetesimal interactions and collisions in our simulations for reasons of computational speed.
The gas disc masses simulated were 1, 1.5 and 2 times the mass of the minimum mass solar nebula (MMSN) (Hayashi, 1981). We also vary the disc metallicity so that the initial solids-to-gas mass ratios are equal to 0.5, 1 and 2 times the solar value for the different models. We smoothly increase the mass of solids exterior to the snow line by a factor of 4, as described in Hellary & Nelson (2012). We track the changes in planetary compositions throughout the simulations, as planets can accrete material originating either interior or exterior to the snow line.
Combining the three different gas disc masses, the three values of metallicity/solids-to-gas mass ratio, and the four different planetesimal/boulder sizes gives a total of 36 parameter variations. We ran two instances of each parameter set, where only the random number seed to generate initial particle positions and velocities was changed, giving a total of 72 simulations. The full set of run parameters are detailed in Table 2.
In order to provide context for our N-body simulations, we begin discussion of the results by describing the general evolution of the disc models, and the orbital evolution of the protoplanets and planetesimals. We then recap the main results obtained in CN14 before describing the results of the new simulations. We divide the results of the new runs into three distinct categories: limited planetary growth (LPG); moderate growth and migration (MGM); giant formation and significant migration (GFSM). For each category, we present the details of one or two representative runs, with Table 2 listing the category for each run. Runs that displayed limited planetary growth resulted in no planet masses growing above during the gas disc life time (and hence the amount of type I migration was also modest), although further growth beyond could occur after dispersal of the gas disc. Runs showing moderate growth and migration formed planets in the mass range during the gas disc life time. Simulations categorised as giant formation and significant migration formed planets with masses during the gas disc life time, and generally displayed multiple bursts of planetary accretion accompanied by large scale migration that ended up with one or more planets migrating into the central star. The planets that are formed in the simulations have different compositions in terms of rocky, icy and gaseous material. We use a classification system for the planets based on their compositions, and these are defined in Table 3.
3.1 Typical behaviour
Disc evolution with an active inner turbulent region
Figure 2 shows the evolution of a disc model. Disc surface density profiles are shown in the left panel, temperature profiles are shown in the middle panel, and profiles are shown in the right panel. The times corresponding to each profile are indicated in the middle panel, expressed as a percentage of the disc lifetime. For a disc this is equal to 4.6 Myr. For a disc the life time is 5.5 Myr, and a disc disperses completely after 6.5 Myr. The inclusion of a turbulent inner region where K causes a dip in surface density due to the higher viscosity there, and it can be seen that as time progresses the location of the transition to the turbulent region moves in towards the star because the reduction in surface density reduces the viscous heating rate and the opacity. The turbulent region disappears when the disc temperature no longer exceeds 1000 K anywhere in the disc, as shown by the yellow line in Figure 2. This happens in all of our disc models when the disc mass falls to approximately 10% of the MMSN, which occurs 0.5 Myr before complete dispersal of the gas disc.
The drop in local surface density caused by the active turbulent region creates a planet trap for low-mass planets (Masset et al., 2006). The trap moves in with the active region until it reaches the inner disc edge located at (assumed in our model to be outer edge of the magnetospheric cavity). Once at the disc inner edge, the trap created from the active turbulent region disappears due to the temperature in the disc falling below 1000 K. However the outer edge of the magnetospheric cavity acts as a planet trap for low-mass planets, until they can open a gap in the disc and undergo type II migration into the cavity as discussed in Section 2.2.2. It should be noted that the reduction of the temperature below 1000 K at all disc locations arises because of our adoption of a 1D disc model which neglects irradiation heating of the disc along radial lines of sight, as discussed in Section 5.
On longer time scales the removal of gas by the photoevaporative wind causes the disc to disperse. The loss of mass at large radius results in the inner disc emptying viscously onto the star, followed by removal of the remnant outer disc by the wind (Clarke et al., 2001).
Type I migration of planets is controlled by both Lindblad and corotation torques. In our disc models Lindblad torques are negative and corotation torques are generally positive. Strong, positive corotation torques arise in regions where the radial entropy gradient is negative, and this is usually the case in the inner disc regions where viscous heating dominates over stellar irradiation. Corotation torques may saturate when either the viscous or thermal time scale differs significantly from the periods of horseshoe orbits executed by gas in the corotation region. Figure 3 shows contours that illustrate the migration behaviour of planets as a function of their masses and semimajor axes in a disc with solar metallicity where half of the solid material is assumed to be in large bodies that do not provide any opacity. Dark blue regions correspond to strong outward migration, red regions correspond to strong inward migration, and white contours represent regions of parameter space where the corotation and Lindblad torques balance each other. We refer to these as zero-migration zones. The planet trap created by the inner turbulent region is shown by the innermost blue contour in the first three panels in Figure 3. Planets in blue regions migrate outwards until they come to white regions where they stop migrating. These can and do act as planet convergence zones. Planets in red regions migrate inwards, and if their masses are in the appropriate range they stop when they arrive at zero-migration zones. Over time we see that the migration contours evolve as the disc surface density and thermal time scale decrease, and planets sitting in zero-migration zones slowly drift in towards the star on the disc evolution time scale. A planet that grows in mass so that it exceeds will be too massive to sit in a zero-migration zone in the main body of the disc, and will migrate inwards rapidly before being trapped at the transition to the inner turbulent region. As this disappears the planet will drift into the magnetospheric cavity interior to where it will stop if it is below the local gap forming mass. If it exceeds the gap forming mass then it will migrate to the 2:1 resonance location with the cavity outer edge before halting its migration. If another planet enters the cavity then it may push the previous one through the inner boundary of the computational domain interior to . The decrease in values in the inner disc regions (and with time) means that it becomes possible for quite low mass planets to open gaps in the disc and enter type II migration. Similarly, planets that accrete significant gas envelopes can become giant planets and open gaps. The transition to gap formation and type II migration is shown by the boundary between the red and white contours in the top regions of the panels in Figure 3.
Each panel in Figure 4 shows the migration histories of individual planets of mass (top left), (top right), (bottom left) and (bottom right) embedded in discs with masses 1, 3 and . In each panel we plot the migration tracks of planets that have initial semimajor axes 1, 2, 5, 10 and 20 . Note that we only consider disc masses in the range 1 - 2 MMSN in the full N-body simulations described below, but we include larger disc masses in this discussion to illustrate how migration changes in significantly heavier discs. Looking at the migration trajectories it is clear that planets starting with in a disc cannot migrate interior to because of the corotation torques. Even in heavier discs planets cannot migrate very close to the star and become stranded outside the magnetospheric cavity at . The implications of this are clear. The origin of compact, short-period low mass planet systems such as Kepler 444 (Campante et al., 2015) or Kepler 42 (Muirhead et al., 2012) cannot be explained by formation at significantly larger radii than where they are observed today, followed by large scale inward migration. An in situ formation model, perhaps aided by the inward drift of solids in the form of pebbles, boulders or small planetesimals would seem to be more plausible. More generally, in situ models of planet building cannot rely on the delivery of large numbers of low mass protoplanets to inner disc regions through type I migration because they are not able to migrate across the required distances during gas disc life times. Looking at the migration trajectories, we see that these planets are also unable to reach the inner magnetospheric cavity unless orbiting in heavier discs. Guaranteed arrival of planets to the very innermost regions of the disc only occurs when planet masses reach . Periods of rapid migration observed in the lower left and right panels of Figure 4 arise when the planets saturate their corotation torques. Slow drift arises when the planets are sitting in zero-migration zones.
Planetesimal orbital evolution
Aerodynamic drag causes planetesimal eccentricities and inclinations to be damped and their semimajor axes to decrease. The 10 m boulders in our simulations experience rapid migration such that a body located initially at migrates to the inner turbulent region of the disc within aproximately years, and a 10 m boulder located at reaches there in just over years. A 100 m body located initially at reaches the inner turbulent region within Myr, and one located initially at will reach within the disc life time. The larger 1 km and 10 km bodies show very little drag-induced migration during disc life times.
The levels of planetesimal/boulder eccentricity excitation due to gravitational stirring by protoplanets at the beginning of the simulations depends strongly on their sizes. We find that the mean eccentricity for the 10 m bodies is - 5 , for the 100 m bodies - 4 , for the 1 km bodies and for the 10 km planetesimals - 3 . Given the importance of gravitational focussing in determining planetary growth rates, it is clear that we should expect smaller boulders/planetesimals to accrete much more efficiently onto the protoplanets. The mobility of the boulders also means that planetary embryos can grow beyond their nominal isolation masses on short time scales before they start to undergo significant type I migration. For protoplanets whose masses are too small for type I migration, it is the mobility of boulders and small planetesimals in our models that enables growth to occur above the isolation mass.
Recap on results from CN14
The simulations presented in CN14 adopted an inner boundary to the computational domain at an orbital period of 20 days, and so those calculations were unable to explore the formation of short period compact planetary systems. Disc masses between and were considered, with metallicities of and solar. Planetesimal sizes were 1 km and 10 km. Unlike in the present runs, no reduction in opacity was imposed to account for the growth of submicron dust grains into planetesimals/boulders and planetary embryos that populate the disc at the beginning of the simulations. Instead, the inconsistent initial condition that all disc solids were in the form of planetesimals and protoplanets, but with no diminution of the opacity, was adopted. The disc model with and solar metallity in CN14 is therefore equivalent to the model with and solar metallicity in this paper. The main results of CN14 can be summarised as:
For discs with a low to moderate abundance of solids, only limited growth of planets was observed before gas disc dispersal, although growth of planets to masses was observed due to continued mutual collisions after the disc was gone. As a consequence, only modest migration was observed in these runs, such that essentially no material was lost through the inner boundary.
It was commonly observed that numerous super-Earths and Neptune-mass planets formed in discs with intermediate masses. The bodies frequently migrated out of the disc through the inner boundary, such that in some runs no planets were left in the system at all. In others, a few super-Earth and Neptune-mass planets were able to survive.
The highest disc masses considered usually led to multiple bursts of planets forming. Gas giants with masses formed frequently. In all cases, these planets migrated rapidly through the disc via type I and II migration and out through the inner boundary. In some runs the final burst of planet formation and migration led to the formation of a short-period compact system of super-Earths and Neptunes that was able to survive. The highest mass planet to survive in all runs was a gas-rich Neptune.
CN14 determined the conditions under which giant planets could form and avoid migration into the star in their model of planet formation. They showed that the disc mass and orbital radius at which a core starts to undergo runaway gas accretion and type II migration need to be MMSN and , respectively. These same conditions also apply to the models that we present in this paper, except that our adoption of a magnetospheric cavity prevents an individual planet migrating all the way out of the computational domain interior to . It is extremely difficult in our disc model for a giant planet core to form and undergo runaway gas accretion at orbital radii because the zero-migration zones shown by the contours in Figure 3 migrate in too quickly to favour such an outcome. We note that increasing the disc viscosity and changing assumptions about the opacity model can favour the trapping of higher mass planets in zero-migration zones at larger orbital radii early during disc life times, but these zero-migration zones drift downwards and inwards more rapidly in the mass-period plane (as plotted in Figure 3) in these models because of the faster disc evolution (e.g. Bitsch & Kley, 2010).
3.2 Limited Planetary Growth (LPG)
The mass growth of planets is expected to be slow when either the
abundance of solids in the disc is small, and/or when the main feedstock for planet
building is in the form of large planetesimals whose velocity dispersion is damped weakly by the gas disc. Consequently, in the limit of slow growth, no gas accreting
cores with masses will be able to form before dispersal of the gas
The simulations labelled as LPG in Table 2 all displayed this mode of behaviour, and below we describe in detail the results of runs K10.50.01B and K2210B as they have very different disc properties, but result in similar outcomes.
Run K10.50.01B had a disc mass of , solar metallicity, and boulder radii m. The combined mass in protoplanets and boulders was equal to , distributed between , with the mass in protoplanets being initially (52 protoplanets each of mass ).
The evolution of the protoplanet masses, semimajor axes and eccentricities are shown in Figure 5 (note that boulders/planetesimals are not represented in this and similar plots). Accretion of boulders by embryos, and mutual collisions, led to the growth of protoplanets to masses in the range during the first 1 Myr. These embryos migrated towards the zero-migration zone located at and drifted in towards the star on the disc evolution time. Embryos located beyond grew more slowly, and remained near their initial locations throughout the simulation. We note that a couple of embryos at the inner edge of the solids disc experienced a short lived burst of migration by being shepherded inwards by a swarm of migrating boulders at the beginning of the simulation.
Despite the convergence of planets in the zero-migration zone, the frequency of collisions was limited by bodies entering mean motion resonances. Boulder collisions with embryos were scarce after 1 Myr, due to the drag-induced migration of boulders into the inner disc occurring on this time scale. With the maximum mass of a planetary embryo in the system being throughout the life time of the gas disc, migration remained limited to a slow inwards drift. No planets accreted gaseous envelopes.
The disc photoevaporated after 4.6 Myr, allowing embryo eccentricities to grow dramatically through mutual encounters because gas disc damping had been removed. Collisions among the inner group of protoplanets led eventually to the formation of a system of four inner bodies with masses in the range after 10 Myr when the simulation ended. These bodies all accreted significant amounts of material from beyond the snowline, and we class them as either water-rich terrestrials or water-rich super-Earths, orbiting with periods days. There were a significant number of protoplanets orbiting exterior to still undergoing collisional evolution at 10 Myr when the simulation ended, and these would have continued accreting if the run had been extended.
We turn now to run K2210B, for which the disc mass was , the metallicity was , planetesimal radii were 10 km, and the disc life time was 6.6 Myr. The initial mass in embryos and planetesimals was , this being the most solids-rich disc considered in this paper. In spite of this, planetary growth was very limited because of the weakly-damped planetesimals.
The evolution of protoplanet masses, semimajor axes and eccentricities are shown in Figure 6. Protoplanets grew to masses after 1 Myr, and when the disc dispersed the maximum embryo mass was approximately , there having been a couple of planets that accreted rapidly just prior to the final remnants of the gas being removed. Migration was limited, with the innermost body orbiting at at the point of gas disc dispersal. After removal of the gas the system entered a stage of chaotic evolution, with on-going collisions occurring within the embryo swarm when the run ended at 10 Myr. Approximately 20 planets remained at this stage, the most massive being . No planets accreted gaseous envelopes.
3.3 Moderate growth and migration (MGM)
Table 2 shows that a total of 16 out of 72 simulations exhibited moderate growth and migration. MGM runs are characterised by the formation of planets with masses before the end of the gas disc life time, with little or no loss of planets through the inner boundary of the computational domain. These simulations can result in two distinct planetary system architectures. One in which a dominant Neptune-mass body forms and migrates all the way into the magnetospheric cavity, and another where growth and migration of planets is more moderate, resulting in super-Earths and Neptunes orbiting at greater distances from the central star. Giants do not form because the growth of planets is slow enough that gas envelope accretion starts late during the disc life time, such that only moderate envelope masses have time to accrete. We discuss one representative example of an MGM run that led to the formation of a compact system of super-Earths and Neptunes on short period orbits, but no planet orbiting within the magnetospheric cavity.
Run K120.1B had a disc mass of , solar metallicity, and 100 m planetesimals. The total amount of mass in embryos and planetesimals was .
The evolution of embryo masses, semimajor axes, and eccentricities are shown in Figure 7. Several planets grew to masses during the first 0.5 Myr. A common phenomenon during the simulations involving 10 m boulders or 100 m planetesimals was the formation of shepherded rings of boulders/planetesimals while the gas disc was present, and from time to time rapid growth of a planet was observed if it crossed one of these rings through embryo-embryo scattering. At 2 Myr an embryo of mass located at grew to by accreting planetesimals from a shepherded ring, and hence started to accrete a gas envelope. The increase in mass eventually caused the corotation torques to saturate and the planet migrated in towards the star before forming a gap and transitioning to slower type II migration at Myr. Figure 7 shows that the inward migration of this planet created an inward-migrating resonant convoy, with collisions between embryos and with planetesimals leading to embryos growing within the convoy. Initially consisting of 12 protoplanets, the arrival of the convoy to the inner disc was followed by dynamical instability and collisions that left 4 short period planets remaining at the end of the simulation. These consisted of (moving out from the star) a rocky terrestrial planet, an gas-poor Neptune, a mini-Neptune, and a gas-rich Neptune, with orbital periods of 4.7, 8.3, 12.4 and 19.5 days, respectively. As all of the orbital periods are less than 100 days, this inner group constitutes a compact system, within which only one resonant pair exists, that being a 3:2 resonance between the gas-poor Neptune, and its neighbouring mini-Neptune. Other resonances existed in this group of planets and their progenitors, but were broken when strong interactions and collisions occurred. This run provides a clear example of how a short-period compact system can form through concurrent growth and migration of planets.
In the outer disc regions beyond , the dispersal of the gas disc after 4.6 Myr led to dynamical excitation of the embryos orbiting there. Planetesimals rings that had been shephered by the planets were disrupted, and a number of planets grew in mass by accreting these planetesimals. At the end of the simulation the outer region was still undergoing active accretion, and would have led eventually to the formation of long period water-rich terrestrial and super-Earth planets orbiting between if the run had been continued.
3.4 Giant formation and significant migration (GFSM)
Table 2 shows that only simulations with either 10 m boulders or 100 m planetesimals formed giant planets with masses . Out of 72 runs, 14 resulted in the formation of giants.
Gas giant planet formation ensues because a core with forms early enough that a substantial gas envelope can accrete either before the disc disperses or before the planet migrates into the inner magnetospheric cavity. In agreement with the results of CN14, we find that discs capable of forming giant planets undergo multiple bursts of planet formation and migration, with the first generation of giants being lost through the inner boundary. Unlike CN14, however, our model allows for the survival of migrating giants because they can become stranded within the magnetospheric cavity. Indeed, we formed a total of 5 surviving giants in the simulations, the most massive of which had . The most massive planet formed in any simulation had (in model K220.01A), but was lost through the inner boundary because a second generation of planets arrived in the magnetospheric cavity and pushed it through the inner boundary interior to . We discuss one run below that formed giant planets that experienced significant migration.
Simulation K1.520.1A had an initial disc mass of , a solid abundance equal to 2 solar and planetesimal radii 100 m. The mass in embryos and planetesimals was .
The evolution of protoplanet semimajor axes, masses and eccentricities are shown in Figure 8. Two planets grew above and started accreting gas envelopes within the first Myr. The saturation of corotation torques for the most rapidly growing protoplanet caused it to migrate inwards, creating a resonant convoy of co-migrating interior embryos, one of which also accreted gas. The largest mass body that drove the migration of the chain reached (with an envelope fraction of 87%) before the convoy entered the magnetospheric cavity. Gap formation prevented the planet from halting at the transition to the turbulent inner disc. The interior members of the group were pushed through the inner cavity and out of the computational domain, and the outermost planet stopped accreting gas and parked at the location of the 2:1 orbital commensurability with the outer edge of the cavity.
Shortly after 1 Myr another pair of planets exceeded , accreted gas envelopes and started to migrate rapidly when their corotation torques saturated, driving another resonant convoy inwards. These planets halted when they arrived at the transition to the active turbulent region at approximately 3.4 Myr. The outer planet in the convoy grew to , formed a gap and underwent type II migration into the magnetospheric cavity, pushing the resonant convoy and the earlier formed giant planet ahead of it. All the interior planets apart from an adjacent (formed by a collision within the cavity) were pushed through the inner boundary, leaving the and gas-rich Neptunes orbiting at 0.035 and at the end of the simulation, with gas envelope fractions of 77% and 32%, respectively.
In the interval between 2 and 4 Myr a group of Earth-mass protoplanets drifted in towards the star while sitting in a zero-migration zone, and halted their migration when the gas disc dispersed. Subsequent collisions resulted in the formation of two water-rich super-Earths, a mini-Neptune and a water-rich terrestrial planet orbiting between and with masses in the range . At large radii (2 and , respectively) two water-rich terrestial planets are formed by the accretion of plantesimals after gas disc dispersal, reaching masses at the end of the simulation at 10 Myr.
3.5 Summary of LPG, MGM and GFSM results
We now summarise the results obtained in the simulations according to which class of outcome they fall into.
Simulations classified as showing limited planetary growth led to similar outcomes despite diverse initial conditions: (i) discs with low solids abundances containing boulders and small planetesimals; (ii) discs with relatively high abundances of solids in the form of large planetesimals. The final outcomes of these simulations are summarised in the mass versus period diagram shown in the left panel of Figure 9. We see that no very short period planets were formed, and final masses are all below . The inverse correlation between mass and semimajor axis arises because of modest disc driven migration that caused the most massive bodies to drift in. The colour coding of the symbols shows that the final outcomes are similar for all boulder and planetesimal sizes.
The final states of all runs that exhibited moderate growth and migration are shown in the middle panel of Figure 9. Super-Earths and Neptune-mass planets on short period orbits are formed, and these occur almost always in compact systems (see the lower panels in Figure 16 in the appendix which shows the final outcomes of all individual runs that were classified as MGM). We note a strong inverse correlation between mass and orbital period in Figure 9 caused by migration. Low mass planets on short period orbits were shepherded in as members of resonant convoys driven by more massive planets. Within individual systems this often led to a direct correlation between mass and orbital period because migration was driven by more massive bodies at the outer edge of migrating resonant chains.
Figure 9 shows that the most massive survivors have migrated into the magnetospheric cavity. Their migration was rapid enough to send them in this far, and they are often accompanied by short-period planets that are surviving members of a resonant convoy that avoided being pushed through the inner boundary. As mentioned briefly above, runs classified as MGM can be divided into two sub-classes: those that produce objects that migrate quickly enough to reach the magnetospheric cavity, and those which do not, with faster planet growth in more solids-rich discs and/or containing smaller planetesimals/boulders leading to the first sub-class.
The final outcomes of runs classified as showing giant formation and significant migration are presented in the right panel of Figure 9. It is clear that all of the surviving gas giant planets have migrated into the magnetospheric cavity, and some of them are accompanied by interior lower mass planetary companions.
Only models with 10 m boulders and 100 m planetesimals formed giant planets with masses . All of these planets except for two were gas-dominated giants - the two exceptions being core-dominated giants (see Table 3 for definitions). For 10 m boulders the abundance of solids required to build a gas giant is equivalent to a MMSN disc with metallicity the solar value. For 100 m planetesimals a solids abundance equivalent to a MMSN disc with metallicity the solar value is required. Simulations with 1 km and 10 km planetesimals presented in CN14 show that giants would have formed in our runs if we had considered disc models with a total solids abundance equivalent to a MMSN disc with solar metallicity (e.g. a disc with solar metallicity/solids-to-gas ratio.
It is noteworthy that the most massive surviving (and non surviving) planets all formed in models with 10 m boulders. Fewer low mass planets are left at large radii in the 100 m planetesimal runs because planet growth at these radii continues to larger masses in these runs as the planetesimals do not migrate inwards too rapidly. This allows the more massive planets formed there to also migrate inwards during the gas disc life time.
3.6 Evolution as a function of planetesimal radius
The simulation results show a very strong dependence on the planetesimal size adopted, and to highlight this point we have plotted planet evolution tracks in the mass–period plane in Figure 10 for simulations with fixed disc properties (disc mass , metallicity solar) and varying planetesimal/boulder sizes: 10 m - left panel; 100 m - middle panel; 10 km - right panel. Lines ending in a black filled circle represent the formation of a surviving planet. The left panel shows the formation and rapid inward migration of gas giant planets. The middle panel shows the formation and inward migration of super-Earths and Neptune-mass planets. The right panel shows much slower growth of planets up to approximately one Earth mass and very little migration.
3.7 Evolution as a function of solid abundance
The simulation outcomes show strong dependence on the total mass in solids for a fixed planetesimal size. This is illustrated in Figure 11, which shows mass-period evolution tracks for planets in discs of varying mass and metallicity for 100 m planetesimals. The left panel shows results obtained from an anaemic disc with a mass and metallicity solar. Moderate growth and migration is observed in the middle panel for a disc mass of and metallicity solar. The right panel shows the dramatic change in evolution when the solids abundance is raised, leading to the formation of numerous Neptune-mass and gas giant planets in successive bursts, with a gas-rich Neptune remaining in the magnetospheric cavity at the end of the simulation.
4 Comparison with observations
It is important to re-emphasise that our simulation set does not constitute an attempt at population synthesis. The aim is much simpler: to examine whether or not the model of planet formation and migration presented here is able to form planetary systems similar to those that have been observed within the context of plausible disc models. We have not used a Monte Carlo approach to select initial conditions from observationally derived distribution functions, and so the frequency with which different types of systems arise in our simulations is not relevant when judging whether or not the planet formation model is successful. Comparing with observations allows us to determine whether or not the model is capable of producing planets with properties that match those of the observed population (or at least a sub-set of it), and provides a guide for understanding where model improvements are needed.
4.1 Mass versus period
Figure 12 is a mass versus period diagram for the surviving planets from all simulations, along with all confirmed exoplanets (Han et al., 2014). The vertical dashed line located at days shows the position of the disc inner edge in our simulations (i.e the location of the magnetospheric cavity).
The large number of long-period ( d) low mass planets () produced by the simulations arises because of the large number of runs that displayed limited growth (21 out of 36 disc models). These are located in a part of the mass-period diagram that is poorly sampled by radial velocity and transit surveys which are biased towards finding massive planets on short period orbits. Microlensing surveys sample this region of parameter space and although relatively few planets have been discovered, constraints obtained from statistical analysis of the data suggest that planets should be common in this region of the diagram (Gould et al., 2010).
There is good overlap between the simulation outcomes and the large numbers of observed short period terrestrial/super-Earth/Neptune-mass planets. In our simulations these planets tend to form in compact multi-planet systems, similar to those discovered by Kepler (Fabrycky et al., 2014) and radial velocity surveys (Mayor et al., 2011), as discussed in more detail below. The observational data also indicate that there are numerous systems containing a single planet or which have low multiplicity. The most recent release of Kepler data, for example, contains more than 3000 single transiting planet candidates (Mullally et al., 2015). In general, our simulations only produce systems with a short period planet and few objects (if any) orbiting significantly further out when a dominant object (Neptune or gas giant) forms and migrates through the system to the inner cavity. This scenario can clear other planets from the system, leading to low levels of multiplicity. Examples of where this occurred can be seen in Figure 16, which shows the final outcomes from all runs with short period planets. Forming single planets or low multiplicity systems without a close orbiting dominant body would seem to be difficult in the planet formation scenario presented here, and this may indicate that our choice of inserting 52 planetary embryos at the beginning of the simulations does not match the mode of planet formation occurring most commonly in nature. The prevalence of single or low multiplicity systems may be an indication that planet formation often proceeds by only forming relatively few embryos, in contrast to traditional scenarios of oligarchic and giant impact growth (Ida & Makino, 1993; Chambers & Wetherill, 1998).
The collection of very short period planets ( days) with masses in the range from the simulations all arose because they migrated into the magnetospheric cavity and were pushed closer to the star by an exterior body that was driving a resonant convoy. These outer planets, that stall finally near the 2:1 resonance with the cavity edge, are also apparent in Figure 12 and sit in a region of parameter space where there are very few observed planets. We can ascribe these distinct orbital period features in the simulated planet population as being due to adopting a single location for the cavity edge, whereas in reality it will vary from system to system (and with time) due to differences in stellar magnetic field strengths and accretion rates through protoplanetary discs. This will have the effect of blurring the locations of the planets at the 2:1 resonance location and the interior planets that have been pushed inwards. The group of more massive planets at 2 days have masses that are not commonly observed, and this may be an indication that our model fails because these bodies should have accreted more gas to become part of the hot-Jupiter population (represented by observed planets with masses ), or should experience substantial evaporation of their atmospheres by stellar X-ray irradiation on Gyr time scales (Owen & Jackson, 2012), leaving planets with smaller masses in better agreement with observations. Erosion of the atmosphere through an evaporative wind can also exert a torque on the planet allowing the planet to migrate a few percent of its semimajor axis, if the wind is anisotropic (Teyssandier et al., 2015).
One clear failing in the simulation results is the lack of surviving giant planets with masses . As mentioned earlier, the most massive planet to form in the simulations had , but migrated into the star. The formation of giant planets within our simulation occurred in the inner regions of the disc (orbital radii ), and during times when there were significant amounts of gas remaining. These giants always migrated into the magnetospheric cavity, before getting trapped at the 2:1 resonance with the disc inner edge. Generally, the last planet that migrated into this region survived, along with a less massive companion if the companion migrated in convoy. Earlier arriving planets are pushed through the inner boundary of the disc by these late arrivers. The later formation time of these surviving planets causes their masses to be smaller, as the amount of material available for accretion was reduced, explaining why there are not any genuine hot Jupiters or hot Saturns remaining at the ends of the simulations. Once again, the high multiplicity of our simulated planetary systems may be causing short period giant planets to be removed from the simulations, thus reducing the level of agreement between the models and the observations. In other words, the choice of initial conditions where embryos are equitably distributed throughout the disc may lead to too many planets forming, preventing the survival of early-forming gas giants.
Finally, we note that our models do not even come close to explaining the long period cold-Jupiter population. This is a feature of our simulations that was discussed at length in CN14, where it was shown that for giant planets to have formed and survived type II migration in our simulations, they would have had to have initiated runaway gas accretion at large orbital radii (typically ) and during sufficiently late periods of the disc life time when the total disc mass remaining was less than a few tenths of a minimum mass disc. Forming under these conditions would allow planets to undergo only a moderate amount of type II migration, allowing them to survive at large orbital radii. Trapping giant planet cores at large orbital radii until late times is difficult in our model, however, because the saturation of entropy-related corotation torques leads to rapid inwards type I migration. This point is illustrated by the migration contours shown in Figure 3.
4.2 Comparison with Kepler-like systems
Figure 13 shows a comparison between a selection of compact Kepler systems, Gliese 581 and Wasp 47 and a selection of our simulated systems. A similar figure is presented in the appendix showing all of the simulated planetary systems that arose from runs resulting in either moderate growth and migration or giant formation and significant migration.
Inspection of the simulated planetary systems in Figure 13 (and Figure 16 in the appendix) shows that we obtain two basic architectures, one where either a gas-rich Neptune or a gas giant planet has migrated through the system into the inner cavity, and another where the migration has been more modest as planet masses have not grown so massive. The runs K221B, K20.50.01B, K120.1B and K1.50.50.01B displayed the latter type of behaviour, whereas runs K220.01B, K20.50.01A and K110.01A displayed the former type. We obtain outcomes in which the planets are well separated and not in resonance, such as K221B (for which there was a lot of scattering and growth after the gas disc dispersed) and outcomes such as K1.50.50.01B where the planets are in a chain of resonances at the end of the simulation. Note that Figure 16 shows which pairs of planets in the final systems are in mean motion resonances. We also find a small number of coorbital planets at the end of the runs (3 trojan systems and 1 horseshoe system were found orbiting within 200 days across all runs. These systems are shown as being in 1:1 resonance in Figure 16). All coorbital planets were found in systems where at least one planet underwent rapid and large scale migration, causing bodies to be scattered onto eccentric orbits that quickly damped once the rapid migrator had passed through the system. This concurs with previous studies of coorbital planet formation which showed that these bodies are a direct consequence of violent relaxation in a highly dissipative environment (Cresswell & Nelson, 2006).
While it is difficult to perform a quantitative comparison between the simulated and the observed planets, certain similarities can be noted. For example, Kepler 444 looks similar to the inner four planets of K1.50.50.01B. These four rocky-terrestrial planets were shepherded in by the exterior more massive water-rich terrestrials, and hence formed a resonant convoy. This is one way in which the Kepler 444 planets could have arrived at their observed locations and provides an alternative to in situ formation (but relies on there being a more massive, undetected planet orbiting further from the star). Kepler 169, 186 and 80 look similar to K20.50.01B, and Kepler 11 and 33 have broad similarities with K120.1B. Although the Kepler sample does not contain examples of compact multi-systems with massive, short period planets (perhaps because these are more dynamically disturbed and therefore not transiting or close to resonances such that they are detectable through transit timing variations), Gliese 581 and Wasp 47 provide two examples that have architectures similar to K210.1B and K220.01B.
4.3 Period ratios and orbital spacings
Figure 14 compares the cumulative distributions of period ratios between neighbouring planets with masses and orbital periods less than 100 days obtained from the simulations (upper blue curve) and the Kepler systems (lower red curve). The sample of Kepler planets was defined by choosing bodies with orbital periods days and radii . This lower radius limit was adopted to account for possible incompleteness in the Kepler sample for planets with small radii. It is clear that the simulated systems are generally more closely packed after run times of 10 Myr, and the structure observed in the distribution shows that this is due in part to there being a number of planet pairs in resonance. The step-like features in the plot show that the 7:6, 6:5, 5:4, 4:3, 3:2 and 2:1 resonances are occupied. Whereas just an isolated pair of migrating planets are likely to be trapped in either the 2:1 or 3:2 resonances if they undergo smooth migration (Paardekooper et al., 2013), we find that migration in a crowded system allows diffusion through successive resonances to occur such that high degree resonances can be occupied, in agreement with earlier studies by Cresswell & Nelson (2006, 2008). Although resonant systems are relatively rare in the Kepler data, it is worth noting that Kepler 36 has two planets very close to the 7:6 resonance (Carter et al., 2012; Paardekooper et al., 2013), and some of the planet pairs in Kepler 444 are reported to be in 5:4 (Campante et al., 2015). Other examples of systems in resonance or near resonance, including 3 body resonances and resonant chains, are Kepler 50 (6:5), Kepler 60 (5:4, 4:3) (Steffen et al., 2012), Kepler 221 (3 body resonance where the mean motion combination has been found to librate around 180 degrees) (Fabrycky et al., 2014).
Furthermore, it has been noted in numerous studies (e.g. Fabrycky et al., 2014) that the distribution of planet period ratios contains an excess of planets just outside of 3:2 and 2:1, suggesting that the resonances have been dynamically important during the evolution but may have been broken by stochastic migration in a turbulent disc (Adams et al., 2008; Rein & Papaloizou, 2009), by tidal interaction with the central star (Terquem & Papaloizou, 2007), by orbital repulsion due to damping of nonlinear spiral waves (Baruteau & Papaloizou, 2013), by overstability in librations about resonant centres (Goldreich & Schlichting, 2014), or because of scattering due to interactions with or accretion of residual planetesimals (Chatterjee & Ford, 2015). We observe that in a handful of simulations, planetesimal scattering after full gas disc dispersal does occur, breaking mean-motion resonances between neighbouring planets, in agreement with Chatterjee & Ford (2015). It is noteworthy that a number of the compact systems are orbiting in regions where their nascent protoplanetary discs are expected to have sustained MRI turbulence due to the local temperature being in excess of 1000 K (Umebayashi & Nakano, 1988), and so may have been subjected to stochastic forcing of their orbits while the gas disc was present. To seek evidence for this transition to turbulence we have examined the minimum periods of planets in the compact Kepler multi-systems to see if they correlate with the effective temperature of the host star, but there is no evidence of a correlation. At present there is no clear evidence that the transition to turbulence in the inner regions of the protoplanetary discs that formed the Kepler systems played a decisive role in dynamically shaping these systems.
It is possible that a number of our simulated systems may be dynamically unstable on time scales much longer than the 10 Myr run times, such that subsequent mutual collisions increase separations between adjacent planets. In a recent study, Pu & Wu (2015) used N-body simulations to show that compact Kepler-like multi-planet systems tend to remain stable for Gyr time scales only if the typical mutual separation between neighbouring planets is approximately 12 mutual Hill radii. Figure 15 shows the distribution of separations between neighbouring planets present at the end of the simulations, and while many planet pairs are well separated there are a significant number whose orbital spacings may be too small for long term stability. Running the simulations for long enough to test this goes beyond the scope of this paper, but will be studied in future work as it may be the case that the mean motion resonances discussed above provide protection against instability. We note that the objects with period ratios of unity shown in Figure 15 are the coorbital planets mentioned previously.
5 Discussion and conclusion
We have implemented a model of planet formation based on a scenario in which numerous planetary embryos are distributed across a wide range of semimajor axes, embedded in a sea of boulders or planetesimals that act as the primary feedstock for planetary growth. The model has a comprehensive list of ingredients: planetary embryo growth through boulder/planetesimal accretion and mutual collisions; a 1D viscous gas disc model, subject to irradiation from the central star and a photoevaporative wind; type I migration using the most up-to-date prescriptions for Lindblad and corotation torques; a transition to gap formation and type II migration when gap formation criteria are satisfied; gas accretion onto solid cores. The disc has an increase in viscosity where the temperature K, to mimic unquenched MHD turbulence developing in the inner disc, and a magnetospheric cavity that creates an inner edge in the gas disc at an orbital period of 4 days. The aim of this study is to determine which types of planetary systems emerge from the planet formation model as a function of disc parameters (mass and metallicity) and planetesimal/boulder sizes. The main results from our simulations can be summarised as follows.
(i) System evolution can be categorised into three distinct modes that depend on the total amount of solids present in the disc and the sizes of the boulders/planetesimals.
- When planetesimal/boulder radii are small ( m) limited planetary growth arises when the inventory of solids is small. When planetesimal radii are large ( km), limited growth arises for all discs models considered, except the one that is the most massive and solids-rich. Planets with maximum masses form during the gas disc life times, and show only very modest migration.
- Moderate growth and migration arises in only the most solids-rich disc considered when planetesimal sizes are 1 km, and for disc models with intermediate abundances of solids when the planetesimal/boulder sizes m. Planets are able to grow to super-Earth or Neptune masses during the disc life time, and may undergo large-scale migration.
- Giant formation and significant migration is observed in the most solids-abundant discs when boulder/planetesimal sizes were m, but did not arise in any of the runs with larger planetesimals. Generally, multiple episodes of planet formation occur, and gas giant planets with masses form and undergo large scale migration before stalling in the magnetospheric cavity. The final surviving short period planets are normally the last ones to arrive in the magnetospheric cavity, with the earlier arrivals being pushed through the inner boundary by the planets that arrive there later.
(ii) Considering systems of short-period planets, we can identify two basic architectures that emerge from the simulations. The first normally consists of a combination of terrestrial planets, super-Earths and low mass Neptunes, where no planet managed to migrate into the magnetospheric cavity. The shortest period orbits in these systems are normally 4 -5 days. The second architecture consists of at least one dominant planet (a gas giant or a relatively massive Neptune) that migrated and stalled in the magnetospheric cavity with a period of days. In approximately 50% of cases, this planet has an interior companion (terrestrial planet, super-Earth or Neptune) which is almost never in resonance because of dynamical interactions and collisions with other planets during the evolution. In most cases where a dominant short period planet formed, there are a number of exterior planets orbiting with periods in the range days.
(iii) The planetary systems display a range of heterogeneity in composition versus orbital period. Systems that formed under relatively quiescent conditions, without a rapidly migrating gas giant or Neptune, have rocky bodies orbiting interior to volatile rich bodies. Systems that contained rapidly migrating giants or Neptunes, that end up in 2 day orbits, often experienced significant scattering, and these systems can have rocky bodies in exterior orbits in close proximity to volatile-rich bodies.
(iv) The planetary systems that emerge from the simulations tend to be closer packed than the observed Kepler systems. The most common spacing between neighbouring planets is 10 - 12 mutual Hill radii, and Pu & Wu (2015) have shown that such systems are likely stable over Gyr time scales. There are, however, numerous simulated planet pairs where the ratio of spacing to mutual Hill radius , and these might cause the systems to evolve and change their spacing through collisions if evolved beyond the 10 Myr that we have considered, improving the agreement with observations. We note, however, that mean motion resonances may help stabilise our simulated systems compared with those considered by Pu & Wu (2015).
(v) One reason for the difference in the distributions of observed versus simulated period ratios is that mean motion resonances are common among our final planetary systems. We find examples of 7:6, 6:5, 5:4, 4:3, 3:2 and 2:1, with the latter three resonances being rather common. It is well known that most of the compact Kepler systems do not display mean motion resonances, even though there is evidence for the 2:1 and 3:2 resonances having been dynamically important in the past, and a few individual systems appear to host resonant pairs or triples. One possible explanation for the greater numbers of resonant systems arising in the simulations is the neglect of stochastic forces in the inner disc regions due to MHD turbulence (Nelson & Papaloizou, 2004; Nelson, 2005) which can cause planets to diffuse out of resonance (Adams et al., 2008; Rein & Papaloizou, 2009). It remains to be seen whether or not inclusion of this effect can increase the agreement between observations and theory on the frequency of mean motion resonances. One further point worthy of note is that the frequency of resonances arising in the simulations is higher for those architectures that contain a dominant planet orbiting with a 2 day period. Systems without a dominant short period planet underwent more quiescent evolution during the gas disc life time, but also experienced more scattering after removal of the disc and this leads to systems that contain few resonances (see Figure 16). Thus, it is important to note that there is a mode of planet formation that includes large scale migration but which does not result in systems that are members of resonant chains.
(vi) A number of coorbital planets were formed in our simulations (three trojan systems, and one undergoing mutual horseshoe orbits, were found to orbit with periods days). These all formed in systems where at least one dominant planet underwent migration through the planetary swarm, causing large amounts of scattering. In earlier work Cresswell & Nelson (2006, 2008) have shown that coorbital planets arise as a consequence of violent relaxation in crowded planetary systems with strong eccentricity damping, and our results are in agreement with these earlier findings.
(vii) Numerous gas giant planets were formed in our simulations, and some survived after migrating into the magnetospheric cavity. The most massive planet to form was a gas giant, but this was pushed through the inner boundary of the computational domain by a planet that arrived in the magnetospheric cavity at a later time. The most massive surviving planet was a “hot Saturn” on a 2 day orbit. CN14 undertook a detailed examination of the conditions required for the formation and survival of longer period giant planets against type II migration, and showed that a Jovian mass planet halting its migration at needs to start runaway gas accretion and type II migration at a distance of from the central star. This has not occurred in any of our simulations (this paper, or CN14, or in the many low resolution test simulations that we have run and not published), because of the difficulty of forming a core and keeping it at such large orbital radius. We have concluded that forming and retaining long period giant planets requires a set of disc conditions that are quite different from those that we have considered thus far. A potential solution to the problem will be presented in a forthcoming paper (Coleman & Nelson in prep.).
5.1 Formation of Kepler 444 and 42
The Kepler 444 and 42 systems are examples of short period compact low mass planetary systems. All have radii substantially smaller than the Earth’s. Kepler 444 is a 5-planet system orbiting a M K0V star with , where the innermost orbital period is 3.6 days and the outer planets are close to the 5:4, 4:3, 5:4 and 5:4 mean motion resonances (Campante et al., 2015). Kepler 42 is a 3-planet system orbiting a M M3V star with . Orbital periods are 0.453, 1.214 and 1.865 days (Muirhead et al., 2012), so there are no first-order mean motion resonances. We showed in Sect 3 that planet masses need to be in excess of for migration over large distances to be effective, and given the low metallicities of these systems they are most likely explained by in situ formation after delivery of solids through drag-induced drift into the disc inner regions. Although large scale migration of these planets is implausible, the resonant or near-resonant configuration of the Kepler 444 planets suggests that modest migration may have occurred. The outermost planet being the largest (and presumably most massive) would lead to the necessary convergent migration.
5.2 Formation of short period super-Earths in low metallicity discs
Our simulations demonstrate how difficult it is to grow planets that are massive enough to undergo significant type I migration during the gas disc life time when growth is dominated by the accretion of large ( km) planetesimals in discs with a moderate inventory of solids. This is because growth time scales are slow for large planetesimals. In addition, if a planet approaches its local isolation mass it will not be massive enough to migrate such that it can accrete from undepleted sources of planetesimals. The situation becomes more difficult in a low metallicity environment, and the existence of short period super-Earths around stars such as Kapteyn’s star (Anglada-Escudé et al., 2014), Gliese 581 (Udry et al., 2007), HD 175607 (Mortier et al., 2015) and the numerous low-metallicity hosts of Kepler systems (Buchhave et al., 2014) suggests that these planets did not form via the classical oligarchic growth picture of widely distributed embryos accreting from a swarm of large planetesimals. These systems instead point towards planetary embryos growing into type I migrating super-Earths by accreting from a supply of highly mobile small planetesimals, boulders or pebbles (e.g. Ormel & Klahr, 2010; Lambrechts & Johansen, 2012), as this is the only means available of exceeding local isolation masses. On the other hand, the requirement for the local solids-to-gas ratio to be approximately twice solar in order for the streaming instability to operate and generate large planetesimals that can acts as the seeds of growing planets (Johansen et al., 2009b) suggests that small particles must first concentrate in specific disc regions due to the existence of zonal flows (Johansen et al., 2009a; Bai & Stone, 2014), vortices (Fromang & Nelson, 2005) or dead zone interfaces (Lyra et al., 2009) in order to create local enhancements of solids. Such a collect-and-grow scenario would appear to offer the best hope for explaining the existence of planets in the lowest metallicity environments.
5.3 Future work and directions
The long term aim of this project is two-fold: to construct a simulation tool for modelling planet formation that comprises accurate prescriptions for the important ingredients for planet building and evolution; to determine whether or not it is possible to explain the diversity of known planetary systems using a comprehensive model of planet formation, loosely based on the classical core accretion model, operating under different initial conditions and environments. A particular issue of interest is explaining the known population of gas giant planets, and we will present a study of this in a forthcoming paper. Areas of future improvement to our model include:
(i) Calculation of gas envelope accretion using self-consistent computations that take account of the changing local nebula conditions, rather than using fits to the Movshovitz et al. (2010) models as is done now. The atmosphere models of Papaloizou & Nelson (2005) are being incorporated into the code, and results from these calculations will be presented in a forthcoming publication.
(ii) Improving the disc model so that stellar irradiation of the disc inner regions is treated more accurately. Our 1D treatment of stellar irradiation underestimates the level of heating near the star, and this allows the temperature to fall below 1000 K everywhere in the disc at late times, such that no region of the disc maintains fully developed MRI turbulence. A more realistic treatment would allow the temperature to always be above 1000 K out to a radius , which for our model corresponds to a distance of from the star.
(iii) Include the effects of stochastic migration when planets and planetesimals enter disc regions where K. This will influence the ability of planet pairs to maintain mean motion resonances.
(iv) Improve the migration model by including 3D effects (Fung et al., 2015) and the influence of the planet luminosity Benítez-Llambay et al. (2015).
(v) Small boulders and planetesimals are able to migrate inwards from beyond the snow line, and in principle these should sublimate quite rapidly. We have not included submlimation in our models, and analysis of the results indicates that planets accreting icy planetesimals that have migrated interior of the ice line increase their masses by at most . Nonetheless a model of planetesimal sublimation should be included for self-consistency.
Simulations incorporating these improvements will be presented in future publications in order to determine how these modifications change the simulation outcomes. Once a suitably sophisticated model of planet formation has been constructed it will be used to generate a synthetic planet population based on initial conditions and physical parameters drawn from observational constraints to determine the level of agreement with exoplanet observations.
Appendix A Presentation of all simulated compact systems
Figure 16 shows all of the compact systems that formed and survived in the simulations. The planets shown in the upper panels all formed in simulations classified as giant formation and significant migration, and the rest formed in simulations classified as moderate growth and migration. Pairs of planets that are coorbital or are in first order mean motion resonances are indicated by the integers printed above and between the relevant pair.
- pagerange: On the formation of compact planetary systems via concurrent core accretion and migration–References
- We note that planetary atmospheres may form via outgassing, but this effects goes beyond the range of physical processes considered in our models. Furthermore, H/He rich envelopes can settle onto relatively low mass planets (Lammer et al., 2014), and although we consider the effect of this on planetesimal accretion, we do not report gas envelope masses for planets with in this paper.
- Adachi I., Hayashi C., Nakazawa K., 1976, Progress of Theoretical Physics, 56, 1756
- Adams F. C., Laughlin G., Bloch A. M., 2008, ApJ, 683, 1117
- Alexander R. D., Armitage P. J., 2009, ApJ, 704, 989
- Alibert Y., et al., 2006, AA, 455, L25
- Anglada-Escudé G., et al., 2014, MNRAS, 443, L89
- Bai X.-N., Stone J. M., 2014, ApJ, 796, 31
- Baruteau C., Papaloizou J. C. B., 2013, ApJ, 778, 7
- Benítez-Llambay P., Masset F., Beaugé C., 2011, AA, 528, A2
- Benítez-Llambay P., Masset F., Koenigsberger G., Szulágyi J., 2015, Nature, 520, 63
- Bitsch B., Kley W., 2010, AA, 523, A30
- Buchhave L. A., et al., 2014, Nature, 509, 593
- Campante T. L., et al., 2015, ApJ, 799, 170
- Carter J. A., et al., 2012, Science, 337, 556
- Chambers J. E., 1999, MNRAS, 304, 793
- Chambers J. E., Wetherill G. W., 1998, Icarus, 136, 304
- Chatterjee S., Ford E. B., 2015, ApJ, 803, 33
- Chatterjee S., Tan J. C., 2014, ApJ, 780, 53
- Clarke C. J., Gendrin A., Sotomayor M., 2001, MNRAS, 328, 485
- Coleman G. A. L., Nelson R. P., 2014, MNRAS, 445, 479
- Cossou C., Raymond S. N., Pierens A., 2013, AA, 553, L2
- Cresswell P., Nelson R. P., 2006, AA, 450, 833
- Cresswell P., Nelson R. P., 2008, AA, 482, 677
- Desch S. J., Turner N. J., 2015, ApJ, 811, 156
- Dullemond C. P., Hollenbach D., Kamp I., D’Alessio P., 2007, Protostars and Planets V, p. 555
- Fabrycky D. C., et al., 2014, ApJ, 790, 146
- Fendyke S. M., Nelson R. P., 2014, MNRAS, 437, 96
- Fromang S., Nelson R. P., 2005, MNRAS, 364, L81
- Fung J., Artymowicz P., Wu Y., 2015, ApJ, 811, 101
- Goldreich P., Schlichting H. E., 2014, AJ, 147, 32
- Gould A., et al., 2010, ApJ, 720, 1073
- Han E., Wang S. X., Wright J. T., Feng Y. K., Zhao M., Fakhouri O., Brown J. I., Hancock C., 2014, \pasp, 126, 827
- Hands T. O., Alexander R. D., Dehnen W., 2014, MNRAS, 445, 749
- Hansen B. M. S., Murray N., 2012, ApJ, 751, 158
- Hayashi C., 1981, Progress of Theoretical Physics Supplement, 70, 35
- Hellary P., Nelson R. P., 2012, MNRAS, 419, 2737
- Herbst W., Mundt R., 2005, ApJ, 633, 967
- Ida S., Lin D. N. C., 2010, ApJ, 719, 810
- Ida S., Makino J., 1993, Icarus, 106, 210
- Inaba S., Ikoma M., 2003, AA, 410, 711
- Johansen A., Youdin A., Klahr H., 2009a, ApJ, 697, 1269
- Johansen A., Youdin A., Mac Low M.-M., 2009b, \apjl, 704, L75
- Jontof-Hutter D., Rowe J. F., Lissauer J. J., Fabrycky D. C., Ford E. B., 2015, Nature, 522, 321
- Kasting J. F., Whitmire D. P., Reynolds R. T., 1993, Icarus, 101, 108
- Lambrechts M., Johansen A., 2012, AA, 544, 32
- Lammer H., et al., 2014, MNRAS, 439, 3225
- Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846
- Lin D. N. C., Bodenheimer P., Richardson D. C., 1996, Nature, 380, 606
- Lissauer J. J., et al., 2011, Nature, 470, 53
- Lovis C., et al., 2006, Nature, 441, 305
- Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
- Lyra W., Johansen A., Zsom A., Klahr H., Piskunov N., 2009, AA, 497, 869
- Marcy G. W., Weiss L. M., Petigura E. A., Isaacson H., Howard A. W., Buchhave L. A., 2014, Proceedings of the National Academy of Science, 111, 12655
- Masset F. S., Morbidelli A., Crida A., Ferreira J., 2006, ApJ, 642, 478
- Mayor M., et al., 2011, preprint, (arXiv:1109.2497)
- McNeil D. S., Nelson R. P., 2009, MNRAS, 392, 537
- McNeil D. S., Nelson R. P., 2010, MNRAS, 401, 1691
- Mortier A., et al., 2015, preprint, (arXiv:1511.03941)
- Movshovitz N., Bodenheimer P., Podolak M., Lissauer J. J., 2010, Icarus, 209, 616
- Muirhead P. S., et al., 2012, ApJ, 747, 144
- Mullally F., et al., 2015, \apjs, 217, 31
- Nelson R. P., 2005, AA, 443, 1067
- Nelson R. P., Papaloizou J. C. B., 2004, MNRAS, 350, 849
- Ormel C. W., Klahr H. H., 2010, AA, 520, A43
- Owen J. E., Jackson A. P., 2012, MNRAS, 425, 2931
- Paardekooper S.-J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950
- Paardekooper S.-J., Baruteau C., Kley W., 2011, MNRAS, 410, 293
- Paardekooper S.-J., Rein H., Kley W., 2013, MNRAS, 434, 3018
- Papaloizou J. C. B., Nelson R. P., 2005, AA, 433, 247
- Pierens A., Nelson R. P., 2007, AA, 472, 993
- Podlewska-Gaca E., Papaloizou J. C. B., Szuszkiewicz E., 2012, MNRAS, 421, 1736
- Pu B., Wu Y., 2015, ApJ, 807, 44
- Rein H., 2012, MNRAS, 427, L21
- Rein H., Papaloizou J. C. B., 2009, AA, 497, 595
- Shakura N. I., Sunyaev R. A., 1973, AA, 24, 337
- Steffen J. H., et al., 2012, ApJ, 756, 186
- Terquem C., Papaloizou J. C. B., 2007, ApJ, 654, 1110
- Teyssandier J., Owen J. E., Adams F. C., Quillen A. C., 2015, MNRAS, 452, 1743
- Udry S., et al., 2007, AA, 469, L43
- Umebayashi T., Nakano T., 1988, Progress of Theoretical Physics Supplement, 96, 151
- Weidenschilling S. J., 1977, MNRAS, 180, 57
- Wu Y., Lithwick Y., 2013, ApJ, 772, 74