JSTOR

2.1. Background and Equations

We begin by assuming the major components of a solid exoplanet are limited to an iron core, a silicate mantle, and a water–ice outer layer. In other words, we assume the interior of the planet is differentiated with the denser materials interior to the less dense materials. We further assume each layer to be homogeneous in its composition.

We can then model the interior of a solid exoplanet by using:

1.	the equation for mass of a spherical shell where m(r) is the mass included in radius r and ρ(r) is the density at radius r;
2.	the equation for hydrostatic equilibriumwhere P(r) is the pressure at radius r; and
3.	the equation of state (EOS) that relates P and ρ.

The EOS is different for each different material. We used Fe (ϵ) for the planet core, MgSiO₃ perovskite for the silicate mantle, and water-ice VII, VIII, and X for the water–ice outer layer. See Seager et al. (2007) for a detailed discussion of the EOSs including their source. The temperature has little effect on the EOS especially in the high pressure regime (Seager et al. 2007); we ignore the temperature dependence of the EOS. This simplifies the equations and their solution, while enabling a relatively accurate analysis.

In this problem, we have five variables:

1.	the iron mass fraction (α);
2.	the silicate mass fraction (β);
3.	the water–ice mass fraction (γ);
4.	the total mass of the planet (Mp); and
5.	the mean radius of the planet (Rp).

The variables α, β, and γ are not independent of each other. Based on the assumption that the planet only consists of iron, silicate, and water we have: α + β + γ = 1, which can also be expressed as γ = 1 - α - β. We therefore have four variables: α, β, M_p, and R_p. Given any three of these variables, we can determine the fourth variable uniquely. We can also see that given M_p and R_p, there is a relationship between α and β. There is not a single value of α and β that produces a given M_p and R_p. Instead, there are infinite pairs of α and β that give the same M_p and R_p, and we call this a degeneracy in the interior composition.

2.2. Algorithm for Solving the Differential Equations

Our program integrates from the surface r = R_p inward to the center of the planet. The outer boundary condition is m(R_p) = M_p and P(R_p) = 0. That is, at the surface, the mass is the specified total planet mass and the pressure is approximately zero.

We aim to interpret observations of a planet of a given mass and radius. We therefore choose to integrate inwards instead of outwards based upon the known parameters of the planet (M_p, R_p, and P(R_p)). The independent variable is r, decreasing from r = R_p to r = 0. The interior boundary condition is m(r) = 0 at r = 0. Typically, when integrating m(r), m does not equal zero at r = 0. We therefore must iterate, tuning the mass fraction of each layer until m = 0 and r = 0 is reached.

Given a single M_p and R_p the computer program finds all possible combinations of α, β, and γ( = 1 - α - β) that give the same planet mass and radius. The program begins with a chosen value of α, and then takes a guess of β. Next, the computer program integrates differential equations (1) and (2) using the above boundary conditions to find the value of the planet radius. By comparing this radius to the desired R_p, the computer program tunes the value for β using the bisection method. This process is repeated several times, until β is found to a satisfactory accuracy of 1/1000. By varying α within the range of 0 to 1, we can get all possible combinations of α and β, which produce a specific M_p and R_p.

2.3. Algorithm for Generating a Database of M_p = 0.5–20 M_⊕

M_p and R_p are observed parameters and one usually wants to find the corresponding allowed α and β. For a range of M_p and R_p—corresponding to observational uncertainties—it can be very time consuming to use the first algorithm described in § 2.2. We therefore generate a database that is a discrete representation of the relation R_p = R_p(α,β,M_p). Figure 1 illustrates the 4D database.

This database is a 3D array that contains the data of R_p corresponding to each combination of α and β (0 to 1 with 1% spacing) and M_⊕ (ranging from 0.5 to 20 M_⊕ with 0.25 M_⊕ spacing). The database can be used via linear interpolation to find R_p for any given M_p, α, and β. A conservative estimate of the fractional error in the database interpolation is 1/1000. For the same range of M_p and R_p, interpolation in the database is about 45 times faster than solving the differential equations.

To generate a database of all values of α, β, and R_p for M_p ranging from 0.5 to 20 M_⊕. Given M_p, α, and β, this algorithm also integrates from the surface inward to r = 0 and m(r) = 0 to find R_p. In contrast to the first algorithm (which solves for a given M_p and R_p), a single integration in radius results in the desired solution of R_p, for a given M_p, α, and β. In other words there is no iteration required, making this algorithm much more efficient.

2.4. Instructions for Downloading and Using the Code

The code is based in MATLAB and can be downloaded.3 If using this computer code please cite this paper and also Seager et al. (2007).

We have made two different codes available. The codes have the same output, but the first is based on a differential equation solver (§ 2.2) and the second code is based on interpolation of the large database (§ 2.3). For the codes, the planet mass must be in the range 0.5–20 M_⊕. The inputs to the codes are: the planet mass in Earth masses (M_p), the planet mass uncertainty in Earth masses (σ_Mp), the planet radius in Earth radii (R_p), and the planet radius uncertainty in Earth radii (σ_Rp). The values σ_Mp = 0 and σ_Rp = 0 are allowed. If the combination of input values M_p and R_p are unphysical, the code will return an error.

ExoterDE(M_p,σ_Mp,R_p,σ_Rp). This code solves the two differential equations described in § 2.1. This code consists of three subroutines (each of which must be downloaded) that are automatically called by the above command. The first subroutine is the differential equation solver, which also reads the equations of state. The second subroutine contains the actual differential equations. The third subroutine plots the ternary diagrams; this subroutine calls a ternary diagram plotting routine 4 that plots a single line for each of the 1, 2, and 3 σ contour lines (see § 4.1). An example from this code is shown in Figure 4.

ExoterDB(M_p,σ_Mp,R_p,σ_Rp). This code reads in the database of M_p, R_p, and fractional composition (α and β). The output is a ternary diagram, shaded throughout the 1, 2, and 3 σ contour curves. This subroutine uses the same ternary diagram plotting routine as described above.

The differential equation solver ExoterDE is much slower than the database extracter ExoterDB. In principle, ExoterDE is more accurate than ExoterDB.

3.1. 2D Cartesian Diagram

For a given M_p and R_p we want to know the interior composition of the relative mass fraction of the three components. There are three variables we have solved for α, β, and γ, but only two of them are independent (since γ = 1 - α - β). Therefore points on a 2D diagram can describe all the possible combinations of α, β, and γ for a given M_p and R_p. We show such a solution in Figure 2. We note that 0 ≤ α ≤ 1, 0 ≤ β ≤ 1, and α + β ≤ 1, and therefore, not every point in the 2D plane will correspond to a set of α, β, and γ. Only the points that are in a right-angled triangular region will represent the set of allowed solutions.

3.2. Ternary Diagram

Ternary diagrams to describe the interior composition of exoplanets were introduced by Valencia et al. (2007). In a ternary diagram, α, β, and γ are each one axis of an equilateral triangle. Although γ is extraneous, the ternary diagram is useful because it is more intuitive to see the three components of the planet interior (in a symmetric way) compared to a 2D Cartesian diagram with only two of the components. Figure 3 shows how to read a ternary diagram.

3.3. Relationship Between the 3D and 2D Cartesian Diagrams and the Ternary Diagram

To explain the full origin of a curve on the ternary diagram we start with the 3D Cartesian diagram with all solutions of M_p, R_p, and composition (in terms of α and β), as shown in Figure 1. We take an isoradius and isomass surface as shown in Figure 1. As an example, in Figure 1 the red surface is the isoradius surface of R_p = 1.7 R_⊕ and one of the blue colored planes is the isomass surface of M_p = 4 M_⊕. These two surfaces intersect each other and result in a curve. This curve can be projected vertically to the x–y plane, which is the isomass plane. We therefore have a (isoradius and isomass) curve on the isomass plane. This curve is shown in Figure 2 in a Cartesian diagram.

The Cartesian and ternary diagrams are two different ways to represent the same information. There exists a linear coordinate transformation between the two. That means if a function is a straight line appearing in the 2D Cartesian diagram, it will still be a straight line in the ternary diagram.

The transformation from 2D Cartesian coordinates to the ternary diagram coordinates is Here x and y are the coordinates of a point in a 2D Cartesian diagram, the x_ternary and y_ternary are the coordinates of the point in ternary diagram in the Cartesian grid variables. Figure 3 shows the same R_p = 1.7 R_⊕, M_p = 4 M_⊕ curve represented by a ternary diagram.

4.1. Observational Uncertainties

Real planet mass and radius measurements have uncertainties. The planet mass and radius uncertainties are typically 5 to 10% (e.g., Selsis et al. 2007), and even smaller for the most favorable targets. We now present examples of ternary diagrams that include the mass and radius uncertainties.

We consider uncertainties of 1, 2, and 3 standard deviations (σ) from the measured value. See Figures 4 and 5. In more detail, the uncertainty in composition on the ternary diagram is where comp refers to composition and compM and compR refer to composition uncertainties caused by the planet mass and radius uncertainty, respectively. Here we have assumed that the uncertainties in mass and radius are independent from each other and have assumed the linearity of the superposition of small uncertainties.

Figure 4 shows a planet with M_p = 10 ± 0.5 M_⊕ and R_p = 2 ± 0.1 M_⊕. We see that taking the 3 σ limit, almost the entire ternary diagram is filled. In other words, for a 5% 3 σ (i.e., 15%) uncertainty on the planet mass and radius, the interior composition in terms of fractional composition of iron, silicates, and water cannot be determined. The reason this example fills almost the whole ternary diagram is that a planet with 10 M_⊕ and 2 R_⊕ has an average density in between two extreme cases (purely iron or purely water). Therefore, a large variety of different combinations of iron, silicate, and water can result in a similar M_p and R_p. Even taking a 1 σ uncertainty of the planet mass and radius, the uncertainty in internal composition is large.

We note that an uncertainty in R_p has more of an effect on the uncertainty in the interior composition than an uncertainty in M_p. This is because the planet’s average density ρ ∼ M_p/R³. Considering error propagation, the uncertainty in radius has a 3 times larger effect on the uncertainty in average density than does the mass uncertainty.

We show ternary diagrams for planets with various masses, radii, and 5% fractional uncertainty in Figure 5. Only solutions in part of the ternary diagram are allowed, despite considering the 3 σ range. In Figure 5b, the upper 3 σ boundary is absent because it goes below the lowest allowed density of the 2 M_⊕ and 1.5 R_⊕ planet and is, thus, unphysical. In Figure 5c we see the opposite case, where the lower 3 σ boundary is absent because it goes above the highest planet density allowed and is, thus, unphysical.

A ternary diagram for a fixed planet mass and radius that includes observational uncertainties is one of the primary outcomes of this paper.

4.2. Model Uncertainties

The model and computer code we present assumes a differentiated planet composed of an iron core, a silicate mantle, and a water–ice outer layer. The division into three major materials is based on the point that differences in the densities of iron, silicate, and water from each other are greater than any minor compositional variant of each individual material. The model neglects phase variation and temperatures, which, as argued in Seager et al. (2007), have little effect on the total planet radius (to an uncertainty of about 1 to 3% uncertainty in planet radius, decreasing with increasing planet mass).

Low-pressure phase changes (at < 10 GPa) are not important for a planet’s radius because for plausible planet compositions most of the mass is at high pressure. For high pressure phase changes we expect the associated correction to the equation of state (and hence derived planet radii) to be small because at high pressure the importance of chemical bonding patterns to the equation of state drops.

Regarding temperature, at low pressures (≲10 GPa) in the outer planetary layers, the crystal lattice structure dominates the material’s density and the thermal vibration contribution to the density are small in comparison. At high pressures the thermal pressure contribution to the EOS is small because the close-packed nature of the materials prevents structural changes from thermal pressure contributions.

Although the code can model radii for planets in the mass range 0.5 to 20 M_⊕, the model is more accurate for planets above a few Earth masses (Seager et al. 2007).

The model also neglects variation in composition, such as a light element in the iron core as Earth and Mercury are believed to have. The model also omits other impurities in the mantle and water layer, including iron in the mantle. Molten cores have also been omitted. At the present time, these model uncertainties are expected to have an effect on the planet radius much less than the ∼5% radius observational uncertainty.

For all of the above reasons, we therefore argue that for the present time the observational uncertainties dominate the model uncertainties; the model presented here is adequate for an estimate of planet bulk composition. In any event, the main results of our work are described in the following subsections despite any model uncertainties.

It is possible to rule out parts of the ternary diagram as being physically unplausible (e.g., Valencia et al. 2007). This is based on the initial composition of the protoplanetary nebula and on planet differentiation. For example, a pure iron planet is unlikely to exist, because removing all of the mantle would be difficult. A pure water planet is also unlikely to exist. Where water–ice forms, so do silicate-rich and iron-rich materials, making planet accretion of pure water unlikely. We prefer to leave the omission of parts of the ternary diagram to users of the code, because in exoplanets surprising exceptions to the “rules” of planet characteristics are not uncommon.

4.3. Spacing, Shape, Direction, and Rotation of Curves on the Ternary Diagram

We now turn to a discussion of the spacing, direction, and shape of the curves in the ternary diagrams. A quantitative and qualitative description of these is a main point of this paper. We emphasize that each curve in the diagrams shown in Figure 5 represents a different mass and radius. The curves to the lower right are more dense, as they have a higher mass and lower radius than the curves moving to the upper left.

All behavior results from the equations of state of the materials, and, in some cases, how they behave differently under pressure.

We begin with an equation that we use repeatedly in this section. We consider the simplified case that the planet core has a uniform density, where is the average density of Fe in the core, is the average density of silicate in the mantle, and is the average density of water in the outer water layer. We then have where α + β + γ = 1.

4.4. Shape and Direction of an Isoradius Surface in our 3D Representation

We return to the 3D representation of the relationship between mass, radius, and composition shown in Figure 1. An isoradius surface (red oblique surface) is shown for R_p = 1.7 R_⊕. An isomass surface (one of the blue flat planes) is shown for M_p = 4 M_⊕. At the bottom tip of the isoradius surface the silicate mass fraction = iron mass fraction = 0, and the planet is composed of 100% water. This is the minimum mass for this radius. The isoradius surface also has a maximum mass reached by a composition of 100% iron.

For the same radius, if either the silicate or iron mass fraction is increased, the planet mass must also increase. This is because iron and silicate are denser than water–ice. We further note that for an increase in the iron mass fraction, the mass of the planet must increase more steeply than for an increase in the silicate mass fraction. This is seen by the different length and shape of the “edges” of the isomass surface in the iron-mass plane and silicate-mass plane in Figure 1.

We now show that the shape of the isoradius surface is concave. The isoradius surface can be considered as the isovolume surface, where the volume is the sum of the core, the silicate mantle, and the water crust. To calculate the total mass of a point on the isoradius curve, we use equation (6), and since R_p is constant on this surface we can rewrite this equation (with C a constant) as, The term in brackets is a linear function and its inverse is hyperbolic. For example, if we let α = 0 (no iron), then we have whereBoth terms in the square brackets are positive and 0 ≤ β ≤ 1. This is in the form where and We have C₁ > C₂ > 0. This is the form of a concave hyperbola, because, taking the first derivative, we find The slope increases as β increases because C₁ - C₂β decreases but is always greater than zero.

Although we have made the simplification that the average densities of each layer remain constant, our qualitative description holds because as long as the planet is differentiated into layers of increasing density toward the planet center.

A similar argument shows that the total iron fraction versus the total mass is also a hyperbola, making the whole isoradius surface concave.