DYNAMICAL ACCRETION OF PRIMORDIAL ATMOSPHERES AROUND PLANETS WITH MASSES BETWEEN 0.1 AND 5 M_⊕ IN THE HABITABLE ZONE

For the numerical description of planetary atmospheres, some assumptions and simplifications of the physical scenario and the geometry of the problem can be adopted.

First, the Hill radius R_Hill defines the sphere where the gravitational potential of the planet dominates over the influence of the central host star (e.g., Mizuno et al. 1978):

$$\begin{eqnarray}&&{R}_{\mathrm{Hill}}\;=\;{\left(\displaystyle \frac{{M}_{\mathrm{pl}}}{3{M}_{\star }}\right)}^{1/3}a.\end{eqnarray} \tag{ 1 }$$

Here, M_pl and M_⋆ are the masses of the planet and of the host star, respectively, and a is the semimajor axis of the planet. The Hill radius, by specifying the gravitational sphere of influence, gives the largest possible spatial extension of the planetary atmosphere, and for radii with $r\ll {R}_{\mathrm{Hill}}$ the gravitational field can be well approximated in spherical symmetry.

Another important consideration is the vertical thickness of the protoplanetary disk. The vertical structure of the disk is essentially determined by hydrostatic equilibrium and can be quantified from the pressure scale height H_p of the vertical density stratification. A reasonable approximation of an embedded planetary atmosphere in spherical symmetry thus requires that the vertical pressure scale height of the disk is larger than the spatial extent of the atmosphere.

Combining these two criteria, we assume in this study that $r\leqslant {R}_{\mathrm{Hill}}\lt {H}_{p}$ throughout, justifying a spherical symmetric description. Spherical symmetry also implies that the effect of the angular momentum of the gas flow with respect to the planetary core is neglected.

Apart from the geometry, it is also possible to simplify the physical description of the problem. The full time-dependent physical scenario can be simplified by considering only stationary solutions, skipping the time derivatives in the equation of motion and in the equation of energy. The former corresponds to the assumption of dynamic (i.e., hydrostatic) equilibrium, and the latter corresponds to the adoption of thermal equilibrium. Such stationary models are the simplest way to study the structure of planetary atmospheres. While these simplifications are probably justified during some parts of the planetary evolution, the thermal relaxation timescale (Kelvin–Helmholtz timescale) can become significant for massive atmospheres, and dynamical processes can also have long relaxation timescales compared to evolutionary timescales of the planetary system (Stökl et al. 2015).

In spherical geometry and in the stationary limit, the structure of planetary atmospheres is described by the well-known stellar structure equations, and with the assumption of suitable boundary conditions for the planetary core and the surrounding protoplanetary nebula, these equations can be readily integrated. Such stationary models of disk-accumulated atmospheres around terrestrial planets have been calculated by Hayashi et al. (1979), Nakazawa et al. (1985), and Rafikov (2006) based on analytically prescribed dust opacities, and by Ikoma & Genda (2006) using realistic data for the equation of state and for dust and gas opacities.

For this study, we employ a radiation hydrodynamics scheme to compute time-dependent, that is, nonstationary, models of planetary atmospheres. A similar scheme has been used by Wuchterl (1991a, 1991b) to study the dynamic behavior of a planetary envelope with a critical core mass, i.e., at the verge of self-gravitational instability. However, Wuchterl (1991a, 1991b) did not investigate less-massive cores (as in the present study) and used this scheme for short-duration dynamical phases only and did not attempt to cover the gas-accumulation phase.

1.1.1. Planetary Luminosity and Accretion of Planetesimals

For most models in our grid of stationary models, the time-dependent calculations demonstrate that the envelopes are in fact dynamically stable and return to their initial configuration after a perturbation. In these cases, the boundary conditions directly, though not uniquely (as discussed above), determine the atmospheric structure. Here in particular, the boundary conditions are expressed through ${T}_{\mathrm{surf}}$ and M_core at the inner boundary and ρ and T at the outer boundary. Therefore, a numerical simulation starting from an arbitrary initial configuration (e.g., spatially constant ρ and T) will ultimately converge to the stable solution, or one of the stable solutions in the case of multiplicity, provided it does not enter self-gravitational collapse during that evolution.

