Fizzy Super-Earths: Impacts of Magma Composition on the Bulk Density and Structure of Lava Worlds

Kiersten M. Boley; Wendy R. Panero; Cayman T. Unterborn; Joseph G. Schulze; Romy Rodríguez Martínez; Ji Wang

doi:10.3847/1538-4357/acea85

1. Introduction

As terrestrial planets evolve, they are thought to transition through a magma ocean phase directly following formation (Elkins-Tanton 2012; de Vries et al. 2016; Schaefer & Elkins-Tanton 2018). During this phase, the temperatures are sufficiently high for global melting of the silicate mantle (Elkins-Tanton 2012; Chao et al. 2021). Although global magma oceans would generally be considered a transient phase in the lifetime of a terrestrial planet, it is possible for them to be long-lived. In particular, planets on close-in orbits will experience longer magma ocean lifetimes, maintaining their surface temperatures through the radiation of their host star (Chao et al. 2021). Ultra-short-period (USP) planets orbiting a host star with a strong magnetic field are also subject to an additional form of heating due to induction that could produce magma layers in the upper mantle as well (Kislyakova et al. 2017, 2018; Kislyakova & Noack 2020). These nontransient magma oceans or lava worlds undergo a distinct evolution that is unique to their temperatures and orbital distance.

After formation, many of these planets will initially experience significant atmospheric mass loss due to the high irradiation of their host stars, making them unlikely to maintain a substantial atmosphere (Lopez et al. 2012). There are several works that find USP super-Earths that are consistent with the absence of such atmospheres (Kreidberg et al. 2019; Keles et al. 2022). Although these planets would require rare initial conditions to maintain their primary atmosphere so close to their star, some planets may still have partial atmospheres due to the vaporization of silicate rock on their permanent dayside as a result of tidal locking (Léger et al. 2011; Zieba et al. 2022). During their evolution, they inevitably become tidally locked with their host star, resulting in a dayside and a nightside (Kasting et al. 1993; Barnes 2017). There is evidence that the dayside may possess an atmosphere at low pressures (P ≤ 10⁻⁵ bar, 1.5 Pa), consisting of rock vapor (Léger et al. 2011; Zieba et al. 2022). The proximity to their host star makes short-period planets not only more easily detectable but also optimal for characterization.

The vast majority of likely terrestrial planets with masses ≲10 M_⊕ or radii ≲1.7 R_⊕ discovered are on orbits less than 10 days. These discoveries have naturally prompted the characterization of these planets in order to understand their properties. For the majority of terrestrial planets without atmospheres, the stellar abundances of the host stars are correlated with the bulk density of the planet (Adibekyan et al. 2021; Schulze et al. 2021). Previous studies have specifically focused on solid planets while including the low densities with high equilibrium temperatures of some of the planets in their samples, e.g., 55 Cnc e and WASP-47 (McArthur et al. 2004; Becker et al. 2015; Adibekyan et al. 2021; Schulze et al. 2021). In these works, high-density planets are explained by an excess of iron. However, the description of low-density planets is largely unexplored due to their nonunique relationships between planet composition and density.

There are several reasons for which one might attempt to explain the low-density of these planets using magma oceans. The first of which is that the surface temperatures of many of these planets are above 1350 K, which is approximately the zero-pressure melting temperature of Earth-like rocks (Katz et al. 2003). Therefore, their surfaces may exceed the melting temperature of typical rocks, which range from 900 to 1500 K, making them unlike any planet in our Solar System (Takahashi & Kushiro 1983; Hirose & Kushiro 1993; Liu et al. 2010). Materials also expand when they are heated and melt; however, magmas are also highly compressible, potentially counteracting any expansion at depth (Rivalta & Segall 2008; McCormick et al. 2016). While Earth-like magmas are up to 24% less dense than their solid at the surface, they experience 20% compression compared to their zero-pressure volume at a pressure equivalent to 150 km below the Earth's surface (5 GPa) (Dannberg & Heister 2016), compared to a 6% compression of a solid of the same composition. As a result, it is possible that planets with deep magma ocean are denser than an equivalent-mass solid planet. However, it is not well-understood the degree to which this would impact the bulk density of the planet, and whether this difference is observable.

This paper will demonstrate that global magma oceans have little impact on bulk density of likely lava worlds. We will quantify the degree to which melting can lower the observed densities. For close-in super-Earths with masses ≳6 M_⊕, such as 55 Cnc e (McArthur et al. 2004) and WASP-47 e (Vanderburg et al. 2017), there have been many studies that propose several explanations for their low densities, which range from the presence of an atmosphere to a core-free planet (Gillon et al. 2012; Ito et al. 2015; Demory et al. 2016; Dorn et al. 2019; Dorn & Lichtenberg 2021). However, there has been little investigation on the impact of magma on planet density, independent of an atmosphere. This is particularly important for less massive planets that are too small to maintain a substantial atmosphere, such as TOI-561 b (Lacedelli et al. 2021).

We will address the impact of magma on the mantle structure of a lava world. Given that volatile species readily dissolve within magma, we explore the implications for trapping these volatiles in the interior over planetary lifetimes with implications for outgassing. Previous studies have shown evidence for basal magma oceans to occur in likely lava worlds and during the phases of Earth's evolution (Labrosse et al. 2007; Andrault 2019). H₂O and CO₂ are volatiles thought to be commonly outgassed from the mantle during the thermal evolution of a planet to form the early atmosphere (Lupu et al. 2014; Schaefer & Fegley 2017). By dissolving H₂O or CO₂ within magma, we can explore the amount of volatiles that may be stored on billion-year timescales. Therefore, we consider the mass–surface temperature relationship, exploring the mantle structures that arise and whether volatiles may be sequestered within the mantle.

The outline of our paper is as follows. Section 2 provides a description of our approach and additions to ExoPlex ⁵ (Unterborn et al. 2018, 2023; Unterborn & Panero 2019). In Section 3, we discuss the resulting density–radius relationships, comparing our findings to likely magma ocean planets, and the impact of magma composition on the bulk density and structure of a planet. We discuss the limitations of our model in Section 4, along with comparing our results to previous studies and known exoplanets. Finally, we provide a summary and conclude in Section 5.

2. Method

