© ESO, 2014

1. Introduction

It is well known (e.g., Ikoma et al. 2000) that a reduction of the opacity κ in the gaseous envelope of a forming giant planet leads to a reduction of the formation timescale. In simulations of concurrent core and envelope accretion (e.g., Pollack et al. 1996, hereafter P96) the reason is that a lower opacity during the so-called phase II (where gas must be accreted in order to allow further core growth) leads to a more efficient transport of released potential energy out of the envelope. Phase II is an intermediate phase that occurs for in situ formation between the moment the core has reached the isolation mass, and the moment when rapid gas accretion starts. This happens when the core has reached the crossover mass when envelope and core mass are equal.

For in situ calculations the overall formation timescale of the planet is dominated by the duration τ_II of phase II, which means that at low κ the overall formation timescale is also reduced. For example, P96 find that an (arbitrary) reduction of the grain opacity κ_gr to 2% of the interstellar medium (ISM) value leads to τ_II = 2.2 Myr for Jupiter formation, while with the full ISM value, τ_II = 6.97 Myr, longer than the mean disk lifetime.

Obviously, one is therefore interested in knowing an estimate for the effective grain opacity in the envelope, instead of having to use arbitrary scalings like the 2%. The grain opacity (by opacity we always mean in this work the Rosseland mean) can be calculated from the microphysics of grain growth via coagulation, grain settling, and evaporation at high temperatures. Podolak (2003) presented such a numerical model, finding that grain growth leads to opacities up to three orders of magnitude smaller than in the ISM because grains grow efficiently. The limitation of this work was a lack of self-consistency as the envelope structure in which the grain growth was studied was calculated with pre-specified different opacities.

Only recently Movshovitz et al. (2010, hereafter MBPL10), presented the first self-consistently coupled calculations of grain evolution and giant planet growth. In their work, the envelope structure is used to calculate at each radius the evolution of the grains at each time step. From the grain size distribution the Rosseland mean opacity is calculated, which is then fed back into the envelope calculation. The main result of these calculations is that the grain opacity is greatly reduced. For Jupiter formation, this leads to a duration of phase II of only 0.52 Myr. This corresponds to a reduction by a factor of 13.4 relative to the P96 full opacity case, and still a factor of 4.2 to the P96 “low opacity case” with a 2% ISM opacity. As forming a giant planet within the typical lifetime is a timing issue (at least if migration is not taken into account; see Alibert et al. 2004), these factors matter.

The complex (and computationally heavy) calculations of grain evolution made by MBPL10 are not the scope of this first paper. Instead, here we follow a practical approach, and determine in the first part of this work the reduction factor f_opa by which interstellar opacities must be reduced in order to obtain the same duration of phase II found by MBPL10. As will be shown, a grain opacity of only about 0.3% the ISM value leads to the best reproduction of the MBPL10 results. Using one uniform reduction factor of the ISM opacity can of course not reproduce the complex structure of the opacity which depends on planetary properties like the core or envelope mass as found in grain evolution models (Movshovitz & Podolak 2008). The simple ISM scaling approach therefore has some important limitations (see the discussion in Sect. 3.5). Our interest in still deriving a global, uniform f_opa is to have an intermediate value for the opacity between the two extremes (full ISM opacity vs. grain free) at a low computational cost for population synthesis simulations. The goal of this work is rather to study the global effects of κ on planetary populations. For this, we compare in the second part of the paper important statistical properties of synthetic planets found with different values of f_opa with observational constraints from extrasolar planets.

The opacity controls the rate at which a core of given mass accretes gas during the formation phase, therefore different f_opa lead to different bulk compositions that manifest in the mass–radius relationship that we can observe after a few Gyr of evolution. Population synthesis calculations that are built on self-consistently coupled planet formation (Alibert et al. 2005) and evolution models (Mordasini et al. 2012b) can predict the synthetic mass–radius relationship. While there are other observational constraints on f_opa that can be derived from the mass distribution of exoplanets alone (in particular the frequency of giant planets), transiting exoplanets are therefore of special interest in this study. What is needed for the comparison are transiting extrasolar planets with a well defined mass and radius, and a rather large semi-major axis ≳0.1 AU to minimize the impact of stellar irradiation. The extrasolar planets of this class must also have a mass–radius relationship that implies that they contain significant amounts of primordial H/He, since this is the envelope type considered here. It also means that we avoid the part of the mass–radius relationship that is most degenerate (e.g., Valencia et al. 2010; Rogers & Seager 2010). Instead, for such planets it is possible to infer, at least in a rough way, the bulk composition (global heavy element content) as has been demonstrated by Guillot et al. (2006), Burrows et al. (2007), Guillot (2008), or Miller & Fortney (2011). The results from the latter study are used in the second part of this study. Our simulation therefore also aims at establishing a bridge between physical processes like grain evolution that govern gas accretion during formation, and observable quantities. This information can be used in the ideal case to feed back into the microphysical models for the grains (or into specialized models in general).

The structure of the paper is as follows: in Sect. 2 we show the modifications of our giant planet formation model to take into account reduced grain opacities. The formation model itself was first described in Alibert et al. (2005), while the version used here is described in Mordasini et al. (2012b). In Sect. 3, we determine f_opa by comparison with MBPL10. We then turn to the population synthesis calculations (Sect. 4) and use them to study the envelope mass as a function of core mass (Sect. 5). The associated mass–radius relationship mainly of low-mass planets is addressed in Sect. 6. In Sect. 7 we compare the enrichment of giant planets relative to the host star for different f_opa with the results derived by Miller & Fortney (2011) for actual extrasolar planets. Finally, in Sect. 8 we give a summary and present our conclusions. In Appendix A we derive a semi-analytical model for the core and envelope mass in phase II and determine its parameters in Appendix B.

In the second paper of this series (Mordasini 2014, hereafter Paper II), we present a simple analytical model for the opacity due to grains in protoplanetary atmospheres. It calculates the grain opacity based on the comparison of the timescales of the governing microphysical processes like grain settling, coagulation and evaporation (Rossow 1978). It therefore takes into account that the grain opacity is a dynamically changing quantity that depends on planetary properties like core and envelope mass, or the accretion rates. In a subsequent work, we will couple this model with our upgraded population synthesis code that can now simulate the concurrent formation of many embryos in one disk during the formation phase (Alibert et al. 2013) and includes atmospheric escape during the long term evolution (Sheng et al. 2014).

2. Updated opacity model

The opacity of the gas in the envelope is calculated using two different tables: first, combined ISM grain opacity (multiplied by a factor f_opa) and high-temperature gas opacities from the Bell & Lin (1994) analytical opacity laws. Second, low-temperature molecular and atomic opacities of a grain-free, solar composition ([M/H] = 0) gas from Freedman et al. (2008).

For the scaled grain and high-temperature gas opacities we use the expressions of Bell & Lin (1994) that analytically give the opacity as a function of density and temperature for eight different opacity regions, where the transitions between the different regions are smoothed¹. The opacity regions are in increasing temperature order dominated by: ice grains, ice grain evaporation, metal grains, metal grain evaporation, molecules, H-scattering, Kramer’s law and finally electron scattering. This means that at low temperatures, the opacity is assumed to be caused by grains only. While for a full ISM grain opacity this is justified, this is no more the case if the grain opacity is reduced by a large factor. The code was modified so that the opacity in the regions where grains dominate can be scaled freely with a grain opacity reduction factor f_opa ≤ 1 in a way that the opacities remain continuous also with such a factor. The high-temperature gas opacities remain unchanged.

When the grain opacity is very low, gas opacities begin to dominate also at low temperatures. This effect is seen in the calculations of MBPL10. Therefore it is necessary to include opacity tables for a grain-free gas also at low temperatures. Freedman et al. (2008) provide such tables of molecular and atomic opacities of a grain-free solar composition ([M/H] = 0) gas. These tables were also used by MBPL10. For the present study, the results of Freedman et al. (2008) were converted into tables giving κ(T,ρ), using the SCvH EOS (Saumon et al. 1995) for the conversion of pressure P and temperature T into density ρ.

Finally, the (scaled) opacities from Bell & Lin (1994) and Freedman al. (2008) were combined, taking for simplicity always the maximum of the two values. This underestimates the opacities in regimes where both gas and grains contribute in a similar way, but avoids problems in the overlapping regime where both Bell and Freedman opacities are due to gas. This means that even when the grain opacity is very low, the combined opacity will not fall below a lower limit given by the gas.

3. Determination of f_opa by comparison with Movshovitz et al. (2010)

The question we want to answer in this section is simple: what is the ISM grain opacity reduction factor that leads to the same duration of phase II for Jupiter in-situ formation as in the coupled grain evolution and envelope accretion calculations of MBPL10? Only the duration of the phase II is considered because the opacity is the dominant factor for its length².

3.1. MBPL10 comparison mode

For the MBPL10 comparison calculations, we made a number of modifications to our model in order to emulate as closely as possible the formation code in MBPL10. In particular we assume, in contrast to normal operation, that no migration or disk evolution occurs. The other setting, which are if possible identical as in MBPL10, are listed in Table 1. It is clear that the two codes still differ in some aspects. Examples are the planetesimal-envelope interaction calculation, or the fact that we include for the outer boundary conditions the effect that the planet radiates into the nebula (Papaloizou & Terquem 1999). This is however known not to make a large difference (Bodenheimer & Pollack 1986).

3.2. Simulation setup

As in MBPL10, three initial planetesimal surface densities Σ_s,0 are studied: 10, 6, and 4 g/cm². The corresponding durations of phase II τ_II found by MBPL10 are 0.52, 0.79 and 2.43 Myr, respectively. For each of the three Σ_s,0, the best fitting f_opa was determined individually. In the ideal case the three values are identical or at least similar. It is however clear that the grain growth process depends on Σ_s,0. Different Σ_s,0 directly lead to different core masses. This in turn leads to different dust settling velocities (MBPL10) meaning that the timescales of the various processes of grain evolution will differ depending on Σ_s,0. Therefore it is expected that the best fitting f_opa will not be exactly the same for all three cases. The limitations of one general reduction factor are further discussed in Sect. 3.5.

3.2.1. Initial planetesimal surface density Σ_s,0 = 10 g/cm²

To illustrate the impact of κ_gr, we show in Fig. 1 the result for the calculations for Σ_s,0 = 10 g/cm². The left panel shows the total mass as a function of time for f_opa = 10^-4, 0.001, 0.01, 0.1 and 1, while the right panels shows the total luminosity including internal and shock luminosity. The decrease in τ_II with f_opa is obvious. For f_opa = 10^-4 no real phase II exists any more when looking at the mass evolution. The planet passes almost directly from phase I into phase III, even though that still two local maxima can be distinguished. For all cases the final core mass is close to 36 Earth masses (36.1–36.7 M_⊕), and the final total mass is 317.5–319.4 M_⊕ (we artificially ramp down the gas accretion rate on a short timescale once a total mass of 1 M_♃ is approached). The core mass at crossover is similar in all cases, covering a range from 15.6 to 17.9 M_⊕. MBPL10 found a crossover mass of 16.09 Earth masses.

3.3. Duration of phase II as a function of f_opa