In Figure 1, we have illustrated that, within a certain parameter range, multiple stationary solutions exist for one surface temperature. Figure 2 focuses on the relevant range for the model series with a core mass of 0.6 M_⊕. The results of the time-dependent simulations now show that, from these multiple solutions, only those in the upper and lower branches are stable, while the solutions in the branch with $\tfrac{{{\rm{d}}{L}}_{\mathrm{planet}}}{{{\rm{d}}{T}}_{\mathrm{surf}}}\lt 0$ are unstable. In the simulation, an unstable atmospheric solution (green dot in Figure 2) evolves, depending on the initial numerical perturbation, toward either the upper branch or the lower branch (indicated by dashed lines). Finally, the atmosphere settles to the corresponding stationary solutions on the stable branches for the same core temperature (blue and red dots).

Zoom In Zoom Out Reset image size

The differences between the three stationary solutions in Figure 2 for one and the same surface temperature are illustrated in Figure 3. There are two stable solutions: a thin solution with a low overall density and a high luminosity, and a dense solution characterized by a more compact structure and a low luminosity.

Zoom In Zoom Out Reset image size

Figure 3 also illustrates snapshots from the dynamical evolutions leading from the unstable model to either of the stable configurations. Starting from the unstable state, it first takes some time before the numerical inaccuracies (or initial perturbation) grow to an amplitude that noticeably changes the envelope structure. The more accurately balanced the initial model is, the longer this growth takes. For the calculation illustrated in Figure 3, the growth time is on the order of 10 years. In order to predetermine the direction of evolution, we artificially perturbed the initial model by changing the density in the innermost cell by ±0.1%, which only marginally reduced the growth time.

In the cases where the evolution goes toward the thin configuration, the exponentially growing perturbation leads to a dynamical outflow on the dynamical timescale of about 10⁶ s (12 days). The outflow propagates as a dynamical wave that is partially reflected at the outer boundary (see line for 11.75 years in Figure 3) and subsequently travels back and forth through the envelope with decreasing amplitude until it dies away after about four to five cycles.

If the evolution goes in the direction of the dense stable solution, the atmospheric changes happen at a much slower pace. The growth of the initial perturbation leads to an inflow of disk gas into the atmosphere that is very similar to the accumulation of atmosphere as described in Section 4. In this phase, the timescale of evolution is given by the ability of the atmosphere to assimilate and radiate off the potential energy released in the gas-accretion process. In the example illustrated in Figure 3, it takes about 30 Myr until the envelope finally settles to the dense, stable solution.

This, however, is obviously not a realistic result, as an evolution starting from an unstable stationary configuration is a purely hypothetical scenario. In the next section, we attempt to devise a more realistic setup for the accumulation and evolution of embedded primordial atmospheres.

Figure 4 illustrates the evolutionary tracks of planetary atmospheres around planetary cores between 0.1 and 5 M_⊕. The evolution tracks start at the right margin at T_surf = 10,000 K and zero luminosity; the initial surface pressure is that of the unperturbed ambient nebula, corresponding to habitable zone values of T = 200 K and ρ = 5 × 10⁻¹⁰ g cm⁻³. With the onset of the simulation, very quickly a radiative flux is established, the luminosity peaks at a high level, and the core starts to cool down. The intensive radiation field heats the ambient gas and produces a thin, hot, and expanding bubble around the planetary core. This outflow depletes the gas in the computational domain compared to the unperturbed disk nebula, causing a density inversion where it joins the invariable outer boundary conditions of T = 200 K and ρ = 5 × 10⁻¹⁰ g cm⁻³, equivalent to the situation in the stationary models of luminous atmospheres, as discussed in Section 2.

Zoom In Zoom Out Reset image size

After the initial dynamic phase (lasting less than 1 year for all models), the evolution tracks align along the sequence of stationary models, indicating that, from that point on, the atmospheric structure is always close to the hydrostatic equilibrium.

The high luminosity in this first evolution phase implies that the planetary core cools down quickly; the initial dynamic high-luminosity phase lasts less than 100 years (green crosses) in all simulation runs. The later planetary evolution is therefore independent of the exact value of the initial core temperature. As mentioned above, the high luminosities at this phase are connected to the simple core-temperature model. Consideration of convective energy transport in the interior of the core would probably reduce the maximum luminosities to values on the order of about 10²⁷–10²⁸ erg s⁻¹. The initial high-temperature, high-luminosity phase should therefore not be considered as a realistic scenario but as a numerical approach to correlate the time-dependent evolution tracks with the stationary atmosphere models in a systematic way.

Finally, when the planets have cooled down to a certain temperature—characteristic of the respective core mass—the formation of a substantial atmosphere through accumulation of disk gas sets in: disk gas flows inward over the outer boundary, the atmospheric density and mass increase, the atmosphere becomes optically thick, and consequently the planetary luminosity is much reduced. The cooling of the planetary core virtually stops at this point, as, first, the planetary luminosity is much lower and, second, the luminosity is balanced with the release of gravitational energy from the inflowing gas.

