© ESO 2020

1 Introduction

Over the past few years, the characterization of planet interiors has been the subject of extensive research. The large number and diversity of discovered exoplanets has allowed us to identify multiple planet populations. Among these, there is an increasing interest in super-Earths and sub-Neptunes, which cover the transition from terrestrial planets to gas giants and have no analog in our solar system. Major improvements in observational techniques allow for relatively precise measurements of mass and radius. The precision of the planetary radius is limited by the uncertainty of the stellar size, since the transit depth scales as $R_{\text{P}}^2/R_{*}^2$. Recently Berger et al. (2018) presented revised radii of more than 180 000 Kepler stars, leading to a remarkable improvement of the median radius precision. In some cases, space missions can perform high precision photometry and asteroseismology and can reach relative radius uncertainties of about 3% (e.g., Hatzes 2016). In addition, the current most advanced spectrographs have radial velocity precision of 1 m s⁻¹. This precision was recently improved with instruments such as ESPRESSO (Echelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations; e.g., Pepe et al. 2018 and references therein), which is expected to have an accuracy close to 10 cm s⁻¹. Therefore, a significant improvement in the mass determination is also expected, allowing us to reach a relative uncertainty better than 10%.

The masses and radii can be used to estimate the interior structure and composition of a planet. However, determining the internal structure is extremely challenging as a consequence of intrinsic degeneracy because several compositions can lead to identical mass and radius (e.g., Rogers & Seager 2010; Lopez & Fortney 2014; Dorn et al. 2015, 2017; Lozovsky et al. 2018). Furthermore, for a planet of given mass and composition, the radius depends on several aspects such as the choice of equation of state (EOS), the envelope structure (differentiated, fully mixed, or with a compositional gradient), or the temperature. This degeneracy is critical owing to the large number of free parameters needed to model the interior of an exoplanet and the few observational constraints. In order to determine how well one interior model compares with the other possible models that also fit the data and which structural parameters can be constrained, Dorn et al. (2017) presented a generalized Bayesian inference method to quantify the degeneracy and correlation of the planetary structural parameters.

In this work we explore the limitations of constraining the internal structure of super-Earths and sub-Neptunes (focusing on planets with masses up to 25 M_⊕ and radii up to 3.5 M_⊕). We use a Bayesian inference analysis together with a nested dampling technique (e.g., Skilling 2004) to discuss several aspects that affect interior characterization: observational uncertainties, location in the mass–radius (M–R) diagram, additional observational constraints as bulk abundances, masses derived from published M–R relationships for fixed radii, and the uncertainty related to model assumptions. First, we study the influence of the data uncertainty on the determination of the internal structure. In Dorn et al. (2015), they quantified the information gained by higher data precision. This study aims to show systematically the effect of data uncertainty for a wide range of masses and bulk densities (0.5− 14.7 g cm⁻³). Additional constraints are crucial to reduce the degeneracy, and the assumption that relative abundances of refractory elements (e.g., Fe/Si, Mg/Si) of a planet are similar to that of its host star has been proposed to reduce the existing degeneracy (e.g., Grasset et al. 2009; Dorn et al. 2015). Several solar system and planet formation studies have stated that there is a direct correlation between stellar and planetary relative bulk abundances (e.g., Carter et al. 2012; Lodders 2003; Drake & Righter 2002; McDonough & Sun 1995; Bond et al. 2010; Elser et al. 2012; Johnson et al. 2012; Thiabaud et al. 2015; Wang et al. 2018). We explore under what conditions interior estimates can be improved by constraints on planetary bulk abundances taken from stellar proxies.

Most of the discovered exoplanets do not have measured masses and radii, and published M–R relationships allow us to estimate the mass for a given radius and vice versa. In addition, M–R relationships describe the main properties of various classes of exoplanets. Several studies have been dedicated to the investigation of the M–R relationship of observed exoplanets. The M–R relationships are power laws of the type M = AR^B and are based on exoplanet data (e.g., Weiss et al. 2013; Weiss & Marcy 2014; Wolfgang et al. 2015; Bashi et al. 2017). Recently, Zeng et al. (2016) inferred a semiempirical M–R relationship depending on the core mass fraction, followed by a detailed forecasting model using a probabilistic M–R relationship via Markov chain Monte Carlo (MCMC; e.g., Chen & Kipping 2017). In addition, in Otegi et al. (2020) we presented an updated exoplanet catalog based on reliable, robust, and as much as possible accurate mass and radius measurements of transiting planets up to 120 M_⊕, and we inferred two new empirical M–R relationships corresponding to rocky and volatile-rich populations. We study inferred interior parameters using the mass calculated from above-mentioned M–R relationships.

Finally, we study how variations in the temperature profiles lead to radius uncertainties comparable to observational uncertainties. In this work, we assess the importance of the various sources of uncertainty related to model assumptions and study the effect of variations of the temperature profile on the planetary radius, which has not been explored yet.

2 Method

2.1 Synthetic sample