The base of our planet interior model uses ExoPlex, a thermodynamically self-consistent planet interior software (Unterborn et al. 2018, 2023; Unterborn & Panero 2019). ExoPlex solves five coupled differential equations: the mass within a sphere, hydrostatic equilibrium, adiabatic temperature profile, Gauss's law of gravity in one dimension, and the thermally dependent equation of state for solid planets. We construct a magma module within this software to include the liquid phase of the mantle when the temperature exceeds the melting point. As melting temperatures are compositionally dependent, we constrain this module to Earth-like mantle compositions. In this study, we assume a spherically symmetric planet for simplicity. However, tidal locking would cause lava worlds to have drastic temperature differences between the substellar and antistellar point (Léger et al. 2011). Therefore, this module considers the upper limit that magma may have on a planet's observed properties.

2.1. Mantle

We adopt an Earth-like pyrolitic molar composition of 0.5Na₂O–2CaO–1.5Al₂O₃–4FeO–30MgO–24SiO₂, which is within 1 wt% of the bulk silicate composition of the Earth (McDonough & Sun 1995; Solomatova & Caracas 2021). For a more realistic model, we incorporate both the solid and melt phases of pyrolite. To simplify, we neglect chemical fractionation, assuming the melt phase to have a liquid pyrolite composition at all temperatures above the melt curve. We also include the effects of volatiles within the melt phase using three separate model scenarios: Anhydrous, Hydrous, and Carbonated models.

In the multicomponent systems of planetary mantles, melting occurs over a range of temperatures. The solidus is the temperature at which the melting begins with a very small volume fraction of melt. As the temperature increases, the volume fraction of melt also increases until all components of the rock are liquid. The temperature at which the rock becomes completely liquid is referred to as the liquidus. For most silicate rocks, the solidus is hundreds of degrees lower than the liquidus. When a rock undergoes a small degree of melting just above the solidus, the trapped volatile species tend to enter the melt phase (Hirschmann 2000). As one of the objectives of this study is to determine whether a magma ocean is observable via the bulk density of a planet, use of the solidus as the temperature at which the planet melts will overestimate the density effects of melting. However, it is solid layers, forming a mechanical barrier, that will sequester volatiles deep in the interior. Therefore, we choose to focus on the solidus temperature, acknowledging that the effects on density are lower bounds.

Many previous studies that incorporate the effects of melt in exoplanets have used single-component melt curves (e.g., Stixrude 2014; Dorn et al. 2019; Dorn & Lichtenberg 2021). However, multicomponent systems are more likely to occur, given the diversity of rock-forming refractory materials available within the protoplanetary disk during formation (Helling et al. 2014). Therefore, we adopt the solidus temperature at which significant changes in physical properties occur. Specifically, we model the solidus using a method similar to that of Boukaré et al. (2022). Boukaré et al. (2022) rely on solidus and liquidus curves from Zhang & Herzberg (1994) and Fiquet et al. (2010). Adopting a similar method, we use a solidus that is lower than Zhang & Herzberg (1994), using methods that overcome experimental limitations of detecting small volumes of melts (Nomura et al. 2014).

To determine the solid–melt transition within the mantle, we use a peridotite solidus melt curve from Katz et al. (2003) for pressures less than 10 GPa for the anhydrous and hydrous compositions. This gives us the following equation for the solidus melt within the mantle:

$\begin{eqnarray}&&{T}_{\mathrm{Melt},\mathrm{Anhydrous}}={a}_{1}{P}^{2}+{a}_{2}P+{a}_{3},\end{eqnarray} \tag{ 1 }$

where P is pressure, a₁ = −5.1° C GPa⁻², a₂ = 132.9° C GPa⁻¹, and a₃ = 1358.85 K.

For the hydrous melt, we use the Katz et al. (2003) solidus melt curve for hydrous melts, which has an additional term to account for H₂O within the melt:

$\begin{eqnarray}&&{T}_{\mathrm{Melt},\mathrm{Hydrous}}={a}_{1}{P}^{2}+{a}_{2}P+{a}_{3}-{{KX}}_{{{\rm{H}}}_{2}{\rm{O}}}^{\gamma },\end{eqnarray} \tag{ 2 }$

where K = 43° C wt%^−γ and γ = 0.75. We assume the wt% of H₂O in the melt is ${X}_{{{\rm{H}}}_{2}{\rm{O}}}=2.2$ wt%, which corresponds to the wt% used in Solomatova & Caracas (2021) hydrous magma equation of state.

To extend each solidus to higher pressures for all melt compositions, we fit a power law similar to the Simon melting law (Ross 1969; Stixrude 2014):

$\begin{eqnarray}&&{T}_{\mathrm{Melt}}={T}_{\mathrm{ref}}{\left(\tfrac{P}{{P}_{\mathrm{ref}}}\right)}^{b},\end{eqnarray} \tag{ 3 }$

where T_ref and P_ref are the reference temperature and pressure at higher pressures. We adopt the anhydrous melt temperature and pressure from Nomura et al. (2014) of T_ref = 3570 K at P_ref = 136 GPa. We use these values to provide an upper bound for the transition between solid and melt at higher pressures. We then fit for b so that the computed melt curve intersects the peridotite solidus temperature of T_ref = 3570 K at P_ref = 136 GPa along with the solidus melt curves at lower pressures. We find the b for the anhydrous melts to be b_anhydrous = 0.188 above 10 GPa. For the hydrous melt curve, we similarly use the Nomura et al. (2014) constraint as an upper bound for the hydrous melting between 10 and 136 GPa. As melting temperatures are closely related to density, differences in the volatile content of the melt for pyrolitic compositions become less prominent at higher pressures. Therefore, anhydrous, carbonated, and anhydrous melts converge at higher pressures. For the hydrous model above 10 GPa, we find b_hydrous = 0.201.

For the carbonated melt, we use an approach similar to the above, fitting 5.2 wt% carbonated peridotite solidus data from 0 to 10 GPa Dasgupta & Hirschmann (2006) and the Nomura et al. (2014) high-pressure point (i.e., T_ref = 3570 K at P_ref = 136) to the Simon melting law in Equation (3) to derive b for the carbonated melt to be b_carbonated = 0.265.

