IN SITU AND EX SITU FORMATION MODELS OF KEPLER 11 PLANETS

The calculation of the structure and evolution of the planetary gaseous envelope is based on the assumption that the envelope is spherically symmetric around its center and evolves through states of hydrostatic equilibrium (e.g., Kippenhahn et al. 2013). The envelope lies on a core of condensed matter, whose total mass M_c and radius R_c are both functions of time. The core radius is determined as explained below. The envelope structure is calculated by solving the equations for mass conservation, hydrostatic equilibrium, energy conservation, and radiation diffusion (see Bodenheimer & Pollack 1986). The energy equation includes heating produced by in-falling planetesimals, the work done by gravity, cooling from the release of internal heat, and heating by stellar radiation (when applicable). In convective unstable shells, where the radiative temperature gradient exceeds in magnitude the adiabatic gradient (Kippenhahn et al. 2013), the actual gradient of temperature is set equal to the adiabatic gradient.

The chemical composition of the envelope gas is assumed to be uniform, with hydrogen and helium mass fractions X = 0.74 and Y = 0.24, respectively. The equation of state for this gas mixture is that computed by Saumon et al. (1995), which accounts for the partial degeneracy of electrons and for non-ideal effects in the gas. (Strictly speaking, this equation of state neglects heavier elements and uses $Y=1-X$.)

The envelope opacity arises from the combined contributions of dust, atoms, and molecules. Dust opacities are calculated as described in D'Angelo & Bodenheimer (2013), assuming the presence of a number of different grain species, and grain size distributions with a minimum radius of $0.005\,\mu {\rm{m}}$ and a maximum radius of 1 or $10\,\mathrm{mm}$. The number density of dust grains is proportional to the −3 power of the grain radius. Low-temperature gas opacities are taken from Ferguson et al. (2005), whereas high-temperature gas opacities are taken from the OPAL tables (Iglesias & Rogers 1996). Figure 1 illustrates the Rosseland mean opacity for the two grain size distributions considered here, along with gas opacity. Dust grains are supposed to be present in the envelope only if solids are supplied via gas and/or planetesimal accretion. In the absence of a steady supply, dust grains quickly sediment to deeper layers and evaporate. When this happens, i.e., during phases of zero gas and solids' accretion, low-temperature opacities are replaced with the molecular opacities of Freedman et al. (2008).

Zoom In Zoom Out Reset image size

The planetary evolution code is largely the same as that used by Pollack et al. (1996), Bodenheimer et al. (2000), Hubickyj et al. (2005), and Lissauer et al. (2009). The structure equations are solved by means of the Henyey method (e.g., Bodenheimer et al. 2006), supplemented with boundary conditions at the core-envelope interface, R_c, and at the top of the envelope, R_p. At R_c, the mass is set equal to M_c and the luminosity is set to zero (i.e., there is no energy flux through the core-envelope interface, see the discussion in Bodenheimer & Lissauer 2014). The planet radius, R_p, is assumed to have an upper bound ${R}_{{\rm{A}}}$, the accretion radius, defined by

$$\begin{eqnarray}&&\displaystyle \frac{1}{{R}_{{\rm{A}}}}\equiv \displaystyle \frac{4}{{R}_{{\rm{H}}}}+\displaystyle \frac{1}{{R}_{{\rm{B}}}},\end{eqnarray} \tag{ 1 }$$

where ${R}_{{\rm{H}}}$ and ${R}_{{\rm{B}}}$ are the Hill and Bondi radius, respectively. The limiting envelope radii, ${R}_{{\rm{A}}}\approx {R}_{{\rm{H}}}/4$ and ${R}_{{\rm{A}}}\approx {R}_{{\rm{B}}}$, were estimated from first principles by means of three-dimensional (3D) calculations of the flow dynamics around planets in disks (Lissauer et al. 2009; D'Angelo & Bodenheimer 2013).

At early stages of evolution, the planet is "in contact" with the disk (the so-called nebular stage), and the density and temperature at the top of the envelope are taken as the local disk values. During the nebular stage, gas is added to the envelope in order to restore the condition ${R}_{p}\approx {R}_{{\rm{A}}}$ (see Pollack et al. 1996). Once this condition can no longer be maintained, the envelope contracts inside ${R}_{{\rm{A}}}$, entering a transition stage, which coincides with the run-away gas accretion phase if the gas accretion rate is sufficiently large. In fact, the gas accretion rate of the envelope, ${\dot{M}}_{e}$, is limited by the maximum rate at which the disk can deliver gas to the planet's vicinity. If tidal perturbations by the planet are negligible, the disk-limited accretion rate can be described in terms of simple analytical arguments (D'Angelo & Lubow 2008). If tidal perturbations are not negligible, the problem becomes highly complex and depends on the interplay among viscous torques, tidal torques, and close-range flow dynamics around the planet. Disk-limited accretion rates were derived from 3D hydrodynamics high-resolution calculations, as described in Bodenheimer et al. (2013), extending the parameter space covered by the fitting functions reported therein to include a larger disk viscosity range. It is important to stress that ${\dot{M}}_{e}$ cannot exceed the disk-limited accretion rate, and there are instances in which this limit sets in during the nebular stage of evolution. The boundary conditions applied at R_p during the transition stage are discussed in Bodenheimer et al. (2000).

Eventually, the disk's gas around the planet's orbit is dispersed, typically via photo-evaporation by stellar radiation, gas supply ceases, and the planet enters the isolation stage. Note that the process of gap formation in the disk by tidal torques alone, under typical disk conditions of viscosity and temperature and for planetary masses up to several times Jupiter's mass (${{\rm{M}}}_{{\rm{J}}}$), does not lead to isolation (e.g., D'Angelo & Lubow 2008; Bodenheimer et al. 2013, and references therein). During the isolation stage, standard photospheric boundary conditions are applied at R_p (e.g., Cox 1968)

$$\begin{eqnarray}&&{L}_{p}=4\pi {\sigma }_{\mathrm{SB}}{R}_{p}^{2}{T}_{\mathrm{eff}}^{4}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\kappa }_{{\rm{R}}}P=\displaystyle \frac{2}{3}\displaystyle \frac{{{GM}}_{p}}{{R}_{p}^{2}},\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{\mathrm{SB}}$ is the Stefan–Boltzmann constant, G is the gravitational constant, and ${\kappa }_{{\rm{R}}}$ and P are, respectively, the photospheric values of the Rosseland mean opacity and pressure. In Equation (2), the luminosity on the left-hand side comprises both the internal power generated by the planet and the re-radiated power that arises from the absorption of stellar radiation. Hence, the planet's effective temperature is given by (Bodenheimer & Lissauer 2014)

$$\begin{eqnarray}&&{T}_{\mathrm{eff}}^{4}={T}_{\mathrm{int}}^{4}+{T}_{\mathrm{eq}}^{4},\end{eqnarray} \tag{ 4 }$$

where ${T}_{\mathrm{int}}$ depends on R_p and on the luminosity internally generated by the planet. The equilibrium temperature is such that (e.g., Guillot 2010)

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}^{4}=\displaystyle \frac{1}{4}{T}_{\star }^{4}{\left(\displaystyle \frac{{R}_{\star }}{a}\right)}^{2}(1-{A}_{b}),\end{eqnarray} \tag{ 5 }$$

which assumes full redistribution of the incident radiation. The albedo A_b is taken as a constant equal to 0.1. In Equation (5), T_⋆ and R_⋆ are the effective temperature and the photospheric radius of the star, respectively.

In calculations that allow for orbital migration, a distinction can be made between isolation from the disk's gas and from the disk's solids. Isolation from the planetesimals' disk occurs when the migration speed ${da}/{dt}\approx 0$ and ${\dot{M}}_{p}={\dot{M}}_{c}+{\dot{M}}_{e}\lesssim 0$, or when the solids' surface density at the planet's location becomes very small. Isolation from the solid disk occurs prior to isolation from the disk's gas or shortly thereafter. However, the consequences of this delay are negligible. Thus, when a planet enters the isolation phase, it is basically isolated from both the disk's gas and solids.