Figure 2 shows the duration of the phase II as a function of f_opa for the three surface densities. For Σ_s,0 = 10 g/cm² one sees that the length varies from just about 12 000 years to about 5 Myr. Beginning at the smallest opacities with f_opa ≲ 5 × 10^-5, τ_II is first approximately independent of f_opa. This is because here, the grain opacity is so small that only the gas opacities matter that are independent of f_opa. These cases thus correspond to a grain-free regime that was studied by Hori & Ikoma (2010). For 10^-4 ≤ f_opa ≤ 10^-2, τ_II increases with the grain opacity reduction factor, following the same slope. For f_opa ≥ 10^-2, τ_II still increases, but following another less steep slope.

In Fig. 2, the short black horizontal curve labeled M10 marks the duration of phase II found by MBPL10, 0.517 Myr. The two curves cross slightly below f_opa = 0.005, which is thus the best fitting value. The value is remarkably low. Previous works (P96, Hubickyj et al. 2005; Lissauer et al. 2009) used as (arbitrary) “low opacity case” f_opa = 0.02. The fitting value found here is a factor of four smaller. This difference has a substantial influence on the formation time, as τ_II(f_opa = 0.02) = 1.33 Myr, i.e., almost a factor of three longer.

The dotted red curve in Fig. 2 again shows the duration of phase II, but now for Σ_s,0 = 6 g/cm². For this case, MBPL10 found τ_II = 0.787 Myr. As expected (e.g., P96) the formation timescales are longer because of the smaller core mass at a lower planetesimal surface density. The duration of phase II we find is now between 62 000 years and 64.3 Myr. The dependence of τ_II on f_opa follows a similar pattern as before. The curve meets with the value of MBPL10 at f_opa = 0.002, which is thus the best fitting value. This is again a very low value. The value is smaller than the one for 10 g/cm² by a factor of 2.5, but still on the order of a few times 10^-3. As discussed earlier, we cannot expect to find an identical value for the three simulations.

The dashed blue line in Fig. 2 finally shows the dependence of τ_II on f_opa for the lowest surface density case of Σ_s,0 = 4 g/cm². Here, the formation time with full ISM opacity would be much in excess of 100 Myr. Comparison with the result of MBPL10 (2.428 Myr) leads to a best fitting value for this case of 0.002, i.e., again a few times 10^-3.

The dashed-dotted black lines in the top left and bottom right corner represent lines with linear slope. The comparison of these lines with τ_II(f_opa) shows that over an intermediate domain of f_opa which itself depends on Σ_s,0 (i.e., M_core), there is roughly a linear correlation between the opacity and the duration of phase II. Such a behavior is theoretically expected, see Appendix B.

3.4. Final result for f_opa

Table 2 summarizes the results for the best fitting values f_opa,bestfit. In the table, we also give the mean of the three values. It seems reasonable for a fitting parameter to take simply the mean. This is especially the case as the three best fitting values do not depend in a strictly monotonic way on Σ_s,0. The mean value f_opa,bestfit = 0.003 is thus the nominal value for the population synthesis calculations. The fact that the results for the three Σ_s,0 differ by a few 0.001 indicates that the mean f_opa,bestfit = 0.003 is accurate only to factors of a few for a specific surface density. This is an additional limitation besides the fact that using a constant reduction factor throughout the envelope is itself an oversimplification relative to actual grain evolution models, as discussed in the next section.

3.5. Limitations of a general scaling factor

The numerical grain evolution models of Movshovitz & Podolak (2008) and MBPL10 find a characteristic radial structure of κ_gr in the atmosphere (outer radiative zone) of protoplanets: near the outer boundary of the planet, the opacity is relatively high and comparable to the ISM opacity (κ_gr ~ 1 cm²/g) but it becomes very low (κ_gr ~ 0.001 cm²/g) in the deeper parts near the boundary to the convective zone. This is clearly a consequence of the growth of the grains. Such a radial structure of κ_gr does not result when simply scaling the ISM opacity, which rather leads to a roughly speaking radially constant κ_gr. This is shown in the comparisons in Paper II. On the other hand, it is also found that for f_opa = 0.003, the total optical depth of the outer radiative zone is similar as predicted by the grain evolution models, at least for the cases that were studied and, importantly, that were used to calibrate f_opa. The comparisons with MBPL10 presented here indicate that this is sufficient to lead to a similar gas accretion rate. While getting a physically motivated radial structure of κ_gr is an important task in its own right (and therefore addressed in Paper II), in the context of this first paper, we are only interested in the gas accretion rate, as this eventually determines the planetary bulk composition.

However, even if the optical depth is similar for the cases for which f_opa was calibrated, there is no a priori reason to assume that for planets which differ substantially (e.g., in core or envelope mass) from the calibration cases, the same f_opa is applicable. This is because the microphysical processes that control κ_gr depend on the gravitational acceleration, the density and temperature in the envelope, or the grain accretion rate. This is a major limitation of the ISM scaling approach used in this paper, and means that no definitive conclusions about the value of κ_gr can be derived from the population syntheses presented further down, but only hints that the opacity might be much smaller than in the ISM. It is also a motivation to derive in Paper II an analytical model that dynamically calculates κ_gr based on the planet’s properties. This is because a numerical model as in MBPL10 can currently not be used in population syntheses due to computational time limitations.

In this context it is also important to clarify that the calibrated f_opa is unlikely to apply to a (for example) 1 M_♃ planet due to its much higher mass than the planets used for the calibration. On the other hand, for the subject of this work, the final bulk composition, the value of f_opa is relevant mainly for ~one order of magnitude in core mass, namely from ~1 to ~15 M_⊕ because of the following: for cores with an even lower mass, the envelope is not influencing much the planet’s formation (it is given by the accretion of solids only) due to the tiny envelope mass. Additionally, only planets with M ≥ 2M_⊕ are considered in the syntheses. Cores more massive than ~15 M_⊕ in contrast trigger gas runaway accretion on a timescale that is much shorter than typical disk lifetimes for all reasonable κ_gr (see, e.g., Ikoma et al. 2000). Once their total mass has reached ~40 to 80 M_⊕ (e.g., MBPL10), they pass into the disk limited gas accretion regime. In this phase, the gas accretion rate is no more dependent on the envelope opacity, but controlled by external factors only (gap formation, accretion rate in the disk, etc.). Interestingly, the opacity in the disk limited accretion phase does, however, affect the post-formation entropy and thus luminosity at young ages, see Mordasini (2013), but this is not discussed here.

We have derived f_opa,bestfit based on three simulations of MBPL10 that cover only a small part of the possible parameter space in which planets form in terms of core mass, luminosity, semi-major axis, and boundary conditions. Using the same f_opa,bestfit for all parameters is therefore an oversimplification. Therefore, we use this reduction value in the population synthesis calculation below as the nominal case, but consider it simply as an intermediate value between the extreme cases with f_opa = 0 and 1.

The simulations of MBPL10 cover a range in core mass of a factor of 4 (about 4 to 16 M_⊕ at crossover). The fact that f_opa,bestfit varies over this mass range by a factor of a few indicates that there is as expected a dependency, but potentially not a very strong one. A hint towards such a behavior is that the radial opacity structures in MBPL10 for the different core masses are all, at least to order of magnitude, similar. To investigate this quantitatively, we will couple in a forthcoming work the analytical model of Paper II to our updated population synthesis model. This will also allow to investigate the impact of different outer boundary conditions.

To investigate this further, we have compared the envelope mass as a function of time for different f_opa in our model with another simulation of combined grain evolution and core accretion (run IIa in Rogers et al. 2011). This study uses the same core accretion and grain evolution model as MBPL10, but additionally simulates the formation of a very low-mass planet at an orbital distance of 4 AU instead of 5.2 AU as in the three cases of MBPL10. The initial surface density is 6 g/cm² leading to a low isolation mass of about 2.4 M_⊕. The simulation is continued to a time t = 2 Myr when Rogers et al. (2011) find that the core has grown to 2.65 M_⊕, while the envelope mass is 0.54 M_⊕.

Figure 3 shows the envelope mass as a function of time as found by Rogers et al. (2011) and in our simulations for f_opa = 0.02, 0.003, 0.001, and 0 (grain free). The plot shows that for these initial conditions, f_opa ≈ 0.001 would lead to the best reproduction of the results of Rogers et al. (2011), while the nominal value of 0.003 calibrated with MBPL10 underestimates the actual envelope mass by a factor of 1.7. The f_opa = 0.02 simulation leads to an M_XY that is too low by a factor of 5, while the f_opa = 0 case results in an envelope that is too massive by a factor of 2.6. At a full ISM grain opacity, the envelope is 25 times less massive than in Rogers et al. (2011). This indicates that the mean calibrated value of f_opa = 0.003 is as expected only accurate to factors of a few, but that it still leads to differences that are significantly smaller than with arbitrary (or no) reduction factors. We must however recognize that even this fourth independent comparison case lies relatively close to the simulations of MBPL10 in the parameter space of initial conditions.

3.6. Effect on long term evolution

We have recently extended our formation model to a combined formation and evolution model (Mordasini et al. 2012b). We can therefore study the effect of different grain opacities also on the long term evolution. To do so we followed the evolution of a 1 M_♃ planet at 5.2 AU for different f_opa up to an age of 4.6 Gyr. We assume that f_opa is constant in time, also after the gas accretion has finished. In reality the grain opacity would probably quite rapidly tend to vanish after the planet has reached its final mass. The reason is that grains quickly rain out and get vaporized once the accretion of gas ceases (Rogers et al. 2011).

One way to quantify the impact of f_opa on the long term evolution is to study the planet’s radius and luminosity at 4.6 Gyr. We find that the radius R is constant for f_opa ≲ 0.004, and equal to about 0.99 R_♃. This means that for f_opa ≲ 0.004, the contribution of the residual grain opacity is so small that is does not significantly affect the evolution, but that it is determined by the gas opacities only. For f_opa larger than 0.004, we see as expected a gradual increase in the radius with f_opa due to the less efficient cooling, reaching about 1.15 R_♃ at the ISM opacity. For the luminosity we see an identical pattern with L ≈ 1.1L_♃ for f_opa ≲ 0.004, rising to about 1.8 L_♃ at f_opa = 1. This is a consequence of the fact that at a higher opacity, the planet’s interior remains longer at a higher entropy due to the stronger insulation. At early times, the luminosity of the high opacity case is in contrast lower, but the more efficient cooling of the low opacity case makes that its luminosity eventually falls below the one of the high opacity case (the post-formation entropies are nearly identical). This result indicates that using one constant f_opa = 0.003 during both the formation and evolution phase is an acceptable simplification. This is confirmed in Sect. 6 also in population synthesis calculations.

3.7. Relationship of M_Z and M_XY in phase II

The simulations with different f_opa show that the envelope mass as function of time during phase II depends on the opacity. It is however also interesting to investigate the envelope mass as a function of core mass for different opacities during phase II, because in this phase, there is a special relationship between core, envelope and total mass: if a planet of mass M at a semi-major axis a from a star of mass M_∗ has accreted all planetesimals (surface density Σ_P) within its feeding zone of half width B_L times the Hill sphere radius, it must hold (1)If the total mass of the planet is equal to the core mass (M = M_Z), this leads to the well known isolation mass (Lissauer 1993) (2)In a more general case, M = M_Z + M_XY. From the above equations we can then write (3)which leads to (4)This equation shows that for a given M_iso, there is only one M_XY for some core mass during phase II. The equation in particular means that in phase II, the envelope mass for a given core and isolation mass is independent of κ. At crossover for example, independently of opacity, as already noted by P96, again if there was sufficient time to accrete all planetesimals in the feeding zone³.

To illustrate this, Fig. 4 shows M_XY as function of M_Z for in situ simulations as in Sect. 3. Planets evolve in time along the lines, starting on the left. During phase I, the envelope mass is small and depends on f_opa, but once planets have reached the isolation mass and phase II begins, the envelope mass becomes independent of f_opa, so that the lines merge. The relationship of core and envelope mass is now given by Eq. (4) which is also shown in the plot using the corresponding M_iso. In these simulations without ejection, it remains valid also in the gas runaway accretion phase, which sets in after the crossover point.

However, how does the temporal evolution come in this picture, i.e., what determines how quickly the planet goes through these different stages if M_Z(M_XY) is independent of κ? We know that the gas accretion timescale depends on the opacity. The mechanism can be understood when considering the mass of the envelope as a function of luminosity L and opacity. Both the fully radiative envelope of Stevenson (1982), but also the generalization of Rafikov (2006) to envelopes that are convective in the interior and radiative at the surface (which is the kind of envelope we are dealing with here) have the property that M_XY ∝ 1/(Lκ).

This means that if L is larger, then κ must be smaller to get the same M_XY, as required for a given isolation and core mass in phase II. However, the luminosity L (in phase II ≈ L_Z, the core luminosity) is proportional to the accretion rate of planetesimals Ṁ_Z (which is in turn proportional also to Ṁ_XY, see P96 and Appendix A), so that a smaller κ goes along with a higher gas and solid accretion rate, and this means a shorter growth timescale.

The special relation of core and envelope mass together with a parametrization of the envelope growth timescale makes it possible to derive semi-analytical solutions for the core and envelope mass in phase II. This is presented in Appendix A while in Appendix B these semi-analytical solutions are used to derive the parameterization by comparison with numerical results.

4. Observational consequences of different f_opa

Up to this point, we have determined in this paper a nominal grain opacity reduction factor. We now turn to the question if different f_opa lead to potentially observable consequences.

Observational evidence of low envelope opacities could come from the bulk composition of planets with primordial H/He, and in particular from transiting low-mass low-density planets like the planets around Kepler-11 (Lissauer et al. 2011). The amount of primordial H/He these low-mass planets can accrete during the finite lifetime of a protoplanetary disk is a function of the grain opacity. A complication is that a significant spread of primordial H/He envelope mass for a given planet mass is expected, as the envelope masses depend besides the opacity also on, for example, the core luminosity, and therefore the individual formation history of a planet. This means that a statistically large enough sample of exoplanets is necessary for the comparisons.

In order to establish possible statistical consequences of grain opacity on observable quantities (e.g., the mass–radius relationship) we have conducted planetary population synthesis calculations with different values of f_opa. Regarding the conclusions that can be drawn from the comparisons, one should keep in mind the limitations discussed in Sect. 3.5, in particular the likely lack of generality of one global f_opa.

4.1. Population synthesis calculations

We calculated four populations of synthetic planets around solar-like stars varying the value of f_opa. The main assumptions for the populations are listed in Table 3. Otherwise, the model and probability distributions are identical as in Mordasini et al. (2012c). The reader is referred to this paper for details. The first three populations were calculated with the updated prescription for non-isothermal migration introduced in Mordasini et al. (2012b). These three populations only differ by f_opa which was set to 0 (grain-free gas, only molecular and atomic opacities from Freedman et al. 2008), 1 (full ISM grain opacity from Bell & Lin 1994), and 0.003, the value best fitting the MBPL10 results as derived in Sect. 3.4. The f_opa = 0.003 synthesis is the nominal case, while the other two syntheses represent limiting cases. An additional fourth population was calculated again with f_opa = 0.003, but this population differs in two other aspects: First, orbital migration was completely neglected, so that planets form in situ. This allows to see the differential impact of a migration as a mechanism that is currently not fully understood (e.g., Dittkrist et al. 2014). It also allows to study a situation that is closer to the semi-analytical solutions of Appendix A.1. Second, the accretion of planetesimals by the protoplanet in the disk limited gas accretion phase was set to zero (Ṁ_Z,run = 0). In the nominal case where Ṁ_Z,run ≠ 0 the core continues to accrete planetesimals also after the gas accretion rate due to the contraction of the envelope has exceeded the limiting rate set by the disk (see Mordasini et al. 2012b). The amount of planetesimals that get accreted in this detached phase could be significant due to the rapid expansion of the feeding zone. Hubickyj et al. (2005) or Lissauer et al. (2009) in contrast assume that the core accretion rate falls to zero once the planet is detached. With this synthesis, we study the global effects of this assumption.

5. Envelope mass as a function of core mass

In Appendix A.4, we discuss the envelope mass as a function of core mass as predicted by simple semi-analytical models. Figure 5 shows the envelope mass as a function of core mass as found numerically in the four populations. Planets with a semi-major axis of 0.11 ≤ a/AU ≤ 5 are included. For comparison, we plot for the subcritical planets a dotted line showing , corresponding to the simple approximation of Eq. (A.15) which is strictly speaking not applicable as it assumes Ṁ_Z = 0 (and M_XY ≪ M_Z). The exponent was set to p = 1.70 which corresponds to the mean for a core mass of roughly 6 M_⊕ as given in Table B.1, while b was set to 7.87. This corresponds to the approximate result for f_opa = 0.003, again for M_Z ~ 6 M_⊕. The time t was set to 2 Myr because the mean lifetime of the synthetic disks is 2.2 Myr. The dashed lines simply indicate scalings of M_XY as to guide the eye, where q is 2, 3, and 4. The absolute value was adjusted to fit the M_Z − M_XY relation of low mass planets in the in situ population.

5.1. Low-mass planets

We see that f_opa leaves a clear imprint in the M_Z − M_XY relation. For low-mass planets, we see that for a fixed core mass, the typical envelope mass increases with decreasing opacity, as expected. When going from f_opa = 1 to f_opa = 0.003, the envelope masses approximately increase by a factor of 4. For a 10 M_⊕ core for example, the typical envelope mass in the syntheses including orbital migration is roughly 0.5, 2, and >10 M_⊕ for f_opa = 1, 0.003, and 0, respectively. In the f_opa = 0 population, at a core mass of 10 M_⊕ there are both subcritical planets and giants. There is a large spread in possible envelope masses for a given core mass that covers one to two and sometimes even more orders of magnitudes. These are imprints of different formation histories, which lead in particular to different core luminosities and therefore envelope masses. Also the semi-major axes and the associated boundary conditions for the envelope given by the nebula affect the envelope mass for some planets (Rafikov 2006).

In the syntheses with orbital migration, there are more sub-structures in the M_Z − M_XY relation than in the in situ simulation. This is due to the more complex formation tracks with orbital migration. The tracks depend in particular on special convergence zones in the disk into which planets can get captured (e.g., Lyra et al. 2010; Dittkrist et al. 2014). A low-mass planet captured in a convergence zone migrates at a lower rate, leading to a lower planetesimal accretion rate and core luminosity relative to a migrating planet that is not in a convergence zone. However, in general, in simulations with orbital migration, the core luminosity is higher than in the in situ case for the prescription of planetesimal accretion used in the model (taken from P96) as it leads to high Ṁ_Z. This higher core luminosity results in lower envelope masses. We indeed see that for the synthesis with f_opa = 0.003 and migration, the envelope masses are roughly a factor of five lower than in the population neglecting orbital migration that is also calculated with f_opa = 0.003. A number of low-mass planets in these populations have clearly lower envelope masses ≲0.1M_⊕. This are planets in phase I at the moment when the disk disappears that therefore have high core luminosities.

5.2. Giant planets and critical core mass

It is known from hydrostatic calculations (e.g., Mizuno et al. 1978) that the critical core mass decreases with decreasing opacity (but not in P96 type evolutionary calculations for a given isolation mass⁴, see Sect. 3.7). The syntheses make it possible to study for a large number of evolutionary calculations how the minimal core mass M_Z,min required to form a giant planet depends on f_opa. From Hori & Ikoma (2010) it is known that for a grain-free envelope (f_opa = 0) and Ṁ_Z = 0, a core of just 1.7 M_⊕ can capture gas on a timescale that allows it to become a giant planet during the typical lifetime of a protoplanetary nebula.

Adopting as criterion for classification as giant planet an envelope mass of at least 100 M_⊕, we find the following M_Z,min in the syntheses with migration: M_Z,min = 6, 16 and 29 M_⊕ for f_opa = 0, 0.003, and 1, respectively. These values are the total final amount of heavy elements, including the planetesimals accreted during the gas runaway accretion phase which likely will end up mixed into the gaseous envelope. The actual critical core mass when gas runaway accretion is triggered is therefore smaller. In the synthesis without migration and f_opa = 0.003, M_Z,min is about 7M_⊕. An alternative definition of a gaseous planet is that the crossover point was reached during the presence of the disk, i.e., that M_XY = M_Z. The minimal core masses that fulfill this definition are 2, 9 and 25 M_⊕ for f_opa = 0, 0.003, and 1, respectively, for the simulations with migration, and 4 M_⊕ for the simulation without. This shows that at a low opacity, planets of a very small total mass can be gas dominated.

The value of M_Z,min = 6M_⊕ for the f_opa = 0 case is higher than the result of Hori & Ikoma (2010). This is expected because they considered isolated cores that do not accrete planetesimals, so that the luminosity is minimal. In the syntheses, the core luminosities are higher than in Hori & Ikoma (2010). Due to the lower typical core luminosity without migration, the M_Z,min = 7 M_⊕ found in the in situ population that assumes f_opa = 0.003 approaches the result for f_opa = 0, but including migration (6 M_⊕). We thus see that there is a degeneracy between luminosity and opacity in controlling the envelope gas masses as mentioned by Hori & Ikoma (2010).

5.3. Maximal core mass

Figure 5 can be used to study the maximal core mass M_Z,max. In the case where the core does not grow in the detached phase (bottom right panel), M_Z,max ≈ 60 M_⊕. This value depends on the largest semi-major axis that is considered. If synthetic planets at all semi-major axis are considered instead of only those inside of 5 AU, M_Z,max ≈ 100 M_⊕. This increase comes from the increase in the isolation mass with semi-major axis. In another synthesis with Ṁ_Z,run = 0, but normal orbital migration, M_Z,max ≈ 100 M_⊕. With migration, planets reaching such a heavy element mass are also found inside of 5 AU, showing that orbital migration moves planets with very massive cores closer to the star.

