A water budget dichotomy of rocky protoplanets from $^{26}$Al-heating

TL;DR

This study uses numerical models to show that initial $^{26}$Al levels influence the water content and size of rocky planets, explaining the diversity of planetary compositions observed.

Contribution

It introduces a new model linking $^{26}$Al heating to planetary water content and size, providing insights into planetary formation diversity.

Findings

High $^{26}$Al leads to water-depleted, smaller planets.

Low $^{26}$Al results in larger, ocean-rich worlds.

The model explains the lack of water trend in TRAPPIST-1 planets.

Abstract

In contrast to the water-poor inner solar system planets, stochasticity during planetary formation and order of magnitude deviations in exoplanet volatile contents suggest that rocky worlds engulfed in thick volatile ice layers are the dominant family of terrestrial analogues among the extrasolar planet population. However, the distribution of compositionally Earth-like planets remains insufficiently constrained, and it is not clear whether the solar system is a statistical outlier or can be explained by more general planetary formation processes. Here we employ numerical models of planet formation, evolution, and interior structure, to show that a planet's bulk water fraction and radius are anti-correlated with initial $^{26}$Al levels in the planetesimal-based accretion framework. The heat generated by this short-lived radionuclide rapidly dehydrates planetesimals prior to accretion onto larger protoplanets and yields a system-wide correlation of planet bulk abundances, which, for instance, can explain the lack of a clear orbital trend in the water budgets of the TRAPPIST-1 planets. Qualitatively, our models suggest two main scenarios of planetary systems' formation: high-$^{26}$Al systems, like our solar system, form small, water-depleted planets, whereas those devoid of $^{26}$Al predominantly form ocean worlds, where the mean planet radii between both scenarios deviate by up to about 10%.

Figures

Equations

$$H_{{}^{26}\mathrm{Al}}(t)=f_{\mathrm{Al}}\cdot {(}^{26}\mathrm{Al}{/}^{27}\mathrm{Al}{)}_0\cdot E_{{}^{26}\mathrm{Al}}\cdot \exp (-t/{\tau }_{{}^{26}\mathrm{Al}})/{\tau }_{{}^{26}\mathrm{Al}},$$

$$k_{\mathrm{eff}}=(q/0.89{)}^{3/2}\cdot {\alpha }_{\mathrm{liq}}g{c}_{\mathrm{p}}/(\Delta T^2{\rho }_{\mathrm{s}}{\eta }_{\mathrm{num}}),$$

$$f_{{\mathrm{H}}_2\mathrm{O}}(t)/f_{{\mathrm{H}}_2\mathrm{O},\mathrm{init}}=1-X_{\mathrm{dry}}(t)/X_{\mathrm{plts}},$$

$$\frac{\mathrm{d}\Sigma }{\mathrm{d}t}=\frac{3}{r}\frac{\partial }{\partial r}\left[r^{1/2}\frac{\partial }{\partial r}\tilde{\nu }\Sigma r^{1/2}\right]+{\dot{\Sigma }}_{\mathrm{w}}+{\dot{Q}}_{\mathrm{planet}},$$

$$R_{\mathrm{H}}=a_{\mathrm{planet}}{\left[M_{\mathrm{planet}}/(3M_{\mathrm{star}})\right]}^{1/3},$$

$${\overline{\rho }}_{\mathrm{plts}-\mathrm{ice}}(t)={\rho }_{{\mathrm{H}}_2\mathrm{O}}{f}_{{\mathrm{H}}_2\mathrm{O}}(t)+{\rho }_{\mathrm{rock}}[1-f_{{\mathrm{H}}_2\mathrm{O},\mathrm{init}}].$$

$${\Sigma }_{\mathrm{plts}-\mathrm{ice}}(t)={\Sigma }_{\mathrm{plts}-\mathrm{ice},\mathrm{init}}\cdot \left(f_{{\mathrm{H}}_2\mathrm{O}}(t)+[1-f_{{\mathrm{H}}_2\mathrm{O},\mathrm{init}}]\right).$$

$$R_{\mathrm{star}}=(M_{\mathrm{star}}/M_{\odot }{)}^{0.945}{R}_{\odot },$$

$$L_{\mathrm{star}}=0.628(M_{\mathrm{star}}/M_{\odot }{)}^{2.62}{L}_{\odot },$$

$$T_{\mathrm{star}}=\sqrt[4]{L_{\mathrm{star}}/(4\pi R_{\mathrm{star}}^2\sigma )},$$

$$M_{\mathrm{disk}}\propto (M_{\mathrm{star}}/M_{\odot }{)}^{1.2},$$

$$\frac{\partial m}{\partial r_{\mathrm{p}}}$$

$$\frac{\partial P}{\partial r_{\mathrm{p}}}$$

$$\frac{\partial T}{\partial r_{\mathrm{p}}}$$

Full text

11institutetext: Institute of Geophysics, ETH Zürich, Sonneggstrasse 5, Zürich 8092, Switzerland

22institutetext: Bayerisches Geoinstitut, University of Bayreuth, Universitätsstrasse 30, Bayreuth 95440, Germany

33institutetext: Physikalisches Institut, University of Bern, Sidlerstrasse 5, Bern 3012, Switzerland

44institutetext: Department of Astronomy, University of Michigan, 1085 S. University Avenue, Ann Arbor, MI 48109, United States

55institutetext: Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, Bern 3012, Switzerland

A water budget dichotomy of rocky protoplanets from 26Al-heating

Tim Lichtenberg

  

Gregor J. Golabek A water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heating tim.lichtenberg@physics.ox.ac.uk Corresponding author, present address: Atmospheric, Oceanic and Planetary Physics, University of Oxford, Parks Rd, Oxford OX1 3PU, United Kingdom, Corresponding author, present address: Atmospheric, Oceanic and Planetary Physics, University of Oxford, Parks Rd, Oxford OX1 3PU, United Kingdom,   

Remo Burn A water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heating   

Michael R. Meyer A water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heating   

Yann Alibert A water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heating   

Taras V. Gerya A water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heating   

Christoph Mordasini A water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heatingA water budget dichotomy of rocky protoplanets from 26Al-heating

(February 11, 2019; preprint differing from journal version (free-to-read pdf))

In contrast to the water-poor inner solar system planets, stochasticity during planetary formation (Alibert & Benz 2017; Raymond & Izidoro 2017) and order of magnitude deviations in exoplanet volatile contents (Kaltenegger 2017) suggest that rocky worlds engulfed in thick volatile ice layers (Kuchner 2003; Léger et al. 2004) are the dominant family of terrestrial analogues (Tian & Ida 2015; Ramirez & Levi 2018) among the extrasolar planet population. However, the distribution of compositionally Earth-like planets remains insufficiently constrained, and it is not clear whether the solar system is a statistical outlier or can be explained by more general planetary formation processes. Here we employ numerical models of planet formation, evolution, and interior structure, to show that a planet’s bulk water fraction and radius are anti-correlated with initial 26Al levels in the planetesimal-based accretion framework. The heat generated by this short-lived radionuclide rapidly dehydrates planetesimals (Grimm & McSween 1993) prior to accretion onto larger protoplanets and yields a system-wide correlation (Millholland et al. 2017; Weiss et al. 2018) of planet bulk abundances, which, for instance, can explain the lack of a clear orbital trend in the water budgets of the TRAPPIST-1 planets (Dorn et al. 2018). Qualitatively, our models suggest two main scenarios of planetary systems’ formation: high-26Al systems, like our solar system, form small, water-depleted planets, whereas those devoid of 26Al predominantly form ocean worlds, where the mean planet radii between both scenarios deviate by up to $\approx$10%.

Main Text

In the early solar system, the decay heat from the short-lived radionuclide 26Al ($t_{\mathrm{1/2,^{26}Al}}$ $\approx$ 0.72 Myr) powered the interior evolution of planetesimals, the seeds and building blocks of the rocky planets, and led to silicate melting (Fu & Elkins-Tanton 2014; Lichtenberg et al. 2016a) and degassing of primordial water abundances (Grimm & McSween 1993; Monteux et al. 2018). Here, we explore the systematic effects of 26Al on rocky planetary systems using numerical models of planetary formation (Benz et al. 2014) with 26Al-induced water loss from planetesimals during the main accretion phase (Monteux et al. 2018). We generate synthetic planet populations with internal structures defined by the planets’ composition, which result in statistical variations of planet water abundance and (transit) radius.

In the models presented, initially Moon-sized protoplanets grow from the accretion of 1–100 km-sized planetesimals and gas, and migrate within the protoplanetary disk of G or M-type systems (Methods). The initial location of the embryos and the starting disk structures and boundaries are randomized to reflect the diversity found in observed young planetary systems (Ansdell et al. 2016). Planetesimals are set to be dry within the snowline and icy outside, with a decreasing water mass fraction over time, calculated from planetesimal interior models that account for the dehydration from internal radiogenic heating of 26Al. Here, we account for accretion of planetesimals only, and ignore the potential contribution from smaller particles, such as pebbles (Johansen et al. 2015; Schiller et al. 2018; Alibert et al. 2018). The heating rate in the planetesimal interior is controlled by the amount of 26Al incorporated upon planetesimal formation, which may vary substantially between planetary systems (Lichtenberg et al. 2016b; Lugaro et al. 2018). We account for this variability by generating synthetic planet populations with different planetesimal radii, $r_{\mathrm{plts}}$ = 3, 10, 50 km, and initial 26Al abundances of 26Al${}_{\mathrm{0}}\in$ [0.1, 10] $\times$ 26Al*⊙, with 26Al⊙* the solar system’s ‘canonical’ (26Al/27Al)0 at CAI formation, and compare them to a nominal case without 26Al-heating (for further details on the models see Methods). For each combination of $r_{\mathrm{plts}}$, 26Al0 and stellar type (G or M) we performed 30’000 single planet simulations (Alibert & Benz 2017; Benz et al. 2014), resulting in a statistically representative set of 540’000 individual simulations over 18 parameter sets (cf. Fig. 1).

The control of 26Al0 and $r_{\mathrm{plts}}$ on the retention of water within planetesimals and resulting planet populations from a given set of initial conditions are shown in Fig. 1. Planetesimals with larger $r_{\mathrm{plts}}$ and higher 26Al0 dehydrate faster and up to 100% for extreme values. Rooted in our conservative choice for dehydration (Methods), the total water loss divides the parameter range into two distinct regimes. The first consists of almost pristine water-rock ratios for small planetesimals with low 26Al0. However, for 26Al0 $\gtrsim$ 26Al*⊙* and planetesimals with $r_{\mathrm{plts}}$ $\gtrsim$ 10 km, water loss is nearly complete (Fig. 1a,b).

For distinct combinations of 26Al0 and $r_{\mathrm{plts}}$, we simulate the influence on the expected planet population for planet masses $M_{\mathrm{P}}$ $\in$ [0.1, 10] $\times$ $M_{\mathrm{Earth}}$ (Fig. 1c). Because the timescale for water-loss caused by 26Al-heating is significantly shorter than the accretion timescale, sufficiently 26Al-enriched planetesimals are mostly dry when they accrete onto protoplanets. Therefore, the final planet water mass fractions are correlated with the retained water fraction in planetesimals due to 26Al-heating. The planet desiccation caused by the accretion of ever-more dehydrated planetesimals reduces the inherent scatter and range in $f_{\mathrm{H_{2}O}}$ in the synthetic planet populations (cf. Fig. 1c and 2a,b). For fixed planetesimal radius and increasing 26Al0, accreting planets receive more relative mass contribution from dry objects and end up water-depleted relative to nominal conditions.

In the solar system, the initial planetesimal size frequency distribution is expected to have been dominated by bodies with $r_{\mathrm{plts}}$ $\gtrsim$ 30–50 km (Delbo et al. 2017). For such bodies, the equilibrium between radiogenic heating and surface cooling stabilizes internal temperatures for an extended timespan at spatially isothermal conditions (Lichtenberg et al. 2016a). Therefore, the fractional dehydration in Fig. 1b flattens above $r_{\mathrm{plts}}$ $\gtrsim$ 50 km and becomes nearly independent of planetesimal size. For $r_{\mathrm{plts}}$ = 50 km, dehydration is dominantly controlled by 26Al0 and generates a dichotomy between planets in 26Al-enriched (top black and red histograms, Fig. 2) versus non-enriched systems (top blue histogram, Fig. 2). M and G stars overall display a similar trend, but M stars form smaller planets on average, due to their lower initial budget of planet-forming material.

The emerging trend from our simulations is illustrated in Fig. 3, with a clear distinction between planetary systems that are significantly enriched (26Al0 $\gtrsim$ 26Al*⊙*), and those that are not. In general, 26Al is expected to be abundant but inhomogeneously distributed within young star-forming regions (Lichtenberg et al. 2016b; Kuffmeier et al. 2016). According to our simulations, planets in enriched systems grow from ever-more dehydrated planetesimals and form desiccated planets in their terrestrial planet zone. Depending on the initial planetesimal sizes, final planet water fractions are up to two orders of magnitude below the initial planetesimal water mass fractions, and are strongly correlated with the efficency of dehydration during accretion (Fig. 1).

The bulk volatile mass fraction has the greatest influence on the structure and mass-radius relation of a rocky planet (Noack et al. 2017; Unterborn et al. 2018). Therefore, we anticipate the resulting smaller radii (from lower water mass fraction) for higher 26Al levels to be reflected in the galactic exoplanet population. For deviations in planet bulk water fractions predicted here, the thickness of the volatile layer on top of the silicate mantle constitutes several per cent of the radius (Kuchner 2003; Léger et al. 2004; Alibert 2014; Noack et al. 2017). We calculate this deviation in our synthetic populations by translating the derived planetary masses and compositions into a mean radius in a given mass bin (Fig. 4) using interior structure models that are sensitive to the total planet mass, its water and (captured) hydrogen/helium mass fraction, and the surface pressure (Methods). For the entire populations of planets among G and M stars, the radius deviation reaches up to 2% for 1 $\times$ 26Al*⊙, and can go up to $\approx$4% deviation for Mars-sized planets for 10 $\times$ 26Al⊙* or $r_{\mathrm{plts}}$ = 50 km, respectively. If we only consider planets that accrete a minimum amount of water (planets that receive some mass contribution from beyond the iceline), with 10 $\times$ 26Al*⊙, or 1 $\times$ 26Al⊙* with $r_{\mathrm{plts}}$ = 50 km, the mean-radius shift can reach up to $\approx$10%. Planetary systems with high 26Al0 (26Al0 $\gtrsim$ 1–10 $\times$ 26Al*⊙*) form water-depleted planets and with system-wide smaller radii compared to the non-enriched population.

Such deviations are expected to be measurable by the planned PLATO mission (Rauer et al. 2014), which will aim to characterize a statistical ensemble of planetary radii in the rocky planet regime. The intrinsic compositional scatter in the inferred mean densities from known exoplanets suggests a large stochastic component in the planet formation process. Yet, recent analyses of data based on Kepler multi-planet systems provided strong evidence for intra-system correlation between planetary radii (Millholland et al. 2017; Weiss et al. 2018). Therefore, in the exoplanet census probed so far, the fate and long-term structure of planets seems to be dominated by physical and chemical effects on a system-to-system level, rather than emerging from intra-system stochasticity during accretion, such as impact stripping (Marcus et al. 2010; Inamdar & Schlichting 2016). With future access to a statistical ensemble of low-mass planet radii from exoplanet-focused missions, the highly 26Al-enriched systems, such as the solar system, where planetary radii deviate by several per cent from the population of 26Al-poor systems, may stick out from the mean of the population and provide clues about the underlying 26Al distribution of planetary systems.

For example, the system-wide water depletion of the TRAPPIST-1 planets (Grimm et al. 2018; Dorn et al. 2018) is consistent with 26Al $\gtrsim$ 26Al*⊙*-induced desiccation (Fig. 2). The atmospheres of the TRAPPIST-1 planets seem to be secondary (de Wit et al. 2018), and may have lost several Earth ocean equivalents of water (Bourrier et al. 2017). However, to account for the consistency of especially the outermost planets e-h with near Earth-like volatile abundances and lack of an orbital trend in water budget (Dorn et al. 2018), an order of magnitude depletion mechanism, such as suggested here, must affect all of the planets. Therefore, the retrieved low water mass fractions of the TRAPPIST-1 (Grimm et al. 2018; Dorn et al. 2018) planets are unexpected from formation and evolution models (Tian & Ida 2015; Ciesla et al. 2015; Ormel et al. 2017; Alibert & Benz 2017; Unterborn et al. 2018), and present a severe challenge for current planet formation scenarios. The 26Al desiccation mechanism we put forward achieves system-wide water depletion for G and M stars without the need to fit specific accretion dynamics, as it has been proposed (Ormel et al. 2017; Unterborn et al. 2018). However, because the TRAPPIST-1 system is just the first of perhaps many such systems, coordinated observational efforts will be required to establish population characteristics for similar systems, in order to distinguish between 26Al desiccation and migration-driven mechanisms (Ormel et al. 2017; Unterborn et al. 2018) as the origin of the TRAPPIST-1 planet compositions.

For accretion scenarios where planetesimals represent the primary carrier of water, our models suggest that planetary systems with 26Al abundances similar to or higher than the solar system generically form terrestrial planets with low water mass fractions, $f_{\mathrm{H_{2}O}}$ $\lesssim$ $\mathcal{O}$(wt%). This effect is more pronounced for planets further out from their host star, as embryos in these regions grow preferentially from water-rich solids. For a non-uniform distribution of 26Al in Milky Way star-forming regions, the systematic water depletion in 26Al-enriched systems suggests the existence of two qualitatively distinct classes of planetary systems: water-poor (26Al-rich) and water-rich (26Al-poor) systems, with a systematic mean-radius deviation for sub-Earth terrestrial planets between these classes. The resulting shape of the distribution of dry and wet planetary systems depends on the genuine, but unknown, distribution of 26Al levels (Lugaro et al. 2018) among planet-forming systems and the nature and timing of protoplanet accretion. If rocky planets grow primarily from the accumulation of planetesimals, then the suggested deviation between planetary systems should be clearly distinguishable among the rocky exoplanet census. If, however, the main growth of rocky planets may proceed from the accumulation of small particles, such as pebbles, then the deviation between 26Al-rich and 26Al-poor systems may become less clear, and the composition of the accreting pebbles needs to be taken into account. Therefore, in future work, models of water delivery and planet growth need to synchronize the timing of earliest planetesimal formation (Dra̧żkowska & Dullemond 2018), the mutual influence of collisions (Lichtenberg et al. 2018) and 26Al dehydration, the potential growth by pebble accretion (Johansen et al. 2015; Schiller et al. 2018; Alibert et al. 2018), and the partitioning of volatile species between the interior and atmosphere of growing protoplanets (Ikoma et al. 2018) in order to further constrain the perspectives for rocky exoplanet evolution (Kite & Ford 2018).

Methods

Planetesimal dehydration

We model water loss from instantaneously-formed planetesimals composed of a rock-ice mixture using numerical models that employ a conservative finite-differences, fully-staggered grid method coupled to a marker-in-cell approach (Gerya & Yuen 2007; Golabek et al. 2014). The thermo-chemical evolution of planetesimals is computed in a two-dimensional infinite cylinder geometry on a Cartesian grid, solving the Poisson, continuity, Stokes and energy conservation equations. We assume the planetesimals to be accreted with the temperature of the protoplanetary disk beyond the water snowline, $T_{0}=150$ K, which is kept constant during the evolution of the planetesimal utilizing the free-surface ‘sticky-air’ method (Crameri et al. 2012). Heating is provided by the decay of 26Al, which defines the radiogenic heat source term over time

[TABLE]

with the chondritic abundance of aluminum, $f_{\mathrm{Al}}$ (Lodders 2003), the ratio of 26Al to stable 27Al at the time of planetesimal formation, 26Al0 = (26Al/27Al)0, the decay energy, $E_{\mathrm{{}^{26}\mathrm{Al}}}=3.12$ MeV (Castillo-Rogez et al. 2009), and the mean lifetime, $\tau_{{}^{26}\mathrm{Al}}=1.03$ Myr. We ignore any potential heat contribution from 60Fe, which may further boost radiogenic heating rates in extrasolar systems (Lichtenberg et al. 2016b; Nicholson & Parker 2017). If the planetesimal interior reaches temperatures beyond the rock disaggregation threshold (Costa et al. 2009) at a silicate melt fraction of $\phi\gtrsim 0.4$, where the rock viscosity drops by more than ten orders of magnitude, we approximate the thermal conductivity in the soft turbulence limit (Siggia 1994) with

[TABLE]

with the convective heat flux, $q$, the temperature difference across nodes, $\Delta T$, silicate density, $\rho_{\mathrm{s}}$, thermal expansivity of molten silicates, $\alpha_{\mathrm{liq}}$, silicate heat capacity, $c_{\mathrm{p}}$, local gravity, $g(x,y)$, and lower cut-off viscosity, $\eta_{\mathrm{num}}$. For numerical values used and further details and references on the code see Lichtenberg et al. (2016a). The initial planetesimal water-to-rock ratio beyond the snowline is expected to be between $\approx$0.05 (Ciesla et al. 2015; O’Brien et al. 2018), the water content of carbonaceous chondrites, and $\approx$0.5, as suggested by equilibrium condensation calculations (Lodders 2003). Here, we adopt a value closer to the upper estimate, $f_{\mathrm{H_{2}O,init}}$ = 0.3, but our calculations only marginally depend on the adopted value.

In general, during heat-up of a primordial water ice-rock mixture, ices melt and react with the ambient rock. The liquid water undergoes pore water convection and escapes quickly once the gas phase is reached (Grimm & McSween 1993), but a small fraction of water may be trapped in hydrous silicate phases. Therefore, we numerically account for dehydration of parts of the planetesimal interior at a conservative upper limit of $T\geq T_{\mathrm{dry}}$ = 1223 K, the upper limit of the amphibolite stability field, when any possibly remaining hydrous silicate phases break down. At these high temperatures exsolved water vapor is lost quasi-instantaneously because planetesimals of this size cannot preserve an outgassed atmosphere. We do not resolve potential earlier water-loss from degassing (Castillo-Rogez & Young 2017), residual volatiles above $T_{\mathrm{dry}}$ (Fu et al. 2017), or ice sublimation during late and optically thin disk stages (Machida & Abe 2010). Using these assumptions, we compute the expected ratio of dehydrated to primordial water-rock mixture at time $t$ due to degassing,

[TABLE]

with the dry fraction, $X_{\mathrm{dry}}(t)$, of the total planetesimal interior, $X_{\mathrm{plts}}$, and the initial water-to-rock ratio, $f_{\mathrm{H_{2}O,init}}=0.3$. Under these conditions, a planetary system in the planet formation model is represented by an initial 26Al0 that corresponds to the time of planetesimal formation. If the 26Al content may vary spatially within the disk, as it was suggested (Larsen et al. 2011; Schiller et al. 2015), the solar system itself would be represented by a sub-canonical (26Al0 $\leq$ 26Al*⊙*) value, similar to the effects of delayed planetesimal formation (Fig. 1b, cf. Lichtenberg et al. (2018) for a discussion of the effects on planetesimal evolution).

Planet formation

We compute the formation of planets and generate our synthetic planet populations using an updated version of the model of Alibert et al. (2005). The computer code numerically treats the structure and evolution of the protoplanetary disk, the dynamical properties and accretion rate of planetesimals onto accreting protoplanets, the planetary envelope structure, and disk-planet interactions (Mordasini et al. 2009a, b, 2015; Benz et al. 2014). Here, we provide a brief summary of the most important code modules used in this work.

The protoplanetary disk model relies on the Shakura-Sunyaev (Shakura & Sunyaev 1973) disk viscosity approximation ($\alpha_{\mathrm{disk}}=2\times 10^{-3}$) and computes the surface density evolution over time by solving the radial diffusion equation,

[TABLE]

with the surface density, $\Sigma$, orbital radius, $r$, effective viscosity, $\tilde{\nu}$, and gas accretion onto embryos, $\dot{Q}_{\mathrm{planet}}$, calculated from removing gas in an annulus centered on the embryo with a width of one Hill radius,

[TABLE]

with the planet semi-major axis, $a_{\mathrm{planet}}$, planet mass, $M_{\mathrm{planet}}$, and star mass, $M_{\mathrm{star}}$. Mass loss due to internal (EUV) photoevaporation $\dot{\Sigma}_{\mathrm{w}}$ (Clarke et al. 2001) is set $\propto r^{-5/2}$ outside of a gravitational radius of $\approx 5$ au and external (FUV) photoevaporation (Matsuyama et al. 2003) is constant outside of $\approx 140$ au, with the total mass loss being a free model parameter. The model used to represent the planetesimal disk relies on the initial central temperature and pressure from the disk model to compute the location of the water snowline, thereby neglecting radial drift of planetesimals (Sasselov & Lecar 2000) and that of the snowline. Drift timescales for planetesimals larger than 1 km exceed the disk lifetime by orders of magnitude (Weidenschilling 1977).

We consider rocky planetesimals ($\overline{\rho}_{\mathrm{plts-dry}}=\rho_{\mathrm{rock}}$ = 3200 kg/m3) inside, and rock-ice aggregates ($\overline{\rho}_{\mathrm{plts-ice,init}}=\rho_{\mathrm{H_{2}O}}f_{\mathrm{H_{2}O,init}}+\rho_{\mathrm{rock}}[1-f_{\mathrm{H_{2}O,init}}]$) beyond the snowline, which are fixed in radius and accrete onto the planetary embryo that is embedded in the disk in a single simulation. The residual water mass fraction, $f_{\mathrm{H_{2}O}}(t)$, of the accreting planetesimals is computed from the internal evolution (Planetesimal dehydration) and is translated to a decreasing planetesimal density, $\overline{\rho}_{\mathrm{plts-ice}}(t)$, and disk solid surface density, $\Sigma_{\mathrm{plts-ice}}(t)$, by reducing the planetesimal density as

[TABLE]

In this formulation, lost water is assumed to be replaced by pore space and the planetesimal radius stays constant. The solid surface density available for embryos to accrete beyond the iceline thus changes with

[TABLE]

In our nominal model, a single embryo of initially lunar mass, $M=0.0123$ $M_{\mathrm{Earth}}$, is placed randomly between specific inner and outer bounds within the protoplanetary disk (Parameter space), with a dry composition inside the snowline, and wet outside. It starts accreting solids (planetesimals) and gas, and may migrate in the type I and II regime, depending on the embryo mass and physical structure of the disk at a given orbit (Dittkrist et al. 2014). The solid accretion rate (Inaba et al. 2001; Inaba & Ikoma 2003) takes into account the captured atmosphere. Planetesimal excitation and damping is computed by taking into account self interactions and damping by gas drag (Fortier et al. 2013). We ignore water loss due to collisions, which may further reduce the water inventory (Genda & Abe 2003; Schlichting et al. 2015; Inamdar & Schlichting 2016; Burger et al. 2018) dependent on the frequency of such interactions, and accretion of smaller solid particles (Johansen et al. 2015; Bollard et al. 2017; Schiller et al. 2018) (’pebbles’) that may shift the ratio of dry to wet accreted primitive materials (Connelly & Bizzarro 2016). Gas accretion due to planetary contraction is considered using a dust opacity reduction factor of $0.01$ compared to interstellar values (Pollack et al. 1996; Alibert et al. 2013).

G star settings are identical to the Sun’s values. The properties for the M star runs are scaled down. We choose a fixed mass of $M_{\mathrm{star}}=0.2\;M_{\odot}$ for the M stars. The radius of the star is set to

[TABLE]

with luminosity (Demircan & Kahraman 1991)

[TABLE]

and temperature

[TABLE]

with stellar radius, $R_{\mathrm{star}}$, stellar mass, $M_{\mathrm{star}}$, and Stefan-Boltzmann constant, $\sigma$. The disk dimensions, exponential cut-off radius and the embryo placement boundaries (Parameter space) are reduced to account for the lower masses and sizes of M star disks. Thus, initially all embryos form closer to the star compared to the G star populations. The initial disk mass follows the scaling law (Alibert et al. 2011)

[TABLE]

with the internal photoevaporation rate adapted to match similar mean lifetimes compared to the G star simulations. In reality, these could be anti-correlated with stellar mass (Kennedy & Kenyon 2009), which would increase the efficiency of the 26Al-dehydration mechanism for M stars due to longer accretion timescales.

Interior structure & evolution

The interior structure and the long-term evolution of the planets is calculated as described in Mordasini et al. (2012) by solving the classical one-dimensional radially symmetric interior structure equations of mass conservation, hydrostatic equilibrium, and energy transport (Bodenheimer & Pollack 1986)

[TABLE]

where $r_{\mathrm{p}}$ is the radial distance from the planet’s center, $m$ the enclosed mass, $P$ the pressure, $\rho$ the density, and $G$ the gravitational constant. The intrinsic luminosity of a planet is assumed to be constant as a function of the planet radius. The gradient $\nabla$ depends on the process by which the energy is transported (radiative diffusion or convection). These calculations yield the radii of the planets given their mass and bulk composition, namely the mass fractions of iron, silicates, water, and H/He, as an output from the planet formation and planetesimal dehydration models. For the H/He envelope, the equation of state of Saumon et al. (1995) is used to solve the structure equations, while for the solid part of the planet, including the water content, the modified polytropic equations of state of Seager et al. (2007) are employed. The transit radius is estimated as in Guillot (2010).

The loss of the primordial H/He envelope by atmospheric escape is considered in the energy- and radiation-recombination-limited approximation as described in Jin et al. (2014) and results in the loss of the primary atmosphere for low-mass planets at smaller orbital distances. Because of the limited water solubility in silicate mantles, the radius of planets without primordial H/He envelopes depends strongly on the water mass fraction (Jin & Mordasini 2018), and thus reveals the dehydration pattern caused by different contents of 26Al. Here, we do not treat interior-atmosphere exchange during early magma ocean phases that may further fractionate the volatile distribution within the body, in particular for close-in planets (Ikoma et al. 2018).

Parameter space

From $\gamma$-ray observations there is evidence for a widespread and heterogeneous distribution of 26Al in the galaxy (Lugaro et al. 2018). Observational evidence from young star-forming regions and theoretical work suggest a non-uniform enrichment pattern among planetary systems (Adams et al. 2014; Gounelle 2015; Pfalzner et al. 2015; Lichtenberg et al. 2016b; Kuffmeier et al. 2016; Parker et al. 2017; Nicholson & Parker 2017; Dwarkadas et al. 2017; Lugaro et al. 2018) with order of magnitude deviations from the solar system’s ‘canonical’ 26Al value (Kita et al. 2013) of 26Al*⊙* = (26Al/27Al)0 $=5.25\times 10^{-5}$. To account for these variations, we consider values in the range 26Al0 $\in$ [0.1, 10] $\times$ 26Al*⊙*. In addition to initial 26Al abundance, the radii of planetesimals during accretion yield different thermal evolutionary sequences and thus dehydration patterns (Lichtenberg et al. 2016a; Monteux et al. 2018). Here, we test values in the range $r_{\mathrm{plts}}$ $\in$ [1, 100] km. However, we note that from asteroid-belt inferences and numerical studies of the streaming instability mechanism, radii larger than $\gtrsim$ 30–50 km (Johansen et al. 2015; Klahr & Schreiber 2016; Simon et al. 2017; Delbo et al. 2017; Tsirvoulis et al. 2018) are expected. All parameter models not listed in the Methods are identical to those used in Lichtenberg et al. (2016a) and Monteux et al. (2018). In the planet formation model, the innermost disk radius is of the order 0.1 au and can vary over time. Disk lifetimes are distributed around 5 Myr, which is controlled via the photoevaporation rate (Fortier et al. 2013) and in agreement with current disk surveys (Meng et al. 2017; Kral et al. 2017). The initial embryos are placed within the boundaries of [0.05, 40] au for G stars, and [0.086, 23.4] au for M stars. We vary in a Monte Carlo fashion (Mordasini et al. 2009a) the disk mass, lifetime, dust-to-gas ratio and the exponential cut-off radius (Mordasini et al. 2012) to represent the diversity found in nature (Andrews et al. 2010; Ansdell et al. 2016).

Code and data availability

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Acknowledgements.

The authors thank Richard Parker for helping to initiate this project, Edwin Kite, Thibaut Roger, Caroline Dorn for comments and discussions, and Martin Bizzarro and two anonymous referees for comments that helped to considerably improve the manuscript. T.L. was supported by ETH Zürich Research Grant ETH-17 13-1 and acknowledges partial financial support from the Swiss Society for Astrophysics and Astronomy trough a MERAC travel grant. Y.A. acknowledges support from the Swiss National Science Foundation (SNSF). C.M. acknowledges support from the SNSF under grant BSSGI0$\_$155816 ‘PlanetsInTime’. The numerical simulations in this work were partially performed on the EULER computing cluster of ETH Zürich. Parts of this work have been carried out within the framework of the National Center for Competence in Research PlanetS supported by the SNSF. We acknowledge the software usage of matplotlib (Hunter 2007), scipy (Jones et al. 2001–), numpy (Van Der Walt et al. 2011), pandas (McKinney 2010), and seaborn (Waskom et al. 2018).