INTERIOR PHASE TRANSFORMATIONS AND MASS–RADIUS RELATIONSHIPS OF SILICON–CARBON PLANETS

The chemical diversity of stars in the universe has been suggested to result in an even greater chemical diversity among the planets which they host. According to condensation models (Bond et al. 2010; Kuchner & Seager 2005), a key variable determining the chemical makeup of planets is the carbon-to-oxygen ratio of the disk. In stellar nebulae with carbon–oxygen ratios in excess of 0.8 (Bond et al. 2010), condensation models propose that solid bodies within the ice line would consist primarily of silicon carbide and carbon rather than the silicate materials of our own solar system, leading to the formation of solid planets consisting of silicon, carbon, and possibly iron with minimal oxygen (Bond et al. 2010). Recent work on carbon-based planets has been particularly motivated by the detection of planets in the 55 Cancri system, whose C/O ratio is particularly high. In particular, the planet 55 Cancri e, whose mass–radius relationship suggests that it may be insufficiently dense to have a silicate composition (Demory et al. 2011), has been modeled as consisting of layers of carbon, SiC, and iron (Madhusudhan et al. 2012). Recent interpretations of spectroscopic data Nissen (2013), however, have called into question whether stars with C/O ratios over 0.8 indeed exist. With this in mind, it is important to compute accurate interior models and mass–radius relationships for carbide planets to assist in resolving the controversy of the composition of planets such as 55 Cancri e.

The accuracy of interior models of carbon-based planets has been hampered by a lack of experimental or theoretical data on the behavior of silicon–carbon materials at high pressures. The best available equations of state (EOSs) for high-pressure silicon carbide are based on the extrapolation of experimental data (Aleksandrov et al. 1989; Seager et al. 2007). Although the high-pressure phase diagram of carbon and the EOS have been extensively studied up to extremely high pressures of 100 Mbar using ab initio theory (Martinez-Canales et al. 2012), the phase diagrams of silicon carbide and pure silicon have not been studied at multi-megabar pressures. Furthermore, as we will show, the assumption that SiC remains the sole stable stoichiometry of the silicon–carbon binary system at extreme pressures is not justified, analogous to the anomalous stoichiometries seen at high pressures in materials such as MgSiO₃ (Umemoto et al. 2006) and H₂O (Zhang et al. 2013; Pickard et al. 2013). In order to build accurate interior models of silicon–carbon-rich planets, and hence to determine the expected mass–radius relationships for silicon–carbon planets as a function of composition, it is thus necessary to model more accurately the high-pressure behavior of these materials.

In this work, we compute the high-pressure phase diagram and EOSs of the silicon–carbon binary system up to pressures of 40 Mbar. We use an ab initio random structure search algorithm to find relevant high-pressure phases of silicon, carbon, silicon carbide, and Si_nC_m compounds of alternative stoichiometries. We predict new high-pressure phases of silicon carbide and silicon, and new stable high-pressure compounds with Si₂C and SiC₂ stoichiometry. We compute EOSs for each relevant material and derive mass–radius relationships for silicon–carbon planets with interior pressures up to 40 Mbar. Our EOS for SiC is 12% denser than that used by Seager et al. (2007) leading to substantially smaller radii for silicon–carbon planets than had previously been suggested. We derive layered models for planets which take the novel Si₂C and SiC₂ phases into account. In addition, we revise earlier EOSs for high-pressure iron. Our new results allow a refinement of the interior models of 55 Cancri e derived by Madhusudhan et al. (2012), which allows us to place stronger constraints on possible compositional models for this planet; in particular, we eliminate the possibilities of pure SiC or SiC–iron compositions.

For pressures higher than 4 Mbar that cannot be reached with diamond anvil cell experiments, our understanding of planetary materials relies primarily on theoretical methods. In order to compute the EOSs of silicon and carbon containing materials at high pressure, it is first necessary to know the phase diagram of the silicon–carbon binary system including the onset of phase transitions and the stoichiometric relationships, which requires a search through the space of possible structures to find the ground-state phase and stoichiometry as a function of pressure. The Ab Initio Random Structural Search (AIRSS) algorithm (Pickard & Needs 2006) has emerged as a successful method for finding stable crystal structures of materials, particularly at high pressures, with an efficiency comparable to more algorithmically complex methods such as genetic algorithms. In the AIRSS methodology, randomly generated cell geometries are filled with randomly positioned atoms. Efficient geometry relaxation procedures are used to find the nearest local minimum. Although most randomly generated structures do not lead to the absolute ground-state structure, it is found that a sufficiently large portion do so to allow the identification of ground state. Consistently and repeatedly finding a single structure to be lower in enthalpy than all competing structures may be considered to be reasonable evidence that it is the true lowest-enthalpy structure at that pressure. Although no optimization scheme can be guaranteed to find the lowest-enthalpy structure, the past success of AIRSS (Pickard & Needs 2011) gives us reasonable confidence in its ability to produce ground-state structures.

Here, ab initio random structure searches were undertaken at 10, 20, and 40 Mbar. For each set of pressure and stoichiometry (C, Si, SiC, and Si_mC_n structures detailed below), we begin by generating at least 800 random structures of between 1 and 4 stoichiometric units. All density functional theory (DFT) simulations in this paper used the VASP code (Kresse & Furthmüller 1996), pseudopotentials of the projector augmented wave type (Blöchl 1994), and the exchange-correlation functional of Perdew et al. (1996). The positions are initially optimized using a conjugate gradient algorithm until the difference between successive energies is less than 0.001 eV, using a plane wave cutoff of 500 eV and a grid of 4 × 4 × 4 k points to sample the Brillouin zone. Following the initial runs, the 50 structures lowest in enthalpy are subjected to a second, more accurate minimization which uses a denser 12 × 12 × 12 k-point grid and a 1200 eV cutoff energy for the plane wave expansion. Finally, to obtain accurate enthalpy–pressure curves, we recompute the stablest few structures using a 1200 eV cutoff and 32 × 32 × 32 k points to ensure accurate and comparable enthalpies between structures with different cells. As a sole exception, a 24 × 24 × 24 k-point grid was used for the SiC rocksalt structure due to memory constraints; however, this should not have an appreciable effect on our results.

High-pressure phases of elemental carbon have been the study of a multiplicity of previous studies (Yin 1984; Grumbach & Martin 1996; Correa et al. 2008), and most recently an AIRSS study (Martinez-Canales et al. 2012) which predicted zero-temperature phases of carbon up to 100 Mbar. We find an identical progression of ground-state structures, shown in Figure 1. Our transition pressures match those of previous calculations, with carbon in the diamond phase transitioning to the BC8 phase at 10.0 Mbar and then to the simple cubic structure at 28.9 Mbar, compared to 9.9 and 29.0 Mbar in Martinez-Canales et al. (2012).

Zoom In Zoom Out Reset image size

Pure silicon has been the study of fewer high-pressure theoretical studies. A transition from the hexagonal close packed (hcp) to the face-centered cubic (fcc) structure was predicted and experimentally confirmed to occur at 0.78 Mbar by Duclos et al. (1987). Using AIRSS, we confirm that the fcc structure remains the ground-state structure at 10 and 20 Mbar. However, at 40 Mbar, we find the body centered cubic (bcc) structure to be more stable. From a plot of enthalpy versus pressure (Figure 1), we find a transition from bcc to fcc at 27 Mbar.

Although silicon carbide exhibits a complex phase diagram at lower pressures, the high-pressure phase diagram is relatively simple. At 10.5 kbar, silicon carbide is known experimentally to transform into the rocksalt phase (Sekine & Kobayashi 1997). No other structure has been predicted at higher pressures. Our calculations find the rocksalt phase to remain stable at pressures as high as 10 and 20 Mbar. At 40 Mbar, however, we find a new SiC structure with Cmcm symmetry to be the stablest geometry. This structure is found to be structurally identical to the B33 structure of CrB, and is a layered structure in which each C or Si atom has five equally near neighbors of the opposite species within one layer unit, and two slightly further neighbors of the opposite species in the next layer. The structural parameters of the Cmcm SiC structure are given in Table 1. The transition from the rocksalt to the Cmcm structure is found to occur at 27 Mbar.

We now turn our attention to the search for silicon–carbon structures with alternative stoichiometries. Here, we restrict our attention to structures with simple integer ratios of atomic species—C:Si = 1:2, 1:3, 1:4, 2:3, 2:1, 3:1, 4:1, and 3:2—since these ratios account for the overwhelming majority of known binary compounds. The space of compounds with each of these stoichiometries up to four functional units per cell was searched at 40 Mbar, and the enthalpy of the stablest compounds at each stoichiometry plotted as a function of elemental ratios in a convex hull diagram as seen in Figure 2. The convex hull diagram depicts the stability of material phases as a function of their chemical composition. If a point lies below the line joining two adjacent compositions, then this structure is stable relative to an (unmixed) combination of the adjacent phases. At 40 Mbar, we find that two new stoichiometries have become stable relative to the combination of other compounds: SiC₂ (stabler than SiC + C) and Si₂C (stabler than SiC + Si).

Zoom In Zoom Out Reset image size

The SiC₂ structure found at 40 Mbar is a structure with Cmmm symmetry. The Cmmm structure consists of alternating rows of silicon atoms which are bonded to four silicon atoms in a planar configuration, and a second class of carbon atoms which are bonded to six silicon atoms. We are not aware of any other compound displaying this crystal structure. The structural parameters of the Cmmm structure are given in Table 2. Subsequent AIRSS searches at 20 and 10 Mbar found this structure to remain the ground-state stoichiometry at these lower pressures.

For Si₂C, a structure with I4/mcm symmetry was found to be the ground state at 40 Mbar. Examination of this structure reveals it to be identical to the C16 structure of Al₂Cu. In this structure, carbon atoms form close-packed linear chains with each carbon atom equidistantly spaced from eight Si atoms. The structural parameters of this structure are given in Table 3. Similar to the SiC₂ case, AIRSS searches at 20 and 10 Mbar found this structure to remain the ground-state stoichiometry at these lower pressures.

Next, we computed the EOSs of each of the SiC₂ and Si₂C structures and the formation energy relative to SiC + Si or SiC + C. Figure 3 shows the enthalpy of the Si₂C and SiC₂ phases relative to separate phases of SiC and C/Si. We find that SiC + Si will form the I4/mcm Si₂C phase at 13 Mbar, with the formation enthalpy continuing to increase with pressure. The Cmmm SiC₂ structure becomes stable at a higher pressure of 23 Mbar relative to SiC + C. The phase transition to the simple cubic structure of carbon, however, results in the formation enthalpy of SiC₂ decreasing above 29 Mbar. Although SiC₂ remains stable relative to SiC + C at 40 Mbar, the highest pressure studied here, a transition back to SiC + C stability is possible at higher pressures but lies beyond the scope of this work.

Zoom In Zoom Out Reset image size

Having established the phase diagram of the silicon–carbon system, we now compute the volume and enthalpy results as a function of pressure. Since our goal is to determine mass–radius relationships on a planetary scale, for the purpose of this calculation we ignore low-pressure phases such as graphite and the many phases of SiC existing below 10 kbar, which affect only the first few tens of kilometers of the planet; this leads to an underestimate of the planetary volume on the order of 10 km for an Earth-sized planet. EOS data for all silicon-carbon structures are shown in Table 4.

In Figure 4, we compare our DFT-based SiC EOS calculations with two EOSs that were previously constructed for SiC. First, we show a Birch–Murnaghan fit that was constructed from the high-pressure diamond anvil cell experiments by (Aleksandrov et al. 1989). The density of SiC was determined with X-ray diffraction measurements up to a pressure of 0.425 Mbar. Results were fit to a third order Birch–Murnaghan equation (Birch 1947). Even though this equation was constructed to describe the compression of materials, significant uncertainties are introduced when one extrapolates this EOS fit by two orders in pressure. It is therefore not too surprising that our DFT calculations predict densities for SiC that are between 20 and 35% higher. Seager et al. (2007) combined the experimental results by Aleksandrov et al. with the predictions from the Thomas–Fermi–Dirac theory (Salpeter & Zapolsky 1967) in order to construct a modified polytrope EOS for SiC for the purpose of planetary interior modeling, which was later also used by Madhusudhan et al. (2012). The densities derived from the modified polytrope EOS are approximately 12% lower than we obtained with DFT calculations. This correction directly implies that the radii of SiC planets have previously been significantly overestimated (Seager et al. 2007; Madhusudhan et al. 2012). Since DFT has been validated for a wide range of materials and thermodynamic conditions (Tuckermann 2002; Kirchner et al. 2012; Parrinello 1997), and the Thomas–Fermi–Dirac theory only becomes valid at ultra-high pressures where chemical bonds can no longer exist, we consider our DFT results significantly more reliable in the megabar regime under consideration.

Zoom In Zoom Out Reset image size

Following Seager et al. (2007), we solve the equations of hydrostatic equilibrium and mass conservation to derive the mass–radius relationship of different planets:

$$\begin{eqnarray} \frac{dP}{dr} &=& {-}\frac{G m \rho }{r^2}, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{dm}{dr} &=& 4 \pi r^2 \rho.\quad \end{eqnarray} \tag{ 2 }$$

We start the integration in the planet's center with r = 0, m = 0, and a central pressure, P = P_c, and then integrate outward until the pressure decreases to zero. m(r) indicates the mass that is enclosed in radius r. The second equation describes how m(r) changes when a new mass shell is added. The first equation characterizes the change in pressure that balances the difference in gravitational potential that an additional layer introduces. For efficiency reasons, we solve these equations with a fourth-order Runge–Kutta method using a fixed step size of dr = 50 km. Alternatively a simple Euler integration with dr = 1 km may also be used. The EOS only enters through ρ = ρ(P). We neglect temperature effects, which were estimated by Seager et al. (2007) to increase the radius of planets by only $\tilde{1}.2\%$. For each material's phase, we construct a separate spline function, ρ = ρ(P), to interpolate our DFT results. When a phase transition occurs at a certain pressure, we switch discontinuously from one spline function to the next.

We integrated Equations (1) and (2) for pure SiC planets to learn how differences in the EOS affect the mass–radius relationship. Figure 4 shows that our DFT results imply that SiC planets are approximately 5% smaller than predicted by Seager et al. (2007) and Madhusudhan et al. (2012). Because of this correction, it is no longer valid to model 55 Cancri e as a pure SiC planet, which was one possible scenario that was recently proposed by Madhusudhan et al. (2012) among other interior models. Our DFT results instead predict 55 Cancri e to have another light outer layer in addition to the SiC core. A likely candidate would be a carbon layer.

Following Madhusudhan et al. (2012), we constructed a suite of ternary interior models with an iron core, an SiC mantle, and an outer carbon layer. The central pressure, P_c, and the pressures where we switch from iron to SiC, P₁, and from SiC to carbon, P₂, are free parameters. We constructed a fine three-dimensional grid ranging from 6.5 ⩽ P_c ⩽ 24.5 Mbar and 0 ⩽ P₁/P_c ⩽ 1 and 0 ⩽ P₂/P₁ ⩽ 1. We selected models where the sum of the χ² deviations in mass and radius from the observed values were less than two. We adopted a radius value of 2.173$^{+0.097}_{-0.098}\, R_E$ that combined Spitzer and Microvariability and Oscillations of STars observations (Gillon et al. 2012) rather than constructing interior models for each radius measurement separately. A mass of M = 8.39 ± 0.38 M_E (Endl et al. 2012) was assumed. All valid models are summarized in the compositional diagram in Figure 5. When we performed this analysis with a modified polytrope EOS for SiC, our results are consistent with those by Madhusudhan et al. (2012) and a planet of pure SiC would be consistent with observations. However, when we switch to using our more accurate DFT EOS for SiC, maximum SiC fraction drops to only 52%. Because we predict SiC to be a denser material, a thick outer carbon layer must compensate for this change.

Zoom In Zoom Out Reset image size

In Figure 5, we compare the iron and SiC mass fractions from our Fe–SiC–C models that match the observed mass and radius. In the inset, we display the same information in a ternary composition diagram, as in Valencia et al. (2007), where each corner corresponds to a planet of one material only. When we compared the effects of using our DFT EOS for SiC with the modified polytrope EOS from Seager et al. (2007), all other model parameters were kept constant. Both EOSs consistently predict that 55 Cancri e can only contain up to 18% iron. This limit appears to be insensitive to changes in the SiC and Fe EOSs because the density of iron is so much higher than that of the other materials. Based on our DFT results, we predict 55 Cancri e to be consistent with a carbon-rich planet with a carbon fraction of 48% or more. All permitted models fall into a triangle in composition space that is spanned by three limiting cases: (1) pure carbon planet, (2) an iron-free SiC–C planet with 48% carbon, and (3) an SiC–free iron–carbon planet with 82% carbon. It is not possible to resolve this degeneracy with the existing constraints on mass and radius. In principle, additional information can be obtained from in situ simulations (Bond et al. 2010).

In this section, we explore how the formation of the novel compounds, SiC₂ and Si₂C, will affect the interior structure of Si–C planets. It is our goal to determine under which conditions these compounds form, what layers emerge in the planet's interior, and how much of each compound is produced over time. As it turns out, the mass–radius relation is not affected very much because SiC₂ and Si₂C form at relatively high pressures, but we address these interior questions here so that a planetary model can readily be constructed for other classes of materials. From our DFT simulations, we derived the compositional phase diagram in Figure 3. This diagram describes which minerals would form from a Si–C mixture at a certain pressure. For example, a carbon-rich Si–C mixture would split into SiC and carbon for pressures up to 25 Mbar. For higher pressures, this mixture would split into SiC₂ and SiC or, for carbon fractions larger than 67%, would form SiC₂ and carbon. For silicon-rich assemblages, a similar change is triggered by the formation of Si₂C at 13.6 Mbar.

First, we will discuss the formation of a carbon-rich planet with a carbon atom fraction of 95%. We assume homogeneous accretion of SiC and carbon in fixed proportions. Materials in the planetary interior differentiate and form separate layers that are sorted by density Turcotte & Schubert (2014). Each layer is assumed to be homogeneous but have distinct chemical composition and be fully convective.

For an accreting carbon-rich planet, these assumptions and the phase diagram in Figure 3 imply the formation of a two-layer planet with a SiC core and carbon mantle until the planet exceeds a critical size. For a carbon fraction of 95%, this critical size corresponds to a planet with a radius of 2.886 R_E and a mass of 24.92 M_E as illustrated in Figure 6(a). The SiC core would be comprised of only 2.77 M_E and have a radius of 1.04 R_E. If the pressure at the core–mantle boundary exceeds 25 Mbar, then SiC₂ will form from SiC and carbon. Since SiC₂ has an intermediate density, this new layer will form between the SiC core and carbon layer. If additional SiC–C material is accreted onto the surface of the planet, then the increased gravitational force will temporarily increase the pressure at the SiC₂–C boundary beyond 25 Mbar. When the sinking SiC material arrives at the bottom of the carbon layer, it will react with the carbon present to form additional SiC₂. Assuming that plenty of SiC is available, this implies the existence of a feedback mechanism that stabilizes the pressure at the SiC₂–C boundary at 25 Mbar during accretion.

Zoom In Zoom Out Reset image size

However, determining whether sufficient SiC is available at the SiC₂–C boundary is not straightforward. If one grows the planet assuming a constant total composition and that the SiC₂–C boundary remains at the critical pressure of 25 Mbar, then some reactant SiC material has to come from the SiC core. Assuming that the core provides sufficient SiC, the planet will assume the state of chemical equilibrium shown in Figure 6(b). As the planet mass increases from 27.68 to 29.00 M_E during accretion, the SiC core would shrink from 2.77 to only 1.32 M_E in such a model. This would require a significant amount of gravitational energy and it is not obvious which mechanism could provide that. However, planetary interiors are complex and the equilibrium model is certainly one that needs to be considered. The scenario of a shrinking SiC core would share some similarities with the core erosion that has been proposed to occur in giant planets (Guillot et al. 2004; Wilson & Militzer 2012a, 2012b; Wahl et al. 2013).

In the absence of any obvious energy source that would be needed to shrink the SiC and bring the whole planet into chemical equilibrium, we also wish to discuss an alternative dynamic scenario, where no mass is removed from the SiC core. We assume that the SiC core would grow to the maximum size that is reached when the pressure at the SiC₂–C boundary attains 25 Mbar. From that point on, the SiC₂ layer would only grow from newly accreted SiC material sinking through the carbon layer. As the planet mass increases from 27.68 to 29.00 M_E, the SiC core would only be more compressed and its radius would shrink from 1.04 to 1.02 R_E. A comparatively thin SiC₂ layer would form that comprises only 0.20 M_E compared to 2.37 M_E in the equilibrium scenario. The SiC₂ layer would be starved for SiC and the pressure at the SiC₂–C boundary would reach 26.2 Mbar (Figure 6(c)). This exceeds the 25 Mbar needed for SiC₂ formation but there is insufficient SiC available. This condition is also consistent with the phase diagram in Figure 3.

For a fixed planet mass and carbon fraction, the deviations in the predicted radii between the equilibrium and the dynamic scenario are very small. For a planet of 29 M_E, we obtained 2.892 and 2.895 R_E, respectively. We conclude that either model can be use to compare with observations in the future.

For a silicon-rich planet with a carbon atom fraction of 45%, the differences between the equilibrium and the dynamic scenario are a bit more pronounced because the Si₂C already forms 13.6 Mbar. A planet may accrete up to 11.62 M_E and reach a radius of 2.091 R_E before the Si₂C layer forms (Figure 7(a)). According to the equilibrium model, the SiC core would then be completely absorbed into the growing Si₂C layer (Figure 7(b)) as the planet reaches a total mass of 16.15 M_E. An amount of 4.84 M_E of Si₂C would be formed according to the equilibrium picture, while in the dynamic scenario one would predict an intermediate Si₂C layer of only 1.36 M_E to form. Figure 7 shows that the predicted planet radii are again very similar.

Zoom In Zoom Out Reset image size

In Figure 8, we summarize the mass–radius relationships of different types of planets. For this purpose, we also revisited the accuracy of the iron EOS used in Seager et al. (2007) by performing DFT calculations of the relevant bcc and hcp phases. Our computed EOS data for iron are shown in Table 5. We identified a modest correction. At 10 and 40 Mbar, we predict densities 2.2% and 4.7% higher. The predicted radii of pure iron planets with 5 and 10 Earth masses would shrink by 1% and 2%, respectively.

Zoom In Zoom Out Reset image size

We added the different types of Si–C planets to Figure 8. As expected, Si–C planets with 45% and 95% carbon closely track the results of the pure SiC and pure carbon planets, respectively. The formation of intermediate Si₂C and SiC₂ layers does not significantly affect the mass–radius relationship. More surprising is that the radii of pure SiC planets are very similar to those of rocky planets made of 100% silicates. Our revision of the SiC EOS put this material much closer to silicate rocks. Therefore, SiC cannot serve as a low-density material to explain the interior structure when observations suggest a radius larger than that of pure silicate planets.

We have conducted extensive simulations of the phase diagram and EOS of the silicon–carbon system at pressures up to 40 Mbar. Using ab initio random structure search methods, we predicted a new phase of silicon carbide and a bcc to fcc transition in silicon. In addition, we find two phases, SiC₂ and Si₂C, which are formed at high pressures for carbon-rich or silicon-rich stoichiometries, respectively. Our newly calculated EOS for silicon carbide is approximately 5% denser at high pressures than the extrapolated EOS used in previous works, leading to a significant downwards revision of predicted mass–radius relationships for SiC–C planets and eliminating the possibility of a 55 Cancri e model made of pure SiC. We also present a revised EOS for iron.

At present, 55 Cancri e remains the sole identified candidate for a silicon–carbon dominated exoplanet, however, its C/O ratio remains the subject of controversy. Future work to identify and characterize additional carbon planet candidates may be able to resolve the question of the existence and composition of carbon-based planets.