We studied the internal structure of a sample of 20 synthetic planets with different masses and radii up to 25 M_⊕ and 3.5 R_⊕ (listed in Table 1). We aim to better understand the transition between rocky and volatile-rich exoplanets. Fulton et al. (2017) found that there is a lack of planets with radii between 1.5 R_⊕ and 2 R_⊕, known as “evaporation valley” (e.g., Owen & Wu 2013; Jin et al. 2014; Lopez & Fortney 2014), suggesting a transition between the super-Earth and sub-Neptune populations. In addition, in Otegi et al. (2020) we did a careful analysis to build an exoplanet catalog as reliable as possible and we found a transition region from rocky to volatile-rich exoplanets, which corresponds to a mass in the range 5–25 M_⊕, and a radius in the range 2–3 R_⊕. The M–R range covered by the synthetic planets includes the transition between these two populations. Furthermore, more massive planets may have massive atmospheres in which the effects of electron degeneracy pressure are significant, and these are not accounted for in our atmospheric model.

The 20 synthetic planets for our sample were chosen to represent the properties of the exoplanet catalog presented in Otegi et al. (2020), in which we presented a revisited exoplanet catalog based on robust, reliable, and precise mass and radius measurements for transiting exoplanets (with σ_M∕M = 25% and σ_R∕R = 8%). This exoplanet catalog is dominated by exoplanets for which the masses have been measured through radial velocity, so it is dominated by relativelyshort-period exoplanets. Exoplanets orbiting close to their host stars are expected to have smaller atmospheres (lost through evaporation) and therefore the sample used in this work may be biased toward higher densities. Figure 1 shows histograms of the relative mass and radius error for the observed exoplanets from the NASA Exoplanet Archive¹. Currently, the distributions are peaked at ~8% in radius and ~35% in mass, butupcoming space and ground-based missions are expected to improve these values to a few percent. We used uncertainties of 3, 5, and 10% in radius and 5, 15, and 30% in mass for our synthetic planets to cover the range from the current most common uncertainties to the smallest uncertainties. Figure 2 shows the synthetic planets, the observed population exoplanets (from our revisited catalog in Otegi et al. 2020), and M–R curves for idealized compositions of iron, Earth-like, and water ice. The bulk densities of the synthetic planets cover the range from the physically motivated high limit (pure iron) to the minimum density of an observed exoplanet with a mass up to 20 M_⊕.

The irradiation has a significant effect on the interior structure determination of exoplanets with thick volatile envelopes. The exoplanets belonging to the volatile-rich population in the revisited exoplanet catalog of Otegi et al. (2020) are typically irradiated with fluxes of few hundreds times the Earth flux. We therefore used a default irradiation of 100 F_⊕ for the synthetic planets. The impact of the irradiation on the internal structure determination is studied further in Sect. 3.5.

Fig. 1 Distributions of the relative mass (left) and radius (right) uncertainties for the 115 observed exoplanets less massive than 25 M⊕ from the NASA Exoplanet Archive.

Fig. 2 M–R diagram with the synthetic planets used in the study (black dots). The black, blue, and light blue bars correspond to uncertainties of 5, 10, and 20% in mass and 3, 5, and 10% in radius, respectively. The crosses represent observed planets from the revisited catalog in Otegi et al. (2020). The synthetic planets are plotted against the composition lines of pure-iron (brown), Earth-like (light brown), and water ice (blue).

2.2 Structure model

We used the structure model presented in Dorn et al. (2017), which assumes a pure iron core, silicate mantle, pure water layer, and H–He atmosphere. This model uses self-consistent thermodynamics in the core, mantle, high-pressure ice, and water ocean. However in this study we used EOS for hexagonal close packed iron for super-Earth conditions presented in Hakim et al. (2018) for the core density profile. These are based on density functional theory results up to 137 TPa. UnlikeEarth’s core (e.g., Badro et al. 2007), we did not consider the presence of light elements in the core.

The silicate mantle is assumed to be made of oxides Na₂ O – CaO – FeO – MgO – Al₂ O₃ – SiO₂. Equilibrium mineralogy and density were computed as a function of pressure, temperature, and bulk composition by minimizing Gibbs free energy (e.g., Connolly 2009). For the water layers, we followed the approach presented in Sotin et al. (2010), which uses a temperate Birch-Murnaghan EOS including thermal corrections. The water can be in the solid, liquid, or super-critical phase depending on the pressure and temperature. The surface temperature of the water layer is set to be equal to the temperatureof the bottom of the gas layer.

For the gas layer, the equations of hydrostatic equilibrium, mass conservation, and energy transport are solved. We assumed an envelope with an elemental composition of H, He, C, and O, which are fundamental for the formation of key atmospheric molecules such as H₂, CO, CO₂, and CH₄ (e.g., Madhusudhan 2012; Lodders & Fegley 2002; Visscher & Moses 2011; Heng & Lyons 2016). We used the Chemical Equilibrium with Applications (CEA) package (e.g., Gordon et al. 1984) for the EOS, which performs chemical equilibrium calculations for an arbitrary gaseous mixture. For the energy transport, we used the model presented in Jin et al. (2014), where an irradiated atmosphere is assumed at the top of the gaseous envelope. Within the envelope, the usual Schwarzschild criterion is used to distinguish between convective and radiative layers. More details on the structural model are found in Dorn et al. (2017).