The consequences of extrapolating these melt curves beyond Nomura et al. (2014) are minor as the difference in compression at such high pressures becomes modest, and our structural model is insensitive to the location of the melt curve at higher pressures. Our computed melt curves are shown in Figure 1.

**Figure 1.** Solidus melt curves for the three magma compositions: Anhydrous (purple), Hydrous (teal), Carbonated (magenta). We include the pressure–temperature profile of a planet at a surface temperature, T_surf = 2000 K and 1 Earth mass (orange). These are compared to the melt curves within the work of Dorn & Lichtenberg (2021), who use a piecewise solidus melt curve (green).
Download figure:
Standard image High-resolution image

To calculate the EOS of the solid mantle, we use ExoPlex. ExoPlex relies on a fine-mesh grid approach to calculate the stable mantle mineral assemblage for a given pressure and temperature from PerpleX (Connolly 2009) grids. We construct a grid consistent with the pyrolitic composition used in Solomatova & Caracas (2021), as we also use their data for the melt within the mantle. Within ExoPlex, we set the mantle composition to reproduce the pyrolytic composition as in Solomatova & Caracas (2021), using the following molar fractions: Ca/Mg = 0.067, Si/Mg = 0.8, Al/Mg = 0.05, and Fe/Mg = 0.9. Given that the molar fraction of iron within ExoPlex includes the core and mantle, we specify a mass fraction of 0.079 wt% FeO within the mantle, to remain consistent with Solomatova & Caracas (2021).

For the melt, we use ab initio molecular dynamics data for pyrolite from Solomatova & Caracas (2021) shown in Table 1. Within ExoPlex, we recast the Solomatova and Caracas series of isothermal equations of state between 2000 and 5000 K to a single pressure (P)–density (ρ)–temperature (T) equation of state model for each composition by fitting their results as follows:

$\begin{eqnarray}&&P(\rho ,T)=P({\rho }_{0},2000\,{\rm{K}})+{P}_{\mathrm{th}},\end{eqnarray} \tag{ 4 }$

where P(2000 K) is calculated by the BM3 equation and the 2000 K ρ₀ is as in Table 1. The thermal pressure (P_th) is expressed by the Mie–Grüneisen relation:

$\begin{eqnarray}&&{P}_{\mathrm{th}}=\displaystyle \frac{\gamma }{V}\left[E(T,{\theta }_{{\rm{D}}})-E(2000\,{\rm{K}},{\theta }_{{\rm{D}}})\right],\end{eqnarray} \tag{ 5 }$

where γ is the Grüneisen parameter, $\gamma ={\gamma }_{0}{({\rho }_{0}/\rho )}^{q}$ , and the Debye temperature ${\theta }_{{\rm{D}}}={\theta }_{{\rm{D}},0}\exp [({\gamma }_{0}-\gamma )/q]$ . The harmonic internal energy E(T, θ_D) is calculated from the Debye model (e.g., Fei et al. 1992).

Table 1. ExoPlex Model Parameters Recalculated from Solomatova & Caracas (2021) for Magma with Earth-like Compositions, and with Added 2.2 wt% H₂O and 5.2 wt% CO₂

Melt Composition	ρ₀ (g cm⁻³)	K₀ (GPa)	${K}_{0}^{{\prime} }$	wt%
Pyrolite	2.49	24	4.7	⋯
Pyrolite+ 4 H₂O	2.32	15	5.5	2.2
Pyrolite+ 4 CO₂	2.37	18	5	5.2

Notes. Each value in the table is at a reference temperature of 2000 K fit using the third-order Birch–Murnaghan EOS. For all compositions, we find a best-fit thermal equation of state of γ₀ = 1.7, q = 0.93, and θ_D,0 = 1000 K.

Download table as: ASCII Typeset image

To determine the best-fit thermal parameters, we use a Markov Chain Monte Carlo (MCMC). For all compositions, we use the higher-temperature calculations (up to 5000 K) to constrain γ and q. We assume θ_D,0 = 1000 K, as fluctuation in θ_D,0 accounts for minimal differences in γ and q. Within each step of the MCMC, a γ and a q are proposed, we then use the Solomatova & Caracas (2021) 2000 K data to determine the pressures and densities for the higher temperature (up to 5000 K). The densities at 135 GPa are compared to Solomatova & Caracas (2021) data at higher temperatures using the two-sample Anderson–Darling statistic (Anderson & Darling 1952; Pettitt 1976). We find the best-fit thermal parameters for all compositions to be: γ₀ = 1.7, q = 0.93, and θ_D = 1000 K.

Using the best-fit thermal parameters, we self-consistently calculate the interior temperature along an adiabat that corresponds to the surface temperature of the planet. Therefore, we fit the suite of pressure (P)–density (ρ)–temperature (T) data to the BM3 Mie—Grüneisen equation of state at a reference temperature of 2000 K for dry, hydrous, and carbonated pyrolite from Solomatova & Caracas (2021) as shown in Table 1. We then determine the density along the adiabat, drawing from the melt or solid EOS depending on the relative position of the adiabatic temperature and the solidus.

2.2. Metallic Core

We explore two core scenarios as an upper and lower bound on the density variation that may arise from the core instead of mantle melting. We assume an Earth-like core composition. Given the surface temperatures required to produce a global magma ocean, our model assumes a liquid metallic core, as the core–mantle temperatures will exceed that of the melting temperature of iron.

As an upper bound on the bulk density, we assume a pure liquid iron core. For this scenario, we rely on the liquid Fe pressure and temperature grid within ExoPlex. This grid uses the temperature-dependent equation of state of Anderson & Ahrens (1994) and is calculated up to 15 TPa and 10,000 K (Unterborn et al. 2023).

We include a lower bound on the bulk density with the addition of lighter elements in the core, which results in a decrease in the bulk density of a planet. We use FeO within our models with realistic mass fractions within the core: Fe = 84% and O = 16%. We choose these values to represent a lower bound on the density, as the highest estimate for light alloys in the core of Earth is ∼16% by mass (Hirose et al. 2021). ExoPlex assumes that, to first order, the presence of light elements lowers only the molar weight of liquid Fe while not affecting its compressibility at pressures indicative of exoplanetary cores. This assumption means that, for an Fe core containing 16 wt% O, the density of the planet would be ∼7% lower than a planet with a pure Fe core for a 1 M_⊕ solid planet.