With respect to the sequence of stationary models in Figure 4, the position of the atmosphere has moved to the left along the upper branch until the model sequence bows in to the unstable branch (with $\tfrac{{{dL}}_{\mathrm{planet}}}{{{dT}}_{\mathrm{surf}}}\lt 0$). At this point, the model moves, essentially at constant core temperature, toward the lower stable branch, which, however, is above the limit of self-gravitational collapse only for the least-massive cores.

In this phase, the release of gravitational energy from the inflow of gas and the consequent contraction of the atmosphere dominates the energy balance of the atmosphere. As illustrated in Figure 5, the energy generation is broadly distributed over the atmospheric structure. In the example in Figure 5, about one-half of the total energy gained from gas accretion (∼11 × 10²³ erg s⁻¹) ends up as thermal energy in the atmosphere, while the remaining part (∼6 × 10²³ erg s⁻¹) constitutes the planetary luminosity during this stage. The energy split between atmosphere heating and luminosity is variable and depends on the model and the evolution time.

Zoom In Zoom Out Reset image size

As there is only little energy generation from gas accretion in deep atmosphere layers, the temperature stratification in this part of the atmosphere is essentially isothermal or even exhibits a small temperature inversion. Accordingly, there is no convective and, in an optically thick environment, only little radiative transport. (In Figure 5, there is a tiny inward-directed radiative flux near the surface.) In effect, these layers thus form a thermal insulation for the planetary core, and the core temperature is almost constant for large parts of the planetary evolution in the time-dependent simulations. As the phase of gas inflow and atmospheric growth lasts much longer than the initial evolution, the atmosphere remains in a narrow range of core surface temperatures for most of the time during the nebula-embedded atmospheric evolution.

In the gas-accumulation phase, the atmospheric evolution timescale becomes comparable to or larger than the evolution timescale of the protoplanetary disk. Therefore, in a realistic scenario, the further development of the planetary atmosphere would be determined by the changing nebula properties, the eventual evaporation of the disk, and irradiation from the host star. If we neglect these effects and assume an indefinitely constant disk environment, the hypothetical long-term outcomes of the evolution of the two smallest cores in Figure 4 are stationary solutions on the lower stable branch. These solutions are characterized by the luminosity provided from radiogenic energy production (which is also assumed to be indefinitely constant).

For more massive cores, there are no such stationary solutions because of the onset of self-gravitational instability for low luminosities. Therefore, when the evolution track of a growing atmosphere crosses that limit of self-gravitational stability, the atmosphere goes into a runaway collapse that would eventually lead to the formation of a gas planet.

In our simulations, we never obtained a dynamical pulsation-like behavior as described by Wuchterl (1991b) and Pečnik & Wuchterl (2007). The reason for this discrepancy is unclear.

In our calculations, the limit of runaway collapse is at about ∼0.5 M_⊕, but in general this limit depends on disk properties and the physical parameterization. More important, however, are the timescales of evolution: if we consider a typical disk lifetime of a few megayears, only the most massive cores in our sample would undergo a runaway atmospheric collapse; for the less massive cores, the collapse is a purely hypothetical result because the evaporation of the protoplanetary disk would cut off the supply of gas.

The release of gravitational energy in the collapse produces a steeply rising planetary luminosity. In the evolution tracks in Figure 4, the runaway collapse appears as a marked gain in core temperature as the excessive heating of the atmosphere by the inrush of gas and the corresponding atmospheric compression are transmitted downward to the core. This heating of the core from above is a scenario where the simple, constant heat capacity model of the core particularly shows its limitations, leading to clearly unrealistic results. Generally, the numerical results for the more advanced collapse stages have to be viewed with reservations. While the numerical setup, as discussed in Section 1.1, allows a consistent, time-dependent simulation of the atmospheric evolution and gas accumulation right into the onset of the self-gravitational collapse, it does not allow an adequate description of the more advanced collapse processes. In the course of the collapse, some fundamental modeling assumptions are violated. In particular, spherical symmetry is no longer a good approximation when the atmosphere becomes large compared to the thickness of the protoplanetary disk, and for high radial velocities when the angular momentum of the flow with respect to the planetary core becomes important. The high gas accretion rate in the runaway collapse would also quickly exhaust the gas supply in the circumstellar disk (possibly opening an annular gap), and the further evolution would be determined by mass (i.e., angular momentum) transport mechanisms in the disk.