2.3 Bayesian inference based on nested sampling scheme

For the interior characterization we used the Bayesian inference analysis based on a nested sampling scheme using the PyMultiNest package (e.g., Buchner et al. 2014). Bayesian inference computes the posterior probability according to Bayes’ theorem, which states that the probability for a fixed model parameter x given a set of data d is given by $$P(x|d)=\frac{P(x)P(d|x)}{P(d)},$$(1)

where P(x) is the prior probability of x before the data is observed, P(d) the Bayesian evidence, and P(d|x) the likelihood function. The likelihood function represents the probability of observing the data given the model parameter x, and is given by $$P(d|x)=\frac{1}{{(2\pi )}^{N/2}{({\displaystyle \prod _{i=1}^N{\sigma }_i^2})}^{1/2}}\text{exp}\left(-\frac{1}{2}{\displaystyle \sum _{i=1}^N\frac{{(g_i(x)-d_i)}^2}{{\sigma }_i^2}}\right),$$(2)

where N is the number of data points, σ_i the uncertainties of the ith datum, and g(x) the operator linking the model parameters with the data, that is, d = g(x). Our posterior probability distribution cannot be derived analytically, so we used a nested sampling scheme (e.g., Skilling 2004). The aim of Monte Carlo technique is to efficiently evaluate of the Bayesian evidence, but also to produce posterior probability distributions. The main strengths of nested sampling with respect to other sampling methods are the small amount of problem-specific tuning required and high efficiency.

In short, the nested sampling scheme works as follows. The algorithm samples some number of live points randomly from the prior P(x). The likelihood P(d|x) is evaluatedat each of these points. At each iteration the point with the lowest likelihood L_min is replaced by a new point sampled from the region of prior with likelihood P(d|x) > L_min, keeping the number of live points constant. This process is continued until Bayesian evidence reaches some specified value (typically 0.5 in log-evidence). It produces a list of samples that can be used to produce marginalized posterior distributions. The sampling scheme is controlled by two main parameters: the number of live points and the maximum efficiency parameter. The number of live points has to be large enough to adequately sample the parameter space, for which we used a recommended value of 1000. The maximum efficiency controls the sampling volume at each iteration, which is equal to the sum of the volumes enclosing the active point set. We set it to 1 to obtain the maximum efficiency.

2.4 Parameters and priors

The composition of the planets are given by the masses of the iron core, mantle layer, water layer, and H–He envelope. As already mentioned we assume a pure iron core, the mantle composition is given by the Mg/Si_mantle, Fe/Si_mantle, Al/Si_mantle, Ca/Si_mantle, and Na/Si_mantle fractions. The water layer is considered to be pure H₂O and for the atmosphere the metallicity Z, luminosity L, and the irradiation temperature T_irr are used to model the atmospheric structure. If not stated otherwise we fix the Mg/Si_mantle, Fe/Si_mantle, Al/Si_mantle, Ca/Si_mantle, and Na/Si_mantle fractions, Z, and T_irr while only the layer masses are varied. The luminosity is scaled as L ∝ M^2.76, which fits the Jovian planets (e.g., LaViolette 2006).

For the nested sampling scheme we assumed a uniform prior distributions on the layer masses from 0 to the target planets mass. We would like to emphasize that for a given target planet, depending on the scientific question, a different choice of priors might be necessary. The set of prior presented here is one possible way of setting up the analysis but not the only one². The impact of different priors should be subject of further studies.

As data variables we chose the total mass and total radius of the planet, except in Sect. 3.4. where we also use the bulk Mg/Si and bulk Fe/Si fractions as data; in that case the Mg/Si_mantle and Fe/Si_mantle fractions arealso varied, assuming a uniform prior. The total mass is calculated as the sum of all layer masses. The total radius and the bulk Mg/Si and Fe/Si fractions are an outcome of the structure model. In Table 2 we summarize the model parameters, their priors, and the data variables.

3 Dependence of internal structure determination on observational uncertainties

In this section we explore how the internal structure determination depends on the observational uncertainties for exoplanets of different masses and radii fixing the irradiation to 100 F∕F_⊕ and solar Fe/Si and Mg/Si abundances. In particular, we explore how the posterior distribution of the internal parameters depend on the observational uncertainty for the synthetic planets in Table 1. An example is shown in Fig. 3, in which the posterior distributions of the internal parameters corresponding to a planet of observed mass and radius of 6 M_⊕ and 1.5 R_⊕ with different observational uncertainties are shown. The first four columns show the posterior distributions of the iron core, mantle, water, and atmospheric masses (output) renormalized by the sum, and the last column shows the mass distribution of the planet (input). The rows correspond to different observed radius uncertainties and the colors to different observed mass uncertainties; black, red, and green for 30, 15, and 5%, respectively. We find that for a planet of 6 M_⊕ and 1.5 R_⊕ a decrease in the uncertainties improves the determination of the core mass, but has a negligible effect on the determination of the other internal parameters.