During all stages, if ${R}_{p}\gt {R}_{{\rm{A}}}$, mass is gradually removed from the envelope. In a more realistic context, during the nebular and transition stages, an inflated planet would lose unbound mass hydrodynamically, carried away by the surrounding disk flow (D'Angelo & Bodenheimer 2013). This mechanism of mass loss is different from those operating during the isolation stage (see Section 2.3). As a result of gas loss, while M_c is a monotonic function of time, M_e and ${M}_{p}$ may not be.

Accretion of solids is treated as in Pollack et al. (1996). All solids accreted by the planet are assumed to sink to the core, in a condensed form, and increment the core mass. As originally derived by Safronov (1969), the accretion rate of solids can be written as

$$\begin{eqnarray}&&{\dot{M}}_{c}={{ \mathcal S }}_{\mathrm{eff}}\,{\sigma }_{Z}(a){\rm{\Omega }}(a){F}_{g},\end{eqnarray} \tag{ 6 }$$

where ${{ \mathcal S }}_{\mathrm{eff}}$ is the effective cross section for planetesimal capture of the planet, ${\sigma }_{Z}(a)$ is the solids' surface density at the planet's orbital radius, a, ${\rm{\Omega }}(a)$ is the planet's orbital frequency, and F_g is the ratio of the gravitational to the geometric cross section (Greenzweig & Lissauer 1990, 1992), known as the gravitational enhancement factor. Further details can be found in Pollack et al. (1996 and references therein). The accretion rate in Equation (6) neglects the contribution of the dust entrained in the accreted gas ($\lesssim 1$% by mass). The planetesimal radius is assumed to be $100\,\mathrm{km}$, although smaller planetesimals were also tested.

A planetary embryo accreting planetesimals at a fixed orbital radius will deplete an annular region around its orbit of full width about equal to $8\,{R}_{{\rm{H}}}$, at which point the condensed core becomes detached from the planetesimals' disk. In these calculations, the secular evolution of planetesimals is neglected and therefore, once detached, the core reaches its final mass

$$\begin{eqnarray}&&{M}_{c}^{f}\approx \sqrt{\displaystyle \frac{{(16\pi {a}^{2}{\sigma }_{Z})}^{3}}{3{M}_{\star }}},\end{eqnarray} \tag{ 7 }$$

which neglects the contribution of the envelope mass to ${R}_{{\rm{H}}}$ and is hence appropriate when ${M}_{p}\approx {M}_{c}$. The situation is more complex for a migrating planet, since the depletion rate of solids in the disk tends to be initially slower than the migration rate through the disk. Hence, the planet cuts a swathe through the solids' disk, which deepens as $| {da}/{dt}|$ reduces, eventually detaching itself. The final mass of the core in this case is more difficult to predict, as it depends on both the accretion and migration history. During the long isolation stage, a planet may be subjected to stochastic impacts, which may alter the core and envelope mass and the planet's orbit if the impactors are sufficiently massive. This possibility is not considered here.

As mentioned above, solids sink to the top of the core. This process releases energy in the envelope, affecting the local energy budget, but not the local chemical composition of the envelope. These calculations consider accretion of hydrated, partly hydrated, and anhydrous planetesimals, depending on the local disk temperature. Rocky planetesimals may reach the core nearly intact if they are large enough (they may be held together by their own gravity) or if the ram pressure does not exceed their compressive strength (D'Angelo & Podolak 2015, and references therein). Ice-rich planetesimals are more easily disrupted or entirely ablated in the envelope because the critical temperature of H₂O is only $\approx 650\,{\rm{K}}$, a value reached in relatively shallow layers of the envelope (their mass is nevertheless assumed to sink to the core).

An important part of the calculation is represented by the capture of planetesimals, which determines self-consistently the cross section ${{ \mathcal S }}_{\mathrm{eff}}$ in Equation (6), and by their interaction with the planet's envelope, which determines depth-dependent mass and energy deposition rates. This part is based on the protocols described in Pollack et al. (1996), enhanced with an improved integration algorithm of the planetesimals' trajectories (D'Angelo et al. 2014). In brief, a number of trajectories with a varying impact parameter ($\leqslant {R}_{p}$) are integrated through the envelope. The largest impact parameter for which the body hits the core surface, breaks up, or is entirely ablated provides the radius for planetesimal capture and hence the cross section ${{ \mathcal S }}_{\mathrm{eff}}$. At this point, an additional series of trajectory integrations is performed, with an impact parameter up to the capture radius, to record the ablation history and the fate of the body as a function of the impact parameter. This collective information provides the mean energy and mass deposition rates in each envelope layer.

The methods outlined above apply to both in situ and ex situ calculations. The basic difference is that the boundary conditions at R_p change over time in ex situ models, whereas in situ calculations allow only for a linear decline with time of the gas surface density at R_p. In ex situ models, the surface density ${\sigma }_{Z}$ in Equation (6) also varies as a function of the distance from the star and of the disk temperature. Additionally, ex situ models use stellar properties from stellar structure models of solar-type stars to compute the equilibrium temperature (Equation (5)), the mass-loss rate during isolation (Equation (8)), and the irradiation temperature of the disk (Equation (16)). The calculations apply a stellar model of a solar-mass and protosolar metallicity ([Fe/H] = 0.0, Asplund et al. 2009; Lodders 2010) star from Siess et al. (2000), whose radius, effective temperature, and luminosity are plotted in Figure 2. For comparison, some calculations are repeated by adopting a Yonsei-Yale model (Spada et al. 2013) for a $0.95\,{M}_{\odot }$, [Fe/H] = 0.0 star, also represented in Figure 2. In contrast, in situ models are based on fixed, solar-type values for radius, effective temperature, and luminosity of the star.

Zoom In Zoom Out Reset image size

In situ formation calculations determine the core radius, R_c, from its current mass, M_c, using tables of results from Rogers et al. (2011). The core is composed of iron and silicates, with Earth-like mass fractions of 30% and 70%, respectively. Applied to an Earth-mass planet, the results predict a radius within 2.9% of the Earth radius, ${R}_{\oplus }$. The radius R_c is then used as an inner boundary condition for the H+He envelope.

Ex situ formation calculations allow for accretion of planetesimals whose composition varies as a function of time and distance from the star. At disk temperatures below $150\,{\rm{K}}$, the planetesimals are ice-rich, 50% by mass (45% silicates and 5% iron). They become progressively ice-poor (and rich in silicates and iron) at higher disk temperatures and are anhydrous at temperatures above $250\,{\rm{K}}$ (D'Angelo & Podolak 2015), where a terrestrial-type composition is adopted (70% silicates and 30% iron by mass). The mass fractions of iron, silicates, and ice are linearly interpolated in temperature between 150 and $250\,{\rm{K}}$.

Since the composition of the condensed core may vary during the evolution of ex situ models, as the composition of accreted solids changes, detailed calculations of the core structure are performed. The core is assumed to be spherically symmetric about its center and described by the equation of mass conservation and hydrostatic equilibrium (see Appendix A). The core is taken to be fully differentiated into an iron nucleus, a silicate mantle, and an outer shell of condensed H₂O. Two-layer (with any combination of the three materials) or one-layer structures are also possible. Each material is characterized by a "cold" equation of state (EoS) relating density and pressure. We experimented with a combination of Birch–Murnaghan, Vinet, and Generalized Rydberg EoS (for a review, see Stacey & Davis 2008), and extend them to very high pressures by means of the Thomas–Fermi–Dirac theory. Details are given in Appendix A. Temperature effects on the EoS are neglected since thermal pressure is expected to provide only a minor correction to R_c for the core masses considered here (see, e.g., Valencia et al. 2006; Seager et al. 2007; Sotin et al. 2007; Sohl et al. 2014b). However, temperature effects are included to account for phase transitions within the layers (see Appendix A for details). Applying this module to an Earth-mass planet with a composition of 70% silicates and 30% iron by mass, the calculated radius is within 0.1% of ${R}_{\oplus }$. For a cold-Earth analog (67.3% silicates and 32.7% iron by mass, Sohl & Schubert 2007, pp. 27–68), the resulting radius is within 0.8% of ${R}_{\oplus }$. More detailed structure calculations, including heat transfer and additional material phases, are also presented in Appendix A. They are compared to isothermal core structures in Appendix C.

The integration of the equations proceeds from the center outward, varying the central pressure in an iterative fashion, until the pressure at the core radius, R_c, matches (within 1%) the pressure at the bottom of the envelope, which is provided by the envelope structure module (Section 2.1). The integration of the core structure equations is performed whenever M_c increases by $\geqslant 1$% or if the pressure or temperature at R_c changes by $\geqslant 10$%.

Removal of envelope gas after the planet becomes isolated, when stellar photons can directly impinge on the planet dayside, is based on the energy-limited hydrodynamics escape driven by stellar X-ray and EUV (XUV) radiation (Watson et al. 1981; Erkaev et al. 2007; Murray-Clay et al. 2009; Lopez et al. 2012). In this limit, the mass-loss rate of the envelope can be approximated as

$$\begin{eqnarray}&&{\dot{M}}_{e}^{\mathrm{iso}}=\displaystyle \frac{\varepsilon \pi {R}_{\mathrm{XUV}}^{2}{F}_{\mathrm{XUV}}}{K(\xi )}\left(\displaystyle \frac{{R}_{p}}{{{GM}}_{p}}\right),\end{eqnarray} \tag{ 8 }$$

where ${R}_{\mathrm{XUV}}\approx 1.1\,{R}_{p}$ is the envelope radius at which the atmosphere becomes optically thick to the incoming stellar XUV radiation and most of the flux ${F}_{\mathrm{XUV}}$ is absorbed (Erkaev et al. 2007; Murray-Clay et al. 2009), $K(\xi )=1\mbox{--}3/(2\xi )+1/(2{\xi }^{3})$ is a reduction factor of the planet's potential energy caused by tidal forces of the star, and $\xi ={R}_{{\rm{H}}}/{R}_{p}\gt 1$ (Erkaev et al. 2007). Escape is assumed to take place at the equipotential surface passing through the collinear Lagrange point L₁. Note that Equation (8) diverges for $\xi \to 1$. During the evolution in isolation, it is assumed that ${\dot{M}}_{e}=-{\dot{M}}_{e}^{\mathrm{iso}}$.

Observations of the histories of EUV and X-ray fluxes of solar-type stars suggest that mass loss is most vigorous at ages of $t\lt 0.1\,\mathrm{Gyr}$ (Ribas et al. 2005). Here, we set ${F}_{\mathrm{XUV}}=3\times {10}^{-4}{L}_{\star }/(4\pi {a}^{2})$ for $t\lesssim 0.1\,\mathrm{Gyr}$ and ${F}_{\mathrm{XUV}}=3\times {10}^{-6}{(5\mathrm{Gyr}/t)}^{1.23}{L}_{\star }/(4\pi {a}^{2})$ at later times (Ribas et al. 2005), where L_⋆ is the stellar bolometric luminosity. The quantity ε is an efficiency factor intended to roughly account for radiative losses from the envelope (Erkaev et al. 2007), so that only a fraction of the incident flux ${F}_{\mathrm{XUV}}$ can effectively drive mass loss. This factor is quite uncertain. These calculations are based on an efficiency of $\varepsilon =0.1$ (e.g., Murray-Clay et al. 2009; Lopez et al. 2012). However, other values of ε are considered for a sensitivity study.

In situ and ex situ calculations handle radiation-induced mass loss of the envelope in similar ways. The only basic difference is that, in Equation (8) and in the function $K(\xi )$, ${R}_{\mathrm{XUV}}$ replaces R_p in the in situ models. Additionally, the stellar luminosity of ex situ models varies in time, according to the applied stellar evolution model (see Figure 2), whereas in situ simulations use ${L}_{\star }={L}_{\odot }$. Since the XUV stellar output is assumed to be proportional to L_⋆, the mass-loss history of the envelope during isolation may differ from that occurring at a constant value of L_⋆, as used by in situ calculations.

Other mass-loss mechanisms during the beginning of the isolation phase have been considered. Ikoma & Hori (2012) and Ginzburg et al. (2016) studied the effect of loss of pressure support at the planet's outer boundary once the disk's gas disperses; the energy for the mass loss is supplied by the planet's cooling luminosity. Owen & Wu (2016) considered a similar mechanism, initiated by the loss of pressure from the disk, in which a "Parker" wind is driven by a combination of the stellar continuum radiation and the gravitational energy released as the planet contracts. The assumptions made in these works are considerably different from those made here, and it is not clear if these processes would be significant in our models. First, the planet radius in the above papers is implicitly assumed to be close to ${R}_{{\rm{B}}}$, while in the models discussed here, just before disk dispersal, R_p is roughly a factor of $\approx 3$ to $\approx 6$ smaller than ${R}_{{\rm{B}}}$ (see also Section 3). Second, as the disk disperses, the possible up-lifting of the outer envelope layers due to loss of disk pressure during the nebular stage is implicitly included in the structure calculation through the boundary conditions at R_p. After disk dispersal, the boundary conditions transition to those of an isolated photosphere (Equations (2) and (3)). The transition occurs on the cooling timescale of the outer envelope layers, which is shorter than the disk dispersal timescale. As a result, the photospheric pressure increases and the planet further contracts inside ${R}_{{\rm{B}}}$. The Owen & Wu mechanism, for example, is not significant for ${R}_{p}\lesssim 0.1\,{R}_{{\rm{B}}}$. Clearly, whether these wind mechanisms are in fact unimportant in calculations like ours needs to be tested in more detail.

These models consider the evolution of an axisymmetric gaseous disk, of surface density Σ, driven by turbulence viscosity ν (of some nature), stellar-induced photo-evaporation, tidal torques due to gravitational interactions with an embedded planet, and accretion of gas on the planet. Indicating with ${{ \mathcal T }}_{\nu }(r)$ the viscous torque acting between two adjacent disk's annuli at a distance of r from the star and with $-{ \mathcal T }(r)$ the tidal torque exerted by the planet on the disk's gas at radius r, conservation of mass and momentum within the disk leads to the following disk's evolution equation (e.g., Lin & Papaloizou 1986)

$$\begin{eqnarray}&&\pi r\displaystyle \frac{\partial }{\partial t}({\rm{\Sigma }}+{{\rm{\Sigma }}}_{\mathrm{pe}}+{{\rm{\Sigma }}}_{\mathrm{ac}})=\displaystyle \frac{\partial }{\partial r}\left\{\displaystyle \frac{1}{r\,{\rm{\Omega }}}\displaystyle \frac{\partial }{\partial r}[{{ \mathcal T }}_{\nu }(r)-{ \mathcal T }(r)]\right\},\end{eqnarray} \tag{ 9 }$$

where Ω is the disk's rotation rate, ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$ is gas mass removed by photo-evaporation per unit disk surface and unit time (see Section 2.6), and ${\dot{{\rm{\Sigma }}}}_{\mathrm{ac}}$ is gas mass removed by accretion on the planet per unit disk surface and unit time. Recalling that ${{ \mathcal T }}_{\nu }=-2\pi {r}^{3}\nu {\rm{\Sigma }}\,\partial {\rm{\Omega }}/\partial r$ (Lynden-Bell & Pringle 1974) and approximating Ω to the Keplerian rotation rate $\sqrt{{{GM}}_{\star }/{r}^{3}}$, the viscous torque becomes ${{ \mathcal T }}_{\nu }=3\pi {r}^{2}\nu {\rm{\Sigma }}{\rm{\Omega }}$ and Equation (9) assumes the more familiar form

$$\begin{eqnarray}&&r\displaystyle \frac{\partial }{\partial t}({\rm{\Sigma }}+{{\rm{\Sigma }}}_{\mathrm{pe}}+{{\rm{\Sigma }}}_{\mathrm{ac}})\\ &&\quad =\displaystyle \frac{\partial }{\partial r}\left[3\sqrt{r}\displaystyle \frac{\partial }{\partial r}(\nu {\rm{\Sigma }}\sqrt{r})-\displaystyle \frac{2\,{\rm{\Sigma }}}{{\rm{\Omega }}}\displaystyle \frac{\partial { \mathcal T }}{\partial m}\right].\end{eqnarray} \tag{ 10 }$$

The quantity $\partial { \mathcal T }/\partial m$ is the torque density distribution, i.e., the gravitational torque per unit disk mass (${dm}=2\pi {\rm{\Sigma }}{rdr}$) arising from tidal interactions with the planet. This function is discussed in Section 2.5. If multiple planets orbit in the disk, $\partial { \mathcal T }/\partial m$ is the sum of all partial torque density distributions. The presence of the tidal torque term in Equation (10) naturally accounts for (planet-induced) gap formation in the density distribution. Notice that, if $\nu {\rm{\Sigma }}$ is constant in radius, a viscously evolving planet-less disk is in a steady state.

The mass removed from the disk, via accretion on the planet, per unit surface and unit time, is written as ${\dot{{\rm{\Sigma }}}}_{\mathrm{ac}}=\delta (r-a){\dot{M}}_{e}/(2\pi a)$ so that $\int 2\pi {\dot{{\rm{\Sigma }}}}_{\mathrm{ac}}{rdr}={\dot{M}}_{e}$, the planet's gas accretion rate, which ensures conservation of the mass transferred between the planet and the disk. Contrary to ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$ (which is zero or positive), ${\dot{{\rm{\Sigma }}}}_{\mathrm{ac}}$ can be positive, null, or negative. In the latter case, mass is transferred from the planet to the disk. The planet's envelope can lose mass (${\dot{M}}_{e}\lt 0$) if its radius exceeds the accretion radius ${R}_{{\rm{A}}}$ (see Section 2.1). Numerically, to avoid discontinuities, the mass added to (removed from) the planet is removed from (added to) a disk region around the planet's orbit of radial width a few times ${R}_{{\rm{H}}}$.

The presence of an accreting planet can change the mass transfer through the disk and thus alter the disk's surface density (Lubow & D'Angelo 2006). This effect cannot be described by Equation (10). To account for it, the approach of Lubow & D'Angelo (2006) is applied. The accretion rate through the disk is ${dm}/{dt}=2/(r{\rm{\Omega }})\partial ({{ \mathcal T }}_{\nu }-{ \mathcal T })/\partial r$ (adopting the convention that ${dm}/{dt}\gt 0$ for an inward transfer of mass). By indicating with $\langle {dm}/{dt}{\rangle }_{\mathrm{ext}}$ and $\langle {dm}/{dt}{\rangle }_{\mathrm{int}}$ accretion rates averaged over narrow rings, respectively, exterior and interior to the planet's orbit (and sufficiently apart from it), the condition is imposed that

$$\begin{eqnarray}&&{\left\langle \displaystyle \frac{{dm}}{{dt}}\right\rangle }_{\mathrm{ext}}={\left\langle \displaystyle \frac{{dm}}{{dt}}\right\rangle }_{\mathrm{int}}+{\dot{M}}_{e}.\end{eqnarray} \tag{ 11 }$$

Equation (11) is used to adjust dm/dt (i.e., Σ) inside and outside of the planet's orbit in a mass-conservative manner. This correction is applied only if $\langle {dm}/{dt}{\rangle }_{\mathrm{ext}}-{\dot{M}}_{e}\gt \langle {dm}/{dt}{\rangle }_{\mathrm{int}}\gt 0$. For stability reasons, mass adjustments are spread over several grid zones.

The thermal structure of the disk is determined by imposing a simple energy balance involving viscous heating, radiative cooling from the surface of the disk, and irradiation heating by the star

$$\begin{eqnarray}&&{Q}_{\nu }-{Q}_{\mathrm{cool}}+{Q}_{\mathrm{irr}}=0.\end{eqnarray} \tag{ 12 }$$

In the case of Keplerian rotation, the energy flux produced by viscous dissipation is (Mihalas & Weibel Mihalas 1999)

$$\begin{eqnarray}&&{Q}_{\nu }=\displaystyle \frac{9}{4}\nu {\rm{\Sigma }}{{\rm{\Omega }}}^{2}.\end{eqnarray} \tag{ 13 }$$

The energy flux escaping from the disk's surface and the heating flux generated by stellar photons are (Hubeny 1990)

$$\begin{eqnarray}&&{Q}_{\mathrm{cool}}=2{\sigma }_{\mathrm{SB}}\,{T}^{4}{\left(\displaystyle \frac{3}{16}{\kappa }_{{\rm{R}}}{\rm{\Sigma }}+\displaystyle \frac{1}{2}+\displaystyle \frac{1}{2{\kappa }_{{\rm{P}}}{\rm{\Sigma }}}\right)}^{-1}\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{Q}_{\mathrm{irr}}=2{\sigma }_{\mathrm{SB}}\,{T}_{\mathrm{irr}}^{4}{\left(\displaystyle \frac{3}{16}{\kappa }_{{\rm{R}}}{\rm{\Sigma }}+\displaystyle \frac{1}{2}+\displaystyle \frac{1}{2{\kappa }_{{\rm{P}}}{\rm{\Sigma }}}\right)}^{-1},\end{eqnarray} \tag{ 15 }$$

respectively. In the equations above, T is the mid-plane temperature of the disk, ${T}_{\mathrm{irr}}$ is the stellar irradiation temperature, and ${\kappa }_{{\rm{R}}}$ and ${\kappa }_{{\rm{P}}}$ are the Rosseland and Planck mean opacities. These opacity coefficients are calculated following the method of D'Angelo & Bodenheimer (2013), for grain size distributions of up to $1\,\mathrm{mm}$ in radius, and connected to the gas opacities of Ferguson et al. (2005). The advantage of using the fluxes in Equations (14) and (15) is that, in the approximation of vertically integrated quantities, they also describe optically thin disks. The irradiation temperature is written as (e.g., Menou & Goodman 2004)

$$\begin{eqnarray}&&{T}_{\mathrm{irr}}^{4}=\displaystyle \frac{1}{2}({T}_{\star }^{4}+{T}_{\mathrm{acc}}^{4}){\left(\displaystyle \frac{{R}_{\star }}{r}\right)}^{2}\left[\displaystyle \frac{2}{5}\left(\displaystyle \frac{{R}_{\star }}{r}\right)+\displaystyle \frac{2}{7}\displaystyle \frac{H}{r}\right].\end{eqnarray} \tag{ 16 }$$

The disk scale height, assuming vertical hydrostatic equilibrium, is given by $H=\sqrt{\gamma \,{k}_{{\rm{B}}}T/(\mu {m}_{{\rm{H}}})}/{\rm{\Omega }}$. The adiabatic index γ is set between 1.4 and 1.6 and the mean molecular weight is $\mu =2.39$ (${k}_{{\rm{B}}}$ is the Boltzmann constant and ${m}_{{\rm{H}}}$ the hydrogen mass). Equation (16) includes the contribution from luminosity released by accretion on the star (Pringle 1981)

$$\begin{eqnarray}&&{T}_{\mathrm{acc}}^{4}=\displaystyle \frac{1}{8\pi }\left(\displaystyle \frac{{{GM}}_{\star }{\dot{M}}_{\star }}{{\sigma }_{\mathrm{SB}}{R}_{\star }^{3}}\right).\end{eqnarray} \tag{ 17 }$$

The stellar accretion rate is computed as ${\dot{M}}_{\star }=2/(r{\rm{\Omega }})\partial {{ \mathcal T }}_{\nu }/\partial r$, from the solution of Equation (10) at the inner boundary of the disk (where ${ \mathcal T }=0$).

The disk also contains a solid component, assumed to be formed of $100\,\mathrm{km}$-radius planetesimals, whose surface density is ${\sigma }_{Z}$. In principle, this planetesimal disk would evolve through gravitational encounters, including collisions, and interactions with any embedded planet. Gas drag would also affect the evolution of this solid component, but over rather long timescales, given the size of the bodies considered here. However, for the sake of simplicity and tractability, it is assumed that ${\sigma }_{Z}$ only varies because of depletion by accretion of solids on the planet (${\dot{M}}_{c}$) and because of scattering by the planet's gravity (Ida & Lin 2004)

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{scat}}={\left(\displaystyle \frac{{{GM}}_{c}}{{R}_{c}}\right)}^{2}{\left(\displaystyle \frac{a}{2{{GM}}_{\star }}\right)}^{2}{\dot{M}}_{c}.\end{eqnarray} \tag{ 18 }$$

The equation above is a very simple approximation, based on energy arguments, and assumes that ${M}_{p}\approx {M}_{c}$, which is appropriate for the planets modeled here. The surface density ${\sigma }_{Z}$ also changes in response to temperature variations in the disk, allowing for the vaporization of ice (see Section 2.2).

Equation (10) is solved by means of a hybrid implicit/explicit numerical scheme, based on the fourth/fifth-order Dormand–Prince method (embedding backward differentiation) with an adaptive step-size control for the global accuracy of the solution (Hairer et al. 1993). Additional details and tests can be found in D'Angelo & Marzari (2012). Toward the end of the disk's life, when the evolution is entirely driven by photo-evaporation and the disk is quickly dispersed, the algorithm transitions from implicit to explicit, with a time step condition that constrains the maximum amount of mass removed from any disk annulus. Details on the solution method of Equation (12) for the energy balance are also given in D'Angelo & Marzari (2012). The disk extends in radius from the larger of $0.01\,{\rm{au}}$ and ${R}_{\star }={R}_{\star }(t)$ to $900\,{\rm{au}}$, and is discretized over 6000 grid points by imposing a constant ratio ${\rm{\Delta }}r/r$. Beyond $900\,{\rm{au}}$, to ensure that the grid boundary does not interfere with viscous spreading, a buffer zone of 400 additional grid points (at a degraded resolution) brings the outer disk edge to $\approx 3.5\times {10}^{5}\,{\rm{au}}$.

In order to describe tidal interactions between the disk and the planet, we apply the formalism of D'Angelo & Lubow (2008, 2010) for local isothermal disks. This is based on the torque density distribution, which is defined by the integral

$$\begin{eqnarray}&&{ \mathcal T }(a)\equiv 2\pi {\int }_{0}^{\infty }\displaystyle \frac{\partial { \mathcal T }}{\partial m}\,{\rm{\Sigma }}(r)r\,{dr},\end{eqnarray} \tag{ 19 }$$

where ${ \mathcal T }(a)$ is the total torque applied to the planet. In actuality, the integral is performed over the disk's radial extent. In a disk whose properties vary smoothly with radius, the theory of disk resonances (e.g., Meyer-Vernet & Sicardy 1987; Ward 1997) suggests that

$$\begin{eqnarray}&&\displaystyle \frac{\partial { \mathcal T }}{\partial m}={ \mathcal F }(x,\beta ,\zeta ){{\rm{\Omega }}}^{2}(a)\,{a}^{2}{\left(\displaystyle \frac{{M}_{p}}{{M}_{\star }}\right)}^{2}{\left(\displaystyle \frac{a}{{{\rm{\Delta }}}_{p}}\right)}^{4},\end{eqnarray} \tag{ 20 }$$

where ${ \mathcal F }$ is a dimensionless parametric function of $x=(r-a)/{{\rm{\Delta }}}_{p}$ with ${{\rm{\Delta }}}_{p}=\max [H(a),{R}_{{\rm{H}}}]$, whose extrema are at $x\approx \mp 1$. The parameters $\beta =-d\mathrm{ln}{\rm{\Sigma }}/d\mathrm{ln}r$ and $\zeta =-d\mathrm{ln}T/d\mathrm{ln}r$ are calculated as averages between $x=-4$ and 4 (the function ${ \mathcal F }$ is practically zero outside of these limits). D'Angelo & Lubow (2010, hereafter DL10) tested the validity of Equation (20) and provided analytic approximations of the function ${ \mathcal F }$, for wide ranges of the parameters β and ζ, based on 3D hydrodynamics calculations of disk-planet interactions.

Gravitational interactions transition from a linear to a nonlinear regime when $| { \mathcal T }(a)| \gtrsim | {{ \mathcal T }}_{\nu }(a)|$, or

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{M}_{p}}{{M}_{\star }}\right)}^{2}\gtrsim 3\pi \left(\displaystyle \frac{\nu }{{a}^{2}{\rm{\Omega }}}\right){\left(\displaystyle \frac{{{\rm{\Delta }}}_{p}}{a}\right)}^{3}.\end{eqnarray} \tag{ 21 }$$

Nonlinear interactions can cause order-of-magnitude variations in the surface density (relative to the unperturbed disk), as the planet mass grows. However, under typical disk conditions, the torque density $\partial { \mathcal T }/\partial m$ varies smoothly across the transition, and the maximum and minimum of the function ${ \mathcal F }$ change only by factors of the order of unity. The variation of $\partial { \mathcal T }/\partial m$ across the transition is implemented as explained in DL10.

The rate of change of the planet's orbital radius is found by imposing conservation of orbital angular momentum, which yields ${da}/{dt}=2{ \mathcal T }(a)/[{M}_{p}\,a\,{\rm{\Omega }}(a)]$. Since ${ \mathcal T }(a)$ is defined through Equation (19), the migration speed becomes

$$\begin{eqnarray}\displaystyle \frac{{da}}{{dt}} & = & 4\pi \,{\rm{\Omega }}(a)\left(\displaystyle \frac{a}{{M}_{\star }}\right)\left(\displaystyle \frac{{M}_{p}}{{M}_{\star }}\right){\left(\displaystyle \frac{a}{{{\rm{\Delta }}}_{p}}\right)}^{4}\\ & & \times \displaystyle \int { \mathcal F }(x,\beta ,\zeta ){\rm{\Sigma }}(r)r\,{dr},\end{eqnarray} \tag{ 22 }$$

where the integration is performed over the entire disk. In the linear regime, one can show that the integral on the right-hand side of Equation (22) is $\propto {H}^{2}{\rm{\Sigma }}(a)$, hence ${da}/{dt}\propto {({M}_{p}/{M}_{\star })(\pi {a}^{2}{\rm{\Sigma }}(a)/{M}_{\star })(a/H)}^{2}$, which is proportional to both the planet mass and the local disk mass ($\pi {a}^{2}{\rm{\Sigma }}$). A comparison between a direct 3D calculation of planet migration and Equation (22) is shown in Figure 9 of DL10. In the nonlinear regime, the integral depends on the planet mass through functions ${ \mathcal F }$ and Σ. As the density gap deepens, the integral has nonzero contributions mostly from regions near the gap edges. D'Angelo et al. (2006) showed that this formalism provides a good agreement with results from hydrodynamics calculations of planet migration also in the nonlinear regime (for non-highly eccentric orbits). It should be noted that there are regimes of fast orbital migration in which $\partial { \mathcal T }/\partial m$ can also depend on da/dt and which may not be fully captured by the formalism applied here (D'Angelo & Lubow 2008). However, the conditions required by these extreme regimes are not met in this study.

The formalism used here for disk–planet tidal interactions relies on the local isothermal approximation of disk's gas. The resulting torques agree well with analytical estimates (Tanaka et al. 2002), when the comparison is possible (DL10; Masset & Casoli 2010). Adiabatic disks can produce torques that may behave differently (see Kley & Nelson 2012; Baruteau et al. 2014, pp. 667–689, for recent reviews). However, while prescriptions are available for the total torque ${ \mathcal T }(a)$ acting on a low-mass planet in the adiabatic limit (Masset & Casoli 2010; Paardekooper et al. 2011), there is no formalism for the description of the torque density $\partial { \mathcal T }/\partial m$ in this limit. It is important to stress that the use of the distribution function $\partial { \mathcal T }/\partial m$, but not of ${ \mathcal T }(a)$, fulfills the action-reaction principle within the disk–planet system, thus accounting for disk–planet tidal interactions. Additionally, a description based on $\partial { \mathcal T }/\partial m$, but not on ${ \mathcal T }(a)$, allows for a continuous transition between different regimes of orbital migration, without the need of relying on some gap formation criterion and imposing different migration rates. In fact, as planet mass and disk thermodynamical conditions change, the tidal interactions (and hence da/dt) adapt consistently to the changing conditions. Finally, inside $\approx 5\,{\rm{au}}$, outward migration in adiabatic disks may occur for planet masses somewhat greater than $\approx 7\,{M}_{\oplus }$ (Baruteau et al. 2014, pp. 667–689), possibly affecting the largest simulated planet. However, by the time this planet attains that mass, the local disk has become radiatively efficient.

Disk–planet tidal interactions also affect orbital eccentricity. In the linear regime, orbits tend to be circularized on timescales shorter than the migration timescales (e.g., Artymowicz 1993; Tanaka & Ward 2004). In the strong nonlinear regime, the outcome of tidal interactions is more complex (e.g., Lubow & Ida 2011, pp. 347–371), though this regime is not relevant in these calculations.

Orbital migration of a planet during formation may also be driven by interactions with planetesimals (e.g., Minton & Levison 2014 and references therein). However, since the secular evolution of the planetesimals' disk is neglected, so is planetesimal-induced migration.

The disk photo-evaporation follows an approach along the lines of Alexander & Armitage (2007), in which the total amount of gas removed from the disk per unit surface and unit time is

$$\begin{eqnarray}&&{\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}={\dot{{\rm{\Sigma }}}}_{\mathrm{dif}}+{\dot{{\rm{\Sigma }}}}_{\mathrm{rim}}.\end{eqnarray} \tag{ 23 }$$