Figure 6 illustrates the total atmospheric mass as a function of time in the planetary evolution simulations with core masses between 0.1 and 5 M_⊕. The upward S-curve in the lines between 10 and 100 years corresponds to the settling of the atmospheric structure into a quasi-hydrostatic configuration, as discussed above. This effect is more pronounced for more massive planetary cores.

Zoom In Zoom Out Reset image size

Then follows a prolonged smooth phase of atmospheric growth and gas accumulation (at essentially constant core temperature). Cores more massive than 0.5 M_⊕ at some point enter a runaway collapse, while less massive cores eventually settle into a stationary solution. The filled circles on the lines indicate the point where the atmospheric mass equals the core mass. This crossover mass can serve as a rough estimate for the onset of gravitational instability (Bodenheimer & Pollack 1986; Ikoma et al. 2000), and in our results these points indeed approximately coincide with the beginning of the runaway collapse. The time at which the collapse occurs depends on the mass of the planetary core. A comparison with the typical lifetime of a protoplanetary disk of 1–10 Myr (gray shaded area in Figure 6) shows that only the atmospheres around the 4 M_⊕ and 5 M_⊕ cores will actually undergo collapse. For the less massive cores, the atmospheric growth would come to a halt once the protoplanetary disk evaporates, producing planets with extended hydrogen envelopes around a rocky core, often referred to as "mini-Neptunes."

For the more massive cores, where a comparison is possible, the gravitational accumulation of disk gas into an envelope and the onset of runaway accretion, as illustrated in Figure 6, agree qualitatively with calculations of core-collapse gas planet formation as in Ikoma et al. (2000), Ikoma & Hori (2012), or Broeg & Benz (2012). A detailed quantitative comparison, however, is not sensible because of the different physical assumptions (Jupiter versus Earth-like orbit) and numerical description, in particular the dust opacity data. As the atmospheric structure is close to a hydrostatic configuration during most of the evolution, the hydrodynamic (versus quasi-static) description does not seem to be essential for the calculation of the atmospheric accumulation timescale. During the accretion process, there is only little physical interaction between the planetary core and the atmosphere. The numerical treatment of the planetary core and its initial temperature are therefore also not crucial for the description of the atmospheric growth phase.

Table 1 compiles the evolution time up to the runaway collapse for different core masses. The growth times agree very well with those in Ikoma et al. (2000), but this might be coincidental considering the differences in the numerical description. The initial dynamic adjustment phase in our simulations does not make an appreciable difference in the long-term accretion evolution. Except for a somewhat different scaling on the time axis, our results are very similar to the curves given in Ikoma et al. (2000) for the case with constant core mass (${\dot{M}}_{\mathrm{core}}\;=\;0$).

In order to verify the dependence of atmospheric growth and evolution on the initial temperature of the planetary core, we performed test runs with a lower initial core temperature of 1000 K.

Figure 7 compares the temporal evolution of atmosphere mass and core temperature for cores with 1 M_⊕ and initial temperatures of 1000 K and 10,000 K (reference model). The thermal insulation of the core once an optically thick atmosphere has formed, as discussed above, suggests that the atmospheric evolution is quite similar in the more evolved stages. In fact, the curves for the atmospheric mass soon become indistinguishable after the initial dynamic phase. In the 10,000 K case, the high initial luminosity causes a steep decline of the core temperature so that it is already down to about 6500 K at the left margin of Figure 7 at t = 1 year. After the short high-luminosity phase, the core gets thermally insulated and from then on essentially remains at its characteristic temperature, in this case about 3200 K (compare Figure 4). Only later in the evolution with an increasing gas accretion rate does heating of the core from overlying atmospheric layers become important, and the core temperature starts to rise.

Zoom In Zoom Out Reset image size

The evolution of the atmosphere in the case of an initial core temperature of 1000 K follows the same lines except that there is no initial high-luminosity phase. Accordingly, the core remains constant at its initial temperature throughout most of the atmospheric growth phase. In this case, the heating effect from gas accretion in the late evolution toward runaway collapse is more pronounced because of the larger temperature difference between the core and the atmosphere.

The different temperature gradients above the core surface are apparent in Figure 8, which compares the structure of the atmosphere in the 10,000 and 1000 K runs at a same simulation time of ∼5.2 × 10⁵ years. This is the same time slice as used in Figure 5 to illustrate the energy fluxes. Apart from the different temperature inversion above the core surface and a corresponding density spike, the atmosphere structures are very similar.

Zoom In Zoom Out Reset image size

For comparison, we have also included in Figure 8 a stationary atmosphere model using the same luminosity as in the time-dependent models. The relevant luminosity at the outer boundary of the time-dependent atmosphere can be read from Figure 5 (L = 4.84 × 10²³ erg s⁻¹). In the stationary model, the luminosity is necessarily assumed to be constant throughout the atmosphere.