2.3. Interior Compositions

To investigate the impacts of volatile compositions on the mass and radius of a planet, we consider four different model scenarios in which each scenario assumes a liquid iron core:

1.
Solid rocky mantle.
2.
Solid rocky mantle and anhydrous melt.
3.
Solid rocky mantle and a carbonated melt.
4.
Solid rocky mantle and a hydrous melt.

Figure 2 illustrates each model scenario with masses equal to four Earth masses at a surface temperature of 2000 K, demonstrating that for low masses generally, the hydrous magma composition has the largest radius while the solid planet has the smallest radius.

3. Results

Incorporating melting temperatures into our model, we find three distinct magma ocean mantle structures when surface temperatures exceed the solidus. Specifically, we find that the mantle of a lava world may be completely molten, only molten at the surface of the planet, or exhibit a layered magma ocean structure hereafter referred to as a MOSMO structure. Here, we define a mantle with a MOSMO structure as having a magma ocean (MO) at the surface, a solid rock layer (S) mid-mantle, and a basal magma ocean (MO) at the base of the mantle (see Figure 3). It is this intermediate solid rock layer that has the potential to serve as a mechanical barrier to volatile transport from a basal magma ocean to the surface.

**Figure 3.** Mantles structures produced by our model. We find that the mantle may be mantle magma ocean, a surface magma ocean and solid rock layer, or a MOSMO structure (i.e., Surface Magma Ocean (MO)–Solid Rock Layer (S)–Basal Magma Ocean (MO)).
Download figure:
Standard image High-resolution image

In Figure 4, we show the "phase diagram" for the mantle structure, including several known planets that are likely lava worlds. The planets are colored by their bulk density ratio, which is their observed density divided by the expected density assuming a solid, Earth-like composition with an Earth-like light-element budget and core-mass fraction. Bulk density ratios <1 indicate low-density planets, whereas those >1 indicate likely super-Mercuries. We construct this planet sample using the Exoplanet archive, making the following cuts. We require all planets within the sample to have mass and radius measurements. We remove any planets with flagged or controversial measurements. Using the period-dependent radius gap as an upper bound on planet radius, we remove planets with larger radii (Van Eylen et al. 2018). We then explore the structure and bulk density of planets with surface temperatures above 1350 K, assuming no atmosphere and a surface albedo of zero. We choose 1350 K as it is the lowest temperature that is predicted by the anhydrous melt curve. Our results suggest that MOSMO structure and surface magma ocean planets are slightly more common than planets that are completely molten (see Figure 4).

We investigate the impact of surface temperature on the depth of the solid–melt transition. Using the median mass of our planet sample, 4 M_⊕, we vary the surface temperature, displaying the lower boundary of the surface magma ocean and the upper boundary of the basal magma ocean (see Figure 5). As the surface temperature increases, the depth of the surface magma increases until the mantle no longer contains a solid layer at depth. We find that the planet becomes completely molten at a surface temperature (T_surf) of 2050 K for Models (2) and (3), whereas Model (4) becomes mantle magma ocean at a T_surf = 2100 K. In the interval where a planet's mantle transitions from solid to mantle magma ocean, the planet will form a MOSMO structure.

Figure 5(A) shows that, for a given mass, a planet will have a surface temperature range in which they exhibit a MOSMO structure that spans T_surf ∼ 2000 K for 4 M_⊕ planets. The formation of a basal magma ocean occurs at T_surf ∼ 1872 K, but varies minimally depending on the magma volatile content. Any differences between the models are not significant, given the assumptions of the melt curve at high pressures. Within this range, the majority of the melt is located within the basal magma, with the surface melt only deepening by a few percent with increasing surface temperature.

The depth of the surface magma ocean on a 4 M_⊕ planet is dependent on the magma composition and the imposed surface temperature (see Figure 5(A)). Model (2), the anhydrous magma, has a surface magma ocean that reaches a depth of 68 km (4.8 GPa) given a surface temperature of 1600 K. Comparing Model (3) and Model (4) to Model (2), we find that Model (2) has a shallower surface magma for a given surface temperature, with the most prominent depth difference occurring at T_surf = 1600 K. At this surface temperature, the surface magma oceans of Model (3) and Model (4) respectively have 10% and 12% greater depth than that of Model (2). The difference in surface magma depth between the models becomes less pronounced with increasing surface temperature. Before the mantle becomes mantle magma ocean for Model (2), the depth of the surface magma ocean reaches 732 km (28 GPa) at a T_surf ∼ 1900 K. At 1900 K, Model (3) and (4) only reach ∼4% greater depth than Model (2). Overall, we find that the surface magma oceans for all model scenarios are shallow, consisting of less than 5% of the radius fraction before the mantle becomes mantle magma ocean.

The formation of a basal magma ocean on a 4 M_⊕ planet is invariant with respect to the volatile composition of the mantle (see Figure 5(A)). We find the basal magma ocean forms with a T_surf ∼ 1870 K), melting only at the base of the mantle, 570 GPa. As the surface temperature increases, the top of the basal magma ocean rises to 118 GPa for a T_surf ∼ 2080 K. At this point, the MOSMO structure transitions to a mantle magma ocean as the top and bottom of the solid, mid-mantle region meet. Therefore, the surface temperature at which the basal magma ocean first forms relies on a significant extrapolation of the melting curve, while the existence of the MOSMO structure only relies on solidus temperatures at pressures not requiring extrapolation. Uncertainties in the melting curve parameters in Equations (1)–(3) will shift the top and bottom of the solid mid-mantle region to higher or lower pressures, but will not eliminate the layer.

Although we show a 4 M_⊕ planet in Figure 5(A), we investigate the impact that mass has on magma depth as well. For planets <4 M_⊕, we find that the MOSMO structure spans a smaller range of surface temperatures, which can be seen in Figure 7. The reduced surface temperature range results in a reduced mantle melt fraction of the basal magma ocean. For a 4 M_⊕ planet at T_surf1900 K, the basal magma ocean would account for 19% of the planet radius whereas a basal magma ocean on a 2.5 M_⊕ planet at T_surf1900 K would only account for 5.2% of the planet's radius. The reduction in the basal magma ocean depth occurs due to the relative location of the solidus with respect to the geotherm, the temperature–pressure gradient of the planet. The opposite is true for planets >4 M_⊕, where MOSMO structure exists over a wider surface temperature range and at greater pressures.