By assumption, photo-evaporation is essentially driven by stellar EUV radiation. Gas removal by FUV and X-ray radiation is not considered (but see the discussion in Gorti & Hollenbach 2009; Gorti et al. 2009). The emission rate of EUV ionizing photons by the star is ${10}^{42}\,{{\rm{s}}}^{-1}$ (Alexander & Armitage 2007). The component ${\dot{{\rm{\Sigma }}}}_{\mathrm{dif}}$ represents the removal rate due to the "diffuse" stellar radiation, whereas the additional component ${\dot{{\rm{\Sigma }}}}_{\mathrm{rim}}$ is activated after the disk becomes radially optically thin to stellar photons inside some radius ${r}_{\mathrm{rim}}$ ("rim" photo-evaporation).

Diffuse photo-evaporation depends on the gravitational radius ${r}_{g}={{GM}}_{\star }/{c}_{s}^{2}$ (Hollenbach et al. 1994), where c_s is the sound speed of an ionized hydrogen/helium mixture at $T\approx {10}^{4}\,{\rm{K}}$, the nearly constant temperature of the upper layers of a disk heated by EUV radiation (e.g., Gorti & Hollenbach 2009). It is assumed that ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}\approx 0$ inside of the critical radius ${r}_{\mathrm{crt}}={r}_{g}/10$ ($\approx 0.7\,{\rm{au}}$, Liffman 2003; Gorti et al. 2009), where gas lies too deeply in the gravitational field of the star to escape. The maximum of ${\dot{{\rm{\Sigma }}}}_{\mathrm{dif}}$ is around the radius $r\approx {r}_{\mathrm{crt}}$.

For most of the disk evolution, ${\dot{{\rm{\Sigma }}}}_{\mathrm{rim}}=0$. At later times, when the mass supply rate operated by viscous stresses cannot keep up with the removal rate caused by ${\dot{{\rm{\Sigma }}}}_{\mathrm{dif}}$, the disk's gas becomes locally depleted (typically, somewhat inward of $1\,{\rm{au}}$). Inside of this density gap induced by photo-evaporation, gas viscously drains toward the star on relatively short timescales, of the order of 10⁴ years for the kinematic viscosity adopted in this study. Once the disk develops an inner cavity, becoming optically thin interior to $r={r}_{\mathrm{rim}}$, ${\dot{{\rm{\Sigma }}}}_{\mathrm{rim}}$ provides an additional contribution to photo-evaporation at and around the rim region. As a result of the enhanced ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$, the rim radius ${r}_{\mathrm{rim}}$ increases as the disk disperses from the inside out. The presence of a sufficiently massive planet, with a semimajor axis of $a\approx {r}_{\mathrm{crt}}$, can aid in the formation of the photo-evaporation induced gap through gas depletion by tidal torques. Beyond the critical radius, a planet accreting gas at high rates can reduce the gas density interior to its orbit (see Equation (11)) and hence facilitate gap formation by photo-evaporation.

Table 2 gives a summary of the final properties of the six simulated Kepler 11 planets, along with the deduced value of ${\sigma }_{Z}$. The final core mass, the final envelope mass, the final core radius, the final planet radius, the final planet luminosity, the composition, and the orbital position are presented.

Figures 3 and 4 show, for various planets, the evolution of the core mass, envelope mass, and outer radius. In general, because of the high solid surface densities required, the core mass increases rapidly, on a timescale of $\lt {10}^{5}$ years. In fact, by using Equation (6), the accretion timescale of the core is

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{c}}{{\dot{M}}_{c}}=\displaystyle \frac{4}{3}\displaystyle \frac{{\rho }_{c}{R}_{c}}{{\sigma }_{Z}{F}_{g}}\sqrt{\displaystyle \frac{{a}^{3}}{{{GM}}_{\star }}},\end{eqnarray} \tag{ 24 }$$

where the cross section for planetesimals' capture is ${{ \mathcal S }}_{\mathrm{eff}}\approx \pi {R}_{c}^{2}$ when the gas bound to the core is very tenuous. If ${F}_{g}\approx 1$, Equation (24) gives initial accretion timescales (see a and ${\sigma }_{Z}$ in Table 2) $\lesssim {10}^{5}\,\mathrm{years}$. Nonetheless, at an initial envelope mass of $\approx {10}^{-3}\,{M}_{\oplus }$, even $100\,\mathrm{km}$-size planetesimals are affected by gas drag in the envelope and ${{ \mathcal S }}_{\mathrm{eff}}$ becomes $\gg \pi {R}_{c}^{2}$ (D'Angelo et al. 2014). In the calculations, the actual timescales are all of the order of ${10}^{4}\,\mathrm{years}$ (see Figures 3 and 4). Thus, any uncertainties in Equation (6) or in the choice of the initial core mass have practically no effect on the final result. The core mass levels off at the isolation mass, as given by Equation (7). The envelope mass increases more slowly, on a timescale of ${10}^{6}\,\mathrm{years}$. The outer radius during the formation phases shows an initial rapid rise corresponding to the rapid core growth, and then a nearly flat section since the outer boundary condition is essentially determined by the nearly constant core mass. Once the transition at ${t}_{\mathrm{iso}}=3.5\,{\rm{Myr}}$ is reached, the radius decreases rapidly as a result of the transition to isolated boundary conditions. Beyond that time, the radius decreases slowly as a result of contraction and cooling. As a further effect, the envelope mass and outer radius decline as a result of mass loss induced by stellar XUV radiation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some properties of the Kepler 11 planets at the time of disk dispersal ($t={t}_{\mathrm{iso}}$), according to our in situ models, are reported in Table 3. Planet Kepler 11b would lose its entire H+He envelope in $4\times {10}^{7}$ years. For planets Kepler 11f through c, proceeding inwards, mass-loss rates are 1 to $5\times {10}^{-9}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$ at $10\,{\rm{Myr}}$, decreasing to 1 to $5\times {10}^{-10}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$ at $100\,{\rm{Myr}}$. At the final age of $8\,{\rm{Gyr}}$, these rates are down to $2.3\times {10}^{-13}$ to $1.3\times {10}^{-12}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$.

Table 3.  Properties of In Situ Formation Models of Kepler 11 Planets at Isolation^a

Planet	${t}_{\mathrm{iso}}$[Myr]	${M}_{c}/{M}_{\oplus }$	${M}_{e}/{M}_{\oplus }$	${R}_{c}/{R}_{\oplus }$	${R}_{p}/{R}_{\oplus }$	(Fe, Si, H+He)%b	${T}_{\mathrm{eq}}$ c [K]	$\langle {\dot{M}}_{e}\rangle$ d $[{M}_{\oplus }\,{{\rm{yr}}}^{-1}]$
b	3.5	1.96	0.05	1.19	5.3	$(29.2,68.3,2.5)$	927	$-4.0\times {10}^{-9}$
c	3.5	5.76	0.53	1.60	9.5	$(27.5,64.1,8.4)$	880	$-1.3\times {10}^{-9}$
d	3.5	5.01	0.74	1.53	11.7	$(26.1,61.0,12.9)$	731	$-1.2\times {10}^{-9}$
e	3.5	6.66	1.65	1.67	14.3	$(24.0,56.1,19.9)$	623	$-1.0\times {10}^{-9}$
f	3.5	2.84	0.16	1.33	8.4	$(28.4,66.3,5.3)$	550	$-2.6\times {10}^{-10}$
g	3.5	5.01	0.75	1.53	9.4	$(25.7,59.9,14.4)$	409	$-1.0\times {10}^{-10}$

Notes.

^aThe isolation time, ${t}_{\mathrm{iso}}$, is the time at which the disk's gas is assumed to disperse. ^bPercentage of the planet mass at time $t={t}_{\mathrm{iso}}$. ^cConstant equilibrium temperature, ${T}_{\mathrm{eq}}$, during the isolation phase. ^dRate of change of the envelope mass averaged over the first $100\,{\rm{Myr}}$ of evolution in isolation. For Kepler 11b, $\langle {\dot{M}}_{e}\rangle$ is an average over $10\,{\rm{Myr}}$.

Download table as:  ASCIITypeset image

Specifically, the main results from our in situ models can be summarized as follows.

Kepler 11b: a low-mass H+He envelope forms around the core mass of $1.96\,{M}_{\oplus }$, but during the isolation phase this envelope is entirely lost. The final mass is consistent with the measured mass ${1.9}_{-1.0}^{+1.4}\,{M}_{\oplus }$ (L13), but the final radius, the core radius of $1.19\,{R}_{\oplus }$, is far below the measured value of ${1.80}_{-0.05}^{+0.03}\,{R}_{\oplus }$, as shown in the top-right panel of Figure 3. Even by taking the upper limit of the measured mass, a 100% silicate core would still have too small a radius, only $1.5\,{R}_{\oplus }$. A steam envelope (not modeled here) is probably required to achieve consistency (Lopez et al. 2012). This possibility, however, is inconsistent with in situ formation inside $0.1\,{\rm{au}}$ of a solar-type star because of the lack of ice in the core. Another possibility is the release of gas sequestered by the core during formation. The result that the entire envelope mass is lost remains valid even if the assumed core mass is increased to $3\,{M}_{\oplus }$.

Kepler 11c: about half of the accreted envelope mass is lost during the isolation phase, but the final radius agrees with the measured radius of ${2.87}_{-0.06}^{+0.05}\,{R}_{\oplus }$ (see Figure 3, center-right). The final total mass of $6.02\,{M}_{\oplus }$ falls just above the one-standard-deviation upper limit for the measured mass of ${2.9}_{-1.6}^{+2.9}\,{M}_{\oplus }$, but is certainly within the uncertainties in the theoretical models.

Kepler 11d: about one-third of the accreted envelope mass is lost during the isolation phase. The final computed radius of $3.24\,{R}_{\oplus }$ is within 4% of the measured value of ${3.12}_{-0.07}^{+0.06}\,{R}_{\oplus }$ (see Figure 3, bottom-right). The final computed total mass of $5.5\,{M}_{\oplus }$ is just below the one-standard-deviation lower limit for the measured mass of ${7.3}_{-1.5}^{+0.8}\,{M}_{\oplus }$. These small discrepancies are well within the uncertainties of the models. To reduce the radius to agree with the measured value would require reducing the mass, increasing the discrepancy with the measured value.

Kepler 11e: at a separation from the star of about $0.2\,{\rm{au}}$, the accreted H+He envelope of this object loses only about 12% of its mass during the isolation phase. The final computed radius of $4.24\,{R}_{\oplus }$ agrees to within 1.5% with the measured value of ${4.19}_{-0.09}^{+0.07}\,{R}_{\oplus }$, as indicated in the top-right panel of Figure 4. The final computed total mass of $8.11\,{M}_{\oplus }$ agrees with the measured mass of ${8.0}_{-2.1}^{+1.5}\,{M}_{\oplus }$. A slight reduction in the assumed mass ($\approx 5$%) would bring the radius into agreement with the observed value and the planet mass would still agree with the observed mass, well within one-standard-deviation uncertainties.

Kepler 11f: even though the planet lies farther from the star ($0.25\,{\rm{au}}$) than Kepler 11e, its lower core mass (2.8 versus $6.7\,{M}_{\oplus }$) results in 31% of the H+He envelope being lost during the isolation phase. The final computed radius of $2.42\,{R}_{\oplus }$ agrees with the measured value of ${2.49}_{-0.07}^{+0.04}\,{R}_{\oplus }$ to within 2.5% (see Figure 4, center-right). The final total mass of $2.95\,{M}_{\oplus }$ is slightly above the one-standard-deviation upper limit for the measured value of ${2.0}_{-0.9}^{+0.8}\,{M}_{\oplus }$. To improve the agreement with the observed radius, the assumed mass would have to increase by about $0.1\,{M}_{\oplus }$, increasing the (small) discrepancy with the observed value. In any case, the agreement is within the uncertainties of the theoretical model.

Kepler 11 g: at a distance of $0.466\,{\rm{au}}$, the planet loses only about 2% of its H+He envelope mass during the isolation phase. The computed final radius of $3.34\,{R}_{\oplus }$ agrees with the measured value of ${3.33}_{-0.08}^{+0.06}\,{R}_{\oplus }$. The corresponding computed total mass is $5.75\,{M}_{\oplus }$. The observed mass in this case is not well constrained; it is less than $25\,{M}_{\oplus }$.

The two phases of planet evolution identified in Section 3 can also be defined in ex situ models. The time t = 0 coincides with the time at which the disk evolution starts from the imposed initial conditions. The gaseous disk evolution depends on several quantities (see Sections 2.4–2.6). The gas surface density at t = 0 is ${\rm{\Sigma }}=1110\sqrt{1\,{\rm{au}}/r}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ with an exponential cutoff beyond some radius, so that a total disk mass of $0.03\,{M}_{\odot }$ is initially confined within $\approx 60\,{\rm{au}}$ of the star (see, e.g., Williams & Cieza 2011). This surface density distribution, ${\rm{\Sigma }}(t=0)$, was determined after a number of attempts aimed at reproducing the observed physical properties of planet Kepler 11b and at tentatively matching those of Kepler 11f (the planet farthest from the star for which both R_p and ${M}_{p}$ are constrained by transit observations). Only later it was realized that, fortuitously, ${\rm{\Sigma }}(t=0)$ inside $\approx 10\,{\rm{au}}$ matches quite closely the minimum-mass solar nebula density of Davis (2005). Although no other slope $d\mathrm{ln}{\rm{\Sigma }}/d\mathrm{ln}r$ was tested for the initial Σ, it is unlikely that the adopted initial surface density provides a unique solution to the problem. In any case, the choice of the initial density values is likely to influence the outcomes of the models more than the choice of the initial slope $d\mathrm{ln}{\rm{\Sigma }}/d\mathrm{ln}r$.

These calculations apply a time-constant kinematic viscosity $\nu ={\nu }_{1}\sqrt{r/1\,{\rm{au}}}$, where ${\nu }_{1}\approx 4\times {10}^{14}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$. In terms of the viscosity prescription of Shakura & Sunyaev (1973), the α-parameter quantifying turbulence varies with time and distance from the star. The value around the starting orbital radius of the simulated planets is between $\approx {10}^{-3}$ and $\approx {10}^{-2}$. The disk provides an initial accretion rate toward the star of the order of ${10}^{-7}\,{M}_{\odot }\,{{\rm{yr}}}^{-1}$.

The lifetime of the gaseous disk is determined by ${\rm{\Sigma }}(t=0)$, ν, the initial disk mass, and the photo-evaporation rate. All of these quantities are the same for all planet models. As mentioned in Section 2.6, the presence of a planet may affect disk dispersal, to a smaller or larger extent depending on the planet mass, its gas accretion rate, and orbital radius. In fact, although the planets end up inside the critical radius ${r}_{\mathrm{crt}}$, they do spend most of their disk-embedded evolution at larger radii and, therefore, they can potentially influence ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$. However, in the models presented here, this effect appears to be marginal and the gas inside $1\,{\rm{au}}$ is dispersed in $\approx 4\,{\rm{Myr}}$, with a time-spread among models of about 2%.

The evolution of Σ and T is illustrated in Figure 5 (see the figure's caption for details). Dust opacity transitions are visible in the temperature profiles, the most prominent of which are represented by the evaporation of the silicate species above $1000\,{\rm{K}}$. The fainter opacity transitions associated with the evaporation of icy grains are also visible (located around $3\,{\rm{au}}$ at $t\approx {10}^{5}\,\mathrm{years}$ and $1\,{\rm{au}}$ at $t\approx {10}^{6}\,\mathrm{years}$). In the forming region of Kepler 11 planets ($r\gt 0.09\,{\rm{au}}$), the gas temperature is $\lesssim 1600\,{\rm{K}}$ at $t\approx {10}^{5}\,\mathrm{years}$ and becomes $\lesssim 1000\,{\rm{K}}$ at times $t\gtrsim 5\times {10}^{5}\,\mathrm{years}$. By the time the simulated planets have settled on their final orbits, the local gas temperature varies between $\approx 160$ and $\approx 500\,{\rm{K}}$.

Zoom In Zoom Out Reset image size

Gas photo-evaporation by stellar irradiation produces a density gap somewhat inward of $1\,{\rm{au}}$, around the radius ${r}_{\mathrm{crt}}\approx 0.7\,{\rm{au}}$, where the ratio between the photo-evaporation timescale and the accretion timescale through the disk is smallest. Viscous diffusion quickly removes gas inward of ${r}_{\mathrm{crt}}$ on a timescale of $\sim {r}_{\mathrm{crt}}^{2}/\nu$ (see Section 2.6), generating a cavity at $t\approx 4.1\,{\rm{Myr}}$. Afterwards, rim photo-evaporation dissipates gas inside-out, pushing the cavity edge outward to $r\approx 40\,{\rm{au}}$ by $t\approx 4.8\,{\rm{Myr}}$ (see Figure 5). Inside the disk cavity, which is virtually devoid of gas, the temperature T is set equal to the irradiation temperature, so that ${T}^{4}={L}_{\star }/(4\pi {\sigma }_{\mathrm{SB}}{r}^{2})$. Beyond the cavity edge, the temperature is set by the gas thermal balance, Equation (12) (hence the large temperature transitions in the bottom panel of Figure 5 for $t\gt 4\,{\rm{Myr}}$). The evolution of Σ and T in Figure 5 is the same for all models, except for variations induced by disk–planet tidal interactions and gas accretion on the planet.

The energy output of the central star can impact both the evolution of the disk and that of the planet. The two stellar models considered here (see Section 2.1) show similar luminosities for $t\gtrsim 10\,{\rm{Myr}}$ (see Figure 2), implying similar evolution of the planets in isolation. However, there are differences at earlier times, specifically in the effective temperature T_⋆ and radius R_⋆, which enter the irradiation temperature in Equation (16). Therefore, the disk thermal budget may be affected and hence the planet migration history may differ somewhat (see Section 2.5). The temperature at the planet surface may also change. These differences are assessed for the cases of Kepler 11b and c (see Section 4.2).

Among the main simplifications of this study are the neglect of planet–planet interactions and the fact that only a single planet evolves in the disk, though the models account for the depletion of the planetesimal disk generated by planets that have already initiated the formation process. It is assumed that by the time $t={t}_{i}$, a solid core of mass ${M}_{c}=0.1\,{M}_{\oplus }$ (corresponding to ${M}_{e}\approx {10}^{-7}\,{M}_{\oplus }$, calculated consistently with M_c and the disk boundary conditions) has formed at a current orbital radius of $r={a}_{i}$. There is no speculation about its previous accretion/migration history and a_i is considered to be the initial orbital radius of the simulated planet. (However, the small initial planet mass implies that orbital migration via disk–planet tidal interactions may be negligible at times $t\lt {t}_{i}$). Both ${t}_{i}$ and a_i are free parameters, constrained by the requirement, among others, that no two orbital paths intersect each other. Strictly speaking, this is not a physical requirement (e.g., Hands et al. 2014) but rather a necessity dictated by the absence of gravitational interactions among planets. The time ${t}_{i}$ ranges from $\approx 6\times {10}^{4}\,\mathrm{years}$ (Kepler 11b) to $\approx {10}^{6}\,\mathrm{years}$ (Kepler 11f). Clearly, the time ${t}_{i}$ is also determined by the choice of using the same initial core mass for all planets. The orbital radius a_i ranges from $\approx 2.1\,{\rm{au}}$ for Kepler 11f to $5.35\,{\rm{au}}$ for Kepler 11e.

Table 4 summarizes the final properties (and a_i) of the six simulated Kepler 11 planets, assumed to have formed ex situ. These are referred to as reference models. The initial orbital radius of each planet, a_i, and the epoch ${t}_{i}$ are found by trial and error so that the final model provides reasonable matches to a, R_p, and ${M}_{p}$ at an age of $8\,{\rm{Gyr}}$, as reported by L13. Additionally, as mentioned above, orbital paths must not intersect. Plots of various quantities versus time from the resulting models are illustrated in Figures 6 and 7. Since all planets start well beyond $1\,{\rm{au}}$, contrary to in situ models, the local density of solids at a_i and ${t}_{i}$ is moderate to low: ${\sigma }_{Z}\approx 9\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ for Kepler 11b, $\approx 6\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ for Kepler 11f, and between $\approx 3$ and $\approx 4\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ for the other planets (in ascending order of a_i).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Summary of Results for Ex Situ Formation of Kepler 11 Planets^a

Planet	${M}_{c}/{M}_{\oplus }$	${M}_{e}/{M}_{\oplus }$	${R}_{c}/{R}_{\oplus }$	${R}_{p}/{R}_{\oplus }$	${L}_{p}/{L}_{\odot }$	(Fe, Si, H2O)%b	(Fe, Si, H2O, H+He)%c	ai [au]	af [au]
b	2.10	0.00	1.47	1.47		$(10.6,50.6,38.8)$	$(10.6,50.6,38.8,0.0)$	2.14	0.091
c	4.56	0.10	1.84	2.84	$2.1\times {10}^{-7}$	$(6.6,46.5,46.9)$	$(6.5,45.5,45.9,2.1)$	3.94	0.109
d	5.58	0.31	1.93	3.14	$1.4\times {10}^{-7}$	$(6.0,46.0,48.0)$	$(5.7,43.6,45.5,5.2)$	4.68	0.156
e	6.90	1.29	2.00	4.14	$1.6\times {10}^{-7}$	$(5.7,45.7,48.6)$	$(4.8,38.5,41.0,15.7)$	5.35	0.194
f	2.74	0.07	1.62	2.49	$3.4\times {10}^{-8}$	$(6.1,46.1,47.8)$	$(6.0,44.9,46.6,2.5)$	2.10	0.248
g	5.57	0.69	1.92	3.40	$1.9\times {10}^{-8}$	$(5.0,45.0,50.0)$	$(4.5,40.1,44.4,11.0)$	4.41	0.469

Notes.

^aValues at time $t=8\,\mathrm{Gyr}$. ^bPercentage of the core mass. "Fe," "Si," and "H₂O" indicate, respectively, the iron nucleus, the silicate mantle, and the H₂O outer shell of the core (see Appendix A for details). ^cPercentage of the planet mass.

Download table as:  ASCIITypeset image

Models for each planet were constructed in ascending order of final orbital radius, a_f. At time t = 0, the gas-to-solid mass ratio is set to about 70 beyond the ice condensation line at around $3\,{\rm{au}}$ ($T\lt 150\,{\rm{K}}$). Planetesimals are anhydrous interior to $\approx 1\,{\rm{au}}$ ($T\gt 250\,{\rm{K}}$), where the gas-to-solid mass ratio becomes approximately 140 (see Section 2.2 for details). As the disk evolves, the ice sublimation line moves inward to $r\approx 1\,{\rm{au}}$ by $t\approx 1\,{\rm{Myr}}$ (see Figure 5). The model for Kepler 11b is constructed from this initial surface density of solids, ${\sigma }_{Z}(t=0)$. The distribution ${\sigma }_{Z}$ is depleted by the passage of Kepler 11b. Indicating with a^c_i the initial and a^c_f the final (i.e., observed) orbital radii of Kepler 11c, the distribution ${\sigma }_{Z}$ for modeling this planet is determined by taking the depleted mass in solids between a^c_f and $\approx 1.25\,{a}_{i}^{c}$, and redistributing the mass over the region according to a $1/\sqrt{r}$ power law. The same procedure is used to determine ${\sigma }_{Z}$ for the construction of models for Kepler 11d and e (each based on the depleted reservoir of solids left by the preceding planet). For the models of Kepler 11f and g, which start inside the orbits of preceding planets at significantly later times, the depleted mass in solids is redistributed interior to $\approx 1.3\,{a}_{i}^{e}$ (down to their observed orbital radii).

Typically, models make a transition from disk to photospheric boundary conditions (see Section 2.1) at the isolation time, ${t}_{\mathrm{iso}}$. Although in some models gas accretion can be disk-limited during late stages of formation, as for in situ models a proper phase of rapid gas accretion (i.e., R_p significantly smaller than ${R}_{{\rm{A}}}$) is never reached during the formation phase. Accretion of solids could in principle (and does on occasion) continue beyond the isolation time, until the feeding zone is emptied (which requires ${da}/{dt}\approx 0$). However, since orbital migration becomes very slow much earlier than ${t}_{\mathrm{iso}}$ (see Figures 6 and 7), M_c plateaus well before isolation is achieved, as can be seen in the left panels of Figures 6 and 7. For all practical purposes, a planet is isolated from both the disk's gas and solids at $t\gt {t}_{\mathrm{iso}}$. Some properties of the reference models of Kepler 11 planets at $t={t}_{\mathrm{iso}}$ are listed in Table 5. The small scatter in isolation times is likely caused by the removal of gas via accretion on the planet (when $a\gt {r}_{\mathrm{crt}}$), which tends to lower the accretion rate through the disk for $r\lesssim a$ (Lubow & D'Angelo 2006, see also Equation (11)) and thus operates in concert with disk photo-evaporation to augment gas depletion inside $r\approx a$. In fact, the time ${t}_{\mathrm{iso}}$ is shorter for planets with larger envelope masses, i.e., with larger $\langle {\dot{M}}_{e}\rangle$. However, the contribution of ${\dot{M}}_{e}$ to ${\dot{{\rm{\Sigma }}}}_{\mathrm{pe}}$ is quite marginal in these calculations. In the more realistic situation in which all planets migrated in the disk, the time ${t}_{\mathrm{iso}}$ would be set by the largest planet, Kepler 11e. However, since the formation phases of all planets are basically complete by that time, no significant consequences would be anticipated.

Table 5.  Properties of Ex Situ Formation Models of Kepler 11 Planets at Isolation^a

Planet	${t}_{\mathrm{iso}}$(Myr)	${M}_{c}/{M}_{\oplus }$	${M}_{e}/{M}_{\oplus }$	${R}_{c}/{R}_{\oplus }$	${R}_{p}/{R}_{\oplus }$	(Fe,Si,H2O,H+He)%b	${T}_{\mathrm{eq}}$ c [K]	$\langle {\dot{M}}_{e}\rangle$ d $[{M}_{\oplus }\,{{\rm{yr}}}^{-1}$]
b	4.11	2.10	0.007	1.47	3.2	$(10.6,50.4,38.6,0.4)$	819	$-3.8\times {10}^{-10}$
c	4.08	4.56	0.235	1.84	8.0	$(6.2,44.3,44.6,4.9)$	749	$-5.8\times {10}^{-10}$
d	4.06	5.58	0.395	1.94	12.7	$(5.6,42.9,44.9,6.6)$	627	$-4.0\times {10}^{-10}$
e	4.04	6.90	1.399	2.02	14.9	$(4.8,38.0,40.4,16.8)$	563	$-4.7\times {10}^{-10}$
f	4.09	2.74	0.098	1.62	6.8	$(5.9,44.5,46.1,3.5)$	492	$-1.1\times {10}^{-10}$
g	4.04	5.57	0.701	1.94	9.8	$(4.5,40.0,44.3,11.2)$	361	$-1.0\times {10}^{-10}$

Notes.

^aThe isolation time, ${t}_{\mathrm{iso}}$, is the time at which the disk's gas at $r\lesssim a$ disperses. ^bPercentage of the planet mass at time $t={t}_{\mathrm{iso}}$. ^cEquilibrium temperature of the planet, Equation (5), at $t={t}_{\mathrm{iso}}$. ^dRate of change of the planet's envelope mass averaged over the first $100\,{\rm{Myr}}$ of evolution in isolation. For Kepler 11b, $\langle {\dot{M}}_{e}\rangle$ is an average over $10\,{\rm{Myr}}$.

Download table as:  ASCIITypeset image