In Fig. 4 the standard deviation (σ) of the posterior distributions of the internal structure parameters are shown versus the mass uncertainty (ΔM), for a radius of 1.5 R_⊕ and a mass of 6 M_⊕. We aim to use this approach to explore what interior parameters are better constrained when the mass precision is higher, keeping the radius fixed with zero uncertainty. In this work we assume a flat uniform distributed uncertainty for the planetary mass since it simplifies the analysis and allows us to consider cases with uncertainties close to 0. The equation ΔM = 0 shows the intrinsic degeneracy of the interior parameters when data uncertainty is zero for a planet of 1.5 R_⊕ and 6 M_⊕. When ΔM increases, σ of the core mass increases, meaning that this is the only tested interior parameter that is sensitive to a change in observational uncertainties. It is in agreement with what is shown in Fig. 3 where we use a more realistic Gaussian distribution for the mass and radius. It is interesting to note that the evolution of σ_M with ΔM is nearly linear. The slope δσ_M∕δΔM contains very valuable information: it indicates how much interior estimates (i.e., layer mass fractions) can be improved by increasing data precision (i.e., mass). Figure 5 shows these slopes δσ_M∕δΔM of the core, mantle, water, and atmospheric masses for all synthetic planets. The results allow us to differentiate three different regimes, described in the following subsection.

Fig. 3 Posterior distributions of core, mantle, water, atmospheric, and total masses of a planet with mass of 6 M⊕ and radius of 1.5 R⊕. The rows correspond to different radius uncertainties and colors in the subplots correspond to different mass uncertainties.

Fig. 4 Standard deviation of the posterior distributions of core, mantle, water, and atmospheric masses vs. width of the total mass distribution for a planet of R = 1.5 R⊕ and M = 6 4 M⊕. For the atmosphere, the standard deviation of the posterior distribution in logarithmic scale is used.

3.1 Planets below Earth-like composition line

We refer to planets C and D (Table 1). For these planets, a decrease of the observational uncertainties leads to a better determination of the core mass. In this density regime, planets do not have a significant gaseous envelope and they require a large amount of iron. Then, variation in the mass or radius distributions is reflected in the core mass, since there is no other interior layer that can account for such high planet bulk densities.

As mentioned above, Fig. 5 is constructed assuming flat uniform error distributions for mass and radius. To check whether the conclusions are also valid for more realistic distributions we compare the results using the Gaussian and flat distributions, where σ of the Gaussian is equal to ΔM of the flat distribution. The obtained results with the flat distributions are consistently slightly higher than with the Gaussian results, but the general shape is nearly equal. Figure 6 shows σ_core versus observational uncertainties in mass and radius using Gaussian distributions. It shows that σ does not decrease uniformly with uncertainty and that the internal structure determination is more sensitive to an improvement of data uncertainties when data uncertainties are large. Furthermore, we find that a decrease of the observed mass uncertainty is much more effective to better determine the core mass than a decrease in the radius uncertainty. This is because the core mass directly affects the planetary mass, while its effect on radius is smaller.

Fig. 5 Model parameters and synthetic planets Δσ∕ΔM (see text for further details). A higher value of Δσ∕ΔM indicates that a decrease of the observational uncertainty provides additional information of the internal parameter. The crosses represent the observed planets from Otegi et al. (2020). Composition lines of pure-iron (brown), Earth-like planet (lightbrown), and water ice (blue) are also plotted.

Fig. 6 Observational uncertainties in mass and radius versus σcore (left) and σatm (right) for the two synthetic planets with the highest and lowest densities.

3.2 Planets above pure-water composition line

We find that when a planet lies above the pure-water composition line a decrease in the mass and radius uncertainties leads to a better determination of the H–He envelope mass. The explanation is similar to that discussed above for high-density planets: for a given mass, low-density planets require thick gaseous envelopes to match the observed radius, and this significantly reduces the degeneracy. Therefore a change in the total mass/radius distributions has a direct impact on the properties of the gaseous envelope (i.e., mass fraction, thickness).

The right panel of Fig. 6 shows the σ_atm versus the observational uncertainties in mass and radius for the low-density synthetic planets with realistic Gaussian mass and radius distributions. Interestingly, in this regime decreasing the radius uncertainty is more effective than decreasing the mass uncertainty to constrain the atmospheric mass. This can be understood by the fact that the atmospheric mass fraction of the low-density planets in the synthetic sample is nearly negligible, but the radius fraction is significant.

3.3 Planets between Earth-like and pure-water composition lines

