1 Introduction

Formation and Structure of Low Density Exo-Neptunes

Abstract

Kepler has found hundreds of Neptune-size (2-6 ) planet candidates within 0.5 AU of their stars. The nature of the vast majority of these planets is not known because their masses have not been measured. Using theoretical models of planet formation, evolution and structure, we explore the range of minimum plausible masses for low-density exo-Neptunes. We focus on highly irradiated planets with  K. We consider two separate formation pathways for low-mass planets with voluminous atmospheres of light gases: core nucleated accretion and outgassing of hydrogen from dissociated ices. We show that Neptune-size planets at  K with masses as small as a few times that of Earth can plausibly be formed core nucleated accretion coupled with subsequent inward migration. We also derive a limiting low-density mass-radius relation for rocky planets with outgassed hydrogen envelopes but no surface water. Rocky planets with outgassed hydrogen envelopes typically have computed radii well below . For both planets with H/He envelopes from core nucleated accretion and planets with outgassed hydrogen envelopes, we employ planet interior models to map the range of planet mass–envelope mass–equilibrium temperature parameter space that is consistent with Neptune-size planet radii. Atmospheric mass loss mediates which corners of this parameter space are populated by actual planets and ultimately governs the minimum plausible mass at a specified transit radius. We find that Kepler’s 2-6 planet candidates at –1000 K could potentially have masses . Although our quantitative results depend on several assumptions, our qualitative finding that warm Neptune-size planets can have masses substantially smaller than those given by interpolating the masses and radii of planets within our Solar System is robust.

Subject headings:
Planetary formation — Extrasolar planets —Accretion

1. Introduction

The first 4.5 months of Kepler data provide evidence for hundreds of “Neptune-size” (2 – 6 , where is the Earth’s radius) planets orbiting within 0.5 AU of their stars (Borucki et al., 2011a, b). The prevalence of planet candidates within this size range raises questions about both planetary growth and migration of Neptune-size planets. Assuming that many of these candidates are true planets, what are they, how did they form, and why are they so numerous?

Kepler measures planetary sizes and orbital periods. In some cases, planet masses can be estimated from dynamical interactions between the planet and its star (Doppler method) or among planets in a multi-planet system (Transit Timing Variations, TTVs, Holman et al., 2010). However, the masses of most of Kepler’s Neptune-size planet candidates will be difficult to measure.

We model herein the growth, physical evolution, and interior structure of Neptune-size planets that possess voluminous atmospheres of light gases. Our focus is obtaining estimates of the minimum plausible masses of Kepler’s planet candidates. The maximum plausible mass of a planet of radius can be estimated from the mass-radius relationship for rocky (Earth-like composition) planets (e.g., Valencia et al., 2006; Seager et al., 2007; Fortney et al., 2007; Marcus et al., 2010). In contrast, estimation of minimum plausible masses requires more detailed modeling of planetary growth, because formation of low-mass planets of solar composition demands complicated and contrived scenarios involving large amounts of mass loss. We consider formation of low-density planets both through core nucleated accretion, and through outgassing of low-molecular weight atmospheres. This work does not consider planet formation via gravitational instability (Boss, 1997) or tidal downsizing (Nayakshin, 2010a, b, 2011).

We present new core nucleated planet accretion calculations following the approach of Pollack et al. (1996) and Movshovitz et al. (2010). Whereas all previous papers with this code emphasize the formation of Uranus mass and larger planets, here we present a new application of the code to a lower mass regime . We push to lower planet masses by modeling formation scenarios where the gas disk dissipates well before rapid gas accretion. We also consider lower solid planetesimal surface densities (4 g cm at 5.2 AU and 6 g cm at 4 AU) than most previous calculations (10 g cm at 5.2 AU) to attain lower heavy element core masses. Until recently, high planetesimal surface densities (about 3 times the minimum mass solar nebula at 5.2 AU) were needed to model Jupiter formation on a reasonable time scale. Advances in the modeling of grain physics (Movshovitz et al., 2010) now allow for a reasonable formation time for Jupiter, even with considered here.

We supplement the detailed core nucleated accretion calculations with equilibrium models of Neptune-size planets having H/He envelopes calculated following the approach of Rogers & Seager (2010a) and Rogers & Seager (2010b). The equilibrium model is less computationally time consuming and allows us to more comprehensively sample the parameter space (heavy element core masses, envelope masses, irradiation levels, and intrinsic planet luminosities) of interest. We focus on low-density planets having equilibrium temperatures of 500–1000 K, since these temperatures are relevant to the planet candidates found in Kepler’s first quarter data (Borucki et al., 2011a).

We also explore outgassing during planet formation as a possible origin pathway for low-density Neptune-size planets, and derive a limiting mass-radius () relation bounding the maximum radius/minimum density for planets with primary de-gassed envelopes. Following Elkins-Tanton & Seager (2008b), we consider outgassing of hydrogen gas produced when water reacts with metallic Fe in accreting materials during planet formation. Our outgassed exoplanet models self-consistently treat the connection between the planets’ interior structure (iron core mass and mantle composition) and the mass of H degassed. We thereby provide the first exoplanet relations that include the effect of an outgassed H gas layer. To derive the limiting relation, we study planets that accreted from a mixture of water and material with chemical composition characteristic of the high-iron enstatite (EH) chondrite meteorite class, corresponding to end member scenarios yielding maximum outgassed H.

We begin by describing our equilibrium model for low-mass planets with gas layers in Section 2. This model is applied in later sections to explore the mass-radius () relationships for low-mass planets with voluminous atmospheres of light gases acquired by core nucleated accretion and outgassing of hydrogen. Section 3 describes the formation and properties of Neptune-size planets that assembled through core accretion, and Section 4 describes the formation and properties of planets that outgassed hydrogen from dissociated ices. We consider mass loss from the envelope in Section 5. We discuss our results and conclusions in Section 6.

2. Models of Planets in Equilibrium: Methods

We use equilibrium models—spherically symmetric planets in hydrostatic equilibrium—for two applications. The first (Section 3) is to explore the mass-radius relationships for low-mass planets formed via core nucleated accretion. The second application is to again study mass-radius relationships for planets that acquired an envelope of light gases through outgassing (Section 4).

Our equilibrium model is based upon the planet interior model from Rogers & Seager (2010a, b). We have, however, included updates to the temperature profile in the radiative regime of the envelope, and to the outer boundary conditions of the planet. We use equilibrium models to study instantaneous states of evolving planets assuming the planets are undergoing quasistatic evolution. Our work does not focus on cases where the envelope dynamics or variations in the interior luminosity profile have an important effect.

We assume spherically symmetric and differentiated planets consisting of up to four layers. From the inside out, these layers are an iron core, a silicate mixture, HO, and a gas envelope. The coupled differential equations describing the mass of a spherical shell in hydrostatic equilibrium:

(1)
(2)

and the differential equation describing the radial optical depth, , in the gas layer

(3)

are integrated inward from the top of the planet’s envelope. Above, is the interior mass coordinate, is the distance from the planet center, is the pressure, is the mass density, is the mean opacity at thermal wavelengths, and is the gravitational constant.

Within each chemical layer, the equation of state (EOS) relates the density to the pressure and temperature . In analogy to the models in Section 3, we define the exterior boundary condition on the planet envelope (, , , ) at radial optical depth . We then determine the corresponding pressure by imposing

(4)

where is the Rosseland mean opacity at the photosphere boundary, and is the gravitational acceleration. While the density varies abruptly between the chemical layers, both and are continuous across layer boundaries. For a specified planet composition, energy budget, and mass, , the planet radius, , is iterated until a self-consistent solution satisfying the inner boundary condition (, ) is achieved.