The relative gas content is largest for Kepler 11e, accounting for $\approx 17$% of the total mass at ${t}_{\mathrm{iso}}\approx 4\,{\rm{Myr}}$, and somewhat less at $8\,\mathrm{Gyr}$ (see Table 4). The light elements (H+He) in the other planets make $\lesssim 10$% of the total mass at the isolation time and, in most cases, only a few to several percent at $8\,\mathrm{Gyr}$. At this age, however, the gaseous envelope always accounts for $\approx 35$% to $\approx 50$% of the planet radius. For all planets, the condensible mass fraction of H₂O is $\gtrsim 39$%, indicating that their cores mostly form behind the ice condensation line. Although the composition of the initial core is dictated by the local disk composition of the solids at $t={t}_{i}$, its mass ($0.1\,{M}_{\oplus }$) is small enough to not affect the final core composition much. Kepler 11b contains the smallest mass fraction of H₂O and the largest mass fractions of silicates and iron, due to its small initial orbital radius ($\approx 2\,{\rm{au}}$) and early start time, ${t}_{i}\approx 6\times {10}^{4}\,\mathrm{years}$. Nonetheless, also in this case the substantial fraction of H₂O ($\approx 39$% by mass) implies that the planet accumulates its condensible inventory mostly behind the ice condensation front. Despite an equally small starting orbit, the core composition of Kepler 11f is instead more similar to that of neighboring planets because of its late start time (${t}_{i}\approx 1.1\,{\rm{Myr}}$) and growth in a colder disk environment.

Except for the varying stellar properties (see Figure 2), the isolation evolution of ex situ models behaves as that of in situ models. The surface temperature of the planet closely follows the equilibrium temperature, ${T}_{\mathrm{eq}}$, which is shown in Figure 8 for the isolation evolution of Kepler 11c, e, and g.

Zoom In Zoom Out Reset image size

As for in situ models, the planet radius R_p decreases on a short timescale once the planet becomes isolated. Afterwards, R_p steadily declines as the planet cools. The radius of planets Kepler 11c and d reaches a minimum at an age between 7 and $7.5\,\mathrm{Gyr}$, after which the envelope begins to slowly expand, following the rise of ${T}_{\mathrm{eq}}$ (see Figure 8). However, the expansion is very modest and by an age of $12\,\mathrm{Gyr}$, R_p increases over its minimum value by $\lt 1$%. Kepler 11f follows a similar trend, achieving a minimum radius around the age of $7.5\,\mathrm{Gyr}$ and then inflating slightly (R_p changing by $\lt 1$%). Possibly due to their more massive envelopes and hence larger internal energy, the simulated planets Kepler 11e and g still contract at an age of $\approx 11\,\mathrm{Gyr}$.

The evaporative loss of envelope gas, caused by the absorption of stellar X-ray and EUV photons and given by Equation (8), starts at $t={t}_{\mathrm{iso}}\approx 4\,{\rm{Myr}}$, i.e., as soon as the gaseous disk interior to the planets' orbit is cleared and the isolation phase begins. As for in situ models, a standard value for the efficiency parameter, $\varepsilon =0.1$, is applied (but see also Section 4.2). The resulting envelope masses versus time are illustrated in Figure 9. The rate $| {\dot{M}}_{e}|$ is especially large at early ages, $t\lesssim 100\,{\rm{Myr}}$, because both R_p and the flux ${F}_{\mathrm{XUV}}$ (see Section 2.3) are large. Values of the gas loss rates, averaged over the first $100\,{\rm{Myr}}$ ($10\,{\rm{Myr}}$ for Kepler 11b) of evolution, are listed in Table 5 and range from 10⁻¹⁰ to around $6\times {10}^{-10}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$, i.e., between $\approx {10}^{10}$ and $\approx {10}^{11}\,{\rm{g}}\,{{\rm{s}}}^{-1}$. Mass loss is significant for the inner planets Kepler 11b, c, and d. By $8\,{\rm{Gyr}}$ these planets lose, respectively, 100%, $\approx 57$%, and $\approx 22$% of their M_e at $t={t}_{\mathrm{iso}}$. Envelope loss is also substantial for Kepler 11f, which loses nearly 29% of its gaseous mass during the isolation phase. In relative terms, Kepler 11e and g are more immune to evaporative mass loss, as M_e reduces by about 8% and 2%, respectively, from $\approx 4\,{\rm{Myr}}$ to $8\,{\rm{Gyr}}$. Of the H+He mass removed during isolation, at least 38%–40% evaporates within $\approx 100\,{\rm{Myr}}$. As planets contract, ${R}_{p}/{R}_{{\rm{H}}}\ll 1$ hence $K(\xi )\approx 1$ in Equation (8), and H+He gas is stripped from the envelope at a rate ${\dot{M}}_{e}^{\mathrm{iso}}\propto {R}_{p}^{3}/({M}_{p}{a}^{2})$. After $8\,{\rm{Gyr}}$, mass-loss rates become quite small for all planets, ranging from around 10⁻¹³ (for Kepler 11f and g) to $\approx {10}^{-12}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$ (for Kepler 11c and d).

Zoom In Zoom Out Reset image size

Comparing the evolution of R_p in Figures 3 and 6 after isolation, one can see that planets formed in situ remain somewhat more inflated than planets formed ex situ. The difference is especially large for the case of Kepler 11b, probably due to the fact that the planet formed in situ acquires a more massive H+He envelope (because of the rapid core growth). Differences tend to vanish at later times. This behavior accounts for the differences in the average rates $\langle {\dot{M}}_{e}\rangle$ in Tables 3 and 5. Although there is a factor of 10 difference in the gas loss rates between in situ and ex situ models of Kepler 11b, in this case, the difference is immaterial since both simulated planets lose their entire envelopes within 30–40 Myr.

The main results from our ex situ reference models can be summarized as follows.

Kepler 11b: at an initial orbital distance of $2.14\,{\rm{au}}$ (${t}_{i}\approx 6.3\times {10}^{4}\,\mathrm{years}$) and a surface density of solids ${\sigma }_{Z}\approx 9\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, according to Equation (7) a non-migrating core would achieve a final mass of $\approx 0.7\,{M}_{\oplus }$, and smaller if depletion via scattering (Equation (18)) was taken into account. About 90% of the core mass is accreted from solids orbiting beyond $0.5\,{\rm{au}}$ from the star, hence the presence of large amounts of H₂O in the planet. The envelope mass is maximum around $2.8\,{\rm{Myr}}$, before part of it becomes unbound and is released back to the disk. After the planet becomes isolated, the remaining H+He gas is removed by stellar X-ray and EUV radiation within $\approx 30\,{\rm{Myr}}$. Reducing the efficiency of the evaporative mass-loss rate, ε in Equation (8), allows the planet to retain an atmosphere for somewhat longer. However, even a value as low as $\varepsilon =0.01$ predicts a complete removal of the primordial H+He gas within a few times $100\,{\rm{Myr}}$ (see Figure 6). Despite the presence of abundant H₂O in the core ($\approx 40$% by mass), R_c is still significantly smaller than the observed radius. Although not included, this model naturally accounts for the formation of a steam atmosphere that was proposed to reconcile simulated and observed radii (Lopez et al. 2012). Both final orbital radius and planet mass agree with measured values (see Figure 10).

Zoom In Zoom Out Reset image size

Kepler 11c: the planet starts at roughly $4\,{\rm{au}}$ (${t}_{i}\approx 1.3\times {10}^{5}\,\mathrm{years}$) and grows about 90% of its condensible final mass at $r\gtrsim 0.8\,{\rm{au}}$. Nearly 46% of the planet's final mass is in H₂O, the second largest (after Kepler 11f) relative fraction of all simulated planets (although Kepler 11d, e, and g contain more H₂O in absolute measure). When the local disk disperses and the planet becomes isolated, the envelope includes about 5% of ${M}_{p}$. Roughly 90% of this H+He mass is accreted from the disk's gas at $170\lesssim T\lesssim 460\,{\rm{K}}$. Over the course of the isolation phase, stellar radiation removes more than half of the envelope mass, leaving only $\approx 2$% of ${M}_{p}$ in primordial H+He gas at $8\,{\rm{Gyr}}$, the smallest relative fraction in all simulated planets that retain an envelope. The evolution of R_p in Figure 6 shows a period, between $t\approx 20$ and $\approx 40\,{\rm{Myr}}$, in which the planet contraction slows down. This feature is likely associated with the rise in ${T}_{\mathrm{eq}}$ at that age (see Figure 8), following the brightening of the star (Figure 2). A similar feature appears in the radius evolution of Kepler 11b and f planets, which can more promptly respond to changes in the external incident flux having (with Kepler 11c) the least massive envelopes. The final planet mass and radius agree with measurements, whereas the final orbital distance is just above the one-standard-deviation upper limit of the measured value ($0.108\,{\rm{au}}$).

Kepler 11d: the starting orbital radius is less than $1\,{\rm{au}}$ larger than that of Kepler 11c and the start time is comparable (${t}_{i}\approx 1.4\times {10}^{5}\,\mathrm{years}$). Therefore, similarly to its inner neighbor, the planet accumulates $\approx 90$% of M_c outside of $\approx 0.9\,{\rm{au}}$. The maximum value of M_e is attained shortly prior to $2.5\,{\rm{Myr}}$, but afterwards some envelope gas becomes unbound and returns to the disk. For $\approx 1\,{\rm{Myr}}$ prior to isolation, R_p remains very close to the accretion radius ${R}_{{\rm{A}}}$, preventing further accretion of gas. During its evolution in isolation, the envelope loses somewhat less than $0.1\,{M}_{\oplus }$, about a fifth of M_e at $t={t}_{\mathrm{iso}}$. Comparing the core radii in Tables 4 and 5, a small difference can be noticed. The core masses at $t={t}_{\mathrm{iso}}$ and at $8\,{\rm{Gyr}}$ are virtually identical, yet the pressure applied at the top of the core at the later epoch is nearly twice as large ($\approx 9\,\mathrm{GPa}$ versus $\approx 4.8\,\mathrm{GPa}$) because of the cooler temperature, which in this case accounts for the $85\,\mathrm{km}$ reduction in R_c. The values of the planet's mass and radius and the orbital distance at $8\,{\rm{Gyr}}$ are all within measurement errors.

Kepler 11e: the most massive of the six, in terms of both core and envelope mass, this planet also has the farthest initial orbit at $5.35\,{\rm{au}}$ (${t}_{i}\approx 1.5\times {10}^{5}\,\mathrm{years}$) from the star. The initial local density of solids is $\approx 3\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. Neglecting scattering, a non-migrating planet at that distance would achieve a final core mass of around $2\,{M}_{\oplus }$ before emptying its feeding zone. The planet grows 90% of its core mass beyond $\approx 1\,{\rm{au}}$, but most of its H+He inventory is accreted inside this radius. Essentially, the planet's core fully forms behind the ice condensation line ($\approx 49$% of M_c is H₂O). The average mass-loss rate at the beginning of the isolation phase ($-\langle {\dot{M}}_{e}\rangle$, see Table 5) is higher than that of Kepler 11d despite the larger mass and orbital radius of Kepler 11e. The reason is the strong dependence of ${\dot{M}}_{e}^{\mathrm{iso}}$ in Equation (8) on planet radius. In the case of Kepler 11d, R_p drops below $6\,{R}_{\oplus }$ by $\approx 10\,{\rm{Myr}}$, whereas the radius of Kepler 11e remains $\gt 6\,{R}_{\oplus }$ well after $100\,{\rm{Myr}}$ (see Figure 7). In absolute terms, the planet loses the second largest amount of primordial H+He during the isolation phase ($0.11\,{M}_{\oplus }$). The difference in R_c between $t={t}_{\mathrm{iso}}$ and $t=8\,{\rm{Gyr}}$ is again caused by the pressure difference at the bottom of the envelope. The final values of a, ${M}_{p}$, and R_p all agree with measurements.

Kepler 11f: the observed mass of the planet is significantly smaller than those of its neighbors, possibly suggesting the formation at a smaller orbital distance. To achieve a correspondingly smaller core mass with our ${\sigma }_{Z}$, the planet's initial orbit is interior to the initial orbits of the other planets. Consequently, the planet requires a late start, ${t}_{i}\approx 1.1\,{\rm{Myr}}$, to avoid crossing other orbital paths.⁴ The assembly of the core takes place for the most part beyond $\approx 0.6\,{\rm{au}}$ and at disk temperatures of $T\lesssim 160\,{\rm{K}}$. As a result, the core composition is very similar to that of fully hydrated planetesimals (50% H₂O and 45% silicates by mass). Though rich in H₂O in relative terms, 47% of ${M}_{p}$ (the richest, in fact), because of its small mass the planet contains an amount of H₂O ($1.3\,{M}_{\oplus }$) greater only than the H₂O mass of Kepler 11b. During the evolution in isolation, stellar radiation strips off about 29% of the H+He gas accreted during the formation phase. At $8\,{\rm{Gyr}}$, the planet is left with a gas content of 2.5% by mass. The final planet radius and orbital distance agree with observations, whereas the final total mass, $2.81\,{M}_{\oplus }$, is close to the one-standard-deviation upper limit of the measured value ($2.8\,{M}_{\oplus }$).

Kepler 11g: similarly to its inner neighbor, the planet starts its simulated evolution at an advanced stage of the disk's life (${t}_{i}\approx 7\times {10}^{5}\,\mathrm{years}$) in between the initial orbits of Kepler 11c and d. Roughly 90% of its condensible inventory is collected beyond $\approx 1\,{\rm{au}}$, and all of the core is accumulated behind the ice condensation front. In fact, its formation occurs at disk temperatures below $\approx 170\,{\rm{K}}$ and its core has virtually the same composition as that of fully hydrated planetesimals. Partially hydrated planetesimals account for only 0.03% of the core mass. Mass loss during isolation is negligible, and the planet basically contracts at a constant mass. Again, the core shrinks somewhat as the envelope cools down and the applied pressure at R_c increases. The planet radius at $8\,{\rm{Gyr}}$ is just above the one-standard-deviation upper limit of the measured value ($3.39\,{R}_{\oplus }$), whereas a_f agrees with measurements. The planet mass is not really constrained by transit observations, but the value provided by this ex situ model is consistent with that from the in situ calculations in Table 2. Both simulations point at a planet somewhat less massive than Kepler 11e and of comparable mass to that of Kepler 11d.