Finally, we find that planets with densities between the Earth-like and pure-water composition lines are the most degenerate. For these planets, decreasing the mass and radius uncertainties does not improve the determination of any internal structure parameter. Even when removing one compositional layer, there is a substantial degeneracy among the other three. Figure 7 shows the posterior distributions the interior parameters of a planet with fixed radius of 2.3 R_⊕ and mass of 15 M_⊕. The different colors represent different widths of the total mass distribution. The posterior distributions of the internal parameters do not coverthe same region in the two-dimensional parameter space. The posterior distributions only overlap when projecting them into a one-dimensional histogram. Therefore, for these planets it would be crucial to have additional information that can further constrain the internal structure, such as as atmospheric metallicity measurements from space missions like ARIEL (Atmospheric Remote-sensing Infrared Exoplanet Large-survey; e.g., Tinetti et al. 2016) and James Webb Space Telescope (JWST; e.g., Beichman & Greene 2018).

Fig. 7 Illustration of the degeneracy in the intermediate density regime previously defined. The posterior distributions the internal parameters of a planet with fixed radius of 2.3 R⊕ and mass of 15 M⊕ are shown The different colors represent different widths of the mass distribution.

3.4 Using stellar abundances as an additional constraint

As shown inDorn et al. (2017), abundance constraints of Fe/Si and Mg/Si from the host star can serve as a proxy for the planet bulk abundance and can reduce model degeneracy. This assumption reproduces the composition of the Earth (e.g., Javoy et al. 2010)and Mars (e.g., Khan & Connolly 2008), but not the one of Mercury, for which a post-formation giant impact scenario is often considered (e.g., Benz et al. 2007; Chau et al. 2018). Thiabaud et al. (2015) used a chemical model to link the composition of 18 synthetic stars with solar mass and luminosity with the composition of the hosted planets, and they found a close relation. Nevertheless, Wang et al. (2018) showed that there are some differences for the Earth compared to the Sun, and that volatilization trends in the bulk composition of exoplanets should be considered (e.g., Dorn et al. 2019).

We next explore how the results are affected when assuming that the stellar relative abundances reflect the planetary relative abundances.We then use the bulk Fe/Si and Mg/Si as additional constraints, and introduce two additional model parameters: the Fe/Si and Mg/Si ratios in the mantle. The goal of this section is to evaluate whether using the bulk abundances allows us to improve the estimation of the internal structure using a more detailed interior model. We find that using the bulk Fe/Si and Mg/Si does not necessarily further constrain the planetary internal structure when adding Fe/Si and Mg/Si ratios in the mantle as free parameters. We assume a small uncertainty on the bulk Fe/Si and Mg/Si of 3%. Figure 8 shows the standard deviation of the core mass and mantle mass for a planet of the intermediate-density regime in terms of the width of the total mass distribution. In this case σ_max −σ is displayed, so that a higher value means that the interior parameter is better constrained with respect to the reference, at which σ = σ_max. In the case of the chosen synthetic planets, when we introduce the constraint of Fe/Si_planet = Fe/Si_star with Fe/Si_star < Fe/Si_⊙, their core mass fraction is better constrained. The lower Fe/Si_star, the better the core mass gets constrained. Regarding their mantle mass fraction, it is only better constrained when Fe/Si_star < 1∕2 Fe/Si_⊙. Considering the bulk stellar abundances mostly affects the core mass fraction determination. For Fe/Si_star >2 Fe/Si_⊙ the abundance constraints do not provide additional information on the internal parameters for the tested planets.

Fig. 8 Standard deviation of the core mass (left) and mantle mass (right) for a planet of mass 15 M⊕ and radius 2.3 R⊕ vs. the width of the mass distribution. The red, green, and blue lines represent different Fe/Si in the planet. A higher value of σmax − σ indicates that the internal parameter is better constrained with respect to the reference, at which σ = σmax.

Fig. 9 Atmospheric mass Δσ∕ΔM for the synthetic planets vs. the bulk density (left) and median atmospheric thickness (right) for different irradiation. The light green and red regions represent the ranges in which decreasing the observational uncertainties lead and do not lead to better constraints on the atmospheric mass, respectively. The blue region represents the density range in which planets with an irradiation of 100 Earth fluxes have better constrained atmospheric mass with decreasing uncertainties and planets with an irradiation of one Earth flux do not.

3.5 Sensitivity of atmosphere-rich planets to irradiation

The calculations presented above correspond to an irradiation of 100 F_⊕. In this subsection we explore the sensitivity of the results to the assumed stellar irradiation. It should be noted that our model does not include atmospheric evaporation and the effect of the insolation is only reflected on the equilibrium temperature of the planet. Figure 9 shows Δσ∕ΔM of the envelope mass for all the simulated planets versus the bulk density and atmosphere radius for different assumed insolations. Clearly, there is a difference when changing the insolation for planets with relatively large atmospheres. A higher insolation increases the temperature of the H–He atmosphere, which leads to an expansion and a decrease in H–He layer density. In general, the atmospheric mass can be better constrained for strongly irradiated planets.

Figure 9 shows whether an improvement in observational uncertainties contributes to constrainthe envelope mass depending on the planetary density and irradiation. Planet O (M = 8 M_⊕, R = 2.9 R_⊕) has a null Δσ∕ΔM for Earth-like irradiation, but it is positive for high irradiation. If we look at this depending on the atmosphere radius (median of the posterior distribution) instead of the density, we see that the breakpoint coincides for the two insolations. The determination of atmospheric mass of planets holding atmospheres up to 0.5 R_atm is insensitive to an improvement of observational uncertainties. When the radius of the H–He envelope is lower than 0.5 R_atm there is a strong degeneracy with the other layers, but otherwise a high amount H–He is needed to fit the observed mass and radius, and therefore, variations on the observational uncertainties affect the determination of the H–He envelope. It is important to note that a more realistic model would include water in the atmosphere. However, this would be another free parameter that would increase the degeneracy and would be difficult to determine.

4 Dependence of internal structure determination on the mass derived from M–R relationships

Since not allplanets have measured masses and radii, (e.g., most Kepler planets in the past), published M–R relationships are often used to estimate the mass of a planet for a given radius and vice versa. Figure 4 of Otegi et al. (2020) compares some of the published mass radius relations. This figure shows significant disagreement in mass of approximately 25% for a given radius. Despite these differences, given the degeneracy when determining the internal structure, different masses inferred from various M–R relations could lead to very similar internal structures. In such a case, even if only the radius is measured, it is still possible to infer information on the planetary structure and bulk composition. In this section we explore whether these differences in mass are reflected in the determination of the internal structure parameters, or whether they get diluted because of the degeneracy.

Figure 10 shows the posterior distribution of the internal structure parameters for planets with the same radii and different masses as computed from the published M–R relationships (the masses and radii are listed in Table 3). We do not include other published M–R relationships as those presented in Weiss et al. (2013) or Zeng et al. (2016) because they also depend on an additional parameter (irradiation and core mass fraction, respectively). Since for a radius of 1.5 R_⊕ the masses corresponding to Chen & Kipping (2017) and Bashi et al. (2017) are almost identical we only show one of them (Chen & Kipping 2017). We assume a very small uncertainty of 1% for mass and radius.

For a radius of 1.5 R_⊕ the masses from different M–R relationships lead to significantly different inferred core mass distributions. The mass from Weiss & Marcy (2014) leads to an inferred core mass of 2.1 ± 0.7 M_⊕ and Chen & Kipping (2017) to 0.8 ± 0.5 M_⊕. At this radius the bulk densities lie between 4 and 7 g cm⁻³, so these planets are expected to have a negligible gaseous envelope. In addition, the water mass is small compared to the mass of refractory materials. Consequently, the degeneracy is reduced and the internal structure is more sensitive to a change in bulk density. As previously discussed, a change in the bulk density of planets dominated by refractory materials is mainly reflected in the core mass, and for this reason the core mass distributions for planets with 1.5 R_⊕ are sensitive to slight mass variations. There are some minor differences between the mantle and water mass distributions, especially between Chen & Kipping (2017) and the other considered studies.

For a radius of 2.5 R_⊕ the mass computed from the M–R relationship for rocky exoplanets in Otegi et al. (2020) is much higher compared to the masses from other published M–R relationships, and consequently the inferred internal structure is drastically different. Nevertheless, the difference between the internal structures inferred with the masses from other M–R relationships are very similar to each other. They lie in the intermediate density regime previously discussed, which is mostly degenerate. Consequently, there are not significant differences in the internal structure independently of the M–R relationship used, exceptfor that corresponding to rocky exoplanets as in Otegi et al. (2020); this work estimates a mass for 2.5 R_⊕ which is approximately three times larger than the other M–R relationships.

Finally, for radius of 3.5 R_⊕ Weiss & Marcy (2014) leads to much lower masses than the other M–R relationships, and therefore it leads to a significantly different estimated internal structure. The main difference at this radius arises in the posterior distribution of the atmospheric mass. We find that the choice of the M–R relationship for internal structure determination only in the solid-dominated and volatile-rich radius regimes is relevant.

We conclude that planetary internal structure should not be inferred from the radius alone. Even if for radii smaller than 2 R_⊕ most of the observed exoplanets closely follow the M–R relationships, relatively small variations in mass have a significant impact on the inferred core mass fraction. For radii in the range 2–3 R_⊕ where the degeneracy is the strongest, instead, small variations on the mass lead to almost identical inferred internal structure. However, this radius regime corresponds to the transition from rocky to volatile-rich exoplanets, and the large diversity of masses for a given radius does not allow us to infer the internal structure by radius alone.

Fig. 10 Posterior distributions of the bulk density and internal parameters of the planet for different radii and masses computed with various published M–R relationships.

Fig. 11 Atmospheric temperature and entropy profiles when artificially shifting the temperature gradient (solid lines) and γ (dashed lines). The temperature and entropy profiles of 25 planets with masses from 3 M⊕ to 20 M⊕ (from bottom to top) with a fixed composition of 0.5% H–He mass fraction, 40% water mass fraction, and 59.5% of a rocky core. On the bottom, the profiles corresponding to a mass of 15 M⊕.