The structure of the stationary model is in surprisingly good agreement with the time-dependent simulations. This can be understood, first, from the insensitivity of the outer envelope to the temperature of the planetary core, and, second, from the convective transport in the inner atmosphere layers, which fixes the temperature gradient close to the adiabatic value. This good agreement, however, relies on the suitable specification of the planetary luminosity for the stationary model. As can be read from Figure 4, the planetary luminosity varies by about three to four orders of magnitude during the embedded evolution phase of the planet. The consistent determination of the properties of a planetary atmosphere therefore requires consideration of the whole time-dependent problem.

Even though envelope growth is a time-dependent process and the final amount of atmosphere depends on the duration of the disk-embedded phase, it only takes a comparatively short period to accumulate a significant primordial atmosphere. Table 2 compiles our results for the amount of atmosphere around cores with various masses after being embedded in the disk for 10⁴, 10⁵, and 10⁶ years (compare Figure 6).

Planetary cores less massive than about 0.5 M_⊕ are probably unable to retain gravitationally a dense planetary envelope after the evaporation of the protoplanetary disk (Stökl et al. 2015). More massive cores are able to keep a significant fraction of their nebula-accreted atmosphere, and the final amount of atmosphere depends on the efficiency of atmospheric escape processes. Hydrodynamic atmospheric escape caused by the high soft X-ray and extreme ultraviolet (XUV) radiation of the young host star can remove on the order of about 10²⁵ g of atmosphere gas (Erkaev et al. 2013; Lammer et al. 2014; Johnstone et al. 2015) in the habitable zone of a solar-like star. Nonthermal atmospheric escape processes, such as pickup of planetary hydrogen ions by the stellar wind, also contribute to atmosphere loss but are less important under these conditions (Kislyakova et al. 2013, 2014).

An estimate of the atmosphere that will remain after the phase of high mass loss driven by the young host stars' XUV radiation is given by the numbers in brackets in Table 2. The values were calculated by Johnstone et al. (2015), who combined a static lower atmosphere model with a hydrodynamic upper atmosphere model to derive scaling laws for the mass loss rate as a function of planetary mass, atmospheric mass, and input stellar XUV flux. Using the stellar XUV evolutionary tracks derived by Tu et al. (2015), they calculated the amount of atmosphere lost by hydrogen atmospheres with a range of planetary masses and initial atmospheric masses. The numbers given in Table 2 correspond to the assumption that the star started out its life as an average rotator.

From this estimate of the atmospheric loss, it appears that terrestrial planets with masses larger than about 1 M_⊕—if they form quick enough to get embedded in the disk—can hardly avoid ending up with permanent hydrogen envelopes. These planets would therefore appear enshrouded with dense and extended hydrogen envelopes, and hence they resemble mini-Neptunes more than the terrestrial planets in the solar system.

This result is an important constraint for explaining the formation of large rocky planets (i.e., without dense hydrogen envelopes) in the habitable zone of solar-like stars. Assuming that planets such as Kepler-62e or Kepler-62f (Borucki et al. 2013) and Kepler-452b (Jenkins et al. 2015) with radii considerable larger than 1 R_⊕ are indeed rocky and do not possess an extended envelope (Kaltenegger et al. 2013; Jenkins et al. 2015), one arrives at core masses of several M_⊕. According to our results, such cores—and also their parent planetary embryos, depending on their mass evolution–would be quick in accumulating substantial primordial envelopes if embedded in the disk gas for any length of time. If they do, it is unlikely that they will lose their primordial atmosphere in the later evolution.

Finally, we want to highlight that this result of our study is in agreement with the observation-based statistical analysis of Rogers (2015). Rogers (2015) focuses on the mass range studied in our work and concludes that most 1.6 R_⊕ exoplanets are not rocky. By analyzing the known exoplanets within this radius–mass domain and by including the sample of Kepler transiting sub-Neptunes with the Keck radial velocity follow-up, they showed that the majority of exoplanets with radii of 1.6 R_⊕ and larger show an average density too low to be consistent with a pure iron and silicate composition without a gaseous envelope. More planet mass–radius measurements by satellites like CHEOPS (Broeg et al. 2013) will produce more accurate data that can be used to constrain the transition between rocky planets and planets with extended gaseous envelopes. As demonstrated in Figure 6 and Table 2, cores of about Earth mass can easily capture a large amount of primordial envelope from the nebula that is unlikely to be lost by thermal escape after the dissipation of the disk (Lammer et al. 2014; Johnstone et al. 2015).