The reference models use a standard value of $\varepsilon =0.1$ in the expression of ${\dot{M}}_{e}^{\mathrm{iso}}$, to account for re-radiation of the EUV and X-ray stellar flux by the envelope (see Section 2.3). Since this parameter is uncertain, the evolution in isolation was repeated by applying values of $\varepsilon =0.05$ and 0.15. Results for the evolution of R_p are plotted in Figure 11. The case of Kepler 11b is omitted for obvious reasons. The solid curves of each pair (upper curve for the smaller efficiency) bracket the excursion of R_p (colored regions). Within the considered range of ε, the impact on the radius at $8\,{\rm{Gyr}}$ is typically small. For Kepler 11g and e, the relative excursion of M_e amounts to about 1% and 9%, respectively, and the percentile variations of R_p are 0.4% and 2%. The excursion of M_e is larger for Kepler 11d (26%) and larger still for Kepler 11f ($\approx 40$%), owing to its small mass. The relative change in R_p, though, is $\lesssim 5$% for either. In all these cases, the simulated planets would still provide reasonable matches to the observed planets, with ${M}_{p}$ and R_p lying within or proximate to measurement ranges. Due to its vulnerability to mass loss, Kepler 11c represents the most extreme case, with M_e changing by a factor of three and R_p by $\approx 14$%, which would place the planet radius well outside of the measurement range.

Zoom In Zoom Out Reset image size