5 Dependence of the radius on the temperature profile

Internal structure models depend on the assumed EOS and on several theoretical assumptions such as the internal structure (mixed, differentiated), the envelope structure and its luminosity, and the composition of the planetary layers. Since observational uncertainties are expected to decrease significantly in the upcoming years, it is desirable to explore whether the theoretical uncertainties related to the model assumptions could dominate compared to the observational uncertainties. In this section we investigate how the assumed temperature profiles affect the inferred planetary radius and study what uncertainty in the temperature profile of a planet is required to match the observational uncertainty.

Relative changes in material density have the largest effects on the radius for low-density materials, for example, variations in the metallicity distribution in the envelope significantly affect the inferred planetary radius (e.g., Lozovsky et al. 2018). We therefore explore how variations in the density, temperature, and entropy profiles in the H-He envelope affect the radius of a planet. A common assumption in the envelope structure models is to consider a fully adiabatic temperature profile. However, more realistic models should account for compositional gradients along the envelope, which leads to non-adiabatic temperatureprofiles. We do not know a priori what composition gradients might exist in most planets. In this section, we artificially perturb the temperature profile from the adiabat and study how it affects the inferred planetary radius.

We perturb the temperature profiles in two different ways: shifting the temperature gradient d T∕dR in the atmosphere and the γ as defined in Hansen (2008), which is the ratio of visible to infrared opacities. The perturbations in d T∕dR and γ are selected to produce a change in the radius comparable to the observational uncertainties. The effect of these perturbations on the temperature and entropy profiles is shown in Fig. 11. The top panels show the temperature and entropy profiles of 25 planets with masses from 3 M_⊕ to 20 M_⊕ and different assumed atmospheric temperaturegradients and γ. The bottom panels illustrate the effect of these variations in the internal temperature and entropy profiles of a planet with a mass of 15 M_⊕. Figure 12 shows the composition lines corresponding to a fixed composition and the perturbations applied to the temperatureprofile as presented in Fig. 11. The gray envelope corresponds to a radius uncertainty of 3%, which is expected to be reached with future missions. We find that small variations in the temperature profile, such as those corresponding to uncertainties of 20% in the temperature gradient or the ratio of visible and infrared opacities lead to radius uncertainties of 3%. This implies that theoretical uncertainties that are associated with the model assumptions could be larger than the observational uncertainties and therefore dominate the uncertainties in internal structure models. It would be necessary to do a systematic and detailed analysis of the uncertainties introduced by the model assumptions(e.g., EOS, composition of the layers, layer boundaries, intrinsic luminosity of the planet, and atmospheric opacities).

Fig. 12 M–R relationships for a planet with fixed composition and the perturbations applied to the temperature profile as presented in Fig. 11. The gray envelope corresponds to a radius uncertainty of 3%, which will soon be reached with future missions.

6 Summary and conclusions

We present new internal structure models based on the work of Dorn et al. (2017) with a Bayesian inference analysis using a nested sampling scheme. We explore several aspects that affect the characterization of exoplanets with masses up to 25 M_⊕ and radii up to 3.5 R_⊕, such as how variations in the observational uncertainties or location in the M–R affect the inferred internal structure, how the choice of the mass using different M–R relationships for a fixed radius influences the inferred internal structure, and how the atmospheric temperature profile affects the inferred radius. We should keep in mind that these results were computed for a given set of priors. The sensitivity of the results given various priors should be a topic for future investigation. Regarding the sensitivity of internal characterization to observed mass and radii, our main findings are summarized as follows:

• A decrease in observational uncertainties for planets below the Earth-like composition line leads to a better determination of the core mass.

• A decrease in observational uncertainties for planets above the pure-water composition line leads to a better determination of the atmospheric mass.

• A decrease in observational uncertainties for planets between Earth-like and pure-water composition lines does not significantly improve the determination of any internal structure parameter.

• The density boundaries listed above slightly depend on the used interior model (e.g., luminosity, insolation).

• The atmospheric mass of strongly irradiated atmosphere-rich planets can be better determined than for weakly irradiated planets.

• Using the stellar Fe/Si and Mg/Si abundances as a proxy for the bulk planetary abundances does not always help to constrain the planetary internal parameters when adding two extra model parameters related to the mantle composition. This depends on the actual value of the measured stellar abundances and their uncertainties and on the data of planetary mass and radius. Low stellar Fe/Si ratios improve the determination of the core mass for planets with densities between 2.3 and 6.9 g cm⁻³, where the planetary internal parameters are most degenerate.

We also find that internal structure must not be estimated using radius alone. For a fixed radius, the inferred planetary internal structure using different masses from various published M–R relationships can vary significantly depending on the measured radius. For radii of nearly 1.5 R_⊕ the choice of the M–R relationship significantly affects the inferred core mass, and for radii above 3 R_⊕ this choice significantly affects the inferred atmospheric mass. For planets with radii of nearly 2.5 R_⊕ small differences in mass lead to similar inferred internal structures, but the wide range of possible masses does not allow us to infer the internal structure with radius alone.

