Tidal Downsizing model. II. Planet-metallicity correlations
Core Accretion (CA), the de-facto accepted theory of planet formation, requires formation of massive solid cores as a prerequisite for assembly of gas giant planets. The observed metallicity correlations of exoplanets are puzzling in the context of CA. While gas giant planets are found preferentially around metal-rich host stars, planets smaller than Neptune orbit hosts with a wide range of metallicities. We propose an alternative interpretation of these observations in the framework of a recently developed planet formation hypothesis called Tidal Downsizing (TD). We perform population synthesis calculations based on TD, and find that the connection between the populations of the gas giant and the smaller solid-core dominated planets is non linear and not even monotonic. While gas giant planets formed in the simulations in the inner few AU region follow a strong positive correlation with the host star metallicity, the smaller planets do not. The simulated population of these smaller planets shows a shallow peak in their formation efficiency at around the Solar metallicity. This result is driven by the fact that at low metallicities the solid core’s growth is damped by the scarcity of metals, whereas at high metallicities the fragments within which the cores grow contract too quickly, cutting the core’s growth time window short. Finally, simulated giant gas planets do not show a strong host star metallicity preference at large separations, which may explain why one of the best known directly imaged gas giant planet systems, HR 8799, is metal poor.
Given the bewildering diversity of exoplanet system architectures (e.g., Batalha et al., 2013; Winn & Fabrycky, 2014), it is clear that the outcome of planet formation is a highly stochastic process. Any statistical trends found in the observed planetary populations have special significance as they testament to planet-forming processes so robust that they rise above the stochasticity. A successful planet formation theory must reproduce such trends.
Metallicity correlations of the observed planets is one such correlation. Giant planets are detected much more frequently around metal-rich stars than around metal-poor ones (Gonzalez, 1999; Fischer & Valenti, 2005). This correlation has been argued to provide a direct support to Core Accretion theory for planet formation (e.g., Pollack et al., 1996) since metal-rich environments assemble massive cores much more readily (e.g., Ida & Lin, 2004, 2008; Mordasini et al., 2009). These cores are the crucial step towards forming a gas giant planet in CA framework, thus a positive giant planet–metallicity correlation ensues. Gravitational disc Instability model for planet formation (e.g., Boss, 1998; Helled et al., 2013) was argued to produce a negative rather than positive giant planet–metallicity correlation since planets’ radiative cooling is the fastest at low dust opacities (Helled & Bodenheimer, 2011), increasing planet’s chances of survival.
However, there are both observational and theoretical reasons to study the issue further. Radial velocity observations show that Neptune-mass planet occurrence does not seem to correlate with the host star’s metallicity for FGK stars (Sousa et al., 2008). Similarly, transit method observations with Kepler, sensitive to the planet’s radius rather than its mass, show (Buchhave et al., 2012) that planets smaller than form around hosts with a wide range of metallicities, with the average close to the Solar metallicity. This is clearly surprising in the context of CA, since these smaller planets are the precursors of the gas giant planets. If CA’s explanation for the positive gas giant planet correlation with metallicity is correct then there should be more cores at higher metallicities, which is not what is observed.
On the theoretical side, several important extensions to the classical variant of the Gravitational Instability (GI) model for planet formation have been recently proposed. In particular, GI model that includes planet migration and pebble accretion was showed to produce a strong positive and not negative correlation with the star’s metallicity (Nayakshin, 2015a, b) for coreless gas giant planets. The goal of this paper is to extend this recent work on giant planets with cores and also on core-dominated (smaller) planets and to compare the results with the observations.
The key reason why GI theory may require a major overhaul is the realisation that GI gas fragments can also migrate from AU all the way into the inner disc (Boley et al., 2010). While planet migration was a standard feature for the CA model since 1996 (Lin et al., 1996), GI planets were somehow thought to not migrate until recently. If this were true than GI could at best account for important but rare giant planet systems like HR 8799 imaged directly (e.g., Marois et al., 2008), and could never explain the Solar System: massive self-gravitating proto-planetary discs can hatch clumps only at many tens of AU (e.g., Rice et al., 2005; Rafikov, 2005).
However, simulations, starting with Vorobyov & Basu (2005, 2006), showed that massive gas fragments do migrate inward. Furthermore, Boley et al. (2010) suggested a new way of forming terrestrial-like planets inside a few Jupiter masses gas clumps by grain growth and sedimentation (see also Williams & Crampin, 1971; Boss, 1997), and then releasing them back into the disc by destroying the gas clumps via tidal forces from the host star. Nayakshin (2010a) used analytical estimates of the relevant processes and arrived at similar ideas, proposing the Tidal Downsizing (TD) hypothesis for formation of all types of planets observed so far.
The key question for a planet forming in the framework of the TD hypothesis is how does the clump’s inward migration time scale, which may be as short as years (Baruteau et al., 2011; Michael et al., 2011), compare with the time scale for its internal evolution? If the fragment contraction time is shorter than the inward migration time, then the fragment collapses, becoming a bona fide ”hot start” protoplanet, which may mature into a gas giant planet. If the fragment contraction time is long, then the host’s tidal forces catch up with the fragment when the latter migrates too close to the star, and the fragment is disrupted. The fragment’s short existence may however not be all in vain as there could be a remnant. If grains within the fragment had a sufficient time to grow and sediment to the centre, getting locked into a self-gravitating massive core (see also Kuiper, 1951; McCrea & Williams, 1965; Williams & Crampin, 1971; Boss, 1997; Helled & Schubert, 2008), then the remnant of the disruption is the core – a practically ready rocky planet. No planetesimal accretion is required for the planet to survive at that stage, although ”veneer accretion” of them or pebbles (Johansen & Lacerda, 2010) is possible.
Making detailed predictions in the TD hypothesis setting worthy of comparison to modern observations of exoplanets is however not trivial. Until recently, it was assumed that collapse of the self-gravitating GI gas fragments must occur due to radiative cooling (e.g., Bodenheimer, 1974; Bodenheimer et al., 1980; Cameron et al., 1982; Helled et al., 2008; Nayakshin, 2011; Nayakshin & Cha, 2013). The formation of giant gas planets in this picture is similar to that of low mass stars (e.g., Larson, 1969), except that accretion of gas onto the pre-collapse planets is terminated very early on (e.g., Nayakshin, 2010b). This formation channel favors survival of massive ( a few ) gas giants (Forgan & Rice, 2013), as more massive planets cool more rapidly (e.g., see Nayakshin, 2015a). Also, radiative collapse of H-dominated fragments produces a negative planet frequency – metallicity correlation (Helled & Bodenheimer, 2011) because radiative cooling is more rapid at low metallicities. That is clearly at odds with observations (Gonzalez, 1999; Fischer & Valenti, 2005).
We (Nayakshin, 2015a), however, considered accretion of cm sized grains (usually called ”pebbles”; Johansen & Lacerda, 2010; Ormel & Klahr, 2010) onto the pre-collapse fragments. While gas has pressure gradient forces that may prevent its accretion from the protoplanetary disc onto the fragments (e.g., Nayakshin & Cha, 2013), pebbles weakly coupled to the gas by aerodynamical forces can still accrete onto the fragment in an analogy to the pebble accretion process onto massive solid cores (Lambrechts & Johansen, 2012) in the context of the CA model.
Surprisingly, pebble accretion on GI fragments was found to accelerate their
contraction despite increasing their metallicity and hence
opacity. Physically, addition of pebbles to the fragment increases its weight
without increasing its thermal energy
An important issue not considered by Nayakshin (2015a, b) is formation of a core within the fragment. If pebbles entering the fragment accrete onto the core, the core’s accretion luminosity may be significant and this may lead to expansion of the fragment instead of its contraction. A similar effect – accretion of planetesimals onto the growing core – is well known in the context of the CA model and has to be taken into account when determining the fate of the planet (e.g., Pollack et al., 1996; Alibert et al., 2005). Nayakshin (2015c, paper I hereafter) has recently presented numerical algorithms that extend the work of Nayakshin (2015b) by adding the core formation and growth inside the pre-collapse fragments.
Here we use that framework to study the planet frequency of survival versus metallicity dependence for all types of planets rather than just the coreless giants. We also go beyond of our previous work by studying not just the outcome of the migration vs disruption competition, but also how the population of planets formed in TD model would look like in the planet mass – separation plane.
The structure of our paper is as following. In §2 we present initial conditions, assumptions and numerical algorithms that are employ to calculate an outcome of a gas-dust fragment being born in an outer region of a massive protoplanetary disc. In §3, two example planet formation tracks are shown, one for a gas giant and another for a hot Super-Earth planet. In §4 we present assumptions and methods for calculating a grid of models in a population synthesis like study. §5 presents the population synthesis study result in the planet mass versus planet-host separation plane. In §6 we marginalise over the results to derive the planet-metallicity correlations for our simulated planets. Discussion of main ideas and results of the paper is given in §7, whereas conclusions are presented in §8.
2 Numerical methods
2.1 Disc treatment
We follow the procedures described in paper I with several changes detailed now. Our goal is to simulate formation of planets from their birth to the end of the disc migration phase. The initial disc surface density profile is given by
where is the disc length-scale, set to AU, AU is the inner boundary radius. The constant is calculated so that the disc contains a given initial mass . The disc is ”live”, that is evolved by solving the time-dependent 1D viscous disc evolution equations that include the planet-disc interactions (section 3 in paper I), and now also include the disc photo-evaporation as a sum of the UV and the X-ray driven terms, using the fits to the photo-evaporation rates from Alexander & Armitage (2007) and Owen et al. (2012), respectively. The ionising photon luminosity of the star is set to photons s, while the X-ray flux of the star is erg s.
Gas fragments in the disc are born by self-gravitational instabilities, hence the disc must be initially gravitationally unstable in its outer region. Figure 1 shows the initial disc properties for disc mass and stellar mass . The top panel shows the disc surface density, , as well as the disc photo-evaporation rate profile (which does not change during the simulation). The quantity plotted with the red dotted curve is , in units of year. The bottom panel shows the Toomre -parameter (black solid curve) and the dimensionless cooling time, . For the disc to be self-gravitating and actually fragmenting, and (Gammie, 2001). Both of these conditions are satisfied at AU.
The disc viscosity is described with the Shakura & Sunyaev (1973) viscosity parameter, but we now add a self-gravity contribution to it:
where is radius and time-independent (free) parameter of the model, whereas describes gravito-turbulence
where and is the local Toomre’s parameter. This form is inspired by Lin & Pringle (1987) suggestion that in the self-gravitating regime, on the one hand, and more recent simulations showing that saturates at about 0.1 (Rice et al., 2005), on the other.
2.2 Planet’s birth and migration
The following considerations are important for both the initial and the final phases of our simulations. Simulations of self-gravitating protoplanetary discs show that fragments rarely form in isolation, hence we believe that every star is likely to have hatched a number of fragments early on. Presumably, most of these fragments migrate rapidly and perish by being disrupted (see, in particular, Vorobyov, 2013) or by being driven all the way into the star. This picture is consistent with the suggestion that disruption and swallowing of planets by their host stars power FU Ori outbursts of young protostars (Vorobyov & Basu, 2006; Boley et al., 2010; Nayakshin & Lodato, 2012), statistics of which suggests that every star may have of such episodes (Hartmann & Kenyon, 1996).
Given that the observed giant gas planets’ frequency of occurrence is % (Winn & Fabrycky, 2014) per star, the survival of a giant gas planet is a very rare event indeed. It is reasonable to assume that giant gas planets that survive to be observed are those that were either hatched by the disc late, or those that hang around at the outer disc for a while before migrating in, because such fragments are more likely to avoid the tidal disruption and also being driven all the way into the star.
This view is not unreasonable. Simulations show that presence of one clump in a disc may induce formation of more clumps in the same disc (e.g., Meru, 2013). Multiple clump formation leads to strong interactions between the clumps, so that some of them migrate inward faster while others are scattered on larger orbits (e.g., Cha & Nayakshin, 2011); some of the clumps may be completely ejected from the system by the clump-clump interactions (Basu & Vorobyov, 2012). Clumps on scattered but still bound orbits may take a while to loose their orbital inclinations with respect to the midplane of the disc, and eccentricities, and only then start migrating in.
Furthermore, presence of clumps in the disc may in fact change the fragmentation properties of the disc since the clump induces additional perturbations (see again Meru, 2013, , where discs were found to fragment at surprisingly small radii in the presence of a pre-existent gas clump at larger radii). Nayakshin & Cha (2013) focused on migration of a single fragment in a marginally self-gravitationally stable disc, thus placing the fragment in a massive but stable disc that would not fragment on its own. They found that in the presence of the clump the disc became more unstable and formed additional fragments. This shows that it is possible to hatch new gas clumps in discs with if there are already fragments born earlier.
Clearly, when survival of the fragments is rare, we should not under-estimate the importance of initial conditions giving the clumps the best chance of surviving. Therefore, our population synthesis model aims to account for this by (a) slowing down the type I migration below its nominal rate by some factor larger than unity (see below); (b) allowing the disc mass to be below the marginally gravitationally unstable configuration when the fragment is born, and (c) removing the disc instantaneously after a time randomly sampled between 0.5 and 5 Million years. These choices are certainly not unique but we feel they are reasonable. An extremely rapid external photo-evaporation of the disc by a close-by massive star (Clarke, 2007) or the presence of a secondary at beyond 100 AU that is later removed by star-star interactions may be important for such a rare population of planets as the gas giants.
The disc exchanges angular momentum with the planet via type I and/or type II migration torques. The only difference from paper I is that the type I migration time scale, , is now given by
The factor , used in paper I on the right hand side of this equation, has been dropped, as it is argued that the saturation of the type I migration rate in the limit is automatically taken into account by passing the torque of the planet to the 1D viscous disc evolution code. We find that when this torque is large then a gap may start opening, thus reducing the torque.
In equation (4) is a free parameter that accounts for a slower migration rate of the planet in a non self-gravitating disc compared to the isothermal type I migration rate. Since our discs are sampling the post self-gravitating phase of the disc evolution by design (see below), we may expect type I migration to be slower than that for discs (Baruteau et al., 2011), for which . In the models below is randomly sampled between and 10.
2.3 Planet’s internal evolution
Planet’s internal evolution is calculated as in paper I, utilising the ”follow the adiabats” approximation to convective/radiative cooling of the planet. The planet is initialised as a polytropic sphere of gas of homogeneous composition, with metallicity, , equal to that of the parent disc, with a given central temperature (see §4). The planet cools by emission of radiation assuming interstellar dust opacity (Zhu et al., 2009) multiplied by the factor to account for grain growth and also by the metallicity of the planet in Solar units, that is, by . Since the planet is surrounded by the disc, its irradiation by the disc (or by the star if the planet is directly exposed to it in the final stages of calculation) reduces the rate at which the planet cools, creating a thermal bath effect (Cameron et al., 1982; Vazan & Helled, 2012).
The planet accretes pebbles from the surrounding disc at the Hill’s accretion rate (Lambrechts & Johansen, 2012)
where is the planet’s Hills radius, , and is the grain surface density at radius . The pebble mass fraction, , is a free parameter independent of the planet’s location, fixed for a given run (see §4 below).
Note that pebble accretion rate depends directly on the surface density of gas around the planet. If the planet manages to open a deep gas in the disc, drops to very low values and hence dives to negligible levels. Pebble accretion is also turned off if and when the fragment undergoes the second (Hydrogen molecule dissociation) collapse. It is argued that after the collapse the planet is so compact that it should be able to accrete gas from the disc as well, so that both gas and pebbles would be accreted. We currently do not include gas accretion onto post collapse planets.
Pebbles entering the fragment have a fixed initial size of cm for this paper, and are deposited in the outer several mass grids of the planet, where they are mixed homogeneously with the grains already in those regions. The grains are allowed to grow by sticking collisions and sediment into the centre of the fragment into the core, as described in paper I. Turbulence within the fragment, parameterised by a coefficient , and convection, however, tend to counteract grain sedimentation. We usually find that grain size must increase to about cm before they sediment efficiently. A further significant obstacle to a rapid grain sedimentation into the core is the fact that grains fragment in collisions rather than stick if collision velocity is too high, as shown by laboratory experiments on dust growth (e.g., Blum & Münch, 1993; Blum & Wurm, 2008; Beitz et al., 2011). As explained in paper I, section 4.3, grain growth turns into grain fragmentation if grain sedimentation velocity with respect to the surrounding gas exceed , set in this paper to 3 m s. Grains are also vaporised if the surrouding gas temperature is too high for the given grain species.
Three grain species are included in the model: water ice, CHON, and silicates (”rocks”). The relative abundances of the species are 0.5, 0.25 and 0.25, respectively. The grain growth, fragmentation, vaporisation, and sedimentation equations are solved for each species separately since these three species have very different material and thermal properties.
Grains reaching the centre are allowed to accrete onto the solid core there, whose initial mass is set to a very ”small” value (). Growing core is expected to radiate some of its gravitational potential energy away, but a self-concistent modelling of energy transfer within the core is fraught with many physical uncertainties (e.g., Stamenković et al., 2012). As in paper I (section 4.4), we parameterise the energy release by the core via the Kelvin-Helmholtz contraction time, , of the solid core, which is fixed for all the runs here at years. The luminosity released by the core is an internal energy source for the contracting fragment, analogous to that of a stellar core, and may even exceed the luminosity of the whole planet, causing its expansion.
3 Example runs
Before the grid of models is discussed, we present two example planet formation tracks that illustrate the sequence of events in the TD model for planet formation.
3.1 A gas giant planet
Figure 2 shows an example selected from the grid of runs in which the end result is a gas giant planet of mass at separation of AU. The initial mass of the planet is , so all the extra mass in the end of the run comes from accretion of pebbles from the disc. The initial mass of the disc is , planet location AU, and the initial planet’s central temperature K. The disc metallicity is Solar, that is, , the pebble mass fraction is , disc viscosity parameter . The disc removal time is Million years, but in this simulation the disc dissipates earlier due to photo-evaporation, as we shall see below. Further, the opacity reduction factor , the planet’s type I migration speed is moderated by factor , which is on a high side for the grid of models, so this run presents an example of a rather slowly migrating clump.
Fig. 2a presents the disc surface density profile at five different times. The initial condition at is shown with the black solid curve, all the other curves are marked by their respective times in the legend. The crosses on the bottom of the panel show the position of the planet at times colour-coded in the same way as the surface density curves.
Next panel, fig. 2b, reports the radial location of the planet (black solid curve), the planet’s Hill (red dotted) and its actual (, blue dashed) radii, all ploted versus time. Boxed text and the corresponding arrows mark the times of three important transitions in this run. In particular, the fragment is initially molecular Hydrogen dominated, collapses at time Million years, at which point Hydrogen becomes atomic and partially ionized. During collapse the fragment contracts to the radius of just a few times that of Jupiter, and then continues to contract further due to a rapid radiative cooling. This fragment is never challenged by tidal forces of the star as the planet’s Hill radius (red dotted curve) is always much larger than .
The next notable transition in the planet-disc system occurs at Million years, when the planet opens a deep gap in the disc. The planet is located at AU at that time. From that time on, the planet migrates in type II regime. Note that at time Million years (the dashed blue curve in fig. 2a), the gap becomes deeper than it was before, and the inner disc has a significant gas deficit at AU. This is not only due to the gap openning; the photo-evaporation of the disc is especially strong at a few AU, and is negligible within 1 AU, as can be seen from fig. 1. The planet-facing side of the disc inside the planet’s orbit thus loose mass rapidly due to photo-evaporation.
At this point in time, the inner disc is drained roughly equally efficiently by both accretion onto the star and photo-evaporation. By time Million years, there is a complete gap between the planet and the star. The outer disc is then evaporated too, and the planet stops migrating completely.
Finally, fig. 2c shows the core mass assembled within the fragment, and its decomposition on the silicates (“rocks”) plus water ice and CHON. Since the fragment heats up quickly when pebble accretion commences, water ice grains vaporise before they could sediment into the core, hence they do not contribute to the core at all. CHON grains (green dashed curve) sediment while the centre of the fragment is cooler than K, and then only the rocks are refractory enough to continue to accrete onto the core. This particular planet did not manage to assemble a massive solid core before it collapses (when all grains not yet locked into the core are vaporised); the final core mass is only .
The final metallicity of the planet is , about ten times that of the initial (Solar) disc metallicity. This example thus shows how a metal-rich gas giant planet with a low mass core can be assembled in the context of the TD hypothesis.
3.2 A hot Super-Earth planet
Figure 3 demonstrates how a hot Super-Earth planet can be assembled in the framework of TD. The initial mass of the planet is again . The initial mass of the disc is , planet’s location is AU, and the initial planet’s central temperature K. The disc metallicity is times Solar, the pebble mass fraction is , disc viscosity parameter . The disc removal time is Million years, but again disc dissipates earlier due to photo-evaporation. The planet’s type I migration speed is barely reduced from the isothermal result, e.g., factor , so the fragment migrates in much sooner than it does in §3.1. The maximum grain sedimentation velocity for this simulation is set to m s, which, combined with a higher value, encourages a faster solid core growth.
It appears that the much more rapid migration is to blame for most of the differences in the fate of the fragment here as compared to the much more slowly migrating one in §3.1. Fig. 3 shows that the fragment finds itself in the inner disc much quicker than the one in Fig. 2 did. The gap in the disc is openned by the fragment at around AU, at time Million years. Significantly, the fragment is still in the pre-collapsestate at that moment, and it also has a massive solid core, at the time of the gap openning. The more massive core is due to the higher value in this simulation.
Contraction of the fragment in this simulation occurs almost entirely due to pebble accretion because radiative cooling is very inefficient. This is because the fragment’s dust opacity is high, especially so due to pebble accretion: by Million years the fragment metallicity is . The fragment is also basked in the radiation field of the disc (which at the outer disc radii is dominated by irradiation from the star). Now, when pebble accretion onto the fragment stops, its contraction stops as well. However, since the fragment is metal-rich, the core of the fragment continues to grow in mass rapidly and its luminosity becomes sufficiently high to start puffing the fragment up. This increase in the fragment’s radis after the disc gap openning is easy to spot in Fig. 3b. Since the fragment continues to migrate in, now in type II regime, the planet’s Hills radius continues to shrink, whilst the planet’s radius, , increases with time. A tidal disruption is thus unavoidable, and indeed occurs when the fragment reaches AU.
The core’s mass is then equal to . As the fragment is disrupted, its mass plumets to (in this paper we neglect a possible tightly bound dense gas “atmosphere” around the core, which is not likely to be massive as our cores are quite bright). With this reduced mass, the planet is no longer able to maintain the deep gap in the disc, which is swiftly closed (see the blue dashed curve in Fig. 3a). The mass of the core is however sufficiently high for it to continue its migration in type I now, on a longer time scale yet fast enough to arrive at AU by the time the disc is dispersed at Million years.
The composition of the core is almost exclusively rock, with CHON and especially water ice grains unable to sediment because the core heats up to temperatures above K quickly (cf. paper I for more detail). The core’s composition is significantly different from that predicted by the Core Accretion paradigm, in which ices are very important (Pollack et al., 1996; Alibert et al., 2005) for assembly of cores as massive as this one. This may be a testable difference between the models if composition of Super-Earths is constrained observationally, although the potential presence of a Hydrogen/He atmosphere on top of the core may complicate differentiation between the models.
4 The population synthesis grid of models
A population synthesis involves investigation of the model’s parameter space, which is usually large, by randomly sampling it (e.g., Mordasini et al., 2009) in order to pull out statistical trends of the model that may not be distinguisheable from just a handful of runs. In following this approach, we fix some parameters of the model in order to reduce the parameter space to a manageable but physically meaningful minimum.
Table 1 lists these fixed parameters, their meanings and values. The latter are chosen to be “reasonable” based on experiment (e.g., ) or on what is commongly used in literature (e.g., ). While these values are not unique, we found by varying them that none of the parameters in Table 1 change the metallicity trends that are the focus of this paper significantly, although quantitatively there are of course changes in the results. For example, the normalisation of the giant planet number per star surviving in the end of the simulations does depend on how the disc is removed in the end of the simulations, and hence on and , but the metallicity trends of the survived planets are the same.
One parameter in Table 1 that is difficult to constrain from first principles and which does influences the results is , the grain opacity reduction due to grain growth. The latter is fixed here on a rather large value, , which essentially precludes a significant opacity reduction. This is done following the arguments presented in paper I: the observed metallicity trends are incompatible with low dust opacity models and presumably imply that dust fragmentation (and not only dust growth) is occuring within the fragments to keep the opacity at relatively large values.
Table 2 lists the set of variable parameters, and their minimum and maximum values. The parameters for a run are then randomly picked in the logarithmic space between these respective minimum and maximum values. For example, the disc viscosity parameter, , is very poorly known for protoplanetary discs (and any other astrophysical discs). We introduce a random variable uniformly distributed between and 1, draw a random value , and define for a given simulation by
Similarly, random uniformly distributed variables are generated for all the other entries in Table 2.
Having defined the parameters for a simulation, we model one gas fragment per disc at a time. The fragment starts with a given initial total (H/He + grains) mass, and has the same metallicity as that of the disc. The initial temperature is randomly picked (as described in equation 6 above), between , given by
and twice . Here the mass dependence of reflects the fact that more massive fragments cool more rapidly (e.g., Helled et al., 2008; Nayakshin, 2010b), and are also hotter at formation due to the adiabatic compression. The choice made in equation 7 does not pre-determines or drives the main results of this paper.
Finally, we wish to investigate how our results depend on the metallicity of the host star/disc and the initial mass of the fragment. Therefore, the random models are repeated for a grid of the metallicity and the initial fragment’s mass values. The metallicity grid is given by , with , … 5. The initial planet mass grid is given by , 0.7, 1 and 2. For each of the combinations of and from this grid, e.g., and , the total of 243 models are run. This gives the total of 4860 runs. A typical run takes about 2 hours of physical time on a single CPU, although some runs (usually with higher disc mass and viscosity, as that leads to smaller time steps) had to be executed for up to 48 hours.
Figure 4 presents the results of the 4860 runs in the planet mass versus planet-star separation plane. In order to facilitate a visual comparison to the exoplanetary data, giant planet masses are shown as , where is inclination of the system to the observer, generated randomly assuming a uniform distribution for . This is done because many of the known giant planets were detected by the radial velocity method and so also have the unknown factor. The gas giants are defined here as fragments that did not have a tidal disruption and hence are all more massive than our minimum starting fragment mass, e.g., . All of these planets are above the horizontal dashed line.
Cores, the remnants of tidal disruption of the initial gas fragments, are the planets below the dashed line. In Fig. 4, all of the survived giant planets are shown, but, for the sake of clarity, only a third of the cores, randomly selected from the full sample, is plotted. The mass of the cores, the scale for which is given in units of on the right vertical axis, is not multiplied by the factor since many of the observed planets of this mass are detected by the transit method which is not affected by the inclination uncertainty.
We emphasise that in this paper no consideration for the presence of a bound atmosphere around the cores after the disruption is made, so their mass may be somewhat under-estimated. Further, gas accretion from the surrounding protoplanetary disc onto the cores and onto the post-collapse giant gas planets may well occur but is also not taken into account, thus again possibly under-estimating the final mass of the planets. Finally, an additional uncertainty which may bring down the mass of the giant gas planets is that here we assumed that H dissociation collapse of a gas giant leads to the instantaneous collapse of the whole fragment. However, 3D hydrodynamical simulations of embedded clumps in protoplanetary discs show that only a fraction of the total mass of the fragment collapses into the post-collapse planet immediately, with the rest (up to % of the fragment’s mass Galvagni et al., 2012) collapsing into a circum-planetary disc. It is not clear whether all of this disc or only a fraction is eventually accreted onto the planet. There is thus an uncertainty of a factor of in the mass of the giant planets shown in Fig. 4.
Despite these uncertainties, there are clear and strong metallicity trends in the planet populations in Fig. 4 which we shall investigate in detail shortly. To aid this, the symbols in Fig. 4 are colour-coded to show the metallicity of the host star/disc. For example, the red colour shows the most metal rich systems, , whereas the green show the most metal-poor systems in the sample, . It is easy to notice that most of the giant planets in the inner few AU of the diagram are metal rich; this is not so for cores, however. This observation essentially summarises the main findings of this paper.
Finally, note the vertical column of giant planets at AU. These are the fragments that managed to collapse before being tidally disrupted but that continued to migrate inward rapidly and reached the inner boundary of our computational domain (we stop the calculation if the planet reaches a short distance away from AU). This column contains 2391 planet, e.g., almost a half of all fragments born in the disc, managed to collapse before being tidally disrupted but reached the innermost disc radius. Most of these would probably be completely devoured by the star if our computational domain were extended all the way to the inner boundary.
An additional point to make on this is that our calculation is optimised for survival of the giant planets because we essentially simulate the end phase of the disc evolution. There is no reason why fragments could not be born before the phase we are simulating. These would migrate even faster, and most likely would all end up either being tidally disrupted or being driven all the way to the inner boundary. Hence the number of giant planets migrating all the way into the star is still an under-estimate. This is consistent with the view presented earlier – that destroyed gas giant planets cause FU Ori outbursts and that statistically speaking, there may be many such events per protostar.
6 Metallicity correlations
6.1 Small separations
Figure 5 shows the number of planets formed in our runs at separations AU AU (this excludes the planets that reached the innermost disc radius, e.g., the column of planets on the left-hand side of Fig. 4) versus host metallicity in Solar units, , for three groups of planets: the gas giants, the Super-Earths and the “Earths”. The division on the three groups is made based purely on the planet’s mass. Gas giants are fragments that collapsed rather than were disrupted, so they are all more massive than . Super-Earths are solid cores more massive than , and Earths are planets with mass .
There is a clear and strong correlation with for gas giant planets, although somewhat less strong than found in Nayakshin (2015b). As explained in the latter publication and in Nayakshin (2015a), physically this correlation arises from gas giants collapsing faster at higher abundance of pebbles in the outer disc, since the fragments then accrete pebbles at higher rates. A larger fraction of gas clumps thus collapses before it is tidally disrupted.
The positive correlation with in Fig. 5 saturates at the highest metallicity bin. It is worth emphasising that we assumed here that the fragment birth locations do not depend on metallicity of the disc for simplicity. Higher mean higher dust opacity, and hence longer cooling times in the disc. One may expect the clumps to be born farther away on average at higher due to higher disc opacity (such effect is actually found in simulations by Meru & Bate, 2010). The farther away the clumps are born, the longer it takes for them to migrate in, and hence the more likely they are to survive. This secondary effect, not included in the present simulations, may strengthen the giant planet frequency versus metallicity trend of the model.
On the other hand, Fig. 5 shows, surprisingly at the first glance, no clear correlation between the host metallicity and the populations of Earths or Super-Earths. Naively speaking, one may expect that a higher abundance of metals within the fragments at higher metallicities would cause more massive cores to be built. When some of the fragments are disrupted, these cores should become visible to the observer, so this would predict more massive cores around stars of higher metallicity.
This simple argument does not actually work because the fragment’s properties also vary with metallicity and that should be also taken into account. In particular, a fragment loaded with pebbles at a higher rate indeed has a higher metallicity, so that the core’s accretion rate is higher, but the time window for grains to sediment down into the core is shorter because the fragment heats up too quickly and the grains are eventually vaporised which clearly stops accretion of the core.
This physics is internal to the fragment and has nothing to do with the fragment-disc interaction. Therefore, to investigate it in greater detail, we perform a series of runs in which we turn off fragment’s migration and keep the fragment’s pebble accretion rate constant, parameterised, as in Nayakshin (2015a), in terms of the ”metal loading time”, , defined by
where is the Solar metallicity. Thus is the time scale on which the planet’s metal content increases by the amount contained in that planet at Solar metallicity. Since the fragment does not migrate it is not tidally disrupted in these isolated runs.
Figure 6 shows the results of this series of runs for three initial planet masses, , 1 and . The top panel shows the mass of the core assembled inside the fragment by the time it collapses versus the metal loading time on the horizontal axis. For , the core mass is essentially independent of . The bottom panel of Fig. 6 explains the reason for this finding to some degree. The panel shows the time it took the fragment to collapse. Unsurprisingly, the faster the metals are added to the planet, the faster it collapses, so that the core indeed has a shorter time window to grow. For the two trends – the faster core growth and the shorter time available for that – nearly cancel each other out, yielding a nearly constant result. For , the dependence is mostly flat, but at the longest the core mass actually decreases. This is correlated with the time it takes for the fragment to collapse (the bottom panel of Fig. 6). Physically, for , the longest planets collapse faster because their contraction is dominated by radiative cooling rather than pebble accretion. Radiative contraction is faster at lower metallicities (Helled & Bodenheimer, 2011; Nayakshin, 2015a), and hence the result. Finally, for the lowest mass planets in Fig. 6, the faster the metals are loaded into the planets (the shorter ), the lower the mass of the cores assembled.
The different behaviour shown by the planets of different initial mass in their response to a variance in is due to a fact that a number of factors, rather than just one, control the core’s assembly rate inside the planet. The most important of these factors are: the central density of the planet, its metal abundance, gas temperature, and the core’s luminosity as this affects the convective flux and the convective grain mixing which opposes grain sedimentation.
The weak dependence of the core’s mass on the metal loading time (and hence on the metallicity of the disc), shown in Fig. 6, explains why the populations of the rocky planets in Fig. 5 are broadly distributed over metallicities. Indeed, if a roughly constant fraction of precollapse fragments behaving like those in the isolated simulations used for 6 is tidally disrupted, then there would be cores of similar masses at a range of metallicities.
However, note that there is a precipitous fall in the number of rocky planets at the highest metallicities in Fig. 5 that would not be explained in this picture. We believe that this decline in the number of rocky cores with at high has to do with a changing fraction of precollapse fragments experiencing tidal disruptions. The higher the metallicity of the disc, the less likely the fragment is to be disrupted. This can be seen from the bottom panel of Fig. 6. The fragment collapse times are quite short ( years) at the low end. The fragment migration times are longer than that for most of our runs in the grid of models, and hence most of these fragments are expected to survive the tidal disruptions.
This is why at the highest metallicities most of our precollapse fragments manage to collapse and either survive as a gas giant at the end of the simulation or reach the inner boundary of the disc. Either way, their rocky cores are never released from under the massive gas envelope, and this is why there are so few rocky planets at the highest metallicities in Fig. 5. This fall is probably too strong compared with observations. We come back to the issue in the Discussion section.
6.2 Large separations
Figure 7 shows the metallicity correlations for planets survived the disc migration phase at the outer AU region of the star. For the Jovian mass planets in the figure, the planets that are stranded far from the host star are those that migrate the slowest. Therefore, these fragments are not likely to be threatened by the tidal disruption, so can take as long as they need to contract. This is why the population of gas giants is now broadly distributed over the host star metallicity, with just a shallow minimum at around the Solar metallicity. This is quite different from Fig. 5 where a strong positive corelation with metallicity was present.
The Super Earth and the Earth-like planets in Fig. 7 are the products of tidal disruptions of precollapse gas fragments, so these tend to come from those fragments that were migrating faster (physically ”released to migrate” earlier, when their discs were more massive) than the gas giant planets in the figure. The metal correlation of the population of the Super-Earth planets in the outer 10 AU is not actually that different from that in Fig. 5, with similar physics driving it.
Earth-like planets are however much less abundant at high separations from the host, and especially at the highest metallicity bins. This is due to the fact that very few high metallicity fragments are disrupted in the outer disc. When they are disrupted, they tend to contain more massive cores because these cores had more time to grow by grain sedimentation than cores disrupted in the inner disc (which would usually migrate in much faster). This is another interesting example where predictions of TD do not follow the simple ”more metals more cores” logic. Hopefully this prediction may be observationally tested, although these relatively low mass planets are not likely to be found by direct imaging surveys in near future due to their exceedingly low luminosities.
In this paper we argued for interpreting the metallicity correlations of observed exoplanets in the framework of the Tidal Downsizing hypothesis for planet formation rather than in the context of the Core Accretion model. In the latter theory, the fact that massive cores form at all metallicities but gas giants prefer metal-rich environments is odd. Once a massive core is assembled, accreting a massive gas envelope on the top of that core is actually easier in metal-poor environments. Lower opacities in the envelope allow the heat of the core’s assembly to escape from the gas envelope faster, accelerating its contraction and the eventual collapse. This can be seen from the well known result that the critical core mass for atmosphere collapse in the CA theory increases (although weakly) with opacity of the envelope (e.g., Stevenson, 1982; Ikoma et al., 2000). So, if massive cores are equally abundant at all metallicities, one would then expect more gas giant planets at low metallicities than at high metallicities, which is observationally not at all the case.
In contrast, formation of gas giants in the TD hypothesis is not preconditioned by the assembly of solid cores. As we found here, there is actually no simple relationship between these two types of planets in general. Both the rocky core dominated planets and the gas giants originate from the same source – the precollapse H dominated fragments born by the gravitational instability of the outer massive disc – but the efficiency with which the different kinds of planets are assembled from this raw material is a function of a number of factors.
Our numerical population synthesis models allowed us to study these issues in some detail in this paper. Focusing first on the planets surviving the disc migration phase in the inner 5 AU of the disc (§6.1), we found that increasing metallicity of the environment leads to a more abundant formation of gas giant planets, as earlier found by Nayakshin (2015b) for coreless gas fragments. This result is driven entirely by the pebble accretion accelerating collapse of the fragments at high metallicites and helping more of them to survive.
However, formation of solid cores in our models does not correlate in a monotonic way with the metal content of the parent discs for either Earth-like or Super Earth planets. There are two primary reasons for this outcome. Firstly, higher pebble accretion rates on a gas fragment do not necessarily make more massive cores. The core’s mass is equal to the integral of the core’s accretion rate over the duration of time period when the core can grow, e.g., . The core’s accretion rate, , increases with increasing pebble accretion rate, but the duration of the core assembly phase, , decreases. For these reasons the core’s mass turns out to not depend sensitively on the pebble accretion rate (which is proxy for the disc metallicity), as could be naively expected (see Fig. 6).
Secondly, in accord with the explanation offered for the observed gas giant planet correlation with the metal abundance in TD, fewer precollapse fragments are disrupted by the stellar tidal forces at higher metallicities. Thus, not only the cores are not necessarily much more massive at higher metallicities, but they are also less likely to be released from inside the overlying massive gas envelopes. This is the reason why the population of the rocky cores nose-dives at the highest metallicity bin in Fig. 5. This particular result is however preliminary and may weaken if we are to consider earlier generation of gas fragments that may migrate more rapidly, so that a sizable fraction of precollapse gas fragments is tidally disrupted even at the highest metal abundances (which is not currently the case in our models). In that case there will probably be more rocky planets at the high metallicity end in Fig. 5.
Finally, in §6.2 we compared the metallicity correlations of the simulated planets in the outer disc ( AU) to that found in the inner disc. Our models predict that gas giant planets of all metallicites may survive at these ”cold” regions of the disc, since tidal disruptions of precollapse fragments are far less likely there. One interesting coincidence here is the fact that the best well known system of giant gas exoplanets probably formed by GI of the disc, HR 8799 (e.g., Marois et al., 2008), is metal poor ([Fe/H] ).
We believe the observed metallicity correlations of planets as a function of their size/radius give us a valuable clue as to how planets form. It is highly significant that gas giant planets are found almost always around metal-rich stars (Gonzalez, 1999; Fischer & Valenti, 2005) but planets smaller than or less massive than Neptune are abundant around stars of any metallicity (Sousa et al., 2008; Buchhave et al., 2012).
If solid cores are assembled first as a prerequisite to gas giant planet formation, then it is hard to see why there are massive cores in low metallicity environments but no gas giants since the lower the opacity of the envelope the easier it is to convert a massive core into a gas giant.
On the other hand, if gas fragments are the nurseries of massive solid cores (Helled & Schubert, 2008; Nayakshin, 2011), and if the fragments are preferentially destroyed in metal-poor environments (Nayakshin, 2015a), then the observed metallicity correlations are to be expected as long as the fragments are able to hatch massive cores before they are tidally disrupted. The population synthesis models presented here appear to be consistent with the observations. They also predict that gas giant planets and Super Earths may be abundant at large separations from the star at all metallicities.
Theoretical astrophysics research in Leicester is supported by an STFC grant. This paper used the DiRAC Complexity system, operated by the University of Leicester, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by a BIS National E-Infrastructure capital grant ST/K000373/1 and DiRAC Operations grant ST/K0003259/1. DiRAC is part of the UK National E-Infrastructure.
- pagerange: Tidal Downsizing model. II. Planet-metallicity correlations–References
- pubyear: 2008
- This is in difference to accretion of planetesimals in the CA framework, in which planetesimals impact the growing planet at a high (up to a few km s) velocity and hence may actually heat it. In contrast, pebbles sediment onto the fragment at a moderate a few m/sec velocity, or else they fragment. Their kinetic energy input into the planet is hence negligible.
- Alexander R. D., Armitage P. J., 2007, MNRAS, 375, 500
- Alibert Y., Mordasini C., Benz W., Winisdoerffer C., 2005, A&A, 434, 343
- Baruteau C., Meru F., Paardekooper S.-J., 2011, MNRAS, 416, 1971
- Basu S., Vorobyov E. I., 2012, ApJ, 750, 30
- Batalha N. M., Rowe J. F., Bryson S. T., et al., 2013, ApJS, 204, 24
- Beitz E., Güttler C., Blum J., Meisner T., Teiser J., Wurm G., 2011, ApJ, 736, 34
- Blum J., Münch M., 1993, ICARUS, 106, 151
- Blum J., Wurm G., 2008, ARA&A, 46, 21
- Bodenheimer P., 1974, Icarus, 23, 319
- Bodenheimer P., Grossman A. S., Decampli W. M., Marcy G., Pollack J. B., 1980, ICARUS, 41, 293
- Boley A. C., Hayfield T., Mayer L., Durisen R. H., 2010, Icarus, 207, 509
- Boss A. P., 1997, Science, 276, 1836
- Boss A. P., 1998, ApJ, 503, 923
- Buchhave L. A., Latham D. W., Johansen A., et al., 2012, Nature, 486, 375
- Cameron A. G. W., Decampli W. M., Bodenheimer P., 1982, Icarus, 49, 298
- Cha S.-H., Nayakshin S., 2011, MNRAS, 415, 3319
- Clarke C. J., 2007, MNRAS, 376, 1350
- Fischer D. A., Valenti J., 2005, ApJ, 622, 1102
- Forgan D., Rice K., 2013, MNRAS, 432, 3168
- Galvagni M., Hayfield T., Boley A., Mayer L., Roškar R., Saha P., 2012, MNRAS, 427, 1725
- Gammie C. F., 2001, ApJ, 553, 174
- Gonzalez G., 1999, MNRAS, 308, 447
- Hartmann L., Kenyon S. J., 1996, ARA&A, 34, 207
- Helled R., Bodenheimer P., 2011, ICARUS, 211, 939
- Helled R., Bodenheimer P., Podolak M., et al., 2013, ArXiv e-prints
- Helled R., Podolak M., Kovetz A., 2008, Icarus, 195, 863
- Helled R., Schubert G., 2008, Icarus, 198, 156
- Ida S., Lin D. N. C., 2004, ApJ, 616, 567
- Ida S., Lin D. N. C., 2008, ApJ, 685, 584
- Ikoma M., Nakazawa K., Emori H., 2000, ApJ, 537, 1013
- Johansen A., Lacerda P., 2010, MNRAS, 404, 475
- Kuiper G. P., 1951, in 50th Anniversary of the Yerkes Observatory and Half a Century of Progress in Astrophysics, edited by J. A. Hynek, 357–+
- Lambrechts M., Johansen A., 2012, A&A, 544, A32
- Larson R. B., 1969, MNRAS, 145, 271
- Lin D. N. C., Bodenheimer P., Richardson D. C., 1996, Nature, 380, 606
- Lin D. N. C., Pringle J. E., 1987, MNRAS, 225, 607
- Marois C., Macintosh B., Barman T., et al., 2008, Science, 322, 1348
- McCrea W. H., Williams I. P., 1965, Royal Society of London Proceedings Series A, 287, 143
- Meru F., 2013, in European Physical Journal Web of Conferences, vol. 46 of European Physical Journal Web of Conferences, 7003
- Meru F., Bate M. R., 2010, MNRAS, 406, 2279
- Michael S., Durisen R. H., Boley A. C., 2011, ApJL, 737, L42+
- Mordasini C., Alibert Y., Benz W., 2009, A&A, 501, 1139
- Nayakshin S., 2010a, MNRAS, 408, L36
- Nayakshin S., 2010b, MNRAS, 408, 2381
- Nayakshin S., 2011, MNRAS, 413, 1462
- Nayakshin S., 2015a, MNRAS, 446, 459
- Nayakshin S., 2015b, MNRAS, 448, L25
- Nayakshin S., 2015c, ArXiv e-prints (1411.5264)
- Nayakshin S., Cha S.-H., 2013, MNRAS, 435, 2099
- Nayakshin S., Lodato G., 2012, MNRAS, 426, 70
- Ormel C. W., Klahr H. H., 2010, A&A, 520, A43
- Owen J. E., Clarke C. J., Ercolano B., 2012, MNRAS, 422, 1880
- Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M., Greenzweig Y., 1996, Icarus, 124, 62
- Rafikov R. R., 2005, ApJL, 621, L69
- Rice W. K. M., Lodato G., Armitage P. J., 2005, MNRAS, 364, L56
- Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
- Sousa S. G., Santos N. C., Mayor M., et al., 2008, A&A, 487, 373
- Stamenković V., Noack L., Breuer D., Spohn T., 2012, ApJ, 748, 41
- Stevenson D. J., 1982, P&SS, 30, 755
- Vazan A., Helled R., 2012, ApJ, 756, 90
- Vorobyov E. I., 2013, A&A, 552, A129
- Vorobyov E. I., Basu S., 2005, ApJL, 633, L137
- Vorobyov E. I., Basu S., 2006, ApJ, 650, 956
- Williams I., Crampin D. J., 1971, MNRAS, 152, 261
- Winn J. N., Fabrycky D. C., 2014, ArXiv e-prints
- Zhu Z., Hartmann L., Gammie C., 2009, ApJ, 694, 1045