The reference models are based on time-dependent stellar properties for a $1\,{M}_{\odot }$ ([Fe/H] = 0.0) star computed by Siess et al. (2000). In order to determine the impact of the stellar evolution model, which enters the calculations through Equations (5), (8), and (16), a simulation of Kepler 11b was performed with the Yonsei-Yale stellar model (Spada et al. 2013) in Figure 2 (see the figure's caption and Section 2.1 for further details), starting from the same disk's initial conditions as in the reference model. The three equations depend on L_⋆ (${T}_{\star }^{4}{R}_{\star }^{2}\propto {L}_{\star }$, although Equation (16) separately depends on R_⋆ as well). The largest differences in stellar luminosity between the two stellar models occur for $t\lesssim 20\,{\rm{Myr}}$, i.e., during the disk-embedded phase and the early isolation phase of the planet. The Yonsei-Yale stellar model predicts lower luminosities, which tend to produce a somewhat cooler disk in regions where stellar irradiation is important in the energy budget of the gas (Equation (12)). The cooler temperature in turn reduces the disk thickness and hence affects the migration rate, since $| {da}/{dt}| \propto {(a/H)}^{2}$ (see Equation (22)). The cooler temperatures during the nebular stage of the planet's evolution also change the Bondi radius ${R}_{{\rm{B}}}$, and hence the planet accretion radius ${R}_{{\rm{A}}}$ in Equation (1). By using the same initial orbital radius, a_i, as in the reference model (see Table 4), the somewhat larger (in magnitude) migration speed requires a later start time, ${t}_{i}\approx 1.4\times {10}^{5}\,\mathrm{years}$ (the difference in luminosity is especially large for $t\lesssim 3\times {10}^{5}\,\mathrm{years}$, see Figure 2). The resulting model provides an isolation time of ${t}_{\mathrm{iso}}=4.1\,{\rm{Myr}}$ and values at $t={t}_{\mathrm{iso}}$ in very good agreement with those in Table 5. During the isolation phase, the surface temperature of the planet differs by $\lesssim 20\,{\rm{K}}$ relative to that of the reference model. The H+He envelope is entirely removed by $t\approx 30\,{\rm{Myr}}$, as in Figure 6. Because of the cooler gas temperatures and later start, the core contains a slightly larger H₂O fraction ($\approx 43$% versus $\approx 39$% by mass).

Since the isolation phase of Kepler 11b does not last long, a calculation with the Yonsei-Yale stellar model was also performed for the isolation phase of Kepler 11c. The properties of the model at ${t}_{\mathrm{iso}}$ are those of the reference model listed in Table 5. For most of the isolation phase ($t\gtrsim 50\,{\rm{Myr}}$), the stellar luminosity of the Yonsei-Yale stellar model is around 9% lower, compared to L_⋆ of the Siess et al. model, and thus ${T}_{\mathrm{eq}}$ is only marginally different. The calculation results in an envelope mass at $8\,\mathrm{Gyr}$ of $0.11\,{M}_{\oplus }$ and a planet radius of $2.88\,{R}_{\oplus }$. Both numbers are a little larger than the values in Table 4, but within the errors of the measured mass and radius of Kepler 11c. Clearly, the stellar evolution model does not affect the simulated planets significantly.

Dust grains, either entrained in the accreted gas and/or produced by ablation of accreted solids, can pollute planetary envelopes at temperature $T\lesssim 1500\,{\rm{K}}$ (assuming silicate grains). In such cases, dust opacity regulates the envelope cooling rate and hence the planet's contraction timescale (e.g., Pollack et al. 1996; Hubickyj et al. 2005). The calculations presented here are based on the opacity plotted in the top panel of Figure 1, which assumes a maximum grain radius of $1\,\mathrm{mm}$. To evaluate the effect of grain opacity in the envelope of the reference models, calculations were also performed by using the table plotted in the bottom panel of Figure 1, which assumes a maximum grain radius of $10\,\mathrm{mm}$. The ratio of the grain-dominated, Rosseland mean opacities in the two tables is about eight. The lower opacity provided by the size distribution with larger dust grains is expected to facilitate cooling and allow for higher gas accretion rates. The planet Kepler 11b, which accumulates the least massive H+He envelope, was simulated with the lower grain opacity. By using the same initial conditions for both disk and planet as for the reference model, the calculation provides values of ${\dot{M}}_{e}$ that are initially a few times as large as those of the reference model. The envelope mass achieves a maximum value of ${M}_{e}\approx 0.035\,{M}_{\oplus }$ at a time $t\approx 2.6\,{\rm{Myr}}$, after which the envelope loses mass to the disk as R_p tends to exceed the accretion radius ${R}_{{\rm{A}}}$. The planet attains isolation at a time ${t}_{\mathrm{iso}}=4.1\,{\rm{Myr}}$. Despite the larger envelope mass during the planet's early accretion history, the core mass at ${t}_{\mathrm{iso}}$ is again ${M}_{c}=2.1\,{M}_{\oplus }$ and ${M}_{e}=0.013\,{M}_{\oplus }$ (about twice as large as that of the reference model), which is entirely removed by stellar radiation by an age of $\approx 33\,{\rm{Myr}}$.

Although the opacity test indicates little impact on the simulated Kepler 11b planet, larger effects are to be expected for planets that accumulate higher H+He mass fractions. Indeed, the same test repeated for the most massive planet, Kepler 11e, results in an entirely different outcome. The planet reaches a crossover mass (${M}_{e}={M}_{c}$) of $\approx 7.2\,{M}_{\oplus }$ at $t\approx 2\,{\rm{Myr}}$ and enters the transition stage (see Section 2.1) of fast—and disk-limited—gas accretion by $t\approx 2.2\,{\rm{Myr}}$. On its track to becoming a Hot Jupiter, the planet reaches a mass of $\approx 0.7\,{{\rm{M}}}_{{\rm{J}}}$ by $t\approx 2.5\,{\rm{Myr}}$. The tendency to evolve into a giant planet, with the lower dust opacity of Figure 1 (bottom panel), is also obtained from in situ calculations. Ex situ models matching the observed radius of Kepler 11e with the lower opacity require a smaller core mass and, hence, a tighter initial orbit (${a}_{i}\approx 4.8\,{\rm{au}}$). A calculation resulting in ${R}_{p}=4.18\,{R}_{\oplus }$ at $8\,{\rm{Gyr}}$ provides a mass ${M}_{p}=5.55\,{M}_{\oplus }$, below the one-standard-deviation lower limit of the measured value of $5.9\,{M}_{\oplus }$. A calculation resulting in ${M}_{p}=5.96\,{M}_{\oplus }$ produces a radius ${R}_{p}\approx 4.5\,{R}_{\oplus }$, well over the one-standard-deviation upper limit of the measured value of $4.26\,{R}_{\oplus }$. Similarly, in situ formation models indicate that R_p can be matched with a mass M_p somewhat smaller than the measured one-standard-deviation lower limit.

Both in situ and ex situ models assume formation in a disk of $100\,\mathrm{km}$ radius planetesimals. A reduction in planetesimal size from $\approx 100$ to $\sim 1\,\mathrm{km}$ was studied by the authors in other contexts. First, the orbital eccentricities and inclinations of these bodies decrease, and second, the cross section for planetesimal capture in the planetary envelope (${{ \mathcal S }}_{\mathrm{eff}}$ in Equation (6)) increases (see D'Angelo et al. 2014, Figure 6). Both effects enhance the accretion rate of solids, ${\dot{M}}_{c}$. Whether the planetesimal size has an effect on the accreted gas mass depends on the ratio of the timescale for the buildup of the core to the disk lifetime. The core buildup timescale of the in situ models is already very short, thus the effect is negligible. In the ex situ models, the core accretion time would decrease (considerably in the case of small planetesimals). The faster core growth would increase the migration rate. The planet would be driven more quickly toward the inner disk regions, where solids' densities are higher but the solids' mass available for accretion is lower. Therefore, the outcome is difficult to predict. We performed tests for Kepler 11d with planetesimals of 10 and $1\,\mathrm{km}$ in radius. Compared to the reference case discussed above, the final core mass increased by only about 10% and 15%, respectively.

The data listed in Tables 4 and 5 indicate that there can be marginal changes in core radii during the isolation phase. As mentioned in Section 4.2, these changes arise from differences in boundary pressures caused by contraction, resulting in core compression. Such small changes, however, may be offset by temperature effects through the core. In the calculations discussed above, core thermodynamics is ignored and ex situ models include it only for phase transitions within the silicate mantle (see Section 2.2 and Appendix A). The assumption of neglecting the core thermal stratification is based on previous studies, which argued that thermal pressure should not significantly affect the core radius at these core masses (e.g., Valencia et al. 2006; Seager et al. 2007; Sotin et al. 2007; Sohl et al. 2014b). Nonetheless, for completeness, improved structure models are presented and discussed in Appendix C. These include temperature stratification, energy transfer, and temperature-dependent EoS in a self-consistent fashion; they also include additional material phases. Thermal effects (in H₂O-rich cores) may be relatively important for the determination of R_c during the early stages of the isolation phase ($t\lesssim 0.5\,{\rm{Gyr}}$, see Figure 12). However, at an age of $\approx 8\,{\rm{Gyr}}$, the differences in R_c between thermal and isothermal cores of the simulated Kepler 11 planets are small, typically $\lesssim 0.5$% in most cases and somewhat less than 2% for Kepler 11c. Assessing the impact of thermodynamics for the case of Kepler 11b is more difficult, because the expected steam envelope is not modeled. If the gas pressure and temperature at the bottom of the envelope were, respectively, $\approx 1\,\mathrm{GPa}$ and $\approx 1000\,{\rm{K}}$, then R_c would be 5% larger than the value listed in Table 4 (see the discussion in Appendix C).

Zoom In Zoom Out Reset image size

During isolation, planets gradually cool down. The pressure P_c at the bottom of the envelope varies by factors of the order of unity, whereas temperature variations are larger. To make a simple assessment of the impact on R_c of the varying conditions at the core-envelope boundary as planets cool, thermal structure calculations of the cores (see details in Appendix A) were performed by applying pressure and temperature at the core surface during the isolation phase. Figure 12 indicates that there can be significant variations in R_c, up to $\approx 0.34\,{R}_{\oplus }$. The symbols in the figure are color-coded according to the core surface temperature, T_c, as obtained from ex situ calculations. However, in the figure, 92%–94% of the total difference ${\rm{\Delta }}{R}_{c}$ along each curve is due to the contraction of the cores' H₂O shell. Additional details are given in Appendix C. Assuming ${\rm{\Delta }}{R}_{p}\approx {\rm{\Delta }}{R}_{c}$, the inflated cores may result in planet radii larger by $\approx 3$% at around the isolation time, which may in turn enhance the evaporative mass-loss rate by $\approx 10$% during early isolation times. Results in Figure 11 suggest that the impact on R_p at $8\,{\rm{Gyr}}$ may be relatively small.

Figure 13 illustrates how some properties of the models for Kepler 11d and e at $8\,{\rm{Gyr}}$ depend on the surface density of solids, ${\sigma }_{Z}$, and the initial orbital radius, a_i (see the figure's caption for further details). Where relevant, in situ models are included as well. In the figure, ${\sigma }_{Z}^{\mathrm{ref}}$ and ${a}_{i}^{\mathrm{ref}}$ are the reference values of the models discussed in Sections 3 and 4. In the calculations of the top panel, the initial gas density Σ is re-scaled, so to keep the ratio ${\rm{\Sigma }}/{\sigma }_{Z}$ fixed. The planet mass increases monotonically as ${\sigma }_{Z}$ and a_i increase. In either case, this is a consequence of the larger core mass that facilitates gas accretion, especially at early times. The final envelope mass of in situ models increases monotonically, as it generally (but not always) does also in ex situ models. In the latter calculations, because of the different accretion history (and initial condition a_i), the final orbital radius a_f varies, affecting the evolution of M_e during isolation. The radius R_p tends to grow with the ratio ${M}_{e}/{M}_{c}$, which is not always a monotonic function of ${\sigma }_{Z}$ or a_i in ex situ models (see the left panels of Figure 13). Because of the small uncertainties on the observed radius, only initial conditions in the neighborhood of those adopted for the reference models can match observations. The final orbital distance varies by a factor of up to $\approx 1.8$ in the ex situ calculations of the top panels and $\approx 1.3$ in those of the bottom panels (relative to the reference values ${a}_{f}=0.156\,{\rm{au}}$).

Zoom In Zoom Out Reset image size

The in situ reference models of Section 3 and those represented in Figure 13 assume an initial disk's gas-to-solid mass ratio of 200. Experiments conducted on Kepler 11e, by using ${\sigma }_{Z}=4600\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ (as in the reference case) and varying the mass ratio from 50 to 400, resulted in mass and radius changes (at $8\,{\rm{Gyr}}$) of $\lesssim 10$%. For a gas-to-solid mass ratio of 50, ${M}_{p}=7.7\,{M}_{\oplus }$ and ${R}_{p}=3.98\,{R}_{\oplus }$ whereas these values increase to $8.46\,{M}_{\oplus }$ and $4.44\,{R}_{\oplus }$, respectively, for a mass ratio of 400. The planet mass is still within measurement errors, while the planet radius lies outside the observed range by at most a few percent.

Both in situ and ex situ simulated planets result in radii, masses, and orbital distances in agreement with measured values at the estimated age of the system. Therefore, it is not possible to distinguish between the two modes of formation from these final properties. The two formation scenarios do, however, provide entirely different perspectives of the environment in which the planets grew, of their compositions, and interior structures.

Figure 14 shows the initial surface density of the gas applied to the in situ (circles) and ex situ models (blue line). The figure also illustrates other reference surface densities (Weidenschilling 1977; Hayashi 1981; Davis 2005; Chiang & Laughlin 2013; Schlichting 2014). As noted above, Σ at t = 0 of the ex situ simulations fortuitously matches the surface density constructed by Davis (2005). It should be pointed out that, among these reference density distributions, only those of Chiang & Laughlin (2013) and Schlichting (2014) were explicitly derived for close-in extrasolar planets. The other two are meant to apply to the solar system and are thus simple extrapolations at the short distances from the star of Kepler 11 planets. The dashed line in the figure indicates the gas density above which the disk would be gravitationally unstable to axisymmetric perturbations according to the Toomre stability criterion (Durisen et al. 2007, pp. 607–622), using a constant gas temperature of $1000\,{\rm{K}}$. Even by applying $T=1000(0.1\,{\rm{au}}/r){\rm{K}}$, appropriate for H/r nearly constant, the initial Σ inferred from in situ models would still be stable (although only marginally stable to non-axisymmetric perturbations).

Zoom In Zoom Out Reset image size

In situ formation requires a very large ${\sigma }_{Z}$ within $\approx 0.5\,{\rm{au}}$ of the star, as illustrated in Figure 14 (see also Bodenheimer & Lissauer 2014). The models discussed here suggest densities of solids between 700 and $1.45\times {10}^{4}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ (see Table 2). If most of these solids had to form locally, the initial gaseous mass of the region had to be accordingly large. Given the short accretion timescales involved in the growth process, $\sim {10}^{4}\,\mathrm{years}$ (see Figures 3 and 4), most of the solids necessary for core assembly had to be available prior to the beginning of this process. Therefore, a gradual replenishment of the solids' reservoir from larger radii might not be a viable alternative to a large ${\sigma }_{Z}$. Moreover, this possibility would likely lead to a hierarchical system of planets with inwardly decreasing masses, which is inconsistent with the masses of Kepler 11f and e (and probably of Kepler 11g as well).

A disk model aimed at mimicking the evolution of the initial in situ Σ shown in Figure 14 is presented in Figure 15. The disk has an initial mass of $\approx 0.18\,{M}_{\odot }$ within $\approx 70\,{\rm{au}}$. Both the kinematic viscosity of the gas and the emission rate of stellar EUV ionizing photons were chosen to induce photo-evaporation of the inner disk at around the age of $3.5\,{\rm{Myr}}$ (see Section 3). The region $r\lesssim 1\,{\rm{au}}$ in the disk in Figure 15 is dispersed between $t\approx 3.5$ and $\approx 3.7\,{\rm{Myr}}$. Gas temperatures are initially high, as can also be realized from simple arguments based on energy balance. At high densities, the disk is optically thick in the vertical direction and the main source of energy is viscous heating. Hence Equation (12) reduces to ${Q}_{\nu }={Q}_{\mathrm{cool}}$, which becomes

$$\begin{eqnarray}&&{\sigma }_{\mathrm{SB}}\,{T}^{4}=\displaystyle \frac{27}{128}\left({\rm{\Sigma }}{\kappa }_{{\rm{R}}}+\displaystyle \frac{8}{3}\right)\nu {\rm{\Sigma }}{{\rm{\Omega }}}^{2},\end{eqnarray} \tag{ 25 }$$

where ${\rm{\Sigma }}{\kappa }_{{\rm{R}}}/2$ is the optical depth of the disk's mid-plane. Since the initial accretion rate through the disk in this case is $3\pi \nu {\rm{\Sigma }}\sim {10}^{-7}\,{M}_{\odot }\,{{\rm{yr}}}^{-1}$, Equation (25) and Figure 15 imply that, at $r\approx 0.5\,{\rm{au}}$, temperatures are initially in excess of $2500\,{\rm{K}}$ and drop below $1000\,{\rm{K}}$ at $t\gtrsim 1\,{\rm{Myr}}$. At $r\approx 0.1\,{\rm{au}}$, temperatures become $\lesssim 1000\,{\rm{K}}$ at $t\gtrsim 2.5\,{\rm{Myr}}$. These high gas temperatures would only allow for the presence in the disk of metals and highly refractory solid species, which would be reflected by the compositions of the cores. The high temperatures also prevent the massive inner disk from becoming gravitationally unstable (Durisen et al. 2007, pp. 607–622). Lowering the kinematic viscosity of the gas would somewhat reduce the gas temperatures ($T\propto {\nu }^{1/4}$) but would also extend the disk lifetime (by halving ν it would take over $5\,{\rm{Myr}}$ to photo-evaporate the inner disk). The effects of the initially high disk temperature and of the disk evolution on the in situ formation process require further calculations.

Zoom In Zoom Out Reset image size

Ex situ formation can occur in low-mass disks, in which the initial ${\sigma }_{Z}$ increases from $\approx 5\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ at around $7\,{\rm{au}}$ to $\approx 10\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ at $0.5\,{\rm{au}}$, the disk region that provides the vast majority of the solids to assemble all the planetary cores. Initial gas densities at these radial distances are in the range from a few times $100\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ to $\approx 1.5\times {10}^{3}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ (see Figure 5). Cores are assembled over timescales of the order of $1\,{\rm{Myr}}$ (see Figures 6 and 7). Planets experience relatively low temperatures during their formation phase and spend a long enough time behind the ice condensation line to become enriched in H₂O and other volatile substances (if available).

The percentile composition of the six planets is represented in the histograms of Figure 16, for in situ as well as ex situ models (see the figure's caption for details). Recall that in situ simulations assume a fixed core composition of iron and silicates of 30% and 70% by mass, respectively. In ex situ models, H₂O is assumed to be all condensed in the core, although part of it should be in gaseous form mixed with H+He gas, released by ablating solids passing through the envelope. Both in situ and ex situ calculations predict that the percentile H+He mass of Kepler 11c is comparable to that of the smaller Kepler 11f, because most of the envelope is lost after formation. Kepler 11b simulated in situ acquires a much larger envelope than does the ex situ simulated planet, because the core mass approaches its final value at very early times (compare the growth of M_e in Figures 3 and 6).

Zoom In Zoom Out Reset image size

In order to reconcile observed and simulated radii of Kepler 11b, the in situ formation scenario requires that gas is sequestered in the core during formation and released afterwards. Outgassing of at least $\approx {10}^{-3}\,{M}_{\oplus }$ of hydrogen would be sufficient to account for the observed radius. However, outgassing would still need to compete against evaporative gas loss, which would operate at rates between $\approx {10}^{-12}$ and $\approx {10}^{-11}\,{M}_{\oplus }\,{{\rm{yr}}}^{-1}$. If outgassing began right after the removal of the primordial H+He envelope (at $t\approx 40\,{\rm{Myr}}$) and continued for the age of the planet, the (average) outgassing rate would only need to be marginally higher than the (average) evaporative gas loss rate. This possibility, however, implies that the amount of gas sequestered in the core during formation ought to be significant in relative terms, between $\approx 0.01$ and $\approx 0.1\,{M}_{\oplus }$. The details regarding the processes of sequestration and outgassing remain to be investigated, but the result could have important implications concerning the in situ versus ex situ formation of Kepler 11b.

The ex situ formation scenario predicts the presence of a steam atmosphere (not modeled here). Based on previous assessments of the planet mass and radius (Lissauer et al. 2011a), Lopez et al. (2012) estimated that the planet radius can be matched by a water-world whose composition is about 40% H₂O by mass, comparable to the H₂O content of the simulated ex situ planet (see Table 4). Assuming a hydrostatic and adiabatic atmosphere with an ideal EOS for H₂O, the mass M_e necessary to match the observed radius is $\approx 7\times {10}^{-3}\,{M}_{\oplus }$, though non-ideal effects in the EOS may reduce this estimate. Because of the larger mean molecular weight, evaporative mass loss of H₂O by stellar radiation should be less significant than loss of lighter elements. Additionally, steam can be replenished by the condensed core. Since hydrogen can mix with H₂O at the pressures and temperatures of planetary interiors (Soubiran & Militzer 2015), hydrogen may be present in the atmosphere of Kepler 11b as well.

Transit observations only place an upper limit of $\approx 25\,{M}_{\oplus }$ on the mass of Kepler 11g (L13). Assuming negligible amounts of gas in the envelope, an ex situ type composition of 5% iron, 45% silicates, and 50% H₂O by mass, would result in a radius of $2.8\,{R}_{\oplus }$, short of the observed value ${R}_{p}=3.33\,{R}_{\oplus }$. An in situ type composition would result in an even smaller radius ($2.3\,{R}_{\oplus }$). Hence, the planet must bear a gaseous envelope.

Both in situ and ex situ models indicate that Kepler 11g is somewhat less massive than Kepler 11e and of comparable mass to Kepler 11d. It is reasonable to imagine a scenario in which a lower core mass may attract a sufficiently large envelope to account for the observed radius. In fact, another simulation (with different start conditions) of Kepler 11g results in a radius only marginally smaller than that observed and ${M}_{p}\approx 4.7\,{M}_{\oplus }$. Even smaller masses might be attained, e.g., by lowering the dust opacity in the planet's envelope. However, the opposite scenario of a much more massive planet and a smaller ${M}_{e}/{M}_{p}$ ratio may be more difficult to realize.

Figure 17 shows outcomes from additional ex situ simulations of Kepler 11g. In these models, the start conditions of the reference model (${t}_{i}$, a_i, and initial ${\sigma }_{Z}$) were changed with the aim of obtaining different planet masses and gas-to-condensible mass ratios. The final orbital radii are within roughly 10% of the observed orbital distance. Circles represent simulations that apply the grain opacity of the top panel of Figure 1, as in reference models. Data points are color-coded by their relative gas content, ${M}_{e}/{M}_{p}$. The solid lines indicate the core radius corresponding to the interior composition of in situ (30% iron and 70% silicates) and ex situ models (see Table 4), assuming a pressure at R_c of $1\,\mathrm{GPa}$ and $T={T}_{\mathrm{eq}}$. The observed radius of the planet, within measurement errors, is indicated by the gray-shaded area. The data points include the in situ simulation of Table 2 as well. This model has a comparatively large gas content for its radius, owing to the smaller value of R_c (i.e., the difference ${R}_{p}-{R}_{c}$ is in fact comparably large).

Zoom In Zoom Out Reset image size

From the data plotted in Figure 17, it appears difficult to reconcile observed and simulated radii for planets more massive than about $7\,{M}_{\oplus }$. In fact, when ${M}_{p}\gtrsim 7\,{M}_{\oplus }$, the H+He mass fraction ${M}_{e}/{M}_{p}$ tends to exceed $\approx 16$%, resulting in too large a radius. In the calculations reported herein, the value of M_e at the end of the formation phase is set by dispersal of the gaseous disk. The model with the largest mass in Figure 17 has ${M}_{e}\gtrsim 2.5\,{M}_{\oplus }$ (${M}_{e}/{M}_{p}\gtrsim 0.28$) and R_p far larger than observed. If under different circumstances the accretion of gas was starved by disk dispersal at an earlier age and R_p at $8\,{\rm{Gyr}}$ was similar to the observed value, the planet mass would not exceed $\approx 8\,{M}_{\oplus }$ (an envelope mass somewhat smaller than $1\,{M}_{\oplus }$ appears necessary to account for the observed radius).

According to L13, an ${M}_{p}\approx 15\,{M}_{\oplus }$ planet with a 6.6% gas content by mass can match the radius of Kepler 11g. It may be possible to build a large core mass without attracting an excessively large envelope by augmenting the opacity in the envelope, which operates to slow contraction and inhibits gas accretion. This possibility was tested by using dust opacities calculated from an interstellar medium-type grain distribution, with dust grains of up to $1\,\mu {\rm{m}}$ in radius (for details, see D'Angelo & Bodenheimer 2013). Results from this second set of models are represented by diamonds in Figure 17. As expected, for a given ${M}_{p}$, the ratio ${M}_{e}/{M}_{p}$ is smaller than in the other set of simulations, yielding smaller planetary radii. Nonetheless, as ${M}_{p}$ grows beyond $\approx 8\,{M}_{\oplus }$, ${M}_{e}/{M}_{p}\gtrsim 0.2$ and R_p becomes too large to match observations. Clearly, other scenarios can be envisaged (e.g., most of M_c is accumulated toward the end of the disk lifetime, delaying envelope growth), although they must be compatible with the formation of the inner planets. Within the assumptions adopted herein, both in situ and ex situ models argue in favor of a mass smaller than or comparable to Kepler 11e's mass.