We find that uncertainties of 20% in the temperature gradient or the ratio of visible and infrared opacities γ can lead to radius uncertainties of 3%. Observational uncertainties are expected to decrease significantly in the near future, thereby decreasing the uncertainties in both mass (with advanced spectrographs like ESPRESSO) and radius (with next generation space missions like TESS, CHEOPS or PLATO). While these improved data will certainly help us to better understand planetary populations around other stars, a detailed characterization of individual planets is expected to remain somewhat limited. While the characterization of volatile-rich and dense exoplanets is expected to improve with decreasing observational uncertainties, a significant degeneracy in the internal structure of most of the super-Earth population is expected to remain. We also emphasize the importance of the theoretical uncertainty related to model assumptions (e.g., envelope structure, and compositionof the planetary layers), which may overcome the observational uncertainties soon. We therefore suggest that, along with the great efforts to improve the data, similar efforts should be made on the theoretical front.

Acknowledgements

We thank the referee for valuable comments which significantly improved our paper. C.D. acknowledges support from the Swiss National Science Foundation under grant PZ00P2_174028. J.H. acknowledges the support from the Swiss National Science Foundation (SNSF) under grant 200020_172746. This work has been carried out within the frame of the National Center for Competence in Research PlanetS supported by the SNSF.

References

All Tables

All Figures

	Fig. 1 Distributions of the relative mass (left) and radius (right) uncertainties for the 115 observed exoplanets less massive than 25 M⊕ from the NASA Exoplanet Archive.
In the text

Fig. 2 M–R diagram with the synthetic planets used in the study (black dots). The black, blue, and light blue bars correspond to uncertainties of 5, 10, and 20% in mass and 3, 5, and 10% in radius, respectively. The crosses represent observed planets from the revisited catalog in Otegi et al. (2020). The synthetic planets are plotted against the composition lines of pure-iron (brown), Earth-like (light brown), and water ice (blue).

In the text

	Fig. 3 Posterior distributions of core, mantle, water, atmospheric, and total masses of a planet with mass of 6 M⊕ and radius of 1.5 R⊕. The rows correspond to different radius uncertainties and colors in the subplots correspond to different mass uncertainties.
In the text

	Fig. 4 Standard deviation of the posterior distributions of core, mantle, water, and atmospheric masses vs. width of the total mass distribution for a planet of R = 1.5 R⊕ and M = 6 4 M⊕. For the atmosphere, the standard deviation of the posterior distribution in logarithmic scale is used.
In the text

Fig. 5 Model parameters and synthetic planets Δσ∕ΔM (see text for further details). A higher value of Δσ∕ΔM indicates that a decrease of the observational uncertainty provides additional information of the internal parameter. The crosses represent the observed planets from Otegi et al. (2020). Composition lines of pure-iron (brown), Earth-like planet (lightbrown), and water ice (blue) are also plotted.

In the text

	Fig. 6 Observational uncertainties in mass and radius versus σcore (left) and σatm (right) for the two synthetic planets with the highest and lowest densities.
In the text

	Fig. 7 Illustration of the degeneracy in the intermediate density regime previously defined. The posterior distributions the internal parameters of a planet with fixed radius of 2.3 R⊕ and mass of 15 M⊕ are shown The different colors represent different widths of the mass distribution.
In the text

	Fig. 8 Standard deviation of the core mass (left) and mantle mass (right) for a planet of mass 15 M⊕ and radius 2.3 R⊕ vs. the width of the mass distribution. The red, green, and blue lines represent different Fe/Si in the planet. A higher value of σmax − σ indicates that the internal parameter is better constrained with respect to the reference, at which σ = σmax.
In the text

Fig. 9 Atmospheric mass Δσ∕ΔM for the synthetic planets vs. the bulk density (left) and median atmospheric thickness (right) for different irradiation. The light green and red regions represent the ranges in which decreasing the observational uncertainties lead and do not lead to better constraints on the atmospheric mass, respectively. The blue region represents the density range in which planets with an irradiation of 100 Earth fluxes have better constrained atmospheric mass with decreasing uncertainties and planets with an irradiation of one Earth flux do not.

In the text

	Fig. 10 Posterior distributions of the bulk density and internal parameters of the planet for different radii and masses computed with various published M–R relationships.
In the text

Fig. 11 Atmospheric temperature and entropy profiles when artificially shifting the temperature gradient (solid lines) and γ (dashed lines). The temperature and entropy profiles of 25 planets with masses from 3 M⊕ to 20 M⊕ (from bottom to top) with a fixed composition of 0.5% H–He mass fraction, 40% water mass fraction, and 59.5% of a rocky core. On the bottom, the profiles corresponding to a mass of 15 M⊕.

In the text

	Fig. 12 M–R relationships for a planet with fixed composition and the perturbations applied to the temperature profile as presented in Fig. 11. The gray envelope corresponds to a radius uncertainty of 3%, which will soon be reached with future missions.
In the text