We also include an analysis of the impact of surface temperature on the bulk density of a planet. Generally, we find that bulk density decreases with increasing surface temperature (Figure 5(B)). The density decrease for Model (1) within Figure 5(B) demonstrates the effects of the thermal expansion of solid rock. The bulk density of a lava world decreases at a higher rate than that of a solid planet at surface temperatures that only produce a surface magma ocean, demonstrating the greater thermal expansion of magmas at low pressures. However, Model (2) and Model (3) exhibit an increase in the bulk density within the surface temperature range where the MOSMO structure is present. This result indicates that the magma is highly compressed under the basal magma ocean pressure conditions, resulting in a density greater than solid rock under equivalent conditions. Therefore, the rapid increase in the basal magma fraction causes a local maximum in the bulk density before decreasing once the planet is above the solidus at all depths. Surface temperature, mass, core composition, and magma composition all impact the bulk density of a planet.

3.1. Impacts of Magma Composition

To quantify the effects of magma and volatile composition on the density of the planet, we calculate the percent difference between the magma models and Model (1). We use Model (1) as the control model to determine the impact that the magma composition has on the mass and radius of the planet. Each subsequent model is compared with Model (1), as shown in Figure 6. We set the mass to 4 M_⊕ and surface temperature of each model to T_surf = 1600 K and T_surf = 2000 K, considering the impact of surface temperature on differences in magma composition.

**Figure 6.** The percent difference between the bulk density of each magma ocean planet and the solid planet with the corresponding core (i.e., light alloy core or pure iron core) at surface temperatures of (A) 1600 K and (B) 2000 K. The top axis shows radius values for Model (2). The dashed lines represent models with light-element alloys in the core, whereas solid lines show models with pure iron cores. Model (2), the anhydrous magma, is shown in purple. Model (3), the hydrous magma, is shown in teal and has the greatest decrease in density. Model (4), the carbonated magma, is shown in magenta.
Download figure:
Standard image High-resolution image

Generally, we find that all four models behave similarly at T_surf = 1600 K compared to T_surf = 2000 K, with any features being more prominent at higher surface temperatures. As planet mass increases, the density difference between a magma ocean and an equivalent-mass solid planet decreases, reflecting the effect of steeper pressure gradients in the planet and greater magma compressibility. At both surface temperatures of both 1600 and 2000 K, Model (4) results in the greatest decrease in density whereas Model (2) results in the least. Broadly speaking, all models find that low-mass planets exhibit the largest fractional inflation due to magma, compared to higher mass planets. We will focus on the T_surf = 2000 K results, as that is the average surface temperature of the known low-density planets with surface temperatures above 1350 K. However, we will briefly address the lower surface temperature model.

As described above, we find that planets at T_surf = 1600 K follow similar bulk density trends as a function of planet mass as those at T_surf = 2000 K. Model (4) results in the lowest density, being 4.3% less dense than an equivalent-mass solid planet at a mass of 0.7 M_⊕. Model (2) results in the least density decrease at 2.9% at a mass of 0.7 M_⊕. As mass increases, we find that the decrease in the bulk density due to the presence of magma becomes less pronounced, with the density difference being approximately 0.2% for 10 M_⊕. With the greatest density decrease for planets with a surface temperature of 1600 K being <4.3%, the inflation due to magma would be indistinguishable from a solid planet with observational uncertainties on mass and radius.

For planets at T_surf = 2000 K, we find distinct differences from the T_surf = 1600 K models. At 2000 K, Model (2) results in the least difference in density when compared to Model (1). It exhibits the most significant decrease in density of 9.8% at lower masses and radii, with the largest decrease at a mass and radius of 0.7 M_⊕ and 0.93 R_⊕, respectively. This reflects a planet with a mantle magma ocean mantle ocean with a core–mantle boundary pressure of 93.6 GPa. As mass and radius increase, we find that the anhydrous magma becomes more compressible. Therefore, the anhydrous magma becomes denser than that of the solid planet (Model (1)). We find that this planet density crossover occurs at a mass and radius of 3.14 M_⊕ and 1.45 R_⊕. For the range of masses that we consider, we find the highest density increase of 2.4% compared to Model (1) at a radius of 1.83 R_⊕.

The carbonated magma composition or Model (3) at T_surf = 2000 K results in a moderate decrease in density of 10.4% at a radius of 0.94 R_⊕ and mass of 0.7 M_⊕. Similarly to Model (2), the magma becomes more compressible with increasing mass. At a radius of 1.51 R_⊕ and 3.85 M_⊕, Model (3) becomes denser than Model (1), with a maximum density increase of 1.3%. Compared to Model (2), the planet density crossover occurs at a higher mass and radius.

Model (4), the hydrous magma, results in the greatest decrease in density at 11.1% compared to all other models at a radius of 0.94 R_⊕, mass of 0.7 M_⊕, and surface temperature of 2000 K. Given the EOS parameters for this magma composition (see Table 1), we expect the hydrous magma to start at a lower density. For the mass range that we consider, it would not become more compressible than the solid planet. However, we are limited to a smaller range of masses, due to the hydrous model becoming dynamically unstable at masses greater than 8.7 M_⊕. The solid rock layer within the mantle becomes denser than the basal magma beneath, requiring dynamical simulations.

3.2. Light Elements in the Core

We explore the impact of light elements in the core, assuming an Earth-like composition. In Figure 6, we show the percent difference between the models with pure iron cores to the models with light alloys in the core. As the mass increases, the difference in density of the magma compared to an equivalent-mass solid planet is reduced as a result of the converging compressibility of the materials with increasing pressure (planet depth). When considering Earth-like compositions as described by Fe/Mg, the effect of light elements in the core is to reduce the planet's density, due to the reduced mass of the core. However, the radius of the core increases to account for the added volume of oxygen when Fe/Mg remains constant. Assuming a pure iron core, a planet with a mantle magma ocean results in a bulk density difference of ∼−9% at 1 M_⊕, and a MOSMO structure results in a ∼2% density difference at 10 M_⊕. The combined effects of a light alloy core and a magma ocean are ∼−11% and ∼1.5% at 1 M_⊕ and 10 M_⊕, respectively (Figure 6(B)). Therefore, the impact on the bulk density is dominated by the effects of the magma ocean at constant Fe/Mg. This result suggests that the bulk density of an Earth-like composition lava planet with an Earth-like core light-element budget and core-mass fraction is more sensitive to the presence of magma than light elements in the core. The effects of the presence of light elements are magnified at greater core-mass fractions, potentially equivalent to the predicted density deficits of a surface magma ocean (Unterborn et al. 2023).

