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 
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). 

Open with DEXTER 
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). 

Open with DEXTER 
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). 

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 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). 

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 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). 

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^{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}. 

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 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.
© ESO, 2014
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.