We assume that the gas envelope is in radiative-convective equilibrium, with an outer radiative zone surrounding a convective layer at greater depths. Within the thin outer edge of the envelope, we adopt the isotropic average temperature profile from Equation (29) in Guillot (2010),

(5)

an analytic solution to the “two-stream” gray equations of radiative transfer for a plane-parallel irradiated atmosphere. The irradiation temperature, , characterizes the short wave energy flux received by the planet from the star and relates through the redistribution factor, , to the equilibrium temperature of the planet in the radiation field of the star ( for full redistribution). The intrinsic temperature parameterizes the total intrinsic luminosity of the planet, ( denotes the Stefan-Boltzmann constant). The total intrinsic planet luminosity, , is the sum total of contributions from envelope contraction, radioactive decay, and secular cooling of the core. The ratio of the short-wave and long-wave optical depths is represented by . We take , which Guillot (2010) found provided a good match to detailed calculations of hot Jupiter atmospheres from Fortney et al. (2008).

In highly irradiated planet atmospheres, the radiative regime of the envelope may extend to depths beyond where the plane parallel approximation (assumed when deriving Equation 5) is valid. In these cases, once all of the incoming stellar radiation is absorbed at optical depths , we transition smoothly to the radiative diffusion equation,

(6)

The onset of convective instabilities determines the depth of the transition to the convective layer of the gas envelope. In the convective regime, we adopt an adiabatic temperature profile.

We use the EOS from Saumon et al. (1995). The effect of uncertainties in the H/He EOS (see e.g., Militzer et al., 2008; Nettelmann et al., 2008) will be small compared to the effect of uncertainties in the heavy element composition and distribution, for the low-mass planets we are considering. While the major uncertainties in the EOS are at Mbar pressures or above, the pressures at the base of our H/He envelopes are typically less than a few tenths of a Mbar. As in the formation and evolution models of Section 3, we use Rosseland mean molecular opacities from Freedman et al. (2008). We neglect grain opacities in our equilibrium models, however, since we are interested in the planet radii at late times, after all the grains have settled.

Under the H/He envelope, the rock-ice interior is modeled with differentiated layers of iron, FeMgSiO silicates, and HO. For these materials, we employ EOS data sets from Seager et al. (2007), which were derived by combining experimental data at with the theoretical Thomas-Fermi-Dirac (TFD) equation of state at high pressures, . The equation of state at intermediate pressures between and 10,000 GPa is not well known, since this pressure range is neither easily accessible by experiments nor by TFD theory. For HO, Seager et al. (2007) used density functional theory calculations to fill in the EOS in this pressure regime, while for all other materials, they bridged the pressure gap by extrapolating the empirical Birch-Murnagham EOS and the theoretical TFD EOS to higher and lower pressures, respectively. Thermal effects are neglected in the Seager et al. (2007) EOSs – at the high pressures found in the interior layers, thermal corrections have only a small effect on the density, . An improvement over our previous models is that we now more consistently take into account the Si/Mg/Fe ratios in the mantle by calculating EOSs for mixtures of MgO, FeO and SiO. Core mass-radius relations calculated following this scheme (but neglecting the small contribution to pressure from the gaseous envelope) were also employed in the planet evolution calculations of Section 3.

A major uncertainty in the validity of the models comes from the assumption that the layers of water and H (or H/He) are not mixed. For the  K planets considered in this work, H and HO are miscible at the pressures and temperatures in the model envelopes. So they could, in principle, be homogeneously mixed. By considering the extreme where the H and HO are fully separated, we set an upper bound on the planet radius; typically if hydrogen is mixed into the water layer one expects the planet’s radius to be smaller (e.g., Nettelmann et al., 2010). Although our aim is to model H/He envelope planets, some of our models do have significant water content, and future work should include the miscibility of H and HO.

The planet radii in both our equilibrium and evolution models underestimate the planet radii measured during transit in a predictable way. This “transit radius effect” (Baraffe et al., 2003; Burrows et al., 2003) is a consequence of our exterior boundary condition (Equation 4), which pegs our model planet radii, at a radial optical depth . In contrast, it is the transverse optical depth for transmission through the planet limb that determines the transit radius. Hansen (2008) derived a correction for the transit radius effect,

(7)

where the transit radius is defined at a transverse optical depth of unity, and represents the atmospheric scale height at the planet limb. Equation (7) applies when , and assumes the outer limb of the planet atmosphere is well described by an ideal gas. For the low-mass () planets with hydrogen-rich envelopes we are considering, the transit radius correction is typically between 1% and 10%. Equation 7 assumes a clear cloud-free planet atmosphere; high level clouds or hazes could further enhance the transit radius effect.

3. Planet Formation by Core-Nucleated Accretion

3.1. Methods

Models for the formation and evolution of a planet consisting of a core of heavy elements and a gaseous envelope of solar composition are calculated according to the procedures described by Pollack et al. (1996) and Movshovitz et al. (2010). The formation calculation consists of three major parts: () the accretion rate of planetesimals onto the planet; () the interaction of the infalling planetesimals with the gaseous envelope; and () the evolution of the gaseous envelope and the determination of the gas accretion rate.

The planetesimal accretion rate onto the planetary embryo is based on the equation originated by Safronov (1969). If is the mass of the embryo, then the fundamental equation for its growth is

(8)

where is the effective geometrical capture cross-section, is the surface density of solid material (planetesimals), is the orbital frequency, and the value of the gravitational enhancement factor, , is obtained from Greenzweig & Lissauer (1992), assuming a planetesimal radius of 100 km. If no gaseous envelope is present, then , the heavy element core radius. As in our previous publications, we take the feeding zone from which the embryo can accrete planetesimals to extend 4 Hill sphere radii on either side of the orbit (Kary & Lissauer, 1994), and assume that the solid surface density is constant within that zone. The value of in the feeding zone is adjusted at each time step to take into account depletion of planetesimals by accretion onto the embryo and expansion of the feeding zone into undepleted regions, as the embryo’s mass increases.

The second element of the code calculates the interactions between planetesimals and the gaseous envelope of the protoplanet as they fall through it (Podolak et al., 1988). The details of how this calculation is performed are described in Pollack et al. (1996), Hubickyj et al. (2005), and Movshovitz et al. (2010). These calculations provide the effective capture radius to be used in Eq. (8), as well as the deposition of mass and energy as a function of radius in the gaseous envelope. The effective capture radius can be several times larger than the actual solid heavy element core radius because of the effects of the gas on slowing down, ablating, and fragmenting the planetesimals. It is assumed that the material from the planetesimals that is deposited in the envelope later sinks to the heavy element core, releasing gravitational energy in the process (Pollack et al., 1996). This assumption is not entirely accurate: Iaroslavitz & Podolak (2007) show that some of the heavy-element material should actually dissolve in the envelope. Thus the ‘core masses’ that we calculate actually represent the total heavy element abundance in the planet in excess of the solar metal abundance; most of these heavy elements (including all of the rock and organic compounds) would be expected to reside in the actual core of the planet.

The third element of the simulation is the solution of the four differential equations of stellar structure for the gaseous envelope, with energy sources from accreting planetesimals, from gravitational contraction, and from cooling. The adiabatic temperature gradient is assumed in convection zones. At the heavy element core boundary, the luminosity is set to the energy deposition rate for the planetesimals that hit the heavy element core. Outside the heavy element core, the energy supplied by ablated and fragmented planetesimals is included as a source term in the energy equation.

At the core-envelope interface, the radius is set to that of the outer edge of the heavy element core. The heavy element core mass-radius relation is calculated using the equilibrium model described in Section 2, assuming 10% Fe, 23% FeMgSiO, and 67% HO, by mass. The heavy element core composition we adopt is motivated by comet compositions, and represents rock with a Fe/Si ratio near solar mixed with ice in a ratio of 1:2 by mass.

At the surface, gaseous mass of solar composition is added at a sufficient rate to maintain an outer radius , which is given by (Bodenheimer et al., 2000) as

(9)

where is the Hill sphere radius, is the sound speed in the disk outside the planet, is the total planet mass, and, nominally, . Note that when is large compared with the Bondi accretion radius, , the expression reduces to the Bondi radius, while in the case of the opposite limit, . In developments after the above expression was formulated, it turned out that had to be modified. Three-dimensional calculations of disk-planet interaction (Lissauer et al., 2009) gave the result that not all the gas passing through the Hill sphere is actually accreted by the planet; some of it simply flows through and rejoins the disk’s azimuthal motion. The 3D simulations provided an estimate of the effective planetary radius, which corresponds to , the value used in this paper.

The density and temperature at the planet’s surface are set to assumed nebular values , , respectively. The value of is constant in time, while decreases linearly to zero with time, over a time scale Myr. In a variation of this boundary condition, is constant in time up to a time comparable to , then it is linearly reduced to zero on a time scale of yr. These assumptions roughly characterize the dissipation of the gaseous disk. is held constant while the planet is accreting; our model incorporates migration only through temperature increases subsequent to the conclusion of the planet’s growth. Modeling simultaneous migration and accretion is beyond the scope of this work.

When approaches zero, the accretion of gas is halted and the evolution is calculated at constant mass over time scales up to 3–4 Gyr. The envelope mass at cutoff in these simulations is always small enough that rapid runaway gas accretion does not occur, and Eq. (9) is always valid for the determination of the gas accretion rate. The accretion rate required to keep remains much lower than the limit imposed by disk physics in supplying material to the Hill sphere of the planet (Lissauer et al., 2009). Once gas accretion is shut off, rapidly falls below , and the planet becomes isolated from the disk. The surface boundary condition changes to that of a hydrostatic atmosphere that radiates from the photosphere:

(10a)
(10b)

where is the Stefan-Boltzmann constant, is the surface temperature, is the total luminosity (energy radiated per second) of the planet, and , , and are, respectively, the Rosseland mean opacity, the pressure, and the acceleration of gravity at the photosphere. There are two contributions to : that from the internal luminosity provided by the planet, and that from the energy absorbed from the central star and re-radiated by the planet. Thus:

(11)

where is the internal contribution (generally small), and is the equilibrium temperature of the planet in the radiation field of the star. The former quantity is determined from the evolutionary calculation, while the latter is a parameter that depends on the assumed distance of the planet from the star and the stellar luminosity.

The equation of state of the gas is taken from Saumon et al. (1995), interpolated to our assumed composition of hydrogen mass fraction X = 0.74, helium mass fraction Y = 0.243, and metal mass fraction . Although the equation of state in the outer, low-density layers is essentially that of an ideal gas, the inner regions near the heavy element core can be significantly non-ideal once the envelope has become sufficiently compressed.

The Rosseland mean opacity calculation has three components. At temperatures above 3000 K, the molecular/atomic opacities of Alexander & Ferguson (1994) are used. In practice the details of the opacities in this region are unimportant because the energy transport is almost always by convection. In the temperature range 100–3000 K the molecular opacities, without grains, of Freedman et al. (2008) are used. Grain opacities are then added in the temperature range 100–1800 K. Two sources of grains are taken into account, first, those provided by the ablating planetesimals as they interact with the envelope, and, second, those that accrete along with the gas at the surface of the planet. At each time step of the evolutionary calculation, and at each depth in the envelope, the grain size distribution is recalculated, taking into account the coagulation and settling of grains. The size distribution is represented by 34 bins, covering the size range 1.26 m to 2.58 mm. The effective cross-sections for absorption and scattering are calculated as a function of grain size and frequency; then an integration over grain size and frequency gives the Rosseland mean opacity as a function of depth. The details of the grain physics are given in Movshovitz & Podolak (2008) and Movshovitz et al. (2010). The grains are composed purely of silicates, with a dust-to-gas ratio of about 0.01 by mass; little error results from this assumption compared to the uncertainties in grain shape, sticking probability, and radiative properties (Movshovitz et al., 2010). Grains are assumed to be completely evaporated above 1800 K. The grains are important during the gas accretion phase. Once accretion is shut off, the grains rapidly settle towards the center and are evaporated. This effect is included in the calculations and indicates that any grains remaining in the atmosphere have a negligible effect upon the evolution. Thus in the final constant-mass evolution phase, the molecular opacities completely dominate.

Figure 1.— Mass of the protoplanet as a function of time for Runs I. For Run Ia (black curves) the solid line denotes the mass of the heavy element core, the dotted line the mass of the H/He envelope, and the short-dash-dot line the total mass. For Run Ib, the same line types are shown in red.
Figure 2.— Mass of the protoplanet as a function of time for Runs II. For Run IIa (black curves) the solid line denotes the mass of the heavy element core, the dotted line the mass of the H/He envelope, and the short-dash-dot line the total mass. For Run IIb, the same line types are shown in red.
Run (AU) (g/cm) (g/cm) (K) (Myr) ()
Ia 5.2 4 280 115 3.5
Ib 5.2 4 280 115 2.5
IIa 4.0 6 420 125 2.0
IIb 4.0 6 420 125 0.9
Table 1Input Parameters for Evolutionary Runs
Run () () () () () () ()
Ia 8.13 4.08 4.05 9.8 8.1 14.8
Ib 5.20 3.52 1.68 8.0 6.6 15.7
IIa 3.19 2.65 0.54 6.0 5.0 17.9 11.7
IIb 2.66 2.50 0.16 3.6 3.3 6.7 6.2

The first subscript on the radius gives the evolutionary time in Gyr. The second subscript gives the assumed equilibrium temperature of the planet.

Table 2Results from Evolutionary Runs: Masses and Radii

3.2. Evolution Input Parameters and Results

The planet initially consists of a heavy element core of 1 and a light element envelope of about . The protoplanet is located at either 5.2 AU or 4.0 AU in a protoplanetary disk, with the solid surface density g cm at 5.2 AU and 6 g cm at 4 AU. The initial evolutionary time is set to yr and yr, respectively, for and 6 g cm, approximately the time needed to assemble a heavy element core of mass .

The quantity is set to 115 K at 5.2 AU and 125 K at 4 AU. Then , where is the surface density of the gas component. As mentioned above, in general declines with time. The scale height of the gas in the disk , where is the orbital distance from the star. Once started, the evolution consists of three main phases. The first involves primarily accretion of solids onto the heavy element core, with a relatively low-mass envelope and a low gas accretion rate. The solids accretion rate slows down significantly near the point where the isolation mass () for the core is reached; for g cm at 5.2 AU this mass is about 2.9 and for g cm at 4 AU, about 2.4 . During the second phase, the gas accretion rate is about 3 times as high as the core accretion rate, , and both are nearly constant in time (Pollack et al., 1996). The envelope mass builds up relative to the heavy element core mass, which grows slowly. The phase of rapid gas accretion, which for giant planets begins once the envelope mass becomes about equal to , does not occur in these calculations. Instead, gas accretion is cut off and the planet evolves through a third phase at constant mass with boundary conditions provided by Eq. (10). During the early part of this phase, the planet is assumed to migrate to a position within 1 AU from the star. Representative cases with K and 1000 K are presented.

The input parameters of the four runs are shown in Table 1, which includes , the gas surface density , the surface boundary temperature , and the isolation mass .

Figure 3.— The protoplanet’s total luminosity, including internal and irradiation contributions, as a function of time during the formation phase and the contraction/cooling phase for Run Ia (solid curve), and Run IIa (dashed curve). The equilibrium temperature is increased to 500 K, after the formation phase, during these runs.

The results of our calculations for Runs Ia, Ib, IIa, and IIb are shown in Figures 1, 2, and 3. The masses and radii that are derived for the four runs are listed in Table 2.

Run Ia is based on a disk with g cm at 5.2 AU. This value is only slightly greater than that of the minimum mass solar nebula. But note that our calculation of (Eq. 8) neglects transport of solids into or away from the planet’s accretion zone. Moreover, our planetesimals are all assumed to have the same radius, 100 km. In fact there must be a range of planetesimal sizes, and the effective planetesimal size is not well known. Smaller planetesimals would result in more rapid accretion (see footnote 3 of Lissauer et al., 2009). The accretion rate that is actually calculated may thus correspond to a value of slightly different from 4 g cm.

The details of the calculation with the parameters of Run Ia are presented in Movshovitz et al. (2010), their Run 4. In that paper the run is continued well into the phase of rapid gas accretion, and is terminated with and . The formation time for a giant planet is found to be 4 Myr. In the present run, the accretion of gas and solids is cut off at a time of 3.5 Myr, consistent with estimated lifetimes of protoplanetary disks (Hillenbrand, 2008). At that time the value of is assumed to decrease to zero on a time scale of 10 yr. The calculation is then continued up to Gyr times with constant values of and . At the beginning of this phase the equilibrium temperature is gradually increased, on a time scale of 4 Myr, to an assumed final value of 500 K. A gradual increase in to 1000 K was accomplished in a total time of  yr. The final values of for these two temperatures and for times of 1 and 4 Gyr are given in Table 2; they are close to Jupiter’s radius R, even though the planet’s mass is only 8.13 .

Run Ib also is based on the run 4 from Movshovitz et al. (2010). In this case the accretion of gas and solids was cut off at 2.5 Myr, at which point and . The evolution was again continued into the phase of cooling and contraction at constant mass, with assumed values of of 500 and 1000 K. In the case with  K, the final radii are again comparable to or larger than R. In the case with K the minimum radius is 6.6 , only slightly smaller than the corresponding value in Run Ia.

Run IIa is an entirely new calculation, with the planet forming at 4 AU in a disk with g cm. During the initial phase of rapid core accretion, the luminosity reaches a maximum of L at a time of yr. The heavy element core mass is 2.2 and the core accretion rate yr at this time. Later, at 1 Myr, has decreased to yr, and has increased to yr. The luminosity has decreased to L. Because of computational time limitations, and to obtain a lower envelope mass than that found for Run Ib, the accretion in this run was cut off at 2 Myr, with and . If the evolution had been continued up to 2.5 Myr, the heavy element core mass would have been practically unchanged, and the would have increased by about 0.25 . At the end of the contraction/cooling phase the radii are in the range 5–6 for the case of  K, and for  K they are larger then , close to the values obtained in Runs I for that temperature.

To investigate the effect of an even smaller value of , Run IIb was calculated with the same parameters as Run IIa, but with an arbitrary accretion cutoff at yr. At that point, and . Final radii turned out to be in the range 3–4 for K and in the range 6–7 for K. The significant reduction in envelope mass resulted in final radii that are about half the values obtained for Run IIa.

We neglect heating from radioactive decay in the core nucleated accretion calculations. Including this additional energy source would delay envelope contraction and planet cooling. Consequently, the planet radii at 1 Gyr and 4 Gyr in Table 2 may be systematically underestimated by a small amount. We estimate that, for the cases in Table 2, the planet luminosity from radioactive decay would be roughly one order of magnitude smaller than luminosity from envelope contraction, assuming bulk Earth abundances of K, U, and Th in the heavy element cores (Van Schmus, 1995). The fractional contribution to the planet energy budget from radioactive heating will be higher for older planets (4 Gyr) and cases where the heavy element core contributes a larger fraction of the planet mass (Run II).

3.3. Equilibrium Model Results

In this section we explore planet radii over a wide range of heavy element core masses, envelope masses, irradiation levels, and intrinsic planet luminosities. The planet formation and evolution model described in the Section 3.1 is computationally intensive. Since it is not feasible to simulate planets under all conditions of interest following that approach alone, we enlist an equilibrium planet structure model (Section 2) to cover a wider range of parameter space.

Our equilibrium model shows good agreement with the planet evolution models in Section 3.2 despite the differences in their treatment of the outer radiative regime, the intrinsic planet luminosity, and the effects of stellar insolation. For each entry in Table 2, we applied the equilibrium model to simulate the same combination of , , and . The radii at  K in the two models agree to better than in every case. The planet radii at  K are more sensitive to model assumptions and exhibit larger discrepancies (up to 14%, with the equilibrium model radii systematically below those in Table 2).

We explored the parameter space of , , , and with our equilibrium model. Figures 4 and 5 present a selection of mass-radius () curves at (a)  K and (b)  K. Figure 4 displays the effect on the radius of varying the envelope mass fraction, while Figure 5 shows the effect of varying the planet’s intrinsic luminosity, . The thick solid line is common between Figure 4 and 5, representing and . corresponds to both the evolution model (Run Ia) at 4 Gyr, and the evolution model (Run IIa) at 1 Gyr (independent of ).

Figure 4.— Equilibrium mass-radius relations for various choices of envelope mass fraction, . All data in this plot have , and (a)  K or (b)  K. Each curve corresponds to a different value of : 0.001 (thin solid), 0.01 (thin dashed), 0.05 (thin dot-dashed), 0.1 (thin dotted), 0.2 (thick solid), 0.3 (thick dashed), 0.4 (thick dot-dashed), and 0.5 (thick dotted). Black lines denote our model radii (defined at a radial optical depth ), while the corresponding blue lines represent radii corrected for the transit radius effect. The thick red line is the mass-radius relation for icy heavy element cores having no envelope (). Red triangles present the subset of Table 2 evolutionary run results that have : Run Ia () at 4 Gyr, and Run IIa () at 1 Gyr. The green curves show the effective planet Roche-lobe radius for four different choices of host-star properties representative of spectral classes M5 V, M0 5V, K0 V, and G2 V (in order of increasing Roche lobe radii). The K0 V and G2 V Roche-lobe radii are beyond the scale of the  K plot.

Figure 5.— Equilibrium mass-radius relations for various choices of intrinsic planet luminosity . All data in this plot have , and (a)  K or (b)  K. Each curve corresponds to a different value of : (thin solid), (thin dashed), (thin dot-dashed), (thin dotted), (thick solid), (thick dashed), and (thick dot-dashed). Black lines denote our model radii (defined at a radial optical depth ), while the corresponding blue lines represent radii corrected for the transit radius effect. The green curves show the effective planet Roche-lobe radius for four different choices of host-star properties representative of spectral classes M5 V, M0 5V, K0 V, and G2 V (in order of increasing Roche lobe radii). The K0 V and G2 V Roche-lobe radii are beyond the scale of the  K plot.

The curves for low-mass planets with voluminous gas layers show several notable features. First, the planet radii (at constant envelope mass fraction, , and ) increase dramatically toward low planet masses. This is due to the low surface gravities, and thus large atmospheric scale heights found at low masses. Second, the radius of planets having identical envelope mass fractions, , are remarkably insensitive to the planet mass when . At these masses, increased compression of the envelope offsets the effect of increasing the planet mass. Third, for planets of identical total mass (within the mass range plotted) the planet radius increases monotonically with the envelope mass fraction. Fourth, and both have a stronger effect on the radius of low mass planets compared to their more massive counterparts. This is understandable, because, given the same envelope mass fraction, in lower mass planets the envelope accounts for a larger fraction of the planet radius.

Planet radii between 2 and are of special interest, because Kepler is finding a large number of planet candidates within this size-range (Borucki et al., 2011a, b). We plot in Figure 6 combinations of and that yield planet radii within this range. Planets at can contain at most 0.08% of their mass in H/He at  K, and at most 0.0015% at  K. Larger planets can support more massive envelopes. A planet at  K requires an envelope accounting for at least a few percent of the planet mass. At  K and , between 0.1% and 23% H/He by mass is possible, depending on the planet mass and intrinsic luminosity.

Figure 6.— Planet mass and envelope mass that are consistent with a particular planet radius, for planets comprised of ice-rock interiors surrounded by H and He in protosolar proportions. These models represent planets that formed beyond the snow line by core nucleated accretion. We plot the envelope mass fraction as a function of total planet mass for planets with radii (a) , (b) , and (c) . Black curves represent planets at  K, while red curves correspond to  K. The line style indicates the planet luminosity: (dashed), (dot-dash), and (solid). The thin dotted lines are contours of constant envelope mass loss timescale, . Each contour is labeled with for , and can easily be scaled for other choices of using Equation (13).

It is important to note that the relations in Figures 4, 5, and 6 are not isochrons, but correspond instead to constant total intrinsic luminosity per unit mass, . The total intrinsic luminosity, , is the sum total of heating from radioactive decay, cooling of the planet core, and contraction of the planet envelope. In the evolution calculations from Table 2, the planet luminosity contribution from envelope contraction alone ranges from to at 1 Gyr and from to at 4 Gyr. Some of the low- curves in Figures 5 and 6 do not extend to higher masses because they encounter unphysically low planet interior entropies. Although is a proxy for the age of the planet, the relationship between and planet age depends on the planet’s mass, composition, abundance of radioactive isotopes, insolation history and dynamical history. Since our equilibrium models are presented at a specified , we have side-stepped the issue of relating to planet age and present the model radii in a way such that they can be applied to many different evolution scenarios. Our aim with the equilibrium models is to broadly explore parameter space; it is beyond the scope of this work to relate and age directly by simulating all possible planet evolution histories.

Simulated planet radii for planets at  K may be in error by up to 20%. The problem is in extrapolating the opacity tables at the high pressure end. This in turn makes the radiative-convective boundary uncertain (a deeper radiative-convective boundary makes for a smaller planet). Planets at  K are less affected by this opacity-caused radius problem ( radius uncertainty for ). This issue affects both our equilibrium and evolution models.

4. Planet Formation by Outgassing of Hydrogen

4.1. Model

Outgassing provides a mechanism for low-mass terrestrial planets to acquire an atmosphere even if they fail to accrete H and He from the protoplanetary nebula. In this section we explore the optimum conditions for a planet to acquire a voluminous gas envelope through outgassing. We base our model approach on Elkins-Tanton & Seager (2008b, a), with the improvements of a more detailed interior structure model and a calculation of the planet radius.

We focus on outgassing of H produced when water reacts with metallic Fe in accreting materials during planet formation (Ringwood, 1979; Wanke & Dreibus, 1994; Elkins-Tanton & Seager, 2008b). Hydrogen gas has the potential to yield the most voluminous outgassed atmospheres, being both of low-molecular weight and (for some planetesimal compositions) degassed in substantial quantities. Although we do not consider these processes in detail here, in general, outgassing may also proceed during accretion as impinging planetesimals are heated and vaporized upon impact; during magma ocean solidification as volatiles are partitioned between the atmosphere and melt; and during volcanic/tectonic activity after the planet has formed.

The reaction between water and metals during planetary accretion and differentiation intrinsically links the planet’s interior structure to its initial atmosphere’s mass and composition. Metallic iron forming the planet will either differentiate to contribute to the planet iron core, or become oxidized and incorporate into the planet mantle. Given an initial composition for the primordial material forming a planet, there are two extremes to the eventual planet outcomes. If none of the available water and metals in the accreting materials react (reducing conditions), the planet will have a maximally massive metallic core, relatively iron-poor mantle, minimal outgassed H, and maximal leftover HO. In contrast, if the water and metals react to the maximal extent possible (oxidizing conditions), the planet will have a minimal iron core mass, iron-rich mantle, maximal outgassed H, and minimal leftover HO. When Fe is the limiting reagent, this extreme will correspond to a coreless planet (Elkins-Tanton & Seager, 2008a).

To bound the radii of outgassed rocky planets, we consider the end-member case of a planet formed purely from high iron enstatite (EH) primordial material. The motivation for this choice is three-fold. First, out of all meteoritic compositions, EH material has the potential to degas the most H per unit mass (up to 3.6%, Elkins-Tanton & Seager, 2008b). Second, the oxygen isotope mixing model (Lodders, 2000) predicts that the Earth accreted from material that was 70% EH chondritic matter by mass. Third, heating of EH material releases a low mean molecular weight atmosphere; Schaefer & Fegley (2010) calculated 44% H, 31% CO, 17% HO, 5% CO, and 3% other molecules by volume under their nominal conditions (1500 K, 100 bars). Thus, complete oxidation of an EH planet should achieve effectively the maximum radius plausible for planets with outgassed atmospheres.

For the EH material we adopt the chemical composition of meteorite ALHA77295 from Jarosewich (1990). We distill the mineralogy in our model to include only the most plentiful and important constituents: metallic Fe, FeS, FeO, FeO, MgO, SiO, HO, and H. Following an approach similar to Sotin et al. (2007), less abundant elements are represented by their most similar neighbors in the periodic table: metallic Ni is added to metallic Fe, Ca is added to Mg, and Al is divided equally (by number) between Si and Mg to preserve charge conservation. Other trace constituents (TiO, CrO, MnO, NaO, KO, PO, Co, which combined account for less than 2.2% by mass) are neglected. The resulting simplified composition adopted for the primordial rocky EH planetesimals consists of (by mass) 38.2% SiO, 25.2% metallic Fe, 14.3% FeS, 20.6% MgO, 1.7% HO. Note that HO included in the EH material is adsorbed to the surface or chemically bound to the minerals.

We consider planets initially formed from a mixture of EH material and HO ice. The HO ice is in addition to the 1.7% HO by mass included in the EH minerals. We compute the planet bulk composition after outgassing from stoichiometry (Table 3), assuming some fraction of the accreted iron reacted with water (Fe + H FeO + H) before sinking to form the planet’s metallic core. We note that although we consider only Fe in our reduced EH chemical composition, Ni can also form oxides and be incorporated in silicates. Nickel accounts for 8% of the generalized metallic Fe in our distilled EH chemical composition – the Ni abundance is abundance in ALHA77295 is 1.83% by mass. We do not vary the S mass fraction of the iron core in our models, effectively assuming metallic Fe and FeS oxidize in equal proportions. We do not follow any S released in the conversion of FeS to FeO.

Our interior models of outgassed planets comprise up to four chemically distinct layers: an Fe/FeS core, silicate mantle, water layer, and hydrogen atmosphere. The bulk chemical composition of the planet after outgassing determines the relative masses of the planet layers and the composition of the silicate mantle. All of the degassed H is included in a gas layer surrounding the planet. We place all of the FeS and metallic iron in the planet core. We model the HO in a differentiated water layer surrounding the mantle, although in practice some water may be sequestered into the silicates (e.g., Elkins-Tanton, 2008, and references therein). All of the remaining species (SiO, MgO, FeO, FeO) make up the mantle. The ratio of MgO/FeO sets the Mg # of the silicates (Mg # = Mg/(Mg+Fe) by number). We adjust the mantle equation of state to reflect the relative abundances of SiO, MgO, FeO and FeO, modeling the silicates as a mixture of (Mg,Fe)O magnesiowustite (Elkins-Tanton, 2008), FeO hematite (Wilburn & Bassett, 1978), and stishovite SiO (Andrault et al., 1998). Outgassed bulk compositions and the corresponding planet properties are reported in Table 3.

% Fe oxidized Core wt % Silicate wt % HO wt % excess H wt % Silicate Composition
MgO wt % FeO wt % FeO wt % SiO wt % Mg #
0.0 39.5 58.8 1.7 0.0 35.1 0.0 0.0 64.9 1.00
15.2 33.8 66.0 0.0 0.2 31.5 10.3 0.0 58.3 0.85
50.0 19.5 79.9 -3.7 0.6 25.5 27.3 0.0 47.2 0.62
100.0 0.0 98.8 -8.6 1.2 20.0 42.9 0.0 37.1 0.45
100.0 0.0 98.3 -13.0 1.7 19.1 0.0 45.5 35.4 0.45
Table 3Bulk compositions of EH-composition planets with outgassed H envelopes. The first column represents the fraction of accreted iron that is oxidized and incorporated in the planet’s mantle. The next four columns give the composition of the EH planet after outgassing, assuming all outgassed H is retained. Negative entries in the HO column indicate a water deficit and represent the proportion of additional water (beyond what is included in the EH material) that needs to be accreted in order to oxidize the specified fraction of iron. The last five columns represent the chemical make-up of the silicate mantle, and determine the mantle equation of state. In rows 2–4 we neglect FeO and assume only FeO is produced when iron is oxidized. Row 5 lists the extreme end-member case where all iron is oxidized to FeO.

4.2. Results

We find that planets accreted from solid bodies that were abundant in our solar nebula can degas at most 1.7% of their mass in H. This limit obtains for a fully-degassed coreless EH composition planet that accreted just enough additional water (13.0% by mass) to fully oxidize all available iron to FeO. EH material alone does not contain sufficient HO on its own to oxidize all the metallic Fe within its bulk (only up to 15.2% of the Fe). The accreted material must include an additional 8.6% HO by mass in order to convert all the metallic Fe into FeO, or an additional 13.0% HO by mass to convert all the metallic Fe into FeO. With any more water than this, the metallic Fe becomes the limiting reagent. The maximal outgassed H atmosphere that we derive here is slightly lower than the value 3.6 wt % H found by Elkins-Tanton & Seager (2008b). Differences in the representative EH chemical compositions assumed account for this disparity.

Mass-radius relations for planets harboring H-envelopes from outgassing are shown in Figure 7 at both  K (Figure 7a) and  K (Figure 7b). The blue dot-dashed curve provides an upper limit on the radius of planets accreted from primordial chondritic material alone (without additional water ice), corresponding to the extreme where the oxidizing reaction proceeds until all of the HO bound to the minerals is expended and 0.2% of the planet mass is released in H. After accreting enough additional water (13% by mass) to convert all available Fe to FeO, the magenta solid line represents planets having the maximal fraction of their mass (1.7%) in a degassed H envelope. This curve may be taken to bound the maximum radius/minimum density relation for planets with de-gassed H envelopes, but no free HO.

Planets that accreted more than 13.0% by mass water with the EH chondrite material would have water left over even if all the metals in the planet iron core were expended in the outgassing reaction. In Figure 7 we show relations of an example with initially 20% by mass water ice in the primordial composition (dotted curves). The fully degassed planets with excess water have, in fact, a lower average density compared to the planets with the highest mass fraction of de-gassed H – the effect of the lower density ice-rock interior offsets the decreased proportion of H. In Figure 7, we model the HO layer as a distinct chemical layer below the outgassed H envelope, but mixing of HO and H is another possibility. If HO and H are mixed in the envelope, the planet radii would be smaller than the model radii in Figure 7 due to the decreased atmospheric scale height compared to the differentiated case.

The radii of the outgassed planets depend on the intrinsic luminosity of the planet. In Figure 7, we show mass-radius relations for planets with . Increasing (decreasing) the planet’s intrinsic luminosity by a factor of 10 affects the planet radii in Figure 7 by at most +16% (–9.5%) at and +4.5% (–3.2%) at . Small planet masses and high H contents both increase the radius dependence on .

We explore in Figure 8 the mass of H required by EH composition planets to reach radii of 2 to 3 . Figure 8 is the outgassing analog to Figure 6 for core nucleated accretion. In Figure 8, we restrict our attention to planets without significant amounts of HO on their surface or in their envelopes. The envelope mass fractions, , at a specified radius are not strongly sensitive to the distribution of Fe within the planet interior (i.e., whether the Fe is differentiated in the metallic core, or included in the mantle as oxides) – we show the case where all Fe is oxidized to FeO. Upper bounds on the H wt % for several of the limiting cases in Table 3 are indicated by colored horizontal lines.

Figure 7.— Mass-radius relations for exoplanets with outgassed H envelopes. The planets are assumed to have formed purely from a combination of EH chondrite material and water ice. Accreting material with 20% water ice by mass (dotted lines), 13% water ice by mass (solid lines) 8.6% water ice by mass (dashed lines), and no additional water ice (dot-dashed lines) are considered. The line color indicates the fraction of accreted iron that reacted with water. Black corresponds to planets with no outgassed H and a maximally massive iron core (0% Fe reacted). Blue corresponds to planets where 15.2% of the Fe reacted - the maximum amount possible for pure EH material without added water. Green represents an end-member case wherein all the metallic Fe that accreted to the planet is converted to FeO. Finally, magenta lines correspond to planets that outgassed the maximum possible H for their initial chemical makeup – 100% of their accreted iron is oxidized to FeO. Both the green and magenta relations represent core-less planets, but they differ in the oxidation state of iron inside the planet (FeO versus FeO) and in the overall proportion of H released. Planet equilibrium temperatures of a)  K and b)  K are shown. A fiducial intrinsic luminosity is assumed in all cases. These curves do not include atmospheric escape of H.
Figure 8.— Planet mass and outgassed H envelope mass that are consistent with a particular planet radius, for EH-composition planets without HO on their surface or in their envelopes. We plot the envelope mass fraction as a function of the total mass of the planet for planets with radii (a) , (b) , and (c) . Horizontal lines indicate the maximal H wt % degassed in three limiting cases: if all HO adsorbed in the EH material reacts with metals (0.2%, blue), if all Fe in the EH material is converted to FeO (1.2%, green), and if all Fe in the EH material is converted to FeO (1.7%, magenta). This figure is the outgassing analog to Figure 6 for core nucleated accretion, and all the red and black lines follow the same naming conventions. Black curves represent planets at  K, while red curves correspond to  K. The line style indicates the planet luminosity: (dashed), (dot-dash), and (solid). The thin dotted lines are contours of constant envelope mass loss timescale, . Each contour is labeled with given .

Our main conclusion from this section is that planets of mass with outgassed H atmospheres typically have radii less than (Figures 7 and 8 ). Larger radii are found at the low-mass extreme of the relations in Figure 7, but correspond to planets with very tenuous, loosely bound, envelopes. Outgassing of H from planets accreted from rocky material alone most likely cannot account for the Kepler planet candidates with radii between 3 and 6 .

5. Mass Loss from Low-Density Envelopes

A major question is whether the high , light element, low gravitational binding energy envelopes modeled above are stable and could be retained over gigayear timescales. It is precisely in the low-mass, low-molecular weight, high regime we are considering in which planets are expected to be most susceptible to mass loss. Below we consider, in turn, the importance of Roche-lobe overflow, and XUV-driven atmospheric escape.

Roche-lobe overflow can limit the radii of low-density planets at close orbital separations from their host stars. Our planet interior model assumes spherical symmetry and neglects tidal forces, but this approximation starts to break down for planets near their star. The effective radius of a planet’s Roche lobe is approximated by

(12)

where (Eggleton, 1983). The Roche-lobe radius sets a firm upper limit on the planet radius; any material outside the planet’s Roche lobe is not gravitationally bound to the planet and can escape. We plot planet Roche-lobe radii in Figures 4 and  5 for a sampling of representative host star properties: G2 (, ), K0 (, ), M0 (, ), M5 (, ) (Carroll & Ostlie, 2007). In computing the Roche-lobe radii, we have assumed a planetary albedo when relating to the semi-major axis, ; taking reflection into account with will result in smaller semi-major axes and smaller . Roche-lobe overflow is not an issue for  K planets surrounding a solar analog star. In contrast, when orbiting an M star many of our low-density low-mass planets do fill their Roche lobes. Our equilibrium planet models are not a priori pegged to a given star spectral type. Tidal effects and the Roche-lobe radius set a lower bound on for which the low tail of our equilibrium models are applicable.

XUV-driven mass loss is expected to be very important for low mass, low density planets. This results from the combined effect of large cross-sections to stellar irradiation, low surface gravities, and low envelope binding energies. Predictions for the exoplanet mass loss rates suffer from unknowns in the stellar XUV fluxes, the conditions at the planet exosphere, and the mass loss efficiency. We consider energy-limited mass loss, (e.g., Lammer et al., 2003; Lecavelier Des Etangs, 2007; Valencia et al., 2010),

(13)

The efficiency represents the fraction of the energy in XUV photons incident on the planet that goes into unbinding particles in the planet atmosphere; we take , but can easily be rescaled to another choice of . represents the flux of photoionizing radiation impinging on the planet. is a correction factor that accounts for tidal effects in the Roche potential of planets in close proximity to their star (given by equation 17 in Erkaev et al., 2007). Finally, reflects the planet radius at which XUV photons are absorbed. We estimate following order-of-magnitude arguments gleaned from Section 2 of Murray-Clay et al. (2009),

(14)

where is roughly the column of neutral hydrogen needed to reach , is the pressure at , and is the pressure scale height at (where for visible light).

We take an illustrative example of planets orbiting a solar analog star to explore the order of magnitude of mass loss rates. Figure 9 shows estimated mass loss rates for the planet models presented in Figure 4. For our assumed solar-twin host star, we compute for by scaling the integrated solar XUV flux measured by Ribas et al. (2005) ( at 1 AU). We find that, for and , planets at the low mass extreme of our relations have implausibly short envelope mass-loss timescales  Gyr.

We use energy limited mass loss (Equation 13) to include contours of constant in Figures 6 and 8. The contour values represent corresponding to , but can easily be scaled to reflect other parameter choices:

(15)
(16)

At a specified , the contours are independent of the host star mass so long as tidal effects can be neglected . For the combinations sampled in Figure 6 and 8, this approximation holds for main sequence host stars that are K0 V or earlier, but breaks down for M stars. We emphasize that gives an instantaneous measure of the time that the planet would take to lose its envelope at the calculated current mass loss rate. is expected to vary over a planet’s lifetime. Stars that are more active (e.g., younger) than our Sun would have higher photoionizing fluxes.

We find that planets at the low mass extremes of Figures 6 and 8 have short envelope mass-loss timescales  Gyr (assuming ). One could conceivably choose a threshold envelope loss timescale and then derive a lower bound on the planet mass at a given radius based on that assumption. We elaborate this possibility further in Section 6.3.

Figure 9.— Energy limited mass loss rates for the planet models in Figure 4. Mass loss rates are estimated for the case where the planets orbit a star with similar properties to our sun (, , and ). A mass loss efficiency of is assumed. The line styles have the same meanings and correspond to the same model planets as in Figure 4. Each curve corresponds to a different value of : 0.001 (thin solid), 0.01 (thin dashed), 0.05 (thin dot-dashed), 0.1 (thin dotted), 0.2 (thick solid), 0.3 (thick dashed), 0.4 (thick dot-dashed), and 0.5 (thick dotted). All data in this plot have , and (a)  K or (b)  K.

6. Discussion

6.1. Formation of Low-Density Neptune-Size Planets

Core Nucleated Accretion

Can core nucleated accretion form low-density planets in the size range of 2-6 ? The answer is yes, given appropriate conditions. The solids surface density in the protoplanetary disk must be appropriate for the accretion of heavy element cores a few times as massive as Earth. These cores must grow early enough to accrete significant gaseous envelopes, but gas accretion must end early enough to avoid runaway gas accretion. Our evolution calculations in Section 3 demonstrate that H/He rich planets can form for plausible choices of and disk lifetimes. The values we chose for are only slightly above that in the minimum-mass solar nebula, but high enough so that Jupiter at 5 AU can form in 4 Myr.

A second related question is whether core nucleated accretion with subsequent migration can lead to Neptune-size planets at high irradiation temperatures  K. Our evolution calculations uncover two factors that complicate achieving 2-6  planets following inward migration. First, high irradiation temperatures lead to very large fluffy planets with radii . Second, very long migration timescales are required to heat a planet to while keeping its envelope intact. We elaborate both of these points below.

The salient feature of our evolution calculations is that, despite , the irradiated planet radii at 1 and 4 Gyr are, in many cases, larger than . Specifically, all cases in Table 2 with  K or have radii in excess of . Lower mass envelopes are required to yield Neptune-size planets at these high irradiation levels (Figure 6). Truncating gas accretion earlier (shorter disk lifetime) and subsequent envelope mass loss are two potential avenues toward . While the model radii at 1000 K are very uncertain due to uncertainties in the opacities near the radiative-convective boundary, for the cases in Table 2 the conclusion that is, nonetheless, robust.

We found that slow planetary migration is needed for the low-mass envelopes to stay bound as the temperature at the planetary surface increases. In our evolution calculations, the planets initially assemble at  K or 125 K and then migrate inward to  K or 1000 K. The long migration timescale ( Myr) taken to reach  K with the envelope intact is in tension with typical disk lifetimes (1–10 Myr). The migration timescale to reach 500 K ( Myr) is more plausible. It is possible that evaporative cooling or increases in the envelope mean molecular weight from preferential loss of hydrogen could help the envelope remain bound. The planet evolution tracks presented do not include mass-loss.

There do exist Neptune-size equilibrium configurations at  K for planets with H/He envelopes from core nucleated accretion. Our equilibrium planet structure models in Section 2 explore and map out the (, , , ) parameter space that yield radii within the range 2–6  (Figure 6). It is important to note, however, that the equilibrium models provide “snap shots” of possible equilibrium configurations of planets undergoing quasi-static evolution. The models do not address how a planet could reach a given state, nor the timescale for the planet to evolve out of the state.

Outgassing of H

The second formation pathway to low-density Neptune-size planets we considered was outgassing of H from rocky planets. Outgassed low-mass planets without substantial HO envelopes, however, can only account for radii up to . Even achieving 3  with an outgassed envelope is a stretch, requiring (concurrently) a near-optimal initial planetesimal composition, full oxidization of accreted metals, and retention of most H released. Realistically, the majority of outgassed planets will be smaller than this radius limit, as we elaborate below.

The metal and HO content of the primordial material forming a planet set a strict limit on the amount of H that can be released via the outgassing reaction, 2Fe + 3H FeO + 3H. In this work, we have adopted a primordial chemical composition representative of EH chondrites mixed with just enough additional HO ice to fully oxidize all the metals. Out of the Solar System chondrites, this composition should be near optimum for outgassing of H due to the high proportion of unoxidized iron (in metal or sulfide form) (Elkins-Tanton & Seager, 2008b). Typically planets forming from a mixture of Solar System chondrite-like material (Jarosewich, 1990, within which the proportion of metallic iron varies from 0.1 to 22 wt. %) would have a lower capacity to outgas H.

Even given a high initial amount of reduced metals, a planet’s eventual outgassed envelope mass is contingent upon the fraction of metals that oxidizes. To bound the radii of outgassed planets, we considered the end-member case of complete oxidation of all Fe to FeO. In this limiting case, the planet is core-less with all its iron incorporated in the mantle as oxides (Elkins-Tanton & Seager, 2008a). Planets retaining a metallic core would degas less H. Ultimately, the overall fraction of Fe that reacts with water is determined by the competition between the rate of oxidation and the rate of sinking of metallic Fe to form the planet iron core. For a more detailed discussion see Elkins-Tanton & Seager (2008a).

Finally, the mass-radius relations for outgassed planets in Section 4 considered 100% retention of all outgassed H. Atmospheric escape leads to less massive H envelopes and smaller planets overall (Section 5). Indeed, while the primary outgassed atmospheres surrounding Earth and Mars during their accretion were likely H-dominated (Schaefer & Fegley, 2010), both planets today harbor secondary atmospheres with higher mean molecular weights.

How close can outgassed-planet radii plausibly get to the limiting outgassing relation? Relaxing our assumptions of optimum outgassing conditions, we investigate an intermediate, incomplete-oxidation case in which 50% of the accreted Fe is converted to FeO. This scenario leads to planets with 19.5% of their (initial) mass in an iron core, 0.6% by mass degassed in H, and a mantle Mg# of 0.62 (Table 3). With no loss of H these planets could have radii up to 2.7  (again considering ). Whereas, with atmospheric mass loss, planets that retain only 1 to 10% of the degassed H would have radii up to at most 2.4–2.5  for  K. Thus, radii up to are more realistically achieved by rocky planets with outgassed H envelopes but no free water. Planets with a water layer between the rocky interior and H envelope could be slightly larger, but only if little or no H was mixed in with the HO.

6.2. Maximum Planet Radius at Specified Mass

We have modeled the internal structure of low-mass, large-radius planets with hydrogen-dominated atmospheres. For planets with outgassed H envelopes, we derived a limiting low density relation by leveraging an upper bound on the amount of H that can be degassed from rocky planetesimals. The limiting low-density relation is less clear cut for planets formed from core nucleated accretion, because the initial reservoir of H/He accreted from the nebula need not be a constraining factor. Our detailed planet formation calculations provide discrete examples of planets at  K with only a few Earth masses yet radii larger than Jupiter.

The low-density limit for planets formed from core nucleated accretion depends on the heavy element core and envelope masses achievable at a given equilibrium temperature. The plausible combinations of (, , ) in turn rely on the protoplanetary disk properties and the migration history of the planet. The heavy element core mass is determined by the isolation mass, given the solid surface density and the distance from the star where the planet forms. The isolation mass (and thus ) can have a wide range of values, from less than to more than . The initial mass of H/He accreted by the planet is determined by the availability of a gas supply from the disk as governed by disk lifetime relative to the time taken for the heavy element core to reach isolation mass. Disk lifetimes range over an order of magnitude – from 1 to 10 Myr, with a characteristic value of a few Myr (Hillenbrand, 2008) – leading to some freedom in the initial expected from core nucleated accretion. Mass loss over the planet’s history would serve to decrease over time. Finally, the current equilibrium temperature depends on the migration history of the planet, and can, in principle, be anywhere from 100 K to 2000 K. Thus, due to the large spread in observed disk properties, a wide range of (, , ) from core nucleated accretion are plausible. We have shown detailed planet formation calculations for four reasonable choices of disk planetesimal densities and lifetimes.

We have succeeded in placing a tighter constraint on the low-density relation for outgassed planets than we have for planets from core nucleated accretion. This is due to the inherent limits on outgassed envelope masses; at very most, only a few percent of the mass of a planet can be outgassed in H. The end-member case of a planet that accreted from an optimum mixture of EH material and HO ice, where all the water and iron reacted, and where all released H was retained, sets an upper bound on the transit radius possible at a given mass for a rocky planet with out-gassed H atmosphere (Figure 7). Typically, rocky planets with out-gassed H atmospheres would have mean densities above this limiting relation. It should be noted, that our limiting relation applies to planets formed from material similar to Solar System chondrites. Planets formed from material with higher metallic Fe content would have the potential to outgas more H.

We have so far considered either core nucleated accretion or outgassing due to water-iron reactions as separate pathways for planets to acquire hydrogen rich envelopes. Core nucleated accretion contributes near solar composition material to the envelope (), while water-iron reactions contribute hydrogen but not helium (). If both processes occur on the same exoplanet, an envelope with intermediate, sub-solar, non-zero Helium content () may result.

The assumed chemical make-up of the planet envelope and heavy element core affect the planet relations for planets formed by core nucleated accretion and by outgassing. H/He envelopes in which He is depleted relative to solar will be more voluminous, for the same envelope masses, temperatures, and heavy element core properties. This is largely due to the influence of the mean molecular weight on the atmospheric scale height. For instance, decreasing Y=0.25 to Y=0.0 in the equilibrium planet models of Section 3.3, increases the radial extent of the envelopes by for . For lower mass planets, the change in the gravitational acceleration between the top and bottom of the envelope can be substantial and can have a larger effect on the envelope thickness. Pure H envelopes can be up to twice as thick as the corresponding H/He envelope, near the low mass extreme of the relations in Section 3.3. In our planet structure models, however, the effect of the envelope He abundance is largely offset by the higher density heavy element core composition in our outgassing models (EH chondrite cores) as compared to our core nucleated accretion models (ice-rock cores).

We have mapped out the contribution of low-mass planets with hydrogen-dominated atmospheres to the limiting low-density relation. Although we have not considered them in detail here, planets may also form with high molecular weight envelopes, for instance, after having accreted large amounts of ices beyond the snow line (e.g., Kuchner, 2003; Léger et al., 2004). Higher molecular weight envelopes are more dense (with smaller atmospheric scale heights) than their hydrogen-dominated counterparts, but may be less affected by atmospheric escape. It is possible that planets with high molecular weight atmospheres could also contribute to the limiting low-density relation for Neptune-size planets.

6.3. Minimum Planet Mass at Specified Radius

Our ideal goal was to determine a lower bound on the plausible planet mass given a planet radius in the range 2 – 6  and equilibrium temperature  K. We note that the relation defining the the maximum radius for a given planet mass does not necessarily translate into a relation for the minimum planet mass at a given radius. Indeed, at low masses, in the iso-composition relations for planets with gas envelopes (e.g., Figures 4, 5 and 7). Thus, in order to bracket the minimum planet mass of a transiting planet candidate, we must assess the survivability of low-mass planets for a range of interior compositions.

Mass loss is a major limiting factor that constrains the minimum for strongly irradiated ( K) Neptune-size planets harboring hydrogen dominated envelopes (Section 5). This is true whether the planet acquired its envelope through core nucleated accretion or through outgassing. If the heavy element core mass is small () and is high (1000 K) then the planet will not be able to hold on to very much gas. With the energy limited mass loss rates from Equation (13), we may roughly assess the survivability of potential planet configurations. By choosing a threshold envelope loss timescale , we can derive a lower bound on the planet mass at a given radius based on the requirement . To illustrate this approach, we adopt  Gyr and explore what this implies for planets with low mean molecular weight envelopes from core nucleated accretion (Figure 6) and from outgassing (Figure 8).

We estimate, using Figure 6, the minimum masses of planets with H/He envelopes formed by core nucleated accretion beyond the snow line. For planets, the least massive planet models that satisfy  Gyr are 1.3 to at  K, and 4.0 to at  K (for between and ). Analogously, at , the  Gyr survivability constraint requires that to at  K, and to at  K. At , almost all possible (, , , ) configurations in Figure 6 have sub-gigayear envelope-loss timescales, due to the small planet and envelope masses (, ). An ice/rock core surrounded by an H/He envelope from core nucleated accretion may not be a plausible