Free Access
Volume 566, June 2014
Article Number A141
Number of page(s) 22
Section Planets and planetary systems
Published online 27 June 2014

Online material

Appendix A: Semi-analytical solutions

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

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

The gas accretion rate is parametrized with the growth timescale τg of a planet of total mass M = MZ + MXY where MZ is the total mass of heavy elements in the planet (equal to the core mass if all planetesimals sink to the core) while MXY 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 = 10b 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 fopa. The explicit dependency of k on fopa can be merged into the parameter, writing (at least over a certain range of fopa), where k0 and q are a priori unknown (in Sect. B we will see that q ≈ 1 over an important domain of fopa). 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.

Table A.1

Literature values of parameters for τg = 10bMp.

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 M1/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 MZ(t) as a function of time t (P96) (A.4)which holds very well during phase II in the numerical simulations above5 This equation is also found from differentiation of Eq. (4). For the envelope mass MXY(t) we have (A.5)As boundary conditions we use MZ(0) = MZ,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 MXY(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 MZ(t) and MXY(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 MZ,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.

thumbnail 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 MZ = MXY, and of runaway τrun when MZ and MXY 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 MZ(t) and MXY(t) obtained in this way. In the numerical simulation, MZ,0 = 4.35 M and fopa = 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.

thumbnail 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, MXY 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 MXY(t) = 0. In the simulation, MXY 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 MXY(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 MZ = 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).

thumbnail 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: , MXYMZ,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 MXY(0) = 0) (A.15)

Appendix A.4: Estimates for MXY as function of MZ

thumbnail 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 MXYMZ,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, MXY calculated with Eqs. (A.11) and (A.15) is shown as function of MZ,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 MXYMZ, MXY is proportional to , while it increases stronger and eventually diverges when MZ,0 approaches the critical core mass for runaway MZ,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 MZMXY relation (see Sect. 5 for a synthetic population).

The curves showing MXY (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 MXYMZ, then MXY 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, MZ,0,crit is (A.16)Figure A.5 shows MZ,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.

thumbnail 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 × 105 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 fopa over the relevant domain. In both panels, the results for the three different planetesimals surface densities Σs,0 = 10, 6, and 4 g/cm2 are given.

thumbnail 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 fopa 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 fopa, bringing it to a value between 2.5 to 3.5 for full grain opacity and Σs,0 = 10 g/cm2. 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 fopa, as expected. The three thin lines approximating the numerical results are proportional to log (fopa), meaning that k = 10b is linearly proportional to fopa. 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 ≤ fopa ≤ 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.

Table B.1

Parameters for τg = 10bMp of planets in phase II for 0.001 ≤ fopa ≤ 0.1.

© ESO, 2014

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.