For the three simulations where the core continues to accrete in the detached phase, the maximal mass is higher. Here, core masses can exceed 200M_⊕, and there are planets with extreme cores (M_Z,max ≈ 400 M_⊕). This value does not increase if we include all semi-major axes instead of a ≤ 5 AU, an effect due to orbital migration. The measurement of mass and radius for some transiting exoplanets indicates that they contain several 100 M_⊕ of solids, like for example ~300–600 (200–500) M_⊕ of water (rock) in HAT-P-2b (Baraffe et al. 2008). These extreme cores can only form in solid rich disk, i.e., when both the disk gas mass and the dust-to-gas ratio (i.e., [Fe/H]) come from the upper end of the probability distributions. The product of the dust-to-gas ratio times the gas surface density (which is proportional to the disk gas mass) sets the surface density of the planetesimals, and this controls together with the initial position of the embryo the maximal mass to which a core can grow.

Rapid core growth the outer parts of the disk requires planetesimal random velocities to stay low as in our simulations that follow the eccentricity and inclination prescription of P96. Models of planetesimal dynamics (like Ida & Makino 1993; Ormel & Kobayashi 2012; Fortier et al. 2013) show that this is difficult to achieve. In this sense, our model has a tendency to predict too high heavy element masses.

5.3.1. Factors affecting : the capture radius

Figure 5 makes clear that the maximal core mass of giant planets depends in a significant way on the amount of solids that can be accreted in the detached phase with possibly observable consequences (Sect. 6). Physically, the magnitude of Ṁ_Z,run depends on the growth timescale of the planet, the timescale of gap opening in the planetesimals disk (and therefore the planetesimal size), as well as the capture radius of the planet and therefore its contraction timescale. In a detailed study, Zhou & Lin (2007) directly integrated the trajectories of a planetesimal swarm under the action of gas and tidal drag due to the protoplanetary disk, and the gravity of the star and a growing protoplanet (see also Shiraishi & Ida 2008). They studied what fraction of the planetesimals in the feeding zone get accreted during the gas runaway phase of the protoplanet (in the synthesis, we consider the limiting cases of 0 and 100%). The mass growth of the protoplanet was parametrized, and the actual capture radius for planetesimals was not calculated. Instead, they simply used mean planetary densities of 1 g/cm³ similar to Jupiter nowadays, and 0.125 g/cm³ as typical for a young Jovian planet that has a radius twice as large as Jupiter. For the 1 and 0.125 g/cm³ cases, about 20% and 28% of the planetesimals got accreted.

While this indicates a rather weak dependence on the capture radius (Shiraishi & Ida 2008 derive a dependency ∝), one should keep in mind that the capture radius can be very large during the formation phase. This is illustrated in Fig. 6. The plot shows the capture radius for 100 km planetesimal as a function of mass for a growing giant planet in a simulation with the same initial conditions a J1 in P96. The planet grows in situ at 5.2 AU, and its final mass is 1 M_♃. The capture radius is found by directly integrating the trajectories of planetesimals in the envelope (Podolak et al. 1988; Mordasini et al. 2006).

At a low mass, the capture radius is small and equal to the core radius. It then increases as the mass of the envelope grows, reaching a maximum of about 9 R_♃ at the moment when the planet detaches from the nebula. At this moment, the mass of the protoplanet is about 110 M_⊕. During the remaining runaway growth phase, it gradually decreases to about 4 R_♃. We thus see that the capture radius might be considerably larger than 2 R_♃. Using the scaling with , we crudely estimate with a mean capture radius of ~6 R_♃ that ~50% percent of the planetesimals contained in the feeding zone might get accreted. Then, the maximal core mass M_Z,max would come to lie about halfway between the two limiting cases of the syntheses. Clearly, these are only rough estimates, and dedicated simulations are necessary that concurrently model the planetesimal swarm and the evolution of the internal structure of the planet.

We see that in principle, the maximal core mass could also constrains the contraction timescale of giant planets formed by core accretion. Helled & Bodenheimer (2011) made similar considerations for the gravitational instability mechanism. They find that if grain growth is included, the timescale until collapse of the first core is much reduced. The amount of planetesimals that can be captured is much reduced, too, as the planet has a large capture radius only while it is extended. One finds an equivalent effect for core accretion: at low f_opa, the accretion of planetesimals in the detached phase is slightly reduced compared to the full opacity case. The reason is the faster decrease in the capture radius. For a 10 M_♃ planet we found for example that the timescale on which the capture radius decreases from 7 R_♃ to 2 R_♃ is roughly 84 000 and 37 000 years at f_opa = 1 and 0.003, respectively. The impact on the final core mass is however very small (difference of a few 0.1 M_⊕). This is because the difference in R_capt becomes only important at a moment when most of the planetesimal accretion has already ended.

6. Mass–radius relationship as a function of f_opa

The last section shows that f_opa has important consequences for the bulk composition of planets. The envelope mass of low-mass planets in particular increases by (very roughly speaking) one order of magnitude if κ decreases from the full ISM to grain-free. Since the radius of solid low-mass planets increases rapidly even if only a tenuous H/He envelope is added (e.g., Adams et al. 2008; Mordasini et al. 2012b), different values of f_opa should leave imprints in the planetary mass–radius relationship.

To test that, we calculated the cooling and contraction of the synthetic planets after the formation phase has finished (i.e., after the dispersal of the protoplanetary disk), using the simple, but self-consistently coupled formation and evolution model of Mordasini et al. (2012b). Deuterium burning in massive companions is included (Mollière & Mordasini 2012), but not important at the age of interest here (5 Gyr). Figure 7 shows the resulting synthetic mass–radius relationships for the same four populations as before. Planets with 0.11 <a/AU < 5 are shown. Colors indicate the fraction of heavy elements in the planet, Z = M_Z/M.

Compared to models of, e.g., Burrows et al. (1997), Fortney et al. (2011), or Baraffe et al. (2008), our evolutionary model is simplified in a number of aspects: simple gray boundary conditions are used, special inflation effects are neglected, and all solids are put into the core, so that M_Z = M_core. This might not be the case in reality since planetesimals can dissolve during impact. Our results for M_Z should therefore be associated with the total mass of heavy elements in a planet and not its core mass. Putting all solids into the core instead of homogeneous mixing with the envelope leads to radii which are somewhat bigger (Baraffe et al. 2008). Another simplification is that the opacity during evolution is identical for all planets in one population, while it is likely that different atmospheric compositions lead to different cooling curves (e.g., Burrows et al. 2011).

In the combined formation and evolutionary calculations, f_opa was for simplicity hold constant at the same value during both phases. In reality, a full ISM grain opacity over Gyr are not expected, because grains quickly rain out after the accretion of gas has ended (Rogers et al. 2011). We also calculated a synthesis where f_opa = 0.003 during the formation phase, but f_opa = 0 during the evolution. For this population we found radii at 5 Gyr which are virtually identical to the case where f_opa is kept at 0.003 also during evolution. This shows that during the evolution, gas opacities control the cooling for such a strongly reduced value of grain opacity as already found in Sect. 3.6.

For the f_opa = 1 case, the contraction is in contrast visibly slowed down for giant planets. This weakens the effect that for a given core mass, M_XY is usually lower at f_opa = 1 (Fig. 5). At the maximum of the mass–radius relationship at about 4 M_♃, we find a maximal radius of about 1.13R_♃ for f_opa = 1, while for f_opa = 0 the maximal radius is about 1.09R_♃. While this is not a large difference, it means that in the giant planet regime (M ≳ 1 M_♃), the delayed cooling with f_opa = 1 dominates over the fact that the Z is higher, as indicated by the colors.

For lower mass planets, and in particular for M ≲ 20 M_⊕, the protracted cooling at f_opa = 1 is in contrast more than offset by the much smaller primordial M_XY for a given total mass in this synthesis. The maximal radii of low-mass planets are much smaller in the f_opa = 1 synthesis than in the nominal f_opa = 0.003 case and even more so in the f_opa = 0 case. At M = 10 M_⊕ for example, the maximal radius in the syntheses with migration is approximately 4, 5.5 and 7.5 R_⊕ for f_opa = 1, 0.003, and 0, respectively. The grain opacity thus leads to a clear imprint in the planetary mass radius relationship via the presence or absence of low-mass, very large planets. This is interesting since otherwise, grain growth models are difficult to compare directly with observations. In the synthesis without migration and f_opa = 0.003, the maximal radius is about 6.5 R_⊕, i.e., between the value for f_opa = 0 and 0.003 in the simulations including migration.

6.1. Comparison with the observed mass–radius relationship

Figure 7 also shows the mass–radius relationship of planets in the solar system and extrasolar planets outside of 0.1 AU that have a well defined radius and mass. The observational data set available for comparison with the synthetic mass–radius relationship is obviously still rather small, mainly because of a> 0.1 AU.

However, already with the currently limited data set we note that for f_opa = 1 (bottom left panel), there are two planets (Kepler-11e and Kepler-18d) that lie above the envelope covered with synthetic planets. For Kepler-11e, the distance from the synthetic population is not very large given the error bars, but Kepler-18d lies several σ in mass and radius away. Since both Kepler-11 and Kepler-18 are old stars with estimated ages of 8 ± 2 and 10 ± 2.3 Gyr (Lissauer et al. 2011; Cochran et al. 2011), the larger radii of the two exoplanets cannot be due to a younger age of the host star relative to the synthetic population at 5 Gyr.

We therefore find an interesting result: at f_opa = 1, gas accretion is so inefficient in the model that no synthetic planets can form containing a sufficiently high M_XY to reproduce the observations. The failure of the f_opa = 1 population to reproduce actual low-mass exoplanets having a high gas mass fraction could therefore be an observational indication of grain growth in protoplanetary atmospheres.

Alternative explanations are also possible, like a delay of the contraction due to an opacity even higher than the full ISM grain opacity (but only during the evolutionary phase), or an additional energy source in the interior. Regarding the latter point it was found (Demory & Seager 2011) that at least for giant planets, no visible inflation seems to be active for orbital distances larger than ~0.07 AU around solar-type stars. Given the small eccentricities also internal heating due to eccentricity damping seems unlikely, at least for Kepler-18d (Cochran et al. 2011).

The population with f_opa = 0 (top left) behaves clearly differently. Here, the planets Kepler-11e and Kepler-18d are well within the envelope covered by synthetic planets. At their mass, there are many synthetic planets that are even significantly larger. The comparison with the observed planets shows that no such extremely large low-mass planets have been found to date. While it is possible that such planets will be found in future, the f_opa = 0 population thus has a large sub-population of low-mass, very large (i.e., very gas rich) planets presently without observational counterparts. This could mean that with grain-free opacities, gas accretion in the model is too efficient, and that a small, but still non-zero contribution of grains to the opacity is acting in protoplanetary atmospheres. An alternative explanation is that the planets lose parts of their gaseous envelope after formation via evaporation (or giant impacts), an effect we do not include in the current model. Since atmospheric escape could be important (e.g., Lopez et al. 2012; Owen & Wu 2013), the result that f_opa = 0 leads to too large planets needs to be revised once this effect is taken into account (Sheng et al. 2014).

The population with f_opa = 0.003 (top right) lies between the two extreme cases. The domain in the M − R plane covered by the synthetic planets seems similar to the observational result. Kepler-18d is at the upper boundary of the synthetic mass–radius relationship, but still in agreement when taking into account the error bars.

In the figure, there is a group of actual extrasolar planets with a mass similar to Saturn. In this mass domain, all three populations with migration show a certain excess of planets with a radius smaller than found in the observational data. This could be an effect related to the semi-major axis and orbital migration: in the figure, synthetic planets out to 5 AU have been included, while the actual exoplanets (but of course not Saturn) typically have smaller orbital distances. As discussed in Mordasini et al. (2012c) we find that despite migration, there remains a positive correlation between the amount of heavy elements in a planet and the orbital distance, as expected from the isolation mass. The excess of small synthetic planets in this mass domain is therefore due to planets at larger semi-major axes and disappears if only planets with a ≲ 1 AU are included. It is also absent in the population without migration, because the feeding zone is in general smaller for in situ formation, so that the core masses are lower (also the setting that Ṁ_Z,run = 0 in this population plays a role).

6.1.1. Population with in situ formation and Ṁ_Z,run = 0

The population with Ṁ_Z,run = 0 and in situ formation is shown in the bottom right panel. The change in a fundamental assumption in the formation model (in situ formation instead of orbital migration) does not completely change the mass–radius relationship, but does leave an imprint. We see that for low-mass planets, this population covers in terms of the largest associated R an envelope in the mass–radius plane that lies between the one of the f_opa = 0.003 and 0 population with migration, as expected. For giant planets the mass–radius relationship becomes vertically very thin, because with Ṁ_Z,run = 0, there are no massive cores that can reduce the radius significantly. This becomes particularly obvious for massive giant planets with M ≳ 1000 M_⊕, where the low maximal core masses of only ~60 M_⊕ in this population (see Sect. 5.3) cannot reduce the total radius in a significant way. This very thin vertical spread is of course also due to the fact that the planets have exactly the same age, and that the envelope opacity is exactly identical for all planets in the population. In reality, this will not be the case, so that a larger vertical spread is expected even if no very massive cores exist.

When comparing with the observational data for giant planets, we see that CoRoT-10b and HD 80606b lie slightly below the envelope of points predicted by this synthetic population. Neither CoRoT-10 nor HD 80606 are much older than 5 Gyr. For CoRoT-10, Bonomo et al. (2010) find weak constraints for the stellar age but even favor values smaller than 3 Gyr. For HD 80606, Pont et al. (2009) estimate an age of 5 ± 3 Gyr. Therefore the smaller radii should not be an age effect, but an indication that these planets might have higher metal contents than found in the model under the assumption of in situ growth and Ṁ_Z,run = 0. Indeed, Bonomo et al. (2010) estimate that CoRoT-10b contains between 120 and 240 M_⊕ of rocks, which is more than the maximal mass found in the in situ population for planets at a< 5 AU (Sect. 5.3). We however also see that the vertical distance of the planets from the synthetic population is small compared to the error bars, so that this is a weak constraint.

6.1.2. Observational hints of a strongly reduced opacity

The comparison of the mass–radius relationship of low-mass planets for f_opa = 0, 1, and 0.003 (including migration) with observations seems in summary to indicate that grain growth occurs in protoplanetary atmospheres, leading to an opacity that is much reduced compares to the ISM value. It is currently not possible to constrain from our simulations the specify value of κ_gr, because first, a range of f_opa between the extreme cases would likely lead to a similar M − R as in the observations, as long as f_opa is clearly smaller than unity. This is due to the currently low number of actual exoplanets and the intrinsic large spread. Second, and more importantly, a κ_gr found by scaling the ISM by one general f_opa cannot reproduce the dependency of the microphysical processes that govern the grain growth on planet properties. This dependency eventually influences the resulting M − R. We will therefore use in a forthcoming work the analytical grain growth model of Paper II for improved population syntheses.

Clearly, our results are tentative as both the observational comparison sample is rather small and the theoretical results are obtained within a framework of a formation and evolution model that is much simplified compared to reality. The theoretical results will thus partially depend on other assumptions in the model like for example the core accretion rate or the migration model.

7. Relative heavy element content of giant planets

We now turn to the second observable quantity which potentially contains an imprint of κ_gr besides the M − R relation of low-mass planet. It is the heavy element content of giant planets. In Fig. 7 we compared the position of observed extrasolar planets relative to the envelope of points covered in the syntheses. The color code in the figure shows that not only the envelope of M − R points changes, but also the bulk composition of planets within the covered envelope of M − R loci. In this section we use this result to compare the metal enrichment of giant synthetic planets with constraints derived from observations. The heavy element content of giant planets was studied in several works. Guillot et al. (2006), Burrows et al. (2007), and Guillot (2008) have in particular found that there is a positive correlation between the stellar [Fe/H] and the planetary heavy element content, a correlation that is reproduced by population synthesis models based on the core accretion paradigm (Mordasini et al. 2009b). This correlation will be re-addressed in future work.

7.1. Relative metal enrichment in observed exoplanets

In this work, we focus on anther aspect. We compare the heavy element enrichment of planets relative to the host star Z_pl/Z_star in synthetic and actual transiting extrasolar planets as derived by Miller & Fortney (2011). The planet’s metal fraction is Z_pl = M_Z/M, while the stellar metal fraction is Z_star = Z_⊙10^[Fe/H]. The heavy element fraction of the Sun is taken as Z_⊙ = 0.0142 (Asplund et al. 2009) while [Fe/H] is measured spectroscopically. We thus assume that iron is a good indicator of the global metal content of a star which is generally the case as long as we are dealing with thin-disk stars (Adibekyan et al. 2012).

The sample of transiting extrasolar planets analyzed by Miller & Fortney (2011) consisting of 14 planets that receive a stellar irradiation flux of less than 2 × 10⁸ erg/(s cm²). Below this threshold, no significantly inflated planets are found. This is in agreement with the analysis of Demory & Seager (2011) of the Kepler sample. The 14 planets have masses greater than 20 M_⊕ and have, with one exception, roughly solar-like host stars with 0.7<M_⋆/M_⊙< 1.3. For each planet, Miller & Fortney (2011) determine the heavy element fraction Z_pl by internal structure modeling (see Fortney et al. 2007). In the model, no extra heating sources are included, and the atmospheres have solar composition. The derived planetary heavy-element masses are therefore likely minimum masses because higher opacities (e.g., Burrows et al. 2007) or additional heating sources would require a higher heavy element fraction to obtain the observed radius. The results are obtained based on a number of general assumptions like a fully convective interior or a specific equation of state (Saumon et al. 1995). Regarding the former point, one must keep in mind that there are also semi-convective interior models for Jupiter that agree with observations (Leconte & Chabrier 2012). These models can have a ~60% higher metal fraction. Concerning the latter point, uncertainties in the EOS of H/He alone cause significant variations for the estimated bulk composition of Jupiter (e.g., Guillot 1999; Militzer et al. 2008; Nettelmann et al. 2012). Due to these reasons, Z_pl is currently only poorly known even for Jupiter (and the other giant planets of the solar system). This means that the results for the heavy element mass M_Z in extrasolar planets come with large systematic uncertainties.

The planetary heavy element fractions derived by Miller & Fortney (2011) are the average of the values for the limiting cases of either all metals in the core or homogeneous mixing. This leads strictly speaking to an inconsistency with the formation model where all solids are in the core. The difference in Z_pl derived for the layered and homogeneous case is however rather small with typically 15% higher Z_pl in the layered case. Within the errors bars due to other uncertainties (e.g., mass and radius measurement) the Z_pl for the two cases agree for all 14 planets.

7.2. Relative enrichment for synthetic planets

The relative enrichment Z_pl/Z_star is straightforward to calculate for the synthetic planets, as one of the four Monte Carlo variables (i.e., initial conditions) in the syntheses is [Fe/H] (Mordasini et al. 2009a). Once a synthetic planet has formed, we thus have its Z_pl and the Z_star of the corresponding initial condition. Z_star represents the bulk metal fraction of the star, and not the dust-to-gas ratio in the innermost ~10 AU that we use to calculate the planetesimal surface density (Mordasini et al. 2009a). The latter value can be higher by a factor of 2–4 (2.8 in the formation model) than Z_star due to the advection of dust from the outer parts of the disk (Kornet et al. 2004, Birnstiel et al. 2012).

Two assumptions in the model are of relevance here: first, we assume that the entire metal content in the disk has been incorporated into planetesimals. Consequently, the nebular gas accreted by the protoplanets is composed of pure H/He. Clearly, this is a limiting case. In Sect. 7.4 we study effects of this assumption. Second, we assume that gap formation does not reduce the gas accretion rate of massive planets due to the eccentric instability (Kley & Dirksen 2006). Again this is a limiting case, because in reality, the eccentric instability can only operate for certain configuration of the location of the Lindblad resonances and the width of the gap. Without a reduction of Ṁ_XY due to gap formation, very massive planets can form with masses up to ~40M_♃. The specific value depends on the largest initial disk mass. They are very rare though (Mordasini et al. 2009b), and only appear numerous in the plots because we are dealing with large populations of ~30 000 planets. These very massive planets were already visible in the plots above, but here they are more obvious because in the observational comparison sample of Miller & Fortney (2011), the maximal mass is in contrast only about 4 M_♃.

7.3. Z_pl/Z_star as a function of planet mass

Figure 8 shows Z_pl/Z_star for the observational comparison sample of Miller & Fortney (2011) and the synthetic populations. For the three populations with migration, planets inside of 1 AU are included (the actual exoplanets have a< 0.41 AU). For the population without migration, planets inside of 5 AU are included because no giant planets form in situ inside of 1 AU. The giant planets of the solar system are also shown. For Jupiter and Saturn, we use the M_Z estimated by Saumon & Guillot (2004) for different EOS, while for Uranus and Neptune the bulk composition given by Guillot & Gautier (2009) is used. The conversion into Z_pl/Z_⊙ is again made with the updated protosolar Z_⊙ = 0.0142 of Asplund et al. (2009). This value is substantially lower than earlier estimates, resulting in correspondingly higher enrichment of the planets.

While all populations (and the observed exoplanets) share the property of a decreasing Z_pl/Z_star with increasing mass, there are substantial differences in the absolute level of the metal enrichment. At 1 M_♃, the relative enrichment at f_opa = 1 ranges for example between 4 and 20, while for f_opa = 0, the enrichment is between 1 and 13. The impact of f_opa becomes particularly clear when one considers the least-squares fit to the observational data (red solid line with red dotted lines indicating approximately the 1σ error bars). The fit is given as (5)The parameters α and β are shown in Table 4. Corresponding fits to the synthetic data are also shown in the figure. These fits were made over the same mass domain as for the observations. The figure shows that the decrease in Z_pl/Z_star with mass follows a similar slope α ≈ −0.7 independently of f_opa for all three populations with migration. This slope is in good agreement with the observational value of −0.71 ± 0.10. The absolute level given by β decreases significantly from 8.5 at f_opa = 1 to 3.5 at f_opa = 0, corresponding to a decrease by a factor of ≈2.4. When we compare these values with the observational result β = 6.3 ± 1.0 we see that the full ISM grain opacity leads to too enriched planets, while the grain-free opacity leads to too low enrichments. At f_opa = 0.003, β = 7.2 which is similar to the observational result. We thus find for the bulk composition of giant planets a similar result as we had found for the mass–radius relation of low-mass planets: grain evolution in protoplanetary atmospheres seems to lead to an opacity that is much smaller than in the ISM, but potentially, there is still a non-zero contribution from the grains.

One can ask why the bulk composition also of giant planets depends on f_opa despite the fact that in the disk limited gas accretion phase (when the planet is detached from the nebula), the gas accretion rate is independent of f_opa (in contrast to phase II). The imprint of f_opa on the bulk composition nevertheless remains also for giant planets because at small f_opa, the critical core mass is lower. Also the amount of solids accreted in runaway is lower due to the smaller surface density of planetesimals in the feeding zone. In the population without migration and Ṁ_Z,run = 0, the decrease in the relative enrichment is as expected steeper (α = −0.88), because in this population, the maximal cores masses are lower (Sect. 5.3). Also the absolute level of Z_pl/Z_star is low (β = 2.4). This is related to the fact that without migration, the feeding zones and luminosities are smaller.

As stated earlier, it is clear that the results for κ_gr in protoplanetary atmospheres have been obtained from models that are strongly simplified compared to reality. We have for example not considered the effect of different chemical compositions of the envelope due to envelope pollution. Hori & Ikoma (2011) have shown that a pollution with heavy elements reduces the critical core mass. The relation between composition and grain opacity on the other hand is complex: Movshovitz & Podolak (2008) have found that an increase in the dust-to-gas ratio in the newly accreted can actually lead to a decrease in κ_gr in the deeper layer. This is extensively studied in Paper II. We see that the results found here must be further investigated with models that combine self-consistently the structure of the envelope, the evolution of the grains, the deposition of mass by planetesimals, and the chemical composition of the gas.

7.4. Effect of gas composition

In the observational data set analyzed by Miller & Fortney (2011) where M ≤ 4 M_♃, all planets are consistent with being enriched relative to the host star, and several seem to be enriched even significantly. Figure 8 shows that in the synthetic populations, there are planets with Z_pl/Z_star< 1 for sufficiently large masses. These cases are a direct consequence of the assumption that the accreted gas is pure H/He so that a massive gas envelope can deplete a planet relative to its host star (see Helled & Bodenheimer 2011). Extrapolating Eq. (5) with the parameters of Miller & Fortney (2011) leads to a mass M₁/M_♃ = β^{− 1 /α} where Z_pl/Z_star = 1 of about 13.9 M_♃. In the nominal f_opa = 0.003 population, the lowest planet mass where Z_pl/Z_star = 1 is about 4 M_♃, while the typical M₁ is 10–15 M_♃. These values are a function of f_opa: at f_opa = 0.0, the lowest M₁ is about 1 M_♃, while the typical value is about 5 M_♃. At f_opa = 1, the lowest M₁ is only about 5 M_♃, while typically it is about 17 M_♃. The relative enrichment Z_pl/Z_star decreases more steeply with mass in the population without planetesimal accretion in the runaway accretion phase. The smallest M₁ is here about 1 M_♃, and the typical value is about 2.5 M_♃. These are clearly lower values. For the nominal populations, M₁ is coincidentally similar to the deuterium burning limit at approximately 13 M_♃ also for companions formed via core accretion (Mollière & Mordasini 2012; Bodenheimer et al. 2013). We see that neither deuterium burning nor the absence of an enrichment relative to the host star necessarily excludes core accretion as the formation mechanism.

The observed infrared excess of many young stars shows that hot micrometer-sized dust remains in protoplanetary disks on timescales of 10^{6 − 7} yr (e.g., Haisch et al. 2001). This shows that in reality, not all dust grains get incorporated into planetesimals leaving the gas essentially dust free as assumed in the formation model. Instead, the gas that is accreted by the protoplanets will itself contain a certain (unknown) fraction of metals. To investigate that, we also calculated Z_pl/Z_star for the synthetic planets under the assumption that the gas has the same composition as the star. This is not self-consistent with the way the solid surface density is derived from [Fe/H] in the model (it means that we artificially generate metals), but it gives an impression of the impact of different gas compositions. An identical metal fraction as in the star is, however, not necessarily the other limiting case, as the gas could also be enriched relative to the star due to photoevaporation of gas only. The enrichment of the synthetic planets under this assumption is also shown in Fig. 8.

The consequence is as expected: at low masses, where the core dominates, Z_pl/Z_star does not change appreciably, while for clearly envelope-dominated massive planets, Z_pl/Z_star approaches unity. In the three populations with migration and accretion of planetesimals in the detached phase, Z_pl/Z_star does not fall below ~1.1 because of the massive cores in the planets and the positive correlation of [Fe/H], heavy element mass (Guillot et al. 2006; Burrows et al. 2007; Guillot 2008; Mordasini et al. 2009b) and total mass for the most massive planets (Mordasini et al. 2012a). In the population with Ṁ_Z,run = 0 the convergence towards unity for the most massive planets is clearer.

In the plot, we also include the planets of the solar system. Even if the uncertainty of their heavy element content is large, we see that they follow the general trend seen for the exoplanets. The nominal values of Jupiter’s and Saturn’s enrichment fall quite close to the mean relation derived from the exoplanets. In the Neptunian mass domain, the enrichment becomes very high with Z_pl/Z_star between 20 and 100 in the synthetic planets. This covers the enrichment estimated for Uranus and Neptune. For these planets, traditional three layer models predict a heavy element mass fraction of about 85–95 %, corresponding to Z_pl/Z_star roughly equal to 60 (Guillot & Gautier 2009).

7.5. Comparison with gravitational instability

In the simplest picture of the core accretion mechanism, first a critical core forms, and afterwards only gas is added “diluting” the planet’s enrichment. This would mean that for a given Z_star, the relative enrichment ∝M_Z/M would decrease with mass as α ≈ −1 for massive giant planets where M ≈ M_XY. In the simplest picture of gravitational instability one would in contrast expect α ≈ 0. The numerical (and observational) result that α ≈ −0.7 indicates that for Ṁ_Z,run ≠ 0, the core mass itself increases weakly with total mass roughly ∝M^0.3. Analytically, without orbital migration and planetesimal ejection, and if the planet accretes all planetesimals in the feeding zone, Eq. (3) shows that M_Z/M falls as α = −2/3. The observational result (α ≈ −0.7) thus lies between the cases that the giant planet accretes all planetesimals in the feeding zone (α = −2/3), and no accretion in runaway (α = −1). In the population where Ṁ_Z,run = 0, α is indeed closer to –1 (–0.88).

Our results indicate that giant planets with masses less than ~10 M_♃ formed by core accretion are enriched due to the accretion of planetesimals roughly as 6(M/M_♃)^-0.7. Figure 8 shows that the spread in Z_pl/Z_star around this relation for a given mass is however large, typically about one order of magnitude. For more massive giant planets the composition of the gas becomes important in determining the bulk composition, and both enriched as well as depleted planets may form depending on the efficiency of dust and planetesimals growth.

The composition of giant planets formed by the gravitational instability mechanism was studied by, e.g., Helled & Schubert (2009), Nayakshin (2010), Boley et al. (2011), and Helled & Bodenheimer (2011). The initial composition of a clump should be similar to the composition of the disk and therefore Z_pl/Z_star ~ 1 independent of mass. Several mechanisms can however lead to enrichment: enrichment at birth, planetesimal capture (before and after the collapse of the first core), and differentiation plus tidal stripping. These processes can lead to a large spread of final bulk compositions also for gravitational instability.

Enrichment via planetesimal capture is efficient as long as the planet is extended (i.e., before the second collapse) because of the large capture radius. This is similar to the behavior for core accretion (Sect. 5.3.1). The duration of the pre-collapse phase depends on opacity. If grain growth is taken into account (Helled & Bodenheimer 2011), the duration of the pre-collapse phase is very short (~ 10³ years) so that enrichment by planetesimal capture is precluded over a wide domain of planetary masses and metallicities. This would result in roughly stellar compositions if this is the dominant enrichment mechanism.

Due to the complexity and the multitude of interacting mechanisms that determine the final composition of giant planets in both formation mechanisms (additional effects are, e.g., giant impacts), it seems likely that for many individual planets one cannot distinguish core accretion and gravitational instability based on the bulk composition only. More promising should be statistical correlations, namely the (mean) enrichment as function of mass (as studied here), stellar [Fe/H], and orbital distance. Regarding the orbital distance, Helled & Bodenheimer (2011) find that for giant planets formed by direct collapse, the metal enrichment strongly decreases with orbital distance. The impact of distance and [Fe/H] on the bulk composition for core accretion will be addressed in future work. An important finding from internal structure modeling is here that the heavy element content of giant planets and the stellar [Fe/H] are correlated, as found by Guillot et al. (2006) and Burrows et al. (2007).

8. Summary and conclusions

In this series of papers we investigate the impact of the opacity due to grains in protoplanetary atmospheres on the predicted bulk composition (H/He fraction) of extrasolar planets. In this first paper, we study the impacts found by scaling the ISM opacity. In Paper II, we present an analytical model to calculate κ_gr. In future work, we will couple the analytical model to our updated population synthesis calculations. The results of this first paper are summarized as follows:

• In the context of the core accretion paradigm, we studied theduration τ_II of phase II for in situ formation ofJupiter as a function of the reduction factor of grain opacity f_oparelative to ISM grain opacity. We found that over an importantrange of f_opa, there is a linear relationship between f_opa and τ_II as expectedfrom theoretical considerations (Sect. 3.3).

• We compared the duration of τ_II as function of f_opa with the corresponding duration found by Movshovitz et al. (2010) who conducted simulations of combined giant planet formation and grain evolution. We found that the ISM grain opacity must on average be reduced by a factor f_opa ≈ 0.003 to reproduce their results (Sect. 3.4). As a caveat one must keep in mind that one uniform reduction factor cannot reproduce the complex radial opacity structure found in grain evolution calculations, and also that the calibration was only made in a small part of the parameters space of possible core masses, luminosities, and outer boundary conditions (Sect. 3.5).

• Without migration and planetesimal drift, there is a unique relation between isolation, core, and total mass during phase II if the planet accretes all planetesimals in the feeding zone. The reason is that the core mass M_Z is given by the size of the feeding zone that depends via the Hill sphere on the total mass M. This means that for a given core and isolation mass, the envelope mass M_XY during this phase is independent of f_opa. The opacity merely determines how quickly a planet evolves through the different M_Z − M_XY states (Sect. 3.7).

• We studied the global consequences of a scaled ISM opacity on planets forming via core accretion with population synthesis calculations (Sect. 4.1). In reality, the grain opacity in protoplanetary atmospheres cannot be obtained as the ISM opacity simply reduced by one general f_opa. This is because the grain dynamics depend on a planet’s properties (Movshovitz & Podolak 2008; Paper II). Therefore, no constraints for the value of κ_gr for a specific planet can be derived from our calculations. Due to this limitation, we considered besides the nominal value f_opa = 0.003 also populations with full ISM grain opacity (f_opa = 1) and grain-free opacity (f_opa = 0) as limiting cases. This allows to see if the opacity leads at all to potentially observable imprints, which is the main goal of the paper.

• Grain opacity leaves clear imprints in the planetary bulk composition (H/He envelope mass as a function of core mass). For sub-critical cores, the envelope mass for a given core mass increases with decreasing κ_gr. In the synthesis with f_opa = 0.003, the envelope masses of sub-critical cores (M_Z ≲ 10M_⊕) are about four times more massive than for full ISM opacity. There is a large spread of at least one order of magnitude in M_XY for a given M_Z (Sect. 5).

• The critical core mass for runaway gas accretion decreases with f_opa. The lowest core mass M_Z,min in giant planets (envelope mass M_XY ≥ 100M_⊕) are 6, 16 and 29 M_⊕ for f_opa = 0, 0.003, and 1, respectively. These results are for populations where orbital migration is included. In a synthesis without orbital migration and f_opa = 0.003, M_Z,min = 7 M_⊕. This difference is due to the lower mean luminosity without migration and the dependency of the critical mass on both the luminosity and opacity (Sect. 5.2).

• If planets continue to accrete planetesimals in the disk limited gas accretion phase, and if the planetesimal random velocities stay low, then the maximal heavy element content in giant planets we find can be up to ~400 M_⊕. If in contrast no planetesimals are accreted in the disk limited phase, then the maximal heavy element content is about 50 to 100 M_⊕ depending on the semi-major axis (Sect. 5.3).

• Grain opacity leaves a clear imprint in the mass–radius relationship of low-mass planets. At a low opacity, planets in the super-Earth and Neptunian mass domain have radii for a given mass that are much larger than for ISM opacity because of higher envelope mass fractions. At M = 10 M_⊕, for example, the maximal radii at 5 Gyr in the syntheses with migration are 4, 5.5 and 7.5 R_⊕ for f_opa = 1, 0.003, and 0, respectively (Sect. 6).

• We compared the envelope covered by actual and synthetic planets in the mass–radius plane assuming an age of 5 Gyr for the synthetic planets. One finds that for f_opa = 1, the synthetic super-Earth and Neptunian planets have too small radii (i.e., too low envelope masses) to cover the same domain as the observations. At f_opa = 0.003, the value calibrated with a grain evolution model, the synthetic and actual planets occupy a similar loci in the mass–radius plane. This is a hint that the opacity in protoplanetary atmospheres is much smaller than in the ISM (Sect. 6.1).

• A second consequence of grain opacity that can be compared with observational data is the bulk composition of giant planets. We compared the bulk composition of synthetic planets as expressed in the enrichment of a planet relative to its host star Z_pl/Z_star (Z_pl = M_Z/M, Z_star = Z_⊙10^[Fe/H]) with the results derived from observations (Miller & Fortney 2011). Miller & Fortney (2011) derived the heavy element content of 14 transiting extrasolar planets obtaining low irradiation fluxes by internal structure modeling (Sect. 7.1).

• We found that there is a clear imprint of opacity on the bulk composition of giant planets. The mean relative enrichment Z_pl/Z_star is about 2.4 times higher for full ISM grain opacity compared to the grain-free case (7.3).

• The mean relative enrichment of giant planets as a function of planet mass M can be approximated as Z_pl/Z_star = β(M/M_♃)^α. The decrease in Z_pl/Z_star with mass follows α ≈ −0.7 independently of f_opa for all three populations with orbital migration. This slope is in good agreement with the value derived from observations (−0.71 ± 0.10). In the simplest picture of core accretion, where first a critical core forms, and afterwards only gas is added, α ≈ −1, while in the simplest picture of gravitational instability α ≈ 0. The α ≈ −0.7 found in the simulations indicates a weak positive correlation of M_Z with total mass. A similar exponent (α ≈ −2/3) is expected if giant planets efficiently accrete all solids in their feeding zone (Sect. 7.3).

• The absolute level of enrichment β decreases significantly from 8.5 at f_opa = 1 to 3.5 at f_opa = 0. When we compare these values with the result derived from observations (β = 6.3 ± 1.0) we see that a full ISM grain opacity leads to too enriched planets, while a grain-free opacity leads to a too low enrichment. At f_opa = 0.003, β = 7.2, which is similar to the observational result. The bulk composition of giant planets thus seems to hint that the opacity in protoplanetary envelopes is much smaller than in the ISM, but that there is possibly still a non-zero contribution from the grains (Sect. 7.3). This result is similar as the one found for the mass–radius relationship of low-mass planets.

• Giant planets with masses less than ~10 M_♃ formed by core accretion are enriched due to the accretion of solids, with the relative enrichment Z_pl/Z_star following roughly 6(M/M_♃)^-0.7 in the observational sample and similarly in the nominal synthetic population. The spread around this relation at a given mass is large, about one order of magnitude. For more massive giant planets the composition of the gas becomes important in determining the enrichment, and both enriched and depleted planets are possible. In the limit that the accreted gas is free of solids, planets more massive than ~10–15 M_♃ are depleted relative to the host star (Sect. 7.4).

• We derived a semi-analytical expression for the core and envelope mass as a function of time in phase II using a two parameter expression for the gas accretion timescale τ_g (Appendix A). We determined the parameters of τ_g for different core masses and f_opa by comparison with numerical results (Appendix B).

We have found in this work that the ISM grain opacity must be reduced by a factor f_opa ≈ 0.003 to reproduce the formation timescales of Movshovitz et al. (2010). This value is almost an order of magnitude lower than the previously considered “low opacity case” with 2% ISM opacity (P96, Lissauer et al. 2009). It is however clear that uniform reduction factors cannot capture the physical effects of grain evolution that in reality lead to complex opacity profiles in the envelope (Movshovitz & Podolak 2008, Paper II) which differ from the ones found by scaling the ISM opacity. The interest in a calibrated f_opa is therefore merely to have an intermediate case for the magnitude of κ_gr in simulations that study the global impact of opacity in a parametrized way between the limiting f_opa = 0 and 1 cases.

With this grain reduction factor, cores with a low mass can trigger gas runaway accretion during the typical lifetime of a protoplanetary nebula, the exact value depending on luminosity. As discussed by Movshovitz et al. (2010) and Hori & Inaba (2010), the lower critical core masses found with more realistic opacities are an important result for the core accretion paradigm, as several studies (e.g., Fortier et al. 2007; Ormel & Kobayashi 2012) indicate that building up a 10 M_⊕ core at 5.2 AU within the typical lifetime of a disk is a critical timing issue. These lower opacities contribute substantially in making this timing issue less stringent, even without taking into account additional mechanisms like migration or envelope pollution (Alibert et al. 2004; Levison et al. 2010; Hori & Ikoma 2011). Also, at such low f_opa, cores of just a few M_⊕ can accrete quite significant H/He envelopes, making them potential progenitors of the recently detected low-mass, low-density planets.

The global consequences of grain opacity on planetary populations are strong and multiple. We have studied these imprints by running population synthesis calculations with a wide range of f_opa. We focused on two possibly observable consequences, namely the mass–radius relationship of low-mass, subcritical planets, and the bulk composition of giant planets. Additional imprints also exist in the planetary mass and radius distributions, and in particular in the frequency of giant planets. This frequency increases in the synthetic populations approximately by a factor of three when reducing f_opa from 1 to 0.

The resulting comparison of synthetic and (statistical) observational results establishes a connection between the opacity in protoplanetary envelopes and observable quantities. In this work, we cannot directly derive constraints on the value of κ_gr because we find it by scaling the ISM opacity instead of calculating it based on planet properties. However, in future, a similar approach should allow to test microphysical models of grain evolution that are otherwise difficult to test observationally. The result of this first paper that the observed mass–radius relationship cannot be reproduced with a full ISM opacity is interesting. It could be an observational hint that the prediction of grain growth models (Podolak 2003; Movshovitz & Podolak 2008, MBPL10, and Paper II) that the opacity in protoplanetary atmospheres is much smaller than in the ISM, is correct.

However, it is clear that our results for the predicted bulk composition of extrasolar planets are preliminary. Too complex and too numerous are the actual physical mechanisms occurring during formation and too simple is the model in comparison that uses, for example, one uniform reduction factor, neglects the impact of the chemical composition of the gas, and only handles one embryo per disk. For a better understanding, it is necessary to couple a physically motivated grain opacity model like the analytical model of Paper II with planet formation codes, and to take into account additional important factors like the pollution of the envelope (Hori & Ikoma 2011), the concurrent formation of many planets in one disk (Alibert et al. 2013), or the loss of the H/He envelope during evolution due to atmospheric escape (Sheng et al. 2014).

On the observational side, a significant extension of the sample of relatively cold, transiting planets with well-constraint mass and radius would be very important, for example with CHEOPS (Broeg et al. 2013) and PLATO (Rauer et al. 2013). Then it should become possible to understand the processes that determine the opacity in protoplanetary atmospheres much better than today. The associated statistical results for the enrichment of the planets will allow to distinguish better different formation mechanisms like core accretion and gravitational instability, or to understand from a theoretical point of view the transition from solid to gas dominated planets (Marcy et al. 2014).

Online material

Appendix A: Semi-analytical solutions

The left panel of Fig. 1 in Sect. 3 shows that the curves giving the mass as a function of time for different opacities are self-similar in the sense that they can be approximately converted in each other by stretching the x-axis. This indicates that there could possibly be (semi-)analytical solutions for the mass as a function of time in phase II. Since such solutions are interesting for comparison with numerical results, we derive a semi-analytical solution for the envelope and core mass as a function of time during phase II.

During this phase, the envelope grows by the contraction of the outer layers (e.g., Ikoma et al. 2000), while the core grows due to the resulting expansion of the feeding zone (P96). This growth mode applies if the accretion rate of solids is limited directly by the availability of planetesimals in the feeding zone, rather than the probability that they get captured by the protoplanet. In other words, the timescale of accretion of planetesimals already within the feeding zone is much smaller than the timescale on which the feeding zone expands. This regimes can arise if the planetesimals are dynamically cold, so that the capture probability is high, as assumed in P96.

The gas accretion rate is parametrized with the growth timescale τ_g of a planet of total mass M = M_Z + M_XY where M_Z is the total mass of heavy elements in the planet (equal to the core mass if all planetesimals sink to the core) while M_XY is the envelope mass. It has become customary to approximate the growth timescale of the envelope in the form (e.g., Ikoma et al. 2000; Bryden et al. 2000; Ida & Lin 2004; Miguel & Brunini 2008) (A.1)where we adopted the same form as Ida & Lin (2004) who use the total mass, while Ikoma et al. (2000) use the core mass only. The core and total mass are at least initially during phase II in any case similar. The two parameters k = 10^b and p (b and p are both parameters >0) are found by comparison with numerical simulations. Clearly, the parameters will depend on the core accretion rate and the associated luminosity, the core mass, and the opacity κ in the envelope that is represented in this study by f_opa. The explicit dependency of k on f_opa can be merged into the parameter, writing (at least over a certain range of f_opa), where k₀ and q are a priori unknown (in Sect. B we will see that q ≈ 1 over an important domain of f_opa). Table A.1 gives an overview of values found in the literature. The units in the fit are years and Earth masses. Note that some fits were derived for planets accreting planetesimals, while others were derived for planets where the luminosity is coming from the accretion of the gas itself only. This naturally explains part of the differences. Still, as already discussed by Miguel & Brunini (2008), the (remaining) significant differences and the high powers lead to highly uncertain results for the predicted gas envelope masses. The comparison of the semi-analytical solution presented below and the numerical simulations makes it possible to determine the parameters in a more systematic way.

Appendix A.1: Growing core mass:

First, we consider the full problem where the core mass increases, assuming that the expansion of the feeding zone is the only limiting factor for core growth, and that the core instantaneously accretes all planetesimals within it. The feeding zone is proportional to the Hill sphere radius of the planet which in turn is proportional to M^1/3, therefore (Sect. 3.7) (A.2)Differentiation with respect to the time gives (A.3)This leads to the following system of coupled differential equation for the core mass M_Z(t) as a function of time t (P96) (A.4)which holds very well during phase II in the numerical simulations above⁵ This equation is also found from differentiation of Eq. (4). For the envelope mass M_XY(t) we have (A.5)As boundary conditions we use M_Z(0) = M_Z,0 which is the core mass at the beginning of phase II (typically the isolation mass), while we assume that the initial mass of the gaseous envelope M_XY(0) = 0. This approximation is valid if the antecedent core luminosity is high, and/or if the core mass is low. The time t is counted relative to the beginning of phase II.

While a closed form for M_Z(t) and M_XY(t) solving the differential equations cannot be found, we still find an implicit equation which can be simply solved by determining numerically the root. One finds (A.6)and (A.7)The second equation is simply the solution to Eq. (A.4) and has as such no explicit temporal dependence.

Two examples of the temporal evolution of the core, envelope and total mass obtained with these equations are shown in Fig. A.1. The initial mass of the core is M_Z,0 = 5 M_⊕, and the parameters of Ida & Lin (2004) and Miguel et al. (2011) are used to illustrate the impact of different parameters. As in numerical calculations, the evolution is characterized by a long plateau phase of slow growth, which accelerates rapidly once the core and the envelope mass become comparable.

The semi-analytical solution can be used to calculate the time of crossover τ_cross when M_Z = M_XY, and of runaway τ_run when M_Z and M_XY approach infinity. In reality, the growth gets then limited due to the finite reservoirs of both gas and solids in the disk and/or gap formation.

The core mass at crossover derived from the above equations is equal to as already found by P96. From that we find that the time when crossover occurs (relative to the beginning of phase II) is (A.8)where is the growth timescale of the initial core, while the rest of the equation corresponds to a numerical factor smaller than unity, as the growth timescale decreases with increasing mass.

The time of runaway accretion when the masses diverge must occur when the right hand side of Eq. (A.6) becomes equal to zero, therefore (A.9)This timescale is again proportional to the growth timescale of the initial core, while the positive numerical factor in the brackets decreases as expected with increasing p.

A problem encountered when fitting numerical results to derive b and p is the degeneracy of the two parameters. It is therefore useful that the ratio of the two characteristic times, τ_run/τ_cross is a function of p only. The ratio is given as (A.10)While we cannot solve algebraically for p, we can invert the equation numerically if τ_run/τ_cross is given from a simulation, and then use Eq. (A.9) to also calculate k.

Figure A.2 shows a comparison of a numerical simulation and the corresponding calibrated semi-analytical solution for M_Z(t) and M_XY(t) obtained in this way. In the numerical simulation, M_Z,0 = 4.35 M_⊕ and f_opa = 0.002, leading to τ_cross = 2.50 Myr and τ_run = 3.57 Myr, both relative to the beginning of phase II. For τ_run, we use the moment when the limiting gas accretion rate of 0.01 M_⊕/yr is reached, but the timescales are at that moment anyway very short, so that the specific choice is not important. These numerical results lead to fit parameters p = 1.269 and b = 7.573.

The comparison shows that the calibrated semi-analytical solution provides a fair general approximation to the numerical result. The times of crossover and runaway agree by construction. Before crossover, M_XY is however underestimated, while it is overestimated afterwards. This seems to be a general property of the semi-analytical solution, and is a sign that the fit parameters are in fact not constant during the entire evolution. The lower envelope mass initially predicted by the semi-analytical solution is partially also due to the initial condition that M_XY(t) = 0. In the simulation, M_XY is in contrast about 0.1 M_⊕ at this moment.

Appendix A.2: Constant core mass:

Second, we consider the simpler case that the core does not grow, for example due to the presence of neighboring competing cores (Hubickyj et al. 2005) or the second passage through a region where the core has previously accreted the planetesimals due to outward migration followed by inward migration (Mordasini et al. 2012b).

In this case, we only have to consider Eq. (A.5). With the boundary conditions that M_XY(0) = 0, one finds (A.11)Figure A.3 shows the temporal evolution for a core mass of 5 M_⊕ and the same parameters p and k as in Fig. A.1. The comparison with this figure shows that while the growth timescales are for identical parameters p and k in principle longer for M_Z = const. than for the case with a growing core, they are in reality shorter since k is roughly one order of magnitude smaller when the core does not grow due to the reduced luminosity without planetesimal accretion (Ikoma et al. 2000).

From the solution we find that crossover occurs at (A.12)and runaway gas accretion at (A.13)The ratio of the time of runaway to the time of crossover is here (A.14)For a fixed p this ratio is larger than for the case with . In contrast to this case, we can here algebraically solve for p if τ_run/τ_cross has been determined in a numerical simulation.

Appendix A.3: Constant core, small envelope: , M_XY ≪ M_Z,0

Finally we consider the simplest case where the core does not grow, and where the envelope mass is much smaller than the core mass, as it is the case during the early phase II long before crossover. Here, we immediately find (again with M_XY(0) = 0) (A.15)

Appendix A.4: Estimates for M_XY as function of M_Z

In view of the numerical results for the envelope mass as a function of core mass for synthetic planetary populations, it is instructive to plot the envelope mass as a function of core mass for a given moment in time using the simple semi-analytical solutions. In Fig. A.4, M_XY calculated with Eqs. (A.11) and (A.15) is shown as function of M_Z,0 for t = 2 Myr, which is roughly the mean lifetime of protoplanetary disks around T Tauri stars (e.g., Haisch et al. 2001; Fedele et al. 2010). For M_XY ≪ M_Z, M_XY is proportional to , while it increases stronger and eventually diverges when M_Z,0 approaches the critical core mass for runaway M_Z,0,crit at t = 2 Myr. In the plot, we again use the parameters from Ida & Lin (2004) and Miguel et al. (2011) for τ_g. These parameters are chosen for illustrative purposes only. However, they make clear that one can gain in principle observational constraints on the parameters of τ_g by measuring the bulk composition of many low-mass, sub-critical planets with primordial envelopes, since this allows to determine the slope and absolute scaling of the M_Z − M_XY relation (see Sect. 5 for a synthetic population).

The curves showing M_XY (Eq. (A.11)) in Fig. A.4 are functionally similar to the relationship between the core and the envelope mass in the classical hydrostatic calculations of Mizuno et al. (1978) and the corresponding approximate analytical solution of Stevenson (1982): at low core masses M_XY ≪ M_Z, then M_XY increases rapidly, and beyond a certain critical core mass, there is no physical solution. It also shears the property that the smaller the opacity (and therefore k), the lower the critical mass. A significant difference is however the explicit temporal dependence in the case considered here. For a given time t, M_Z,0,crit is (A.16)Figure A.5 shows M_Z,0,crit as a function of time using the parameters of τ_g from Ikoma et al. (2000) and Ikoma & Genda (2006). The latter work uses more realistic opacities. In both cases, the cores do not accreted planetesimals, and full a ISM grain opacity is assumed.

The plot shows that without planetesimal accretion, gas runaway accretion sets in for massive cores (~ 15 M_⊕) relatively quickly after 1–2 × 10⁵ years even for full grain opacity. For lower mass cores (~ 5 M_⊕), the moment of runaway however approaches or exceeds the lifetime of the disk, so that lower opacities are needed for giant planet formation even for the case of minimal luminosity (Ikoma & Genda 2006).

The explicit temporal dependence of the critical core mass has been derived for a core that does not accrete planetesimals. In this case, the luminosity is generated by the contraction of the envelope itself, i.e., by the −T∂S/∂t term (T is the temperature, S the entropy) in the energy equation which establishes a direct dependence on the earlier state of the envelope. This is different from the situation considered by Stevenson (1982) where the luminosity is generated by the accretion of planetesimals which is a free parameter.

This is also different from phase II considered in the simulations in Sect. 3 where the core accretes. There, the dominant part of the luminosity again comes from the accretion of planetesimals but the accretion rate of planetesimals in no more a free parameter, but depends via the extension of the feeding zone on the gas accretion timescale, which in turn depends on the opacity.

Appendix B: Parameters of the growth timescale

With the semi-analytical solution of Sect. A.1 and the numerical results for the time of crossover and runaway from Sect. 3, we can systematically study the parameters p and b (or k) of the planetary growth timescale as a function of core mass and opacity for the MBPL10 comparison calculations. We thus only consider the case where the core grows due to expansion of the feeding zone. Figure B.1 shows the results as a function of f_opa over the relevant domain. In both panels, the results for the three different planetesimals surface densities Σ_s,0 = 10, 6, and 4 g/cm² are given.

The upper panel shows the exponent p. One sees that the exponent varies by about 20% for a variation of f_opa over two orders of magnitude, so that there is only a weak, but still non-negligible dependency. There is a trend of an increase in p with f_opa, bringing it to a value between 2.5 to 3.5 for full grain opacity and Σ_s,0 = 10 g/cm². This agrees with previous results considering a setup similar to the one found in Bryden et al. (2000) who fit the result of P96. The plot also shows that p is an increasing function of the mass of the core, indicating that the functional dependence of τ_g on the total mass only (Eq. (A.1)) is insufficient for a full parametrization.

In the lower panel we see that b increases with f_opa, as expected. The three thin lines approximating the numerical results are proportional to log (f_opa), meaning that k = 10^b is linearly proportional to f_opa. A linear dependence of the growth timescale on the opacity has already been found by Ikoma et al. (2000) and is found in the analytical solution of Stevenson (1982) for purely radiative envelopes. Rafikov (2006) generalized this result to envelopes that are radiative at the interface to the nebula but

convective below, as it is the case during phase II for the simulations shown here.

Table B.1 summarizes the results by giving approximate expression for p and b as a function of Σ_s,0 (and therefore the core mass) and 0.001 ≤ f_opa ≤ 0.1. These expression can be used together with the semi-analytical solution for rough estimates of the timescales of crossover and runaway. However, we caution that full evolutionary calculations are necessary to determine them under realistic conditions, where the planetesimal accretion rate varies due to orbital migration or the competition with other cores.

---

Acknowledgments

We thank K.-M. Dittkrist, P. Mollière and especially C. Ormel for enlightening discussions. This work was supported in part by the Swiss National Science Foundation and the European Research Council under grant 239605. Computations were made on the BATCHELOR cluster at MPIA. We thank two anonymous referees for important comments that led to a substantial improvement of the manuscript. C.M. thanks the Max-Planck-Gesellschaft for the Reimar Lüst Fellowship.

References

All Tables

All Figures