4. Discussion

4.1. Basal Magma Oceans

Several studies have proposed the potential for basal magma oceans to exist as a planet cools (Elkins-Tanton et al. 2005; Labrosse et al. 2007; Pachhai et al. 2022). Magma oceans may solidify in various ways: downward from the surface, upward from the core, or outward from the mid-mantle. The behavior of the solidification is a function of the mantle composition, and differences in gradient between the melting temperature and the local geotherm. When a planet solidifies from the mid-mantle, the magma ocean beneath the rock layer is a basal magma ocean. However, the inferred mechanical behavior of the intermediate rock layer changes depending on the melt curve that is used.

We use the solidus melt curve to identify the existence of a solid rock layer separating the basal and surface magma. This differs from a liquidus melt curve, as the rock layer would be a potentially porous crystal mush (Marsh 1989). The resulting solid layer that we find could potentially trap volatiles deep within the mantle, as there would not be a path for them to outgas. Lava worlds that are continuously bombarded with radiation from their host star likely have a dry surface melt (Léger et al. 2011; Kite et al. 2016). Earth-sized planets have been shown to be inefficient at outgassing their volatiles (Miyazaki & Korenaga 2022; Salvador & Samuel 2023) over the timescales of the cooling Earth shortly after formation. However, it is unclear if long-lived magma oceans can retain their volatiles over long timescales. For this reason, the presence of a MOSMO structure is potentially significant for maintaining large volatiles within the mantle (Labrosse et al. 2007; Edmonds & Woods 2018; Caracas et al. 2019; Bower et al. 2022), even if the surface magma ocean eventually degasses. For a 4 M_⊕ mass planet with a volatile concentration similar to those of Model (4) and (3), the presence of a basal magma ocean could potentially sequester approximately 130 times the mass of the water in Earth's present-day oceans and 1000 times the carbon in the Earth's surface and crust, respectively, for the water and carbon contents of the modeled magma.

For simplicity, we do not consider a thermal boundary layer within the solid mid-mantle rock layer of the MOSMO structure. However, a nonconvective solid rock layer with a thickness of a few kilometers would be enough to cause a thermal boundary layer (Monteux et al. 2016; Andrault 2019), resulting in thinning of the mid-mantle solid layer. However, a fully convective system, in which material melts and solidifies as it crosses the solid material depth, will reduce the net thermal boundary layer resulting from producing and consuming the heat of fusion.

4.2. Lowest-density Planet Produced

To discuss the impact of magma on reducing the density of a planet, we consider our planet sample of likely lava worlds in terms of density and radius illustrated in Figure 7. We include density–radius curves for each of our models for planets with surface temperatures of 2000 K. Because magma is highly compressible, the greatest effects of its thermal expansion are seen at lower pressures or in Earth-sized planets as opposed to super-Earths. In particular, when considering larger super-Earths' masses (7–8 M_⊕), our results suggest that magma ocean worlds without an atmosphere would be denser than their cooler counterparts.

**Figure 7.** Density–Radius plot for each model that is equivalent to a planet with a surface temperature of 2000 K. We include likely lava worlds selected using criteria described in Section 3. The planets are colored by their equilibrium temperatures. Planets distinguished with white circles and planet names are discussed in detail in Sections 4.2 and 4.3. Model (4), the hydrous magma, is limited to a radius of 1.7 R_⊕, because the density of the magma becomes greater than the density of the solid in the surface magma ocean, making it dynamically unstable.
Download figure:
Standard image High-resolution image

Within our planet sample, TOI-561 b, CoRoT-7 b, and TOI-500 b have the lowest bulk densities with densities of 3.00 ± 0.80 g cm⁻³, 5 ± 1.5 g cm⁻³, and ${4.9}_{-0.88}^{+1.03}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , respectively (Léger et al. 2011; Lacedelli et al. 2021; Serrano et al. 2022). However, we find that the lowest bulk density produced by our models for the same planet mass is only consistent with the measurements for TOI-500 b. For TOI-561 b and CoRoT-7 b, our model overestimates the bulk density by 79% and 50%. If we instead calculate the density with light elements within the core discussed in Section 3.2, we find a bulk density of 5.1 and 4.9 g cm⁻³ for TOI-561 b and CoRoT-7 b. These densities are still 69% and 50% greater. These results suggest that magma alone cannot account for the decrease in planets' bulk density, assuming they have Earth-like compositions. The impact of magma on the density of a planet is less than typical observational uncertainties on mass and radius, such that planets below the expected density for a rocky planet cannot be explained by the presence of melt.

Assuming an Earth-like composition for low-density planets requires the presence of an atmosphere regardless of the presence of a magma ocean. However, there is little evidence to suggest that a USP planet could sustain a substantial atmosphere that would significantly decrease the bulk density of a planet at such close proximity to its host star. Even assuming the presence of rock vapor atmosphere can only increase the radius of a super-Earth planet by approximately 1%. For this reason, it may be more likely that these planets may have significantly different compositions from Earth, with potentially lower core-mass fractions (Schulze et al. 2021; Unterborn et al. 2023).

4.3. Known Lava Worlds

Seventy-five percent of the super-Earths that have been discovered to date are on orbits of less than 10 days. Given the presence of an atmosphere, ∼50% of these planets have surface temperatures large enough to sustain a magma ocean (Dorn & Lichtenberg 2021). However, without an atmosphere, the majority of these planets would still be able to maintain a magma ocean on their dayside.

