Online material

Appendix A: Semi-analytical solutions

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

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

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

Appendix A.1: Growing core mass:

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

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

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

	Fig. A.1 Semi-analytical 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).
Open with DEXTER

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

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

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

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

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

	Fig. A.2 Numerical result (thick lines) and corresponding calibrated semi-analytical solution (thin lines) for the temporal evolution of the core mass (dashed blue), envelope mass (dotted red), and total mass (solid black lines).
Open with DEXTER

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

Appendix A.2: Constant core mass:

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

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

	Fig. A.3 Semi-analytical 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).
Open with DEXTER

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 MXY(MZ,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 MXY ≪ MZ,0. The thin lines use the parameters from Ida & Lin (2004), while the thick lines use those of Miguel et al. (2011).
Open with DEXTER

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

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

	Fig. A.5 Mass of a core MZ,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).
Open with DEXTER

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

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

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

Appendix B: Parameters of the growth timescale

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

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

Open with DEXTER

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

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

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

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

© ESO, 2014