Issue 
A&A
Volume 566, June 2014



Article Number  A141  
Number of page(s)  22  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201321479  
Published online  27 June 2014 
Grain opacity and the bulk composition of extrasolar planets
I. Results from scaling the ISM opacity^{⋆}
^{1} MaxPlanckInstitut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
email: mordasini@mpia.de
^{2} Center for space and habitability, Physikalisches Institut, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland
^{3} Institut UTINAM, CNRSUMR 6213, Observatoire de Besançon, BP 1615, 25010 Besançon Cedex, France
^{4} Department of Astronomy and Astrophysics, University of California, Santa Cruz, USA
Received: 12 March 2013
Accepted: 10 March 2014
Context. The opacity due to grains in the envelope of a protoplanet κ_{gr} regulates the accretion rate of gas during formation, meaning that the final bulk composition of planets with a primordial H/He envelope is a function of that opacity. Observationally, for extrasolar planets with known mass and radius it is possible to estimate the bulk composition via internal structure models.
Aims. We want to study the global effects of κ_{gr} as a poorly known, but important quantity on synthetic planetary populations.
Methods. We first determine the reduction factor of the interstellar medium (ISM) grain opacity f_{opa} that leads to a gas accretion timescale consistent with grain evolution models for specific cases. In the second part we compare the mass–radius relationship of lowmass planets and the heavy element content of giant planets for different values of the reduction factor with observational constraints.
Results. For f_{opa} = 1 (full ISM opacity) the synthetic superEarth and Neptunian planets have radii that are too small (i.e., envelope masses that are too low) compared to observations because at such high opacity, they cannot efficiently accrete H/He during the formation phase. At f_{opa} = 0.003, the value calibrated with the grain evolution models, the synthetic and actual planets occupy a similar mass–radius domain. Another observable consequence is the metal enrichment of giant planets relative to the host star, Z_{pl}/Z_{star}. We find that the mean enrichment of giant planets as a function of mass M can be approximated as Z_{pl}/Z_{star} = β(M/M_{♃})^{α} both for synthetic and actual planets. The decrease in Z_{pl}/Z_{star} with mass follows α ≈ −0.7 independent of f_{opa} in synthetic populations, in agreement with the value derived from observations (−0.71 ± 0.10). The absolute enrichment level β decreases from 8.5 at f_{opa} = 1 to 3.5 at f_{opa} = 0. At f_{opa} = 0.003, one finds β = 7.2 which is similar to the result derived from observations (6.3 ± 1.0).
Conclusions. We find observational hints that the opacity in protoplanetary atmospheres is much smaller than in the ISM even if the specific value of κ_{gr} cannot be constrained in this first study as κ_{gr} is found by scaling the ISM opacity. Our results for the enrichment of giant planets are also important to distinguish core accretion and gravitational instability. In the simplest picture of core accretion, where first a critical core forms, and afterwards only gas is added, α ≈ −1. If a core accretes all planetesimals inside the feeding zone during runaway gas accretion α ≈ −2/3. The observational result (−0.71 ± 0.10) lies between these two values, pointing to core accretion as the likely formation mechanism.
Key words: planetary systems / planets and satellites: formation / planets and satellites: interiors / methods: numerical / planets and satellites: individual: Jupiter
Appendices are available in electronic form at http://www.aanda.org
© 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 socalled 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 selfconsistency as the envelope structure in which the grain growth was studied was calculated with prespecified different opacities.
Only recently Movshovitz et al. (2010, hereafter MBPL10), presented the first selfconsistently 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 selfconsistently 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 semimajor 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 lowmass 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 semianalytical 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 hightemperature gas opacities from the Bell & Lin (1994) analytical opacity laws. Second, lowtemperature molecular and atomic opacities of a grainfree, solar composition ([M/H] = 0) gas from Freedman et al. (2008).
For the scaled grain and hightemperature 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^{1}. The opacity regions are in increasing temperature order dominated by: ice grains, ice grain evaporation, metal grains, metal grain evaporation, molecules, Hscattering, 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 hightemperature 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 grainfree gas also at low temperatures. Freedman et al. (2008) provide such tables of molecular and atomic opacities of a grainfree 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)
Fig. 1 Simulations of the in situ formation of Jupiter for Σ_{s,0} = 10 g/cm^{2} as a function of the grain opacity reduction factor. The left panel shows the total mass as a function of time. The right panel is the total luminosity (including the intrinsic and the shock luminosity) in units of the intrinsic luminosity of Jupiter today. The curves correspond, from left to right, to f_{opa} = 10^{4} (black solid), 0.001 (shortdashed green), 0.01 (longdashed blue), 0.1 (dashdotted magenta), and 1 (dotted red line). The influence of f_{opa} on the duration of the plateau phase II is obvious. 
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 insitu 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^{2}.
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 planetesimalenvelope 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^{2}. 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^{2}
To illustrate the impact of κ_{gr}, we show in Fig. 1 the result for the calculations for Σ_{s,0} = 10 g/cm^{2}. 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^{2} 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 grainfree 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.
Fig. 2 Solid green curve: duration of the phase II as a function of f_{opa} for an initial planetesimal surface density of 10 g/cm^{2}. The short horizontal black line labeled M10 is the result of MBPL10, 0.517 Myr. Dotted red curve: 6 g/cm^{2} case. The short horizontal black line labeled M6 is the result of MBPL10, 0.787 Myr. Dashed blue curve: 4 g/cm^{2} case. The short horizontal black line labeled M4 is the result of MBPL10, 2.428 Myr. The dashdotted black lines in the top left and bottom right corners represent a linear slope to guide to eye. 
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^{2}. 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^{2} 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^{2}. 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 dasheddotted 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.
Grain opacity reduction factors fitting MBPL10.
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^{2}/g) but it becomes very low (κ_{gr} ~ 0.001 cm^{2}/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 postformation 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, semimajor 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.
Fig. 3 Envelope mass as a function of time for initial conditions as in run IIa in Rogers et al. (2011). The thick dotted line is the result of Rogers et al. (2011) obtained in a combined grain evolution and core accretion simulation. The blue lines are our results for different values of the grain opacity reduction factor. The lines are f_{opa} = 0.02 (solid), 0.003 (dotted), 0.001 (dashed) and 0 (dotdashed). 
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 lowmass 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^{2} 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 postformation 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 semimajor 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^{3}.
Fig. 4 Envelope mass M_{XY} as a function of core mass M_{Z} for in situ simulations. The green, blue, and red lines correspond to isolation masses of about 3, 5.5, and 13 M_{⊕}, respectively. For the latter two, different f_{opa} are shown, with smaller values of f_{opa} leading to larger M_{XY} in phase I. The two blue lines show f_{opa} = 0.003 and 0.1, while the three red lines show f_{opa} = 1, 0.1, and 0.01. In phase II, the f_{opa} lines converge. Thin black lines (partially covered) show Eq. (4). The thick gray line shows the crossover point where M_{XY} = M_{Z}. 
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κ).
Settings for the four synthetic populations.
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 semianalytical solutions for the core and envelope mass in phase II. This is presented in Appendix A while in Appendix B these semianalytical 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 lowmass lowdensity planets like the planets around Kepler11 (Lissauer et al. 2011). The amount of primordial H/He these lowmass 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 solarlike 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 nonisothermal migration introduced in Mordasini et al. (2012b). These three populations only differ by f_{opa} which was set to 0 (grainfree 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 semianalytical 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 semianalytical models. Figure 5 shows the envelope mass as a function of core mass as found numerically in the four populations. Planets with a semimajor 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.
Fig. 5 H/He envelope mass M_{XY} as a function of heavy element mass M_{Z} for synthetic planets with semimajor axes 0.11 ≤ a/AU ≤ 5. The populations in the top panels and the bottom left panel only differ by the value of f_{opa} as indicated in the plots. The populations were calculated assuming full nonisothermal type I migration rates. The population in the bottom right was calculated with f_{opa} = 0.003, but without orbital migration and assuming Ṁ_{Z,run} = 0. The lines show simple analytical estimates and scalings (see text). 
5.1. Lowmass planets
We see that f_{opa} leaves a clear imprint in the M_{Z} − M_{XY} relation. For lowmass 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 semimajor 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 substructures 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 lowmass 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 lowmass 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^{4}, 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 grainfree 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 semimajor axis that is considered. If synthetic planets at all semimajor 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 semimajor 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 semimajor 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 HATP2b (Baraffe et al. 2008). These extreme cores can only form in solid rich disk, i.e., when both the disk gas mass and the dusttogas ratio (i.e., [Fe/H]) come from the upper end of the probability distributions. The product of the dusttogas 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^{3} similar to Jupiter nowadays, and 0.125 g/cm^{3} as typical for a young Jovian planet that has a radius twice as large as Jupiter. For the 1 and 0.125 g/cm^{3} 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).
Fig. 6 Capture radius for 100 km planetesimals as a function of mass for a growing giant planet. The simulation uses the same initial conditions as the case J1 in P96. The capture radius is the thick blue line, while the green line is the total radius. 
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}
Fig. 7 Comparison of the mass–radius relationship of synthetic planets at an age of 5 Gyr, the solar system planets, and extrasolar planets with a> 0.1 AU and well constrained mass and radius. Synthetic planets have 0.1 <a/AU < 5. The colors show the fraction of heavy elements Z = M_{Z}/M in a synthetic planet. Orange: Z ≤ 1%. Red: 1 <Z ≤ 5%. Green: 5 <Z ≤ 20%. Blue: 20 <Z ≤ 40%. Cyan: 40 <Z ≤ 60%. Magenta: 60 <Z ≤ 80%. Yellow: 80 <Z ≤ 95%. Brown: 95 <Z ≤ 99%. Black: Z> 99%. 
The last section shows that f_{opa} has important consequences for the bulk composition of planets. The envelope mass of lowmass planets in particular increases by (very roughly speaking) one order of magnitude if κ decreases from the full ISM to grainfree. Since the radius of solid lowmass 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 selfconsistently 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 lowmass 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 lowmass, 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 (Kepler11e and Kepler18d) that lie above the envelope covered with synthetic planets. For Kepler11e, the distance from the synthetic population is not very large given the error bars, but Kepler18d lies several σ in mass and radius away. Since both Kepler11 and Kepler18 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 lowmass 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 solartype stars. Given the small eccentricities also internal heating due to eccentricity damping seems unlikely, at least for Kepler18d (Cochran et al. 2011).
The population with f_{opa} = 0 (top left) behaves clearly differently. Here, the planets Kepler11e and Kepler18d 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 lowmass 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 subpopulation of lowmass, very large (i.e., very gas rich) planets presently without observational counterparts. This could mean that with grainfree opacities, gas accretion in the model is too efficient, and that a small, but still nonzero 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. Kepler18d 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 semimajor 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 semimajor 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 lowmass 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 CoRoT10b and HD 80606b lie slightly below the envelope of points predicted by this synthetic population. Neither CoRoT10 nor HD 80606 are much older than 5 Gyr. For CoRoT10, 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 CoRoT10b 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 lowmass 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 lowmass 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 readdressed 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 thindisk 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^{8} erg/(s cm^{2}). 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 solarlike 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 heavyelement 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 semiconvective 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 dusttogas 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_{♃}.
Fig. 8 Planetary heavy element fraction Z_{pl} divided by the heavy element fraction of the host star Z_{star}. The relative enrichment is shown as a function of planetary mass M. Red dots with error bars are from Miller & Fortney (2011) for transiting exoplanets. The solid red line is a leastsquares fit to the exoplanet data, while the dotted red lines give the 1σ errors in the fit. Giant planets of the solar system are shown with black symbols. Blue dots are synthetic planets with 0.11 ≤ a/AU ≤ 1 (5 AU for the synthesis without migration). The blue line gives the fit to the synthetic data. The four panels show populations calculated with an f_{opa} indicated in the plot. The bottom right panel assumes in situ formation and Ṁ_{Z,run} = 0. The pale green dots are the same synthetic planets, but assuming that the accreted gas has the same metallicity as the host star instead of being metal free. 
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 leastsquares 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 grainfree 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 lowmass 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 nonzero 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 dusttogas 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 selfconsistently 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_{1}/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_{1} is 10–15 M_{♃}. These values are a function of f_{opa}: at f_{opa} = 0.0, the lowest M_{1} is about 1 M_{♃}, while the typical value is about 5 M_{♃}. At f_{opa} = 1, the lowest M_{1} 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_{1} is here about 1 M_{♃}, and the typical value is about 2.5 M_{♃}. These are clearly lower values. For the nominal populations, M_{1} 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 micrometersized 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 selfconsistent 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 envelopedominated 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 precollapse phase depends on opacity. If grain growth is taken into account (Helled & Bodenheimer 2011), the duration of the precollapse phase is very short (~ 10^{3} 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_{opa}relative 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 grainfree 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 subcritical 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 subcritical 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 semimajor axis (Sect. 5.3).

Grain opacity leaves a clear imprint in the mass–radius relationship of lowmass planets. At a low opacity, planets in the superEarth 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 superEarth 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 grainfree 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 grainfree 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 nonzero contribution from the grains (Sect. 7.3). This result is similar as the one found for the mass–radius relationship of lowmass 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 semianalytical 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 lowmass, lowdensity 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 lowmass, 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 wellconstraint 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).
See Cuzzi et al. (2014) for a more recent work on grain opacities in protoplanetary atmospheres.
It is far from clear that phase II actually occurs at all during giant planet formation: low planetesimal accretion rates (Fortier et al. 2007), gravitational shepherding (Zhou & Lin 2007), gap formation in the planetesimal disk (Shiraishi & Ida 2008), and orbital migration (Alibert et al. 2005) all tend to suppress it. In any case, in this work we are interested in it only to calibrate the magnitude of κ_{gr}.
The timescale of gas accretion and the associated expansion of the feeding zone becomes very short at very low f_{opa}, so that some planetesimals in the feeding zone are not accreted. This is the explanation why in the Σ_{s,0} = 10 g/cm^{2} simulation, the crossover mass at f_{opa} = 10^{5} is about 2 M_{⊕} smaller than for f_{opa} = 1.
The equations remain initially also valid in the runaway gas accretion phase as long as the extension of the feeding zone is the limiting factor for Ṁ_{Z}. At some moment, however, the planet’s envelope has collapsed so much that the capture radius for planetesimals becomes small. Concurrently, the ejection of planetesimals becomes important (depending on the semimajor axis of the planet). Then, Ṁ_{Z} diverges from Eq. (A.4) and eventually falls to zero (see Mordasini et al. 2012b).
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 MaxPlanckGesellschaft for the Reimar Lüst Fellowship.
References
 Adams, E. R., Seager, S., & ElkinsTanton, L. 2008, ApJ, 673, 1160 [NASA ADS] [CrossRef] [Google Scholar]
 Adibekyan, V. Z., Santos, N. C., Sousa, S. G., et al. 2012, A&A, 543, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alibert, Y., Mordasini, C., & Benz, W. 2004, A&A, 417, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Baraffe, I., Chabrier, G., & Barman, T. S. 2008, A&A, 482, 315 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bodenheimer, P. H., & Pollack, J. B. 1986, Icarus, 67, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Bodenheimer, P., D’Angelo, G., Lissauer, J. J., Fortney, J. J., & Saumon, D. 2013, ApJ, 770, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Boley, A. C., Helled, R., & Payne, M. J. 2011, ApJ, 735, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Bonomo, A. S., Santerne, A., Alonso, R., et al. 2010, A&A, 520, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, EpJ Web Conf., 47, 3005 [Google Scholar]
 Bryden, G., Lin, D. N. C., & Ida, S. 2000, ApJ, 544, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, A. S., Marley, M. S., Hubbard, W. B., et al. 1997, ApJ, 491, 856 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, A., Hubeny, I., Budaj, J., & Hubbard, W. B. 2007, ApJ, 661, 502 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, A., Heng, K., & Nampaisarn, T. 2011, ApJ, 736, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Cochran, W. D., Fabrycky, D. C., Torres, G., et al. 2011, ApJS, 197, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Cuzzi, J. N., Estrada, P. R., & Davis, S. S. 2014, ApJS, 210, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Demory, B.O., & Seager, S. 2011, ApJS, 197, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Dittkrist, K.M., Mordasini, C., Klahr, H., Alibert, Y., & Henning, T. 2014, A&A, in press, DOI: 10.1051/00046361/201322506 [Google Scholar]
 Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, R., & Oliveira, J. M. 2010, A&A, 510, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortier, A., Benvenuto, O. G., & Brunini, A. 2007, A&A, 473, 311 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortier, A., Alibert, Y., Carron, F., Benz, W., & Dittkrist, K.M. 2013, A&A, 549, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Ikoma, M., Nettelmann, N., Guillot, T., & Marley, M. S. 2011, ApJ, 729, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T. 1999, Science, 286, 72 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Guillot, T. 2008, Phys. Scr., 130, 014023 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T., & Gautier, D. 2009, in Treatise on Geophysics, eds. G. Shubert, & T. Spohn, 10, 439 [Google Scholar]
 Guillot, T., Santos, N. C., Pont, F., et al. 2006, A&A, 453, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haisch, K. E., Jr., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
 Helled, R., & Bodenheimer, P. 2011, Icarus, 211, 939 [NASA ADS] [CrossRef] [Google Scholar]
 Helled, R., & Schubert, G. 2009, ApJ, 697, 1256 [NASA ADS] [CrossRef] [Google Scholar]
 Hori, Y., & Ikoma, M. 2010, ApJ, 714, 1343 [NASA ADS] [CrossRef] [Google Scholar]
 Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419 [NASA ADS] [CrossRef] [Google Scholar]
 Hubickyj, O., Bodenheimer, P. H., & Lissauer, J. J. 2005, Icarus, 179, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Ida, S., & Lin, D.N.C. 2004, ApJ, 604, 388 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Ida, S., & Makino, J. 1993, Icarus, 106, 210 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Ikoma, M., & Genda, H. 2006, ApJ, 648, 696 [NASA ADS] [CrossRef] [Google Scholar]
 Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., & Dirksen, G. 2006, A&A, 447, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kornet, K., Różyczka, M., & Stepinski, T. F. 2004, A&A, 417, 151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leconte, J., & Chabrier, G. 2012, A&A, 540, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levison, H. F., Thommes, E. W., & Duncan, M. J. 2010, AJ, 139, 1297 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J. 1993, ARA&A, 31, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Hubickyj, O., D’Angelo, G., & Bodenheimer, P. H. 2009, Icarus, 199, 338 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011, Nature, 470, 53 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Lyra, W., Paardekooper, S.J., & Mac Low, M.M. 2010, ApJ, 715, L68 [NASA ADS] [CrossRef] [Google Scholar]
 Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Marmier, M., Lovis, C., et al. 2011, A&A, submitted [arXiv:1109.2497] [Google Scholar]
 Miguel, Y., & Brunini, A. 2008, MNRAS, 387, 463 [NASA ADS] [CrossRef] [Google Scholar]
 Miguel, Y., Guilera, O. M., & Brunini, A. 2011, MNRAS, 412, 2113 [NASA ADS] [CrossRef] [Google Scholar]
 Militzer, B., Hubbard, W. B., Vorberger, J., Tamblyn, I., & Bonev, S. A. 2008, ApJ, 688, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, N., & Fortney, J. J. 2011, ApJ, 736, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Mizuno, H., Nakazawa, K., & Hayashi, C. 1978, Prog. Theor. Phys., 60, 699 [Google Scholar]
 Mollière, P., & Mordasini, C. 2012, A&A, 547, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C. 2013, A&A, 558, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C. 2014, A&A, submitted [Google Scholar]
 Mordasini, C., Alibert, Y., & Benz, W. 2005, in Tenth Anniversary of 51 Pegb: Status of and prospects for hot Jupiter studies, eds. L. Arnold, F. Bouchy, & C. Moutou, 84 [Google Scholar]
 Mordasini, C., Alibert, Y., & Benz, W. 2009a, A&A, 501, 1139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C., Alibert, Y., Benz, W., & Naef, D. 2009b, A&A, 501, 1161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012a, A&A, 541, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012b, A&A, 547, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mordasini, C., Alibert, Y., Georgy, C., et al. 2012c, A&A, 547, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Movshovitz, N., & Podolak, M. 2008, Icarus, 194, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Movshovitz, N., Bodenheimer, P. H., Podolak, M., & Lissauer, J. J. 2010, Icarus, 209, 616 (MBPL10) [NASA ADS] [CrossRef] [Google Scholar]
 Nayakshin, S. 2010, MNRAS, 408, 2381 [NASA ADS] [CrossRef] [Google Scholar]
 Nettelmann, N., Becker, A., Holst, B., & Redmer, R. 2012, ApJ, 750, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Kobayashi, H. 2012, ApJ, 747, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105 [Google Scholar]
 Papaloizou, J. C. B., & Terquem, C. 1999, ApJ, 521, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Podolak, M. 2003, Icarus, 165, 428 [NASA ADS] [CrossRef] [Google Scholar]
 Podolak, M., Pollack, J. B., & Reynolds, R. T. 1988, Icarus, 73, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 (P96) [NASA ADS] [CrossRef] [Google Scholar]
 Pont, F., Hébrard, G., Irwin, J. M., et al. 2009, A&A, 502, 695 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rafikov, R. R. 2006, ApJ, 648, 666 [NASA ADS] [CrossRef] [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2013, Exp. Astron., accepted [arXiv:1310.0696] [Google Scholar]
 Rogers, L. A., & Seager, S. 2010, ApJ, 712, 974 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager, S. 2011, ApJ, 738, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Rossow, W. B. 1978, Icarus, 36, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Saumon, D., & Guillot, T. 2004, ApJ, 609, 1170 [NASA ADS] [CrossRef] [Google Scholar]
 Saumon, D., Chabrier, G., & Van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Sheng, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, submitted [Google Scholar]
 Shiraishi, M., & Ida, S. 2008, ApJ, 684, 1416 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, D. J. 1982, Planet. Space Sci., 30, 755 [NASA ADS] [CrossRef] [Google Scholar]
 Valencia, D., Ikoma, M., Guillot, T., & Nettelmann, N. 2010, A&A, 516, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhou, J.L., & Lin, D. N. C. 2007, ApJ, 666, 447 [NASA ADS] [CrossRef] [Google Scholar]
Online material
Appendix A: Semianalytical 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 selfsimilar in the sense that they can be approximately converted in each other by stretching the xaxis. 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 semianalytical 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_{0} 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 semianalytical solution presented below and the numerical simulations makes it possible to determine the parameters in a more systematic way.
Literature values of parameters for τ_{g} = 10^{b}M^{− p}.
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^{5} 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.
Fig. A.1 Semianalytical solution for the temporal evolution of the core mass (dashed blue), envelope mass (dotted red), and total mass (solid black lines). Thin lines use the parameters of Ida & Lin (2004), while thick lines use those of Miguel et al. (2011). 
The semianalytical 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 semianalytical 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.
Fig. A.2 Numerical result (thick lines) and corresponding calibrated semianalytical solution (thin lines) for the temporal evolution of the core mass (dashed blue), envelope mass (dotted red), and total mass (solid black lines). 
The comparison shows that the calibrated semianalytical 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 semianalytical 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 semianalytical 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).
Fig. A.3 Semianalytical solution for the temporal evolution of the envelope mass (dotted red) for the case that the core mass is constant (dashed blue line). The total mass is shown with solid black lines. The thin lines use the parameters from Ida & Lin (2004), while the thick lines use those of Miguel et al. (2011). 
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}
Fig. A.4 Illustration of M_{XY}(M_{Z,0}) at t = 2 Myr. The solid black lines show envelope masses as a function of core mass for the constant core mass case, while the dashed blue lines are the approximation for M_{XY} ≪ M_{Z,0}. The thin lines use the parameters from Ida & Lin (2004), while the thick lines use those of Miguel et al. (2011). 
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 semianalytical 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 lowmass, subcritical 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.
Fig. A.5 Mass of a core M_{Z,0,crit} that undergoes gas runaway accretion after being embedded in the nebula without planetesimal accretion during a time t. The solid black line uses the parameters for τ_{g} from Ikoma et al. (2000) while the dashed blue line is for the updated parameters of Ikoma & Genda (2006). 
The plot shows that without planetesimal accretion, gas runaway accretion sets in for massive cores (~ 15 M_{⊕}) relatively quickly after 1–2 × 10^{5} 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 semianalytical 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^{2} are given.
Fig. B.1 Parameters for the planetary growth timescale (Eq. (A.1)) obtained from the comparison of numerical simulations (Sect. 3) with the semianalytical solution of Sect. A.1. The upper panel shows the exponent p as a function of f_{opa}. The initial surface density of planetesimals is 10 (top, black), 6 (middle, red) and 4 g/cm^{2} (bottom, blue line). The lower panel shows b, again as function of f_{opa} and the three planetesimal surface densities. Thick lines are the numerical results, while thin lines give the approximative linear scaling with f_{opa}. 
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 nonnegligible 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^{2}. 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 semianalytical 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.
Parameters for τ_{g} = 10^{b}M^{− p} of planets in phase II for 0.001 ≤ f_{opa} ≤ 0.1.
All Tables
Parameters for τ_{g} = 10^{b}M^{− p} of planets in phase II for 0.001 ≤ f_{opa} ≤ 0.1.
All Figures
Fig. 1 Simulations of the in situ formation of Jupiter for Σ_{s,0} = 10 g/cm^{2} as a function of the grain opacity reduction factor. The left panel shows the total mass as a function of time. The right panel is the total luminosity (including the intrinsic and the shock luminosity) in units of the intrinsic luminosity of Jupiter today. The curves correspond, from left to right, to f_{opa} = 10^{4} (black solid), 0.001 (shortdashed green), 0.01 (longdashed blue), 0.1 (dashdotted magenta), and 1 (dotted red line). The influence of f_{opa} on the duration of the plateau phase II is obvious. 

In the text 
Fig. 2 Solid green curve: duration of the phase II as a function of f_{opa} for an initial planetesimal surface density of 10 g/cm^{2}. The short horizontal black line labeled M10 is the result of MBPL10, 0.517 Myr. Dotted red curve: 6 g/cm^{2} case. The short horizontal black line labeled M6 is the result of MBPL10, 0.787 Myr. Dashed blue curve: 4 g/cm^{2} case. The short horizontal black line labeled M4 is the result of MBPL10, 2.428 Myr. The dashdotted black lines in the top left and bottom right corners represent a linear slope to guide to eye. 

In the text 
Fig. 3 Envelope mass as a function of time for initial conditions as in run IIa in Rogers et al. (2011). The thick dotted line is the result of Rogers et al. (2011) obtained in a combined grain evolution and core accretion simulation. The blue lines are our results for different values of the grain opacity reduction factor. The lines are f_{opa} = 0.02 (solid), 0.003 (dotted), 0.001 (dashed) and 0 (dotdashed). 

In the text 
Fig. 4 Envelope mass M_{XY} as a function of core mass M_{Z} for in situ simulations. The green, blue, and red lines correspond to isolation masses of about 3, 5.5, and 13 M_{⊕}, respectively. For the latter two, different f_{opa} are shown, with smaller values of f_{opa} leading to larger M_{XY} in phase I. The two blue lines show f_{opa} = 0.003 and 0.1, while the three red lines show f_{opa} = 1, 0.1, and 0.01. In phase II, the f_{opa} lines converge. Thin black lines (partially covered) show Eq. (4). The thick gray line shows the crossover point where M_{XY} = M_{Z}. 

In the text 
Fig. 5 H/He envelope mass M_{XY} as a function of heavy element mass M_{Z} for synthetic planets with semimajor axes 0.11 ≤ a/AU ≤ 5. The populations in the top panels and the bottom left panel only differ by the value of f_{opa} as indicated in the plots. The populations were calculated assuming full nonisothermal type I migration rates. The population in the bottom right was calculated with f_{opa} = 0.003, but without orbital migration and assuming Ṁ_{Z,run} = 0. The lines show simple analytical estimates and scalings (see text). 

In the text 
Fig. 6 Capture radius for 100 km planetesimals as a function of mass for a growing giant planet. The simulation uses the same initial conditions as the case J1 in P96. The capture radius is the thick blue line, while the green line is the total radius. 

In the text 
Fig. 7 Comparison of the mass–radius relationship of synthetic planets at an age of 5 Gyr, the solar system planets, and extrasolar planets with a> 0.1 AU and well constrained mass and radius. Synthetic planets have 0.1 <a/AU < 5. The colors show the fraction of heavy elements Z = M_{Z}/M in a synthetic planet. Orange: Z ≤ 1%. Red: 1 <Z ≤ 5%. Green: 5 <Z ≤ 20%. Blue: 20 <Z ≤ 40%. Cyan: 40 <Z ≤ 60%. Magenta: 60 <Z ≤ 80%. Yellow: 80 <Z ≤ 95%. Brown: 95 <Z ≤ 99%. Black: Z> 99%. 

In the text 
Fig. 8 Planetary heavy element fraction Z_{pl} divided by the heavy element fraction of the host star Z_{star}. The relative enrichment is shown as a function of planetary mass M. Red dots with error bars are from Miller & Fortney (2011) for transiting exoplanets. The solid red line is a leastsquares fit to the exoplanet data, while the dotted red lines give the 1σ errors in the fit. Giant planets of the solar system are shown with black symbols. Blue dots are synthetic planets with 0.11 ≤ a/AU ≤ 1 (5 AU for the synthesis without migration). The blue line gives the fit to the synthetic data. The four panels show populations calculated with an f_{opa} indicated in the plot. The bottom right panel assumes in situ formation and Ṁ_{Z,run} = 0. The pale green dots are the same synthetic planets, but assuming that the accreted gas has the same metallicity as the host star instead of being metal free. 

In the text 
Fig. A.1 Semianalytical solution for the temporal evolution of the core mass (dashed blue), envelope mass (dotted red), and total mass (solid black lines). Thin lines use the parameters of Ida & Lin (2004), while thick lines use those of Miguel et al. (2011). 

In the text 
Fig. A.2 Numerical result (thick lines) and corresponding calibrated semianalytical solution (thin lines) for the temporal evolution of the core mass (dashed blue), envelope mass (dotted red), and total mass (solid black lines). 

In the text 
Fig. A.3 Semianalytical solution for the temporal evolution of the envelope mass (dotted red) for the case that the core mass is constant (dashed blue line). The total mass is shown with solid black lines. The thin lines use the parameters from Ida & Lin (2004), while the thick lines use those of Miguel et al. (2011). 

In the text 
Fig. A.4 Illustration of M_{XY}(M_{Z,0}) at t = 2 Myr. The solid black lines show envelope masses as a function of core mass for the constant core mass case, while the dashed blue lines are the approximation for M_{XY} ≪ M_{Z,0}. The thin lines use the parameters from Ida & Lin (2004), while the thick lines use those of Miguel et al. (2011). 

In the text 
Fig. A.5 Mass of a core M_{Z,0,crit} that undergoes gas runaway accretion after being embedded in the nebula without planetesimal accretion during a time t. The solid black line uses the parameters for τ_{g} from Ikoma et al. (2000) while the dashed blue line is for the updated parameters of Ikoma & Genda (2006). 

In the text 
Fig. B.1 Parameters for the planetary growth timescale (Eq. (A.1)) obtained from the comparison of numerical simulations (Sect. 3) with the semianalytical solution of Sect. A.1. The upper panel shows the exponent p as a function of f_{opa}. The initial surface density of planetesimals is 10 (top, black), 6 (middle, red) and 4 g/cm^{2} (bottom, blue line). The lower panel shows b, again as function of f_{opa} and the three planetesimal surface densities. Thick lines are the numerical results, while thin lines give the approximative linear scaling with f_{opa}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.