Of the hundreds of super-Earths discovered, only six planets⁶ have been discovered with masses and radii that fall below the potential planet density crossover of 3.14 M_⊕ and 1.45 R_⊕ at surface temperatures greater than 1350 K. Therefore, the majority of the hot, relatively low-density planets exist in the regime where the presence of magma would cause them to have a slight increase in bulk density if they possess a magma ocean. For this reason, we cannot attribute the extremely low densities of some likely lava worlds primarily to magma. Instead, models addressing hot, relatively low-density planets should consider an atmosphere or smaller core-mass fraction in addition to magma.

To investigate the impact of magma on known planets, we select planets from our planet sample (Section 3) that have a high probability of being atmosphere-free or having a thin atmosphere that could not result in significant changes in the density of the planet. When available, we select likely lava worlds from our sample evaluated within Schulze et al. (2021) that have a high probability that the inferred compositions from mass and radius alone are statistically indistinguishable from that of their host star at the >1σ level, assuming no atmosphere. Therefore, we include K2-265 b and WASP-47e, which have high probabilities of having no more than a thin atmosphere and a core-mass fraction comparable to their host star (94% and 80%, respectively) (Schulze et al. 2021).

We also include three additional planets that are likely to be atmosphere-free or have thin atmospheres adopting the nominally rocky planet zone (NRPZ) introduced by Unterborn et al. (2023). The NRPZ is the likelihood that the planet's observed mass and radius alone are consistent with the planet having a bulk rocky composition, without the addition of significant surface volatiles (Unterborn et al. 2023). We include HD 80653 b, Kepler-10 b, and K2-141 b, which have respective NRPZ probabilities of 66%, 60%, and 55%.

Below, we discuss these five likely atmosphere-free lava worlds. Due to the current observational uncertainties of planet mass and radius for this population, we cannot distinguish between Models (1)–(4). However, we include a discussion on bulk density, to demonstrate that magma alone cannot cause extremely low densities in planets. We also discuss the required observational uncertainties to distinguish between a solid planet and lava world for each planet. Therefore, we choose Model (4), as it produces the largest differences from Model (1), the solid planet. We also assume a pure iron core.

We provide the following discussion primarily to consider the impacts that high surface temperatures may have on the structure of a planet. By relying on the underlying physics, we can gain insight into the potential structures that may arise assuming an Earth-like composition. We consider the potential melt fraction and structure of the mantle assuming an Earth-like composition.

4.3.1. K2-265 b

K2-265 b is a short-period planet discovered in 2018 on an orbit of 2.37 days with a density of 7.1 ± 1.8 g cm⁻³ (Lam et al. 2018). It has a mass and radius of 6.54 ± 0.84 M_⊕ and 1.71 ± 0.11 R_⊕, respectively (Lam et al. 2018). Given the close proximity its host star, K2-265 b is exposed to strong stellar irradiation, likely resulting in the photoevaporation and loss of its atmosphere (Lopez et al. 2012). Assuming the planet is a blackbody with a zero albedo, the dayside equilibrium temperature is calculated to be T_eq ∼ 1400 K.

Applying Model (4), we find that 1% of the mantle is molten at the surface, due the low equilibrium temperature. We find the bulk density to be 7.94 g cm⁻³, 11% greater than the measured density. The density is consistent with a lava world, due to observation uncertainties; however, it is also indistinguishable from a solid planet, as Model (1) also produces a bulk density of 7.95 g cm⁻³. With such a small percentage of melt, there is minimal impact on the bulk density, and the surface magma ocean is likely entirely degassed. With the low equilibrium temperature, the uncertainties on mass and radius required to differentiate a lava world with 1% melt fraction must be less than 0.01%.

4.3.2. WASP-47 e

WASP-47 e is a USP planet on an orbit of 0.7895 days (Bryant & Bayliss 2022). It has an observed bulk density of 6.29 ± 0.60 g cm⁻³ with a mass and radius of 6.77 ± 0.57 M_⊕ and 1.808 ± 0.026 R_⊕ (Bryant & Bayliss 2022). Assuming a zero albedo, the dayside equilibrium temperature of the planet is 2200 K. As it is a well-characterized super-Earth, previous studies conjecture that the low density of WASP-47 e may be due to a Ca- and Al-rich interior or the presence of a magma ocean and secondary atmosphere (Dorn et al. 2019; Gupta & Schlichting 2021). However, these studies do not account for effects due to magma within their models.

Applying Model (4) to WASP-47 e, we find that the mantle has a MOSMO structure. The total melt within the mantle is 88% by volume. The majority of the melt is within the basal magma ocean within 80%, whereas only 8% of the melt within the mantle is in the surface magma ocean. The basal magma ocean, therefore, has the potential to trap a significant inventory deep in the interior over extended times. Dynamical models of transport across a solid intra-mantle would be required in order to determine if outgassing may still be feeding an atmosphere at the surface.

Considering the density of WASP-47 e, Model (4) produces a density of 7.89 g cm⁻³. Given the uncertainties on the mass and radius, Model (4) produces a radius that is above the uncertainties on the observed density. Model (1) produces similar results, with the bulk density being 7.95 g cm⁻³. Although Model (4) has a slightly lower density, the presence of a magma ocean would not increase the radius enough to account for the observed low density. In order to distinguish between a solid or lava world for WASP-47 e, the uncertainties on mass and radius required are ≲0.14% and ≲0.05%, respectively. Therefore, it is necessary to invoke compositional variables, such as a reduced core-mass fraction, or a thick, unobserved atmosphere in order to account for the low density.

4.3.3. HD 80653 b

Discovered in 2020, HD 80653 b is a USP planet on an orbit of 0.719 day. It has a mass and radius 5.60 ± 0.43 M_⊕ of 1.613 ± 0.071 R_⊕, respectively (Frustagli et al. 2020). The observed bulk density is 7.4 ± 1.1 g cm⁻³. Using a zero albedo, the dayside equilibrium temperature of the planet is 2300 K.

Using Model (4), we find that the mantle of HD 80653 b is a mantle magma ocean structure, where the magma ocean extends from the surface to the core–mantle boundary, and therefore it is unlikely to be able to trap volatiles in the interior over the age of the system. The calculated bulk density is 7.4 g cm⁻³. With Model (1), we find the bulk density to be 7.7 g cm⁻³. Therefore, the density of HD 80653 b is consistent with either a solid or lava world, due to the uncertainties on the observed mass and radius. To differentiate between Model (1) and Model (4), the uncertainties required would be ≲1.8% on the mass and ≲0.6% on the radius of the planet.

4.3.4. Kepler-10 b

As the first rocky planet discovered with Kepler, Kepler-10 b is a well-known USP planet on an orbit of 0.8374 day with a mass of ${3.57}_{-0.53}^{-0.51}$ M_⊕ and radius of ${1.489}_{-0.021}^{-0.023}$ R_⊕ (Batalha et al. 2011; Dai et al. 2019). With an equilibrium temperature of T_eq = ${2130}_{-60}^{120}$ K, it is a low-density lava world that does not deviate far from a solid planet on the density–radius relationship. The observed bulk density is 6.0 ± 1.1 g cm⁻³ (Dai et al. 2019).

We apply our models to Kepler-10 b with Model (4); the mantle consists of a mantle magma ocean structure. Similarly to HD 80653, a magma ocean structure that extends from the surface to core–mantle boundary is unlikely to be able to trap volatiles in the interior over the age of the system. The calculated bulk density is 6.6 g cm⁻³. Applying Model (1), we find the bulk density to be 6.8 g cm⁻³. Therefore, a solid or lava world would be consistent with the observed density, as they both produce densities that are well within the observational uncertainties. The required uncertainties on mass and radius to distinguish between Model (1) and Model (4) are ≲1.9% and ≲0.06%, respectively.

4.3.5. K2-141 b

K2-141 b is another USP planet on an orbit of 0.2803 day with a density of 8.2 ± 1.1 g cm⁻³ (Barragán et al. 2018; Malavolta et al. 2018). It has a mass and radius of 5.08 ± 0.41 M_⊕ and 1.51 ± 0.05 R_⊕ with an equilibrium temperature of T_eq = ${2161}_{-60}^{120}$ K (Malavolta et al. 2018). However, K2-141 is also an active star with strong magnetic activity (Barragán et al. 2018). Therefore, induction heating could also be a potential heat source for K2-141 b. There is evidence showing that K2-141 b is inconsistent with a thick atmosphere (Malavolta et al. 2018; Zieba et al. 2022). However, it could possess a thin rock vapor atmosphere, but this would not significantly lower the bulk density. For this reason, we include this planet assuming that is atmosphere-free, as a thin rock vapor atmosphere would not significantly impact the density.

Applying Model (4) to K2-141 b, we find the planet to have a MOSMO structure. The total melt fraction of the mantle is 82% by volume, with the basal magma ocean containing the majority of the melt (68% by volume). Therefore, the surface melt accounts for only 14%. As with WASP-47e, the MOSMO structure has the potential to trap volatiles in its interior with the basal magma ocean if the solid mid-mantle layer prevents mass transport to the surface.

Considering the bulk density of the K2-141 b, we apply Model (4). We calculate the bulk density to be 7.3 g cm⁻³. Model (1) produces a bulk density of 7.4 g cm⁻³. Due to the uncertainties on the mass and radius, both Model (1) and Model (4) are consistent with the observed density. To distinguish between the solid planet and magma ocean model, the uncertainties on the mass and radius must be ≲0.2% and ≲0.06%, respectively. However, K2-141 b is a likely lava world that is denser than expected for a solid Earth-like planet. Therefore, Model (4) would be more appropriate, but the observational uncertainties required to distinguish between the models are not yet achievable. Given that it is denser than expected for an equivalent-mass solid planet, K2-141 b may also have significant iron enrichment indicative of an Fe-rich super-Mercury, in addition to a magma ocean.

5. Summary and Conclusions

Nearly half of the terrestrial planets discovered to date could maintain magma on their surfaces, in particular on their dayside. Therefore, it is important to understand the way in which magma oceans may impact the structures and observable properties of these likely lava worlds, in order to better characterize them. Given that many of the likely lava worlds on USP have yet to detect substantial atmospheres large enough to significantly decrease its bulk density (e.g., Léger et al. 2011; Kreidberg et al. 2019; Keles et al. 2022; Zieba et al. 2022), compositional factors such as the impact of magma or core-mass fraction must be considered when characterizing these planets.

In this paper, we investigated the impact of magma on the 1D structure and bulk density of atmosphere-free magma ocean planets for the following magma compositions: anhydrous, hydrous, and carbonated magma. The objectives of this study were to determine whether a magma ocean is observable via the bulk density of a planet and to determine if volatiles may be trapped in the interior. Therefore, we constructed our model using a solidus melt curve placing an upper limit on the impact of magma on the bulk density of a lava world and placing a conservative lower limit on the mantle structure. From this study, we present our primary conclusions:

1.
The presence of a magma ocean alone is not sufficient to explain low-density magma ocean planets that are expected to be atmosphere-free or have thin atmospheres (Section 4.2).
2.
For a given mass, there exists a range of surface temperatures in which a planet will have a basal magma ocean, which may sequester a significant amount of dissolved volatiles (Section 3).
3.
The addition of H₂O or CO₂ to the magma does not significantly impact the calculated bulk density of a planet, only resulting in a maximum density difference of ∼1% (Section 3.1).
4.
For magma ocean planets that are atmosphere-free, the presence of magma can impact the bulk density of the planet, causing two distinct regimes where magma ocean planets exhibit a planet density crossover that is dependent on mass and surface temperature. This leads to two regions where a magma ocean planet may be more or less dense than an equivalent-mass solid planet for anhydrous or carbonated magmas (Section 3.1).
5.
For an Earth-like core light-element budget and core-mass fraction, the addition of magma has a greater impact on the bulk density than the addition of lighter elements within the core for planets with masses and radii less than ∼3.14 M_⊕ and ∼1.45 R_⊕ (Section 3.2).

Acknowledgments

K.M.B. acknowledges support from the NSF Graduate Research Fellowship Program under grant No. (DGE-1343012). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. W.R.P. acknowledges support from NSF under grant No. (EAR-1724693). C.T.U. acknowledges support under grant NNX15AD53G. 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. J.W. acknowledges the support by the National Science Foundation under grant No. 2143400.

Software: ExoPlex (Unterborn et al. 2018, 2023; Unterborn & Panero 2019).

Fizzy Super-Earths: Impacts of Magma Composition on the Bulk Density and Structure of Lava Worlds

Article metrics

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction