© ESO, 2016

1. Introduction

The characterization of exoplanet interiors is key to understand planet diversity. Here we focus on characterizing six exoplanets for which mass, radius, and also refractory element abundances of the stellar hosts are known. These exoplanets are HD 219134b, Kepler-10b, Kepler-93b, CoRoT-7b, 55 Cnc e, and HD 97658b. In order to make meaningful statements about their interior structures, it is mandatory not only to find few interior realizations that explain the data, but also to quantify the spread in possible interior structures that are in agreement with data, and data and model uncertainty (Rogers & Seager 2010; Dorn et al. 2015, 2017). This spread is generally large, due to the inherent degeneracy, i.e. different interior realizations can have identical mass and radius. We therefore apply the methodology presented by Dorn et al. (2017) that employs a full probabilistic Bayesian inference analysis using a Markov chain Monte Carlo (McMC) method. This method is validated against Neptune and previously by Dorn et al. (2015) against the terrestrial fitsyno planets. It allows us to determine confidence regions for structural and compositional parameters of core, mantle, ice layer, ocean, and atmosphere.

Planet bulk abundance constraints, in addition to mass and radius, are crucial to reduce the generally high degeneracy in interior structures (Dorn et al. 2015; Grasset et al. 2009). Stellar abundances in relative refractory elements (Fe/Si and Mg/Si) were suggested to serve as a proxy for the planetary abundance (e.g., Bond et al. 2010; Elser et al. 2012; Johnson et al. 2012; Thiabaud et al. 2015). Spectroscopic observations of relative photospheric abundances of host stars may allow us to determine their importance for the interior characterization of planetary companions as outlined by Santos et al. (2015). We proceed along these lines and study the sensitivity of interior parameter predictions to different planet bulk abundances by testing (1) different stellar abundance estimates and (2) different correlations between stellar and planetary abundance ratios.

Over the last few decades, a variety of telescopes and spectrographs have been utilized to measure stellar abundances, which affects the resolution and signal-to-noise (S/N) ratios of the data. In addition, the techniques used to determine the element abundances differ. Then, each research group chooses their own line list and number of ionization stages included, solar atmosphere models, and adopted solar abundances. All of these differences result in element abundance measurements that are fundamentally on different baselines and which do not overlap to within error (Hinkel et al. 2014). Here, we investigate whether these discrepancies significantly affects the determination of interior structure or not.

The six selected exoplanets for which host star abundances are available have bulk densities in the range from 6.8 g/cm³ (Kepler-93b) to 3.4 g/cm³ (HD 97658b) and span across the proposed transition range between mostly rocky and mostly non-rocky exoplanets (e.g., Lopez & Fortney 2013; Rogers 2015). Terrestrial-type planets are generally thought to be differentiated with an iron-rich core, a silicate mantle, and a crust, whereas volatile-rich planets are thought to also contain significant amount of ices and/or atmospheres (e.g., see Howe et al. 2014). Direct atmosphere characterization for the selected planets are not available, except for the only recently characterized atmosphere of 55 Cnc e (Tsiaras et al. 2016; Demory et al. 2016b). For super Earths in general, the diversity of atmospheric structures and compositions are mostly theoretical (e.g., see Burrows 2014, for a review). Studies of mass-radius relations generally consider H₂-He atmospheres (e.g., Rogers et al. 2011; Fortney et al. 2007) and H₂O as liquid and high pressure ice (e.g., Valencia et al. 2007; Seager et al. 2007). However, compositional diversity of atmospheres and ices in planets may exceed the one found in our solar system (e.g., Newsom 1995). But the sparse observational data on exoplanets will only allow us to constrain few structural and compositional parameters. In that light, we assume a general planetary structure of a pure iron core, a silicate mantle, a water layer, and an atmosphere made of H, He, C, and O. Given the data, we constrain core size, mantle thickness and composition, mass of water, and key characteristics of the atmosphere (mass fraction, luminosity, and metallicity). Here, we use the term atmosphere synonymously with gas layer, for which there is commonly a distinction between a convective (i.e., envelope) and a radiative part (i.e., atmosphere).

In order to compute the structure for a planet of a given interior, we employ model I from Dorn et al. (2017). Namely, we employ state-of-the-art structural models that compute irradiated atmospheres and self-consistent thermodynamics of core, mantle, high-pressure ice and ocean.

This paper is organized as follows: in Sect. 2 we review previous interior studies on the planets of interest. In Sect. 3 we introduce the inference strategy, model parameters, data and the physical model that links data and model parameters. In Sect. 4, we apply our method on the six selected exoplanets and test different stellar abundance proxies for constraints on planet bulk abundance. We end with discussion and conclusions.

Fig. 1 Masses and radii of the six exoplanets of interest, including two scenarios for 55Cnc e (see main text). The terrestrial solar system planets are also plotted against idealized mass-radius curves of pure iron, pure magnesium silicate, and a 2-layered composition of 5% magnesium silicate and 95% water.

2. Previous studies of the selected exoplanets

Previous estimates on the interiors of HD 219134b, Kepler-10b, Kepler-93b, CoRoT-7b, 55 Cnc e, and HD 97658b generally do not comprehensively reflect the inherent degeneracy that different interiors can explain the observations. Often a few interior scenarios are discussed but it remains unclear how well specific models compare with the generally large number of possible scenarios. Commonly used mass-radius curves of idealized compositions plotted against the planetary data (Fig. 1) allow a general evaluation whether a planet can either be composed of pure rocks or is volatile-rich.

For Kepler-10b, Batalha et al. (2011) estimate a predominantly rocky composition of iron and silicates, with a mass fraction of water that can be as high as 25%. With the additional constraint of maximum mantle stripping during planet formation (Marcus et al. 2009), they predict a dry iron-rich interior. For both Kepler-10b and CoRoT-7b, Wagner et al. (2012) argue for a large iron core similar to Mercury (~60% iron and ~30% silicates), whereas Dressing et al. (2015), Weiss et al. (2015) suggest an Earth-like interior (17% iron and 83% silicates). The latter conclusion also holds for Kepler-93b (Dressing et al. 2015). For Kepler-10b and CoRoT-7b it remains unclear if they may or may not be remnants of larger volatile-rich planets that have lost their atmospheres (e.g., Kurokawa et al. 2013; Valencia et al. 2010).

We find that HD 219134b is very similar in size and mass to CoRoT-7b. Compared to CoRoT-7b, HD 219134b has higher precision on mass and radius and it receives less energy from its host star. Given a bulk density of 5.76 g/cm³, Motalebi et al. (2015) suggest a rocky composition and describe a best fit two-component model of iron-magnesium silicate mantle (77%) and a core mass fraction of approximately 23%.

The bulk density of 55 Cnc e (R = 2.0R_⊕) is too low to be made of pure rocks and needs to contain significant mass of water, gas, or other light elements (Winn et al. 2011), although Gillon et al. (2012) argue that any kind of atmosphere is unlikely due to the intense stellar irradiation. Previous estimates of water mass fractions that can explain the data are e.g., ~45% (Endl et al. 2012) or ~20% (Gillon et al. 2012). We note that both estimates are based on differently assumed bulk densities and rock/silicate partitions. Recent infrared monitoring revealed a high surface temperature difference between day and night side which implies either an optically thick atmosphere with heat recirculation confined to the planetary dayside or no atmosphere (Demory et al. 2016a). Also, recent spectroscopic observations have suggested a hydrogen-rich atmosphere (Tsiaras et al. 2016). For 55 Cnc e, different radii can be found in the literature (Demory et al. 2016b, and references therein). Here, we simply distinguish two cases 55 Cnc e^L (R = 2.08R_⊕) and 55 Cnc e^S (R = 1.75R_⊕), where the superscripts refer to small (S) and large (L), both of which can be seen in Fig. 1.

The least dense and irradiated exoplanet in this study, HD 97658b, can be explained by a rocky core of at least 60% in mass with an envelope of lighter elements and little H-He components, at most 2% by mass (Van Grootel et al. 2014). But its density is also consistent with a water-dominated planet (Dragomir et al. 2013).

3. Methodology

We use the methodology of Dorn et al. (2017, model I) that determines possible ranges of interiors which can include thin and thick atmospheres. We outline the main aspects of the method in the following.

3.1. Bayesian inference

We employ a Bayesian method to compute the posterior probability density function (pdf) of the model parameter conditional on the data and prior information. The prior probability of the model represents the information known about the interior structure before using data. According to Bayes’ theorem, the posterior distribution for a fixed model parametrization, conditional on data, is proportional to prior information on model parameters and the likelihood function. The likelihood function summarizes the distance (typically a norm of a vector of residuals) between the model simulation and observed data. It is practically impossible to analytically derive the posterior distribution since the addressed problem is high-dimensional and non-linear and data are sparse. We therefore employ a McMC sampling method that iteratively searches the prior parameter space and locally evaluates the likelihood function. Samples of model parameters that are generated with this approach are distributed according to the underlying posterior distribution. As in Dorn et al. (2017), we use the Metropolis-Hastings algorithm to efficiently explore the posterior distribution. We refer to Dorn et al. (2015) and Dorn et al. (2017) for more details.

The large number of models that is needed for the analysis requires very efficient computations. In our work, generating the planet’s internal structure takes on average 40−90 s of CPU time on a four quad-core AMD Opteron 8380 CPU node and 32 GB of RAM. In all, we sampled about 10⁶ models and retained around 5 × 10⁴ models for further analysis.

3.2. Bulk abundance constraints

Theoretical and empirical evidences suggest that relative abundances of refractory elements are correlated between star and planet. For example, there is a high similarity among Fe/Si and Mg/Si for Sun, Earth, Mars, the Moon as well as meteorites (Lodders 2003; Drake & Righter 2002; McDonough & Sun 1995; Khan & Connolly 2008; Kuskov & Kronrod 2001), which are believed to represent building blocks of planets (Morgan & Anders 1980). Planet formation studies have demonstrated that there are direct correlations among relative refractory element abundances (e.g., Fe, Mg, and Si) between star and hosted planets (e.g., Bond et al. 2010; Elser et al. 2012; Johnson et al. 2012; Thiabaud et al. 2015). This is because refractory elements Fe, Si, and Mg condense at similar temperatures corresponding to small distances from the host star (≳0.06 AU).

Thiabaud et al. (2015) verified the direct correlation for a wide range of stellar compositions. For their tested solar-like disk composition (e.g., Fe/Si_Sun = 1.69, Mg/Si_Sun = 0.89), the full distribution of planetary bulk Fe/Si is shown in Fig. 2. For the subset of planets with semi-major axes below 0.1 AU, the distribution is very weakly bimodal with 8% of all planets having a very low iron content (Fe/Si <0.001). For this particular solar-like disk composition, the authors showed that at close distances to the star (0–0.2 AU), the disk’s Fe/Si (10^-6–10^-4) and Mg/Si (0.08–4) mass ratios can significantly differ from solar values (see also Fig. 3 in Thiabaud et al. 2015). This is because the condensation conditions differ among Mg-Si-Fe-bearing species. This seems to indicate that a small fraction of planets that accrete material from respective regions may have very low iron content. Here, we test the sensitivity of our interior predictions on such extremely different stellar abundance proxies for the planetary bulk composition (see Sect. 4.3):

• direct abundance proxy:

  Fe/Si_bulk = Fe/Si_star

  Mg/Si_bulk = Mg/Si_star;

• non-direct abundance proxy:

  Fe/Si_bulk = 10^-5 Fe/Si_star

  Mg/Si_bulk = 0.5 Mg/Si_star.

We note that Thiabaud et al. (2015) does neglect the drift of planetesimals, which may lead to an underestimation of bulk Fe/Si in their simulated planets. Future studies need to improve our understanding of the evolution of element ratios.

Fig. 2 Top: histogram of bulk Fe/Si of ~15 400 simulated planets that accreted from discs of solar composition. Bottom: planet mass over bulk Fe/Si of the same simulations (10 planet embryos were used in each simulation). Note the log-scale. The distribution of planets is a result from Thiabaud et al. (2015).

3.3. Stellar abundance data

We have culled stellar abundances from multiple datasets for HD 219134, 55 Cnc e, and HD 97658b and homogenized them via the Hypatia Catalog (Hinkel et al. 2014). In order to properly compare the datasets with each other, we have renormalized the abundances to the Lodders et al. (2009) solar abundance scale. To maintain coherence between abundance measurements, our analysis did not take into consideration groups who utilized NLTE (non-local thermodynamic equilibrium) approximations. The different stellar abundance estimates of 55 Cnc e, HD 97658, and HD 219134 are listed in Table 1. For Kepler-10 b, Kepler-93 b, and CoRoT-7 b, the stellar abundances have only been determined by one group, that is Santos et al. (2015).

Note, ratios with squared brackets ([X/H], e.g., Table 1) are logarithmic particle ratios in reference to the Sun. The majority of our discussion is with respect to abundance mass ratios without reference to the Sun (X/H, e.g., Table 3).

The abundances for HD 219134, which is located 6.55 pc from the Sun, were measured by a total of 9 groups within the literature. The maximum [Fe/H] determination was by Mishenina et al. (2015) with 0.84 dex, while Ramírez et al. (2007) measured a minimum of 0.04 dex and the median was 0.13 dex. Only three groups measured [Mg/H], where we use the median value of 0.32 dex and spread of 0.28 dex. The [Si/H] ratio spanned 0.04 dex per Valenti & Fischer (2005) to 0.27 dex (Thévenin 1998; Thévenin & Idiart 1999), with a median of 0.12 dex between the five contributing groups. The sodium, aluminum, and calcium abundance ratios were measured by three groups or less. The median for [Na/H] is 0.19 dex and [Al/H] is 0.23 dex, both with a spread of 0.13 dex. The [Ca/H] measurements had a smaller spread with 0.07 dex, while the median was 0.21 dex.

A total of 18 groups measured element abundances for 55 Cnc, which is 12.34 pc from the Sun. The [Fe/H] measurements ranged from 0.29 dex (Mishenina et al. 2004) to 0.67 dex (Takeda 2007), with a median of 0.44 dex. Of the 18 groups, six of them measured [Mg/H], where the minimum determination was 0.34 dex (Delgado Mena et al. 2010), the maximum was 0.80 dex (Allende Prieto et al. 2004), and the median was 0.41 dex. The [Si/H] ratio was determined by 11 measurements and spanned 0.29 dex (Gonzalez et al. 2001) to 0.57 dex (Allende Prieto et al. 2004), with a median of 0.40 dex. The sodium abundance ratio varied from 0.23 dex (Gonzalez et al. 2001) to 0.63 dex per Shi et al. (2004) out of eight datasets, where the median was 0.49 dex. The [Al/H] determinations were only measured by five groups who ranged from 0.43 dex (Kang et al. 2011) to 0.55 dex (Luck & Heiter 2005) and a median of 0.51 dex. Finally, [Ca/H] was measured by 7 groups, with a minimum of 0.09 dex (Zhao et al. 2002) to a maximum of 0.46 dex (Luck & Heiter 2005); the median is 0.24 dex.

HD 97658 is a star 21.11pc away from the Sun and was spectroscopically observed by six different groups. The variation in [Fe/H] goes from −0.05 (Petigura et al. 2011) to −0.28 dex (Mishenina et al. 2008) with a median value of −0.13 dex. The [Mg/H] ratio was only measured by two groups: Mishenina et al. (2008) and Mortier et al. (2013) who determined −0.18 dex and −0.27 dex, respectively. Three groups measured [Si/H] in HD 97658 and had a variation spanning −0.20 dex (Valenti & Fische 2005) to −0.28 dex (Mortier et al. 2013), with a median of −0.21 dex. Both Mortier et al. (2013) and Valenti & Fischer (2005) determined [Na/H], citing −0.24 dex and −0.26 dex, respectively. Finally, only Mortier et al. (2013) measured [Al/H] = −0.23 dex and [Ca/H] = −0.25 dex.

For the host stars HD 219134, 55 Cnc, and HD 97658, we adopt three scenarios for their abundances: a case representing the median estimate from different studies (V0), an iron-rich (V1) and an iron-poor (V2) case. Whenever unspecified, the reference scenario V0 is used. For the three cases, we compile stellar abundance estimates in Table 2 based on Table 1. These stellar abundance estimates serve as proxies for the planet bulk abundance (see Table 3, Cols. 3–7). In Sect. 4.2 we examine how the spread of bulk abundance constraints influence our predictions on the planetary interiors of HD 219134b, 55 Cnc e, and HD 97658b.

3.4. Data

The data dthat we rely on to constrain planetary interiors (see Sect. 3.1) are:

• planetary mass M;

• planetary radius R;

• bulk planetary ratio Fe/Si_bulk;

• bulk planetary ratio Mg/Si_bulk;

• mantle composition of minor elements: CaO, Al₂O₃, Na₂O;

• semi-major axis;

• stellar radius R_star;

• stellar effective temperature T_star.

Table 3 lists all data for the six selected exoplanets. Figure 1 shows their masses and radii plotted against M-R-curves of idealized compositions.

The Fe/Si_bulk mass ratio relates the mass of iron to silicon for the entire planet (core and mantle). Since all magnesium and silicate are in the mantle, Mg/Si_bulk equals Mg/Si_mantle. As stated earlier, we are using two stellar abundance proxies, that is direct and non-direct proxy (Sect. 3.2). Following Dorn et al. (2015), we fix absolute abundances of minor refractories in the mantle (oxides of Ca, Na, and Al) to the stellar proxy. The system of Na₂O–CaO–FeO–MgO–Al₂O₃–SiO₂ is chosen, because it is able to explain 99% of Earth’s mantle. Whenever stellar proxies of these minor compounds are not available, we assume chondritic abundances from Lodders et al. (2009). Semi-major axis, stellar radius, and stellar effective temperature are also fixed parameters.

3.5. Model parametrization

Our exoplanet interior model consists of an iron core surrounded by a silicate mantle, a water layer and a gas layer as illustrated in Fig. 3. The key structural model parameters mthat we aim to constrain are:

• core radius r_c;

• mantle Fe/Si_mantle;

• mantle Mg/Si_mantle;

• mantle radius r_mantle;

• mass of water m_water;

• mass of atmosphere m_env;

• atmosphere Luminosity L;

• atmosphere metallicity Z_env,

where Z_envis defined as the atmosphere mass fraction of elements heavier than H and He (here C and O).

Fig. 3 Illustration of model parametrization. Model parameters are core radius rc, mantle composition c comprising the oxides Na2O–CaO–FeO–MgO–Al2O3–SiO2, mantle radius rmantle, mass of water mwater, mass of atmosphere menv, atmosphere luminosity L, and atmosphere metallicity Zenv.

3.6. Structure model

Data dand model parameters mare linked by a physical model embodied by the forward operator g(·): ${d}\mathrm{=}\mathit{g}\mathrm{\left(}{m}\mathrm{\right)}\mathit{.}$(1)For a given model m, the interior structure (density profile), total mass M, transit radius R, and Fe/Si_bulk are computed for the purpose of comparing with observed data d. The function g(m) is computed combining thermodynamic and Equation-of-State (EoS) modeling which is adopted from Dorn et al. (2017). In the following, we briefly summarize the detailed description given in Dorn et al. (2017, model I).

Iron core.

For the core density profile, we use the EoS fit of iron in the hcp (hexagonal close-packed) structure provided by Bouchet et al. (2013) on ab initio molecular dynamics simulations for pressures up to 1500 GPa. This EoS fit is optimal for the modeling of Earth-like exoplanets up to ten Earth masses. The temperature gradient is assumed to be adiabatic.

Silicate mantle.

We employ a self-consistent thermodynamic method within the Na₂O-CaO-FeO-MgO-Al₂O₃-SiO₂ model system. We assume thermodynamic equilibrium. Equilibrium mineralogy and density are computed as a function of pressure, temperature, and bulk composition by Gibbs energy minimization (Connolly 2009). For these calculations the pressure is obtained by integrating the load from the surface boundary, while temperature is obtained by integrating an Earth-like temperature gradient from the surface temperature of the mantle. This surface temperature equals the temperature at the bottom of the water layer.

Water layer.

In order to compute the density profile of the water layer, we follow Vazan et al. (2013) using a quotidian equation of state (QEOS). This QEOS is in good agreement with widely used ANOES (Thompson & Lauson 1972) and SESAME EoS (Lyon & Johnson 1992). Above 44.3 GPa, we use the tabulated EoS from Seager et al. (2007) that is derived from DFT simulations and predict a gradual transformation from ice VIII to X. We assume an adiabatic temperature profile.

Gas envelope and atmosphere model.

For a given radius and mass of the solid interior, distance to the star a, stellar effective temperature T_star, stellar radius R_star, envelope luminosity L, gas metallicity Z_env, and gas mass m_env, we solve the equations of hydrostatic equilibrium, mass conservation, and energy transport. As in Venturini et al. (2015), we implement the CEA (Chemical Equilibrium with Applications) package (Gordon & McBride 1994) for the EoS, which performs chemical equilibrium calculations for an arbitrary gaseous mixture, including dissociation and ionisation, and assuming ideal gas behavior. We account for gas compositions of H, He, C and O which are key atmospheric elements.

An irradiated atmosphere is assumed at the top of the gaseous envelope, for which the analytic irradiation model of Guillot et al. (2010) is adopted. This irradiation model assumes a semi-gray, global temperature averaged profile. The temperature is defined by envelope luminosity L and stellar irradiation. For the planets of interest, internal temperatures generally range between 1000 and 3000 K, intrinsic heat fluxes from the interior lie between 10 and 350 K. Calibrated opacities for different equilibrium temperatures and solar metallicities from Jin et al. (2014) are used. When non-solar metallicities are considered (Z_env≠ 0.02), computed opacities are thus not self-consistent (see discussion in Sect. 5).

The boundary between the irradiated atmosphere and the envelope is defined where the optical depth in visible wavelengths is much larger than 1 (Jin et al. 2014). Within the envelope, the usual Schwarzschild criterion is used to distinguish between convective and radiative layers. That is, if the adiabatic temperature gradient is larger than the radiative one, the layer is stable against convection. In order to compute the planet radius of a model realization we evaluate where the chord optical depth becomes 2/3 (Guillot et al. 2010). This allows us to compare calculated planet radii with measured radii from primary transits.

3.7. Prior information

Table 4 lists the ranges of prior information. Reasonable choices for the prior distributions are crucial such that we do not arbitrarily favor certain model realizations without justification. We have adopted the prior ranges from Dorn et al. (2017).

Prior bounds on Fe/Si_mantle and Mg/Si_mantle are linked to the host star’s photospheric abundances. Since all Si and Mg are assumed to be in the mantle, Mg/Si_star defines the prior bounds on Mg/Si_mantle, while Mg/Si_star is Gaussian-distributed. The iron content, on the other hand, is distributed between core and mantle. Thus, the bulk constraint Fe/Si_bulk defines only the upper bound of the prior on Fe/Si_mantle. There is also a numerical constraint that the absolute iron oxide abundance in the mantle cannot exceed 70%. For m_envand L, we assume a Jeffrey’s prior, i.e., the logarithm of the parameter is uniformly distributed.

Fig. 4 Sampled two-dimensional (2D) marginal posterior pdfs for the six selected exoplanets (a)–f)) showing the correlation between core and mantle sizes, rc, and rmantleand mantle Fe/Simantle. Bigger dots explain the data within 1σ uncertainty. The straight line describes the lower limit of the mantle size where rmantle= rc. Earth-like parameters are depicted by the red cross.

Fig. 5 Sampled 2D marginal posterior pdfs for the six selected exoplanets (a)–f)) showing the correlation between mantle size rmantleand mass of water mwater. Blue points explain the data within 1σ uncertainty. Earth-like parameters are depicted by the red cross.

Fig. 6 Sampled 2D marginal posterior pdfs for the six selected exoplanets (a)–f)) showing the correlation between core size rcand mass of water mwater. Blue points explain the data within 1σ uncertainty. Earth-like parameters are depicted by the red cross.

Fig. 7 Sampled 2D marginal posterior pdfs for the six selected exoplanets (a)–f)) showing the correlation between gas metallicity Zenvand atmospheric radius fraction renv/R. Blue points explain the data within 1σ uncertainty. Earth-like parameters are depicted by the red cross.

4. Results

We apply our method to the six selected exoplanets with the data listed in Table 3 and first assume the direct abundance proxy (Fe/Si_bulk = Fe/Si_star and Mg/Si_bulk = Mg/Si_star) and best guess stellar abundance estimates. We calculated the posterior distribution of model parameters and analyzed 2-dimensional (2D) correlations between all 8 model parameters. In the following, we highlight a few (2D) parameter correlations: r_mantle, r_c, and Fe/Si_mantle (Fig. 4), r_mantleand m_water(Fig. 5); r_cand m_water(Fig. 6), as well as Z_envand atmospheric radius fraction r_env/R (Fig. 7).

The correlation between mantle composition and core size (Fig. 4) is introduced by the bulk abundance constraint Fe/Si_bulk. For planets for which the uncertainty on Fe/Si_star is lowest (below 34%), namely Kepler-10b, Kepler-93b, and CoRoT-7b, an increase in core-size leads to lower mantle Fe/Si_mantle in order to fit bulk abundance constraints. This is illustrated by the horizontal shift in color (Figs. 4b–d) from yellow to blue. Whenever Fe/Si_star has a large uncertainty (e.g., about 90%, 80%, and 60% for HD 219134b, 55 Cnc, and HD 97658, respectively), the negative correlation between mantle size and mantle Fe/Si_mantle dominates. In these cases the color shifts from yellow to blue vertically. This is because the dominating datum is mass. More specifically, the core size is less well constrained, because for example a decrease in mantle size can be compensated by increasing Fe/Si_mantle in order to fit mass without affecting r_c. For HD 97658b, actually both trends are visible, yellow roughly shifts along a diagonal.

The correlation between r_mantleand m_wateris generally strong as depicted in Fig. 5. As expected from comparisons to idealized mass-radius-curves (see Fig. 1), it is possible to explain the data for HD 219134b, Kepler-10b, Kepler-93b, CoRoT-7b, and 55 Cnc e^S with a rocky interior and no significant volatile layer. The full posterior range of the water mass fraction, however, reveals that significant volatile compounds cannot be excluded by the data. The possible mass fraction of water, m_water, spans from 0–35% (0–17% to explain the data within 1σ uncertainty) depending on the planet. Thus, their probability to be Earth-like is low, also since posterior distributions are generally large. Only for Kepler-10b the range of possible water mass fractions is restricted to maximum few percents. For the least dense exoplanets 55 Cnc e^L and HD 97658b, m_watercan be as large as 0.7 M and 0.95 M, respectively (0.3 M and 0.7 M, respectively, to fit the data within 1σ uncertainty, blue points in Fig. 5).

Also, the (negative) correlation between mantle size r_mantleand m_waterfor all planets (Fig. 5) is much stronger than between core size r_cand m_water(Fig. 6). More specifically, this is because the parameters that most strongly affect planet mass are r_mantleand m_water. An increase in r_mantlecan be compensated by decreasing m_water, while still fitting planet mass. However, relative core size is mostly determined by Fe/Si_bulk.

Figure 7 indicates that the correlation between gas metallicity and atmospheric radius fraction is generally weakly negative. The upper bound on m_envpositively correlates with Z_env(not shown). This is because the same atmospheric thickness of a more metal-rich gas can only be obtained by increasing atmospheric mass. Among the six exoplanets, the upper bound on m_envdiffers depending on bulk density and data uncertainty. Thus the five potentially rocky exoplanets (Figs. 7a–e) can have an atmosphere of maximum mass of 10^-5–10^-7M_⊕for low metallicities (Z_env<0.2). For high metallicities (Z_env>0.8), the maximum mass can be four magnitudes higher. In comparison the Earth’s atmosphere has similar mass fraction of about 9 × 10^-7M_⊕. For the volatile-rich exoplanets (Figs. 7e–f), the correlation between the atmospheric radius fraction r_env/R and Z_envis strongest and the correlation between the upper bound on m_envand Z_envseems to be negative. Whether the planets could retain the constrained atmospheres with regard to possible photo-evaporative mass loss is not addressed in our model, but needs some post-analysis consideration (Sect. 4.1).

Overall, correlations are generally weak since the problem addressed here is highly degenerate. In other words, an arbitrary change in 2 model parameters can be compensated by the remaining 6 model parameters in order to fit data. We find significant correlations between those model parameters that most strongly influence mass and/or radius and that countervail each other regarding this datum (e.g., r_mantleand m_waterboth most strongly affect planet mass).

Finally, the predicted ranges for structural parameters strongly depend on the specific planet data and associated uncertainties. Even though we see a distinction between the potentially rocky (HD 219134b, Kepler-10b, Kepler-93b, CoRoT-7b, 55 Cnc e^S) and volatile-rich planet types (55 Cnc e^L and HD 97658b), the distributions can vary significantly among planets of the same type. Among all studied planets, Kepler-10b seems to be best constrained and its interior is most similar to an Earth-like structure. It is our contention that a case-by-case analysis is essential for a thorough characterization of planet interiors. A comparison of previous results (Sect. 2) with our analysis shows very good agreement.

4.1. Evaporative mass loss

Highly irradiated planets can lose large amounts of their atmospheres through hydrodynamic mass loss driven by the extreme ultraviolet and X-ray heating from their host stars. In order to investigate whether our data-constrained interior structures are stable in the sense that atmospheres are stable against photo-evaporation, we use simplistic analytical relationships by Lopez & Fortney (2013) and Jin et al. (2014). A thorough analysis involving state-of-the-art mass loss models that couple hydrodynamic evaporation and thermal evolution models (e.g. Valencia et al. 2010) lies outside the scope of this study. Both simplistic scaling relations were fit against a large number of model runs and qualitatively describe the importance of mass loss but do not make detailed predictions. Lopez & Fortney (2013) provide a fit for threshold fluxes at which significant atmospheric loss occurs as a function of core mass and mass-loss efficiency for Earth-like, rocky interiors and H/He atmospheres (of up to 60% in mass fraction). Similarly, Jin et al. (2014) provide a simplistic mass-radius-threshold below which planets have likely lost their entire H/He atmosphere. The application of both relationships to our posterior realizations suggest that all studied planets experience significant mass loss, except for HD 97658b (and mass-loss efficiencies below 0.05 regarding the Lopez & Fortney 2013, relationship). These findings are also in agreement with previous studies, e.g., namely with the suggestions of Valencia et al. (2010) that CoRoT-7b would likely not possess a significant atmosphere as well as Demory et al. (2016a) that 55 Cnc e can be devoid of atmosphere.

In general, surface temperatures of the studied planets are high (~2000 K, depending on structure and planet case) which may imply molten rocks at the surface and outgassing from the planetary interior. Outgassing is a likely source for secondary atmospheres (e.g., Schaefer & Fegley 2007). However, due to the high irradiation, significant mass loss may remove most of such atmospheres. Detailed models are needed to determine how much vapor could remain in equilibrium with possible rates of outgassing and evaporative mass loss.

4.2. Sensitivity of structural parameters to bulk abundance constraints

Here, we examine the sensitivity of interior predictions on the spread of planetary bulk abundance constraints. We do this by testing how interior predictions for HD 219134b, 55 Cnc e, and HD 97658b vary in light of different bulk abundance estimates taken from: (I) best guess stellar abundance estimates, (II) iron-rich (V1), and (III) iron-poor (V2) abundance estimates (Table 3). Results are shown in Fig. 8. The blue curves for the best guess abundance estimate refer to the same posterior pdfs previously shown in Figs. 4–7.

The spread in interior parameters is largest for HD 219134b and smallest for HD 97658b. This is because the difference in bulk abundance constraints between V1 and V2 is a factor of 10 for HD 219134b, whereas for 55 Cnc e and HD 97658b the difference is only a factor of 4.5 and two, respectively. Except for Fe/Si_mantle and r_c, differences between the cumulative distribution functions (cdfs) are generally minor. Compared to V2, the 50th-percentile in r_cand Fe/Si_mantle increases by up to 20% and 200% for V1, respectively, for r_mantleand m_waterthese variations are on the order of a few percents. Differences in m_envare minor.

The differences in core size and mantle composition introduced by the use of different stellar abundance estimates are even more pronounced when reduced uncertainties for mass and radius are available. This will be addressed in future parameter studies.

Fig. 8 Application of V0, V1, and V2 stellar abundance estimates: cdfs of sampled 1D marginal posterior (green to blue) and prior (red) for HD 219134b, 55 Cnc eS, 55 Cnc eL, and HD 97658b: a) mwater, c) rmantle, d) rc, and e) Fe/Simantle. Posterior cdfs refer to the different host star abundance estimates (V1 – iron rich, V2 – iron poor) listed in Table 3. The dashed lines are the 95% lower and upper confidence bounds, which almost completely overlap the solid lines. Depending on the stellar abundance proxy estimates, the upper prior bound in e) differs, which is indicated by the vertical colored lines corresponding to the respective proxy.

Fig. 9 Application of direct and non-direct abundance proxies: cdfs of 1D marginal posterior pdfs and prior pdfs for HD 219134b, Kepler-10b, Kepler-93b, CoRoT-7b, 55 Cnc eS, 55 Cnc eL, and HD 97658b: a) menv, b) Zenv, c mwater, d) rmantle, e) rc. Prior (red), posterior pdf in blue and green for direct and non-direct stellar-planetary abundance correlation, respectively. See discussion in Sect. 4.3.

4.3. Direct versus non-direct stellar abundance proxy

Here, we investigate the influence of different stellar abundance proxies by comparing retrieved interior structures obtained with the direct (Fe/Si_bulk = Fe/Si_star and Mg/Si_bulk = Mg/Si_star) and non-direct abundance proxy (Fe/Si_bulk = 10^-5 Fe/Si_star and Mg/Si_bulk = 0.5 Mg/Si_star; see Sect. 3.2) on best guess stellar abundance estimates. The non-direct abundance proxy implies an extremely low bulk iron (nearly iron-free) content of the planets. With the constraint of a nearly iron-free planet, mass and radius of Kepler-93b and Kepler-10b can not be fit within 1σ uncertainty. All other studied planets can be explained within 1σ uncertainty, because data uncertainty is simply large enough or the bulk density (e.g., 55 Cnc e^L and HD 97658b) can be explained without a large iron core.

Figure 9 illustrates the differences in parameter estimates for the direct (blue curves) and non-direct (green curves) abundance proxy. Again, the blue curve refer to the same posterior pdfs previously shown in Figs. 4–7, and Fig. 8 (blue curve).

The obtained structural parameters differ significantly between the two tested abundance proxies. As expected, using the non-direct abundance proxy compared to the direct abundance proxy leads to extremely low iron content in the mantle and tiny iron cores. In order to still fit mass, the size of rocky interior is larger. The 50th-percentiles of the parameter cdfs in Fig. 9 change as follows for the extreme sub-stellar bulk abundance constraints:

• r_c decreases by 78–83%,

• r_mantle increases by 2–13%,

• Fe/Si_mantle decreases by 100%,

• m_water decreases by 9–50%.

Since structural variability in core and mantle is strongly restricted by the extremely low bulk abundance constraint, core and mantle parameters are generally better constrained (i.e., cdf curves are steeper).

With our methodology and the given data, such extreme sub-stellar abundances cannot be excluded for potentially rocky exoplanets, although extremely low iron content is disfavored for Kepler-93b and Kepler-10b. It is important to note that such low iron contents are very unlikely for volatile-rich planets, because accreted volatiles originate from primordial-disk regions beyond the ice-line at which also Mg-Fe-Si-bearing species condensed and likely reached the stellar ratio (Thiabaud et al. 2015).

In order to quantify how much the direct abundance proxy is preferred over the non-direct abundance proxy, we use the posterior odds, K, that is a multiplication of Bayes factor and prior odds (Jeffreys 1998): $\mathit{K}\mathrm{=}\frac{\mathit{p}\mathrm{\left(}{d}\mathrm{\right|}{\mathit{H}}_{\mathrm{d}}\mathrm{\left)}\mathit{p}\mathrm{\right(}{\mathit{H}}_{\mathrm{d}}\mathrm{)}}{\mathit{p}\mathrm{\left(}{d}\mathrm{\right|}{\mathit{H}}_{\mathrm{n}}\mathrm{\left)}\mathit{p}\mathrm{\right(}{\mathit{H}}_{\mathrm{n}}\mathrm{)}}\mathrm{·}$(2)Here, H_d and H_n denote the hypothesis of validity of the direct and non-direct abundance proxy, respectively, and the probability $\mathit{p}\mathrm{\left(}{d}\mathrm{\right|}{\mathit{H}}_{\mathrm{d}}\mathrm{)}\mathrm{=}\mathrm{\int }\mathit{p}\mathrm{\left(}{m}\mathrm{\right|}{\mathit{H}}_{\mathrm{d}}\mathrm{\left)}\mathit{p}\mathrm{\right(}{d}\mathrm{\left|}{m}\mathit{,}{\mathit{H}}_{\mathrm{d}}\mathrm{\right)}\mathrm{d}{m}\mathit{,}$(3)represents how well some data are explained under the assumption of a hypothesis H_d. Instead of the Bayes factor integral we use the maximum likelihood estimate. We obtain K-values around 10¹, except for Kepler-93b and Kepler-10b with K-values of 10². This means that only for Kepler-93b and Kepler-10b, there is a strong evidence that the direct abundance proxy better explains the planetary mass and radius.

Generally, it will be very difficult to study the validity of a direct correlation of relative bulk abundances between star and hosted exoplanet given usual data uncertainties. Further studies on solar system objects are needed to improve our understanding of such abundance correlations and their dependence on planet formation history.

5. Discussion of the assumptions

The success of an inference method is subject to the limitations and assumptions of the forward function g(·). Firstly, we consider pure iron cores, which may lead to a systematic overestimation of core density and thus an underestimation of core size, if volatile compounds in the core are significant. The trade-off between predicted and independently inferred core sizes can be on the order of few to tens of percents, as shown for the terrestrial solar system planets in Dorn et al. (2015).

Secondly, our approach assumes sub-solidus conditions for core and mantle as well as a perfectly known EoS parameters for all considered compositions. For the six exoplanets considered here, their pressures and temperatures exceed the ranges for which these parameters can be measured in the laboratory. Clearly, extrapolations introduce uncertainty.

Thirdly, we have assumed pure water composition of the ice layer, but other compounds such as CO, CO₂, CH₄, NH₃, etc. are also possible since the additional elements are relatively abundant. We have used water as a proxy for these volatiles in general because (1) oxygen is more abundant than carbon and nitrogen in the universe, and (2) water condenses at higher temperatures than ammonia and methane.

Fourthly, our atmospheric model uses prescribed opacities of Jin et al. (2014), which is suitable for solar abundances. In this aspect, the prescription is therefore not self-consistent with the non-solar metallicity cases we present in this work. We have checked that using different infrared and/or visible opacities can introduce errors in radius of ~5%. Models that compute line-by-line opacities with their corresponding atmospheric abundances should be performed in the future to compute planetary radii in a self-consistent way. Another aspect is our assumption of ideal gas behavior. We have checked that for atmospheric mass fractions (m_env/M) of up to 1%, the difference in the radius between using the Saumon et al. (1995) EOS and ideal gas (for H-He atmospheres) can reach 10%. For the vast majority of our results, atmospheric mass fractions are smaller than 0.1%, and thus the ideal gas assumption does not introduce a remarkable flaw.

Since our interior models are static, evolutionary effects on the interior structure are not accounted for in the forward model. As demonstrated in Sect. 4.1, constraints on evolutionary processes such as thermal mass loss of gaseous layers can be roughly approximated by post-analysis considerations. We have used very approximate analytic scaling laws from Lopez & Fortney (2013) and Jin et al. (2014) to rule out gas-rich planet realizations.

6. Summary

We have presented the application of a generalized Bayesian inference method from Dorn et al. (2017) in order to constrain the general interior structure of six selected super Earths: HD 219134b, Kepler-10b, Kepler-93b, CoRoT-7b, 55 Cnc e, and HD 97658b. Their available data are mass, radius, and bulk abundance constraints on relative refractory elements (Fe, Mg, Si) measured in the photosphere of their host stars. We assume a general structure of core, mantle, ice and gas for which we constrain core size, mantle thickness and composition, mass of water, and key characteristics of the atmosphere (mass, intrinsic luminosity, composition). We have investigated the sensitivity of the constrained interior parameters on the spread of bulk abundance estimates as determined by research groups with different techniques (Hinkel et al. 2014). We have also determined how predictions of interior parameters are affected by a very non-direct correlation of relative abundances between star and planet, compared to a direct correlation. Our findings are summarized as follows:

• The available data of the six studied exoplanets enables us to constrain interior structure and to determine parameter correlations. Depending on data and data uncertainty, the ability of constraining interior structure differs. Kepler-10b appears to be best constrained, since data uncertainties are relatively small.

• In general, correlations between interior parameters are weak. We find strongest 2D correlations between the size of the rocky interior r_mantleand water mass fraction m_water.

• We determine an upper range of possible atmospheric masses that positively correlate with gas metallicity. Whether the planets can retain the constrained atmospheres with regard to the high irradiation of their host stars is evaluated with simplistic analytical relationships that determine the importance of photo-evaporative mass loss.

• The data of most of the studied exoplanets can not exclude massive amounts of water. Even for the potentially rocky exoplanets, the possible range of water mass fraction can reach up to 35%, depending on the planet. Only for Kepler-10b, a massive water mass fraction can be excluded (m_water/M< 1%).

• Discrepancies in the stellar abundance estimates that are used as proxies for planetary bulk abundance mainly affect the prediction of mantle composition and core size. Effects on the other parameters of atmosphere, water layer, and size of rocky interior are small.

• Extremely sub-stellar abundance constraints are less in agreement with mass and radius data for the potentially rocky exoplanets, but cannot be excluded without additional considerations. For instance, potentially volatile-rich exoplanets likely accreted material that originated from regions beyond the ice line, where Fe-Mg-Si bearing species reached stellar ratios.

Available mass and radius data of exoplanets cannot be used to verify whether a correlation between stellar and planetary relative abundances of Fe, Mg, and Si differs from a direct 1:1 correlation. This is because degeneracy in interior structure of exoplanets is too large. More detailed studies of terrestrial solar system planets, e.g. Mars and Venus, would greatly help to improve our understanding of the discussed abundance constraints.

In order to rigorously characterize interior structure for an exoplanet, it is crucial to perform a case-by-case inference method, since constrained interior distributions highly depend on data and data uncertainties. The range of possible exoplanet interiors is large, since the problem is highly degenerate. Commonly used mass-radius-curves of idealized compositions for comparison against observed masses and radii can only serve as a rough planet characterization. Bayesian inference methods offer a powerful way to determine exoplanet interiors sensibly.

Acknowledgments

This work was supported in part by the Swiss National Foundation under grant 15-144. It was in part carried out within the frame of the National Centre for Competence in Research PlanetS. We would like to thank James Connolly for informed discussions, Amaury Thiabaud for making available the data published in Thiabaud et al. (2015). N.R.H. would like to thank CHW3. The results reported herein benefited from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate.

References

Appendix A: Element abundance notations

Typically in astronomy, the total number of absorbing atoms is scaled with respect to hydrogen (or 10¹² atoms), which is the most common element in the universe. Geologists often use 10⁶ silicon atoms as a scalar. In addition, the abundance ratio within the star is referenced to the same element abundance ratio within the Sun, the most well measured star, such that the star can be defined to be rich, poor, or solar-like with respect to that element. The unit of dex is described in mathematical

terms as: $\mathrm{[}\mathit{X}\mathit{/}\mathit{H}\mathrm{]}\mathrm{=}\log \mathrm{(}{\mathit{N}}_{\mathrm{X}}\mathit{/}{\mathit{N}}_{\mathrm{H}}{\mathrm{)}}_{\mathrm{star}}\mathrm{-}\log \mathrm{(}{\mathit{N}}_{\mathrm{X}}\mathit{/}{\mathit{N}}_{\mathrm{H}}{\mathrm{)}}_{\mathrm{Sun}}\mathit{,}$(A.1)where N_X and N_H are the number of element X and hydrogen atoms per unit volume, respectively. When an abundance ratio is given in square brackets, it denotes that it is in reference to the Sun. The majority of our discussion is with respect to abundance mass ratios without reference to the Sun, X/H (see Table 3, Cols. 3 and 4).

All Tables

All Figures

	Fig. 1 Masses and radii of the six exoplanets of interest, including two scenarios for 55Cnc e (see main text). The terrestrial solar system planets are also plotted against idealized mass-radius curves of pure iron, pure magnesium silicate, and a 2-layered composition of 5% magnesium silicate and 95% water.
In the text

	Fig. 2 Top: histogram of bulk Fe/Si of ~15 400 simulated planets that accreted from discs of solar composition. Bottom: planet mass over bulk Fe/Si of the same simulations (10 planet embryos were used in each simulation). Note the log-scale. The distribution of planets is a result from Thiabaud et al. (2015).
In the text

	Fig. 3 Illustration of model parametrization. Model parameters are core radius rc, mantle composition c comprising the oxides Na2O–CaO–FeO–MgO–Al2O3–SiO2, mantle radius rmantle, mass of water mwater, mass of atmosphere menv, atmosphere luminosity L, and atmosphere metallicity Zenv.
In the text

	Fig. 4 Sampled two-dimensional (2D) marginal posterior pdfs for the six selected exoplanets (a)–f)) showing the correlation between core and mantle sizes, rc, and rmantleand mantle Fe/Simantle. Bigger dots explain the data within 1σ uncertainty. The straight line describes the lower limit of the mantle size where rmantle= rc. Earth-like parameters are depicted by the red cross.
In the text

	Fig. 5 Sampled 2D marginal posterior pdfs for the six selected exoplanets (a)–f)) showing the correlation between mantle size rmantleand mass of water mwater. Blue points explain the data within 1σ uncertainty. Earth-like parameters are depicted by the red cross.
In the text

	Fig. 6 Sampled 2D marginal posterior pdfs for the six selected exoplanets (a)–f)) showing the correlation between core size rcand mass of water mwater. Blue points explain the data within 1σ uncertainty. Earth-like parameters are depicted by the red cross.
In the text

	Fig. 7 Sampled 2D marginal posterior pdfs for the six selected exoplanets (a)–f)) showing the correlation between gas metallicity Zenvand atmospheric radius fraction renv/R. Blue points explain the data within 1σ uncertainty. Earth-like parameters are depicted by the red cross.
In the text

Fig. 8 Application of V0, V1, and V2 stellar abundance estimates: cdfs of sampled 1D marginal posterior (green to blue) and prior (red) for HD 219134b, 55 Cnc eS, 55 Cnc eL, and HD 97658b: a) mwater, c) rmantle, d) rc, and e) Fe/Simantle. Posterior cdfs refer to the different host star abundance estimates (V1 – iron rich, V2 – iron poor) listed in Table 3. The dashed lines are the 95% lower and upper confidence bounds, which almost completely overlap the solid lines. Depending on the stellar abundance proxy estimates, the upper prior bound in e) differs, which is indicated by the vertical colored lines corresponding to the respective proxy.

In the text

Fig. 9 Application of direct and non-direct abundance proxies: cdfs of 1D marginal posterior pdfs and prior pdfs for HD 219134b, Kepler-10b, Kepler-93b, CoRoT-7b, 55 Cnc eS, 55 Cnc eL, and HD 97658b: a) menv, b) Zenv, c mwater, d) rmantle, e) rc. Prior (red), posterior pdf in blue and green for direct and non-direct stellar-planetary abundance correlation, respectively. See discussion in Sect. 4.3.

In the text