A HARDCORE model for constraining an exoplanet's core size
MMNRAS , 1–9 (2017) Preprint 6 November 2018 Compiled using MNRAS L A TEX style file v3.0 A HARDCORE model for constraining an exoplanet’s core size
Gabrielle Suissa (cid:63) , Jingjing Chen & David Kipping Dept. of Astronomy, Columbia University, 550 W 120th Street, New York NY 10027
Accepted 2018 February 5. Received 2017 December 29; in original form 2017 September 5
ABSTRACT
The interior structure of an exoplanet is hidden from direct view yet likely plays acrucial role in influencing the habitability of Earth analogs. Inferences of the interiorstructure are impeded by a fundamental degeneracy that exists between any modelcomprising of more than two layers and observations constraining just two bulk pa-rameters: mass and radius. In this work, we show that although the inverse problemis indeed degenerate, there exists two boundary conditions that enables one to in-fer the minimum and maximum core radius fraction,
CRF min & CRF max . These holdtrue even for planets with light volatile envelopes, but require the planet to be fullydifferentiated and that layers denser than iron are forbidden. With both bounds inhand, a marginal CRF can also be inferred by sampling inbetween. After validatingon the Earth, we apply our method to Kepler-36b and measure
CRF min = ( . ± . ), CRF max = ( . ± . ) and CRF marg = ( . ± . ), broadly consistent with the Earth’strue CRF value of 0.55. We apply our method to a suite of hypothetical measurementsof synthetic planets to serve as a sensitivity analysis. We find that CRF min & CRF max have recovered uncertainties proportional to the relative error on the planetary den-sity, but
CRF marg saturates to between 0.03 to 0.16 once ( ∆ ρ / ρ ) drops below 1-2%.This implies that mass and radius alone cannot provide any better constraints oninternal composition once bulk density constraints hit around a percent, providing aclear target for observers. Key words: planets and satellites: interiors — planets and satellites: terrestrialplanets — planets and satellites: composition
Despite recent strides in our ability to characterize exoplan-ets (Kaltenegger 2017), knowledge regarding the internalstructure of distant worlds remains almost entirely lacking(Spiegel et al. 2013; Baraffe et al. 2014). Unlike the searchfor exoplanetary atmospheres (Burrows 2014), moons (Kip-ping 2014) or tomography (McTier & Kipping 2017), ourremote observations do not have direct access to that whichwe seek to infer - the planet’s interior. The habitability ofan Earth-like planet, in particular via the likelihood of platetectonics, is likely strongly influenced by the internal struc-ture (Noak et al. 2014) and thus the community is stronglymotivated to infer what lies beneath, as part of the broadergoal of understanding our own planet’s uniqueness.In general, the only information we have about an exo-planet which is directly affected by internal structure is thebulk mass and radius of the planet . Aside from this, wehighlight that there are some special cases where additional (cid:63) E-mail: [email protected] Quantities such as bulk density and surface gravity are of coursederivative of mass and radius information about the planetary interior can become avail-able. For example, Kaltenegger (2010) argue that volcanismand planetary outgassing could be detectable using atmo-spheric characterization techniques. Certain dynamical con-figurations of planetary systems, such as tidal fixed points(Batygin et al. 2009) for example, can also enable inferenceof the planetary tidal properties, which in turn constrainsinternal composition (see also Kramm et al. 2012). Finally,direct measurements of oblateness may also provide con-straints on internal structure (Seager & Hui 2002; Carter &Winn 2010; Zhu et al. 2014).Whilst there is some hope of identifying outgassingof exoplanets, providing clues to the mantle composition(Kaltenegger 2010), and measuring tidal dissipation con-stants in special cases (Batygin et al. 2009), full structureinference will likely be limited to indirect methods based ontheoretical models. In this approach, one takes the basic ob-servables we do have access to, in particular planetary massand radius, and compares them to theoretical models in aneffort to find families of compatible solutions. Since theo-retical models depend on more than just two parameters,accounting for factors such as chemical composition (Valen-cia et al. 2006; Seager et al. 2007), ultraviolet environment c (cid:13) a r X i v : . [ a s t r o - ph . E P ] A p r Suissa et al. (Lopez & Fortney 2013; Batygin & Stevenson 2013) and age(Fortney et al. 2011), the problem is degenerate, in a generalsense.Since we do not have direct access to the interior ofexoplanets, their interior structure is generally modelled byassuming several key chemical constituents. In the case ofsolid exoplanets, extrapolation from the Solar System im-plies that they should be comprised of three primary chemi-cal ingredients, namely water, H O, enstatite, MgSiO , andiron, Fe (Valencia et al. 2006). If we assume the planet is notyoung and has thus become fully differentiated, the equa-tions of state of these three layers can be solved to providetheoretical estimates of the mass and radius of solid bodies(Zeng & Sasselov 2013). A fourth layer describing a lightvolatile envelope can be placed on top to capture the behav-ior of mini-Neptunes, where the light envelope is assumedto have negligible relative mass and thus only affects thebulk radius and not the mass (Kipping et al. 2013; Lopez &Fortney 2014; Wolfgang et al. 2016).These four constituents can be combined in multipleways to re-create the same mass and radius. Even in the ab-sence of a volatile envelope this degeneracy persists, leadingto the common use of ternary diagrams to illustrate theirsymplectic yet degenerate loci of solutions (e.g. see Seageret al. 2007). This degeneracy is a major barrier to inferringunique solutions for planetary interiors, leading authors toeither switch out to simpler and non-degenerate two-layermodels (e.g. Zeng et al. 2016) or adding a chemical proxyfrom the parent star (e.g. Dorn et al. 2017) to break the de-generacy. Whilst these are both certainly promising avenuesfor tackling interior inference, in this work we focus on athird approach based on boundary conditions.The possibility of exploiting boundary conditions wasfirst highlighted in Kipping et al. (2013), where the authorsfocused on the concept of “minimum atmospheric height”.The method works by first predicting the maximum allowedradius of a planet without any extended envelope given itsmeasured mass. This atmosphere-less planet is typically as-sumed to be a pure water/icy body, for which detailed mod-els are widely available. If the observed radius exceeds thismaximum limit, then some finite volume of atmosphere mustsit on top of the planetary interior, and the difference in radiirepresents the “minimum atmospheric height” (MAH). Theapproach therefore formally describes a key boundary con-dition of a general four-layer exoplanet.In this work, we explore the other extreme, asking thequestion under what conditions would an observed mass andradius definitively prove some finite iron-core must exist,and what is the minimum radius fraction that the core mustcomprise? Going further, we argue that the maximum coreradius fraction is another boundary condition in the problemand thus can be derived to provide a complete bounding boxof a planet’s core size in a general four-layer model frame-work.We introduce the concept of the minimum core size inSection 2, as well as our fast parametric model to interpolatethe Zeng & Sasselov (2013) grid models. Section 3 discussesour approach to inverting the relation to solve for core radiusfraction directly, as well as the much more straightforwardmethod for inferring the maximum core size. In Section 4,we demonstrate the approach on both synthetic and realexoplanets, with special attention to sensitivity. Finally, we M g S i O F e H . . . . . . . . . . . . . . . . . . . . . . ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ CRF=0.43CRF=0.60CRF=0.75 CRF>0.43 ∀ compositions } } CRF=0.77 ⊕ ⊕ CRF<0.77 ∀ compositions and Fe + H/He
Figure 1.
Ternary diagram of a three-layer interior structuremodel for a solid planet. All points along the red line lead to a M = M (cid:12) and R = R (cid:12) planet. Although indistinguishable from eachother with current observations, all points satisfy having a core-radius fraction exceeding 43%, a boundary condition we exploitin this work. The largest iron core size allowed is depicted by thelowest sphere, where the volatile envelope contributes negligiblemass. discuss the anticipated value of this work, as well its limita-tions, in Section 5. Exoplanets are generally expected to display a diverse rangeof physical characteristics, owing to their presumably dis-tinct formation mechanisms, chemical environments and his-tories, as three examples. Even if the mass and radius of anexoplanet were known to infinite precision, and the bodywas known to be definitively solid, these two observed pa-rameters are insufficient to provide a unique solution forthe relative fractions of water, silicate and iron typically as-sumed to represent the major constituents of solid planets(Kipping et al. 2013). In other words, mass and radius alonecannot confidently reveal an exoplanet’s CRF or CMF (coreradius fraction or core mass fraction, respectively).As a concrete example, a planet composed of water andiron can have the same mass and radius as an iron-siliconplanet, and thus have very different CRFs (as illustrated inFigure 1).As touched on in Section 1, four-layer theoretical mod-els of solid planets are degenerate for a single mass-radiusobservation. However, across the suite of loci able to serve asviable solutions, there exists a boundary condition when thecomposition is pure silicate and iron. At this point, the CRFtakes the smallest value out of all possible models, since thesecond-layer (the mantle) is now as heavy as it can be, be-ing pure silicate (the second densest material). Therefore,for any given mass-radius pair, we can solve for the corre-sponding CRF of a pure silicate-iron model (which is not adegenerate problem) and define that this CRF must equalthe minimum CRF,
CRF min , allowed across all models.Similarly, another boundary condition we can exploit isto consider the maximum allowed core size. As depicted inFigure 1, this occurs when all of the mass is located within
MNRAS , 1–9 (2017)
HARDCORE model for constraining an exoplanet’s core size a pure iron core, padded by a second layer of a light volatileenvelope. The core can’t possibly exceed this fractional sizeelse the mass would be incompatible with the observed value.Note that although we refer to CRF rather than CMFhere and in what follows, once armed with a CRF, the massand radius, it is easy to convert back to CMF. Note too thathere and throughout in what follows, we refer to the CRFstrictly in terms of an iron core. Although technically weacknowledge that a water-silicate body could be describedas having a finite sized silicate core, that core would notqualify as being a “core” in this work.Before continuing, we highlight some key assumptionsof our model, for the sake of transparency: (cid:4) The planet is fully differentiated and is not recentlyformed. (cid:4)
The outer volatile envelope has insufficient massand thus gravitational pressure to significantly affect theequation-of-state of the inner layers. (cid:4)
The densest core permitted is that of iron i.e. heavy-element (e.g. uranium) cores are forbidden. (cid:4)
We have accurate models for a planet’s mass and radiusgiven a particular compositional mix.We stress that what has been described thus far includesthe possibility of a light volatile envelope, and is not limitedto some special case of two- or three-layer conditions, asdiscussed earlier.Under these assumptions, the limiting core radii frac-tions should be determinable, although violating any of theassumptions listed above would invalidate our argument. Anobvious one is that the theoretical models used are invalid orinaccurate, for example because their assumed equations ofstate are wrong. In this work, we primarily use the Zeng &Sasselov (2013) model but we point out that should a userbelieve an alternative model to be superior, it is straight-forward to reproduce the methods described in this paperusing the model of their preference. The actual existence ofa boundary condition remains true.Another more serious flaw would be if the planetarybody in question has a significant mineral fraction based onsome heavier element than iron, for example a uranium core.Such a body could feasibly have a significantly smaller corethan that derived using our approach. If evidence for suchcores emerges in the coming years, then we advice againstusers employing the model described in this work.The remaining two assumptions, a differentiated, non-young planet and a light volatile envelope, mean that youngsystems are not suitable and gas giants would not be either.In general then, the model described is expected to be validfor most planets smaller than mini-Neptunes.
CRF min
For any combination of mass and radius, we need to beable to predict what the corresponding CRF would be for asilicate-iron two-layer model, in order to determine
CRF min .Inferring
CRF max is far simpler and is briefly explained laterin Section 2.3. In what follows, we use the models of Zeng& Sasselov (2013), which are made available as a regulargrid of theoretical points. Whilst we could simply performa nearest neighbor look-up, this is unsatisfactory since ourprecision will be limited to the grid spacing and resulting posteriors would be rasterized to the same grid resolution.Instead, we seek a means to perform an interpolation of thegrid.The first successful literature interpolation of the Zeng& Sasselov (2013) models comes from Kipping et al. (2013),who found that for a specific fixed CRF, each of the varioustwo-layer models of Zeng & Sasselov (2013) are very well-approximated by a seventh-order polynomial of radius withrespect to the logarithm of mass, given by RR ⊕ = ∑ i = a i ( CRF ) × log (cid:16) MM ⊕ (cid:17) i (1)Temperature does not feature in this expression, as Zeng& Sasselov (2013) find its effects on the density profile aresecondary compared to pressure and can be safely ignored.Equation 1 is attractive since it is parametric, linear (andthus can be trained using linear least squares) and extremelyfast to execute as an interpolative model once the coefficientshave been assigned. Accordingly, we are motivated to pursuea similar strategy in what follows.Consider first the mass-radius relation for a 100% sili-cate planet. Plotting the radius against the logarithm of themass indeed reveals a series of points that are well describedby a seventh-order polynomial, as shown in Figure 2. Wefound that the range for which this interpolation works bestis M > . M ⊕ and thus we set this as a truncation pointduring training. Note that any planet that lies beneath thispolynomial curve must have an iron core.However, our goal is to not only determine whether ornot a planet has an iron core, but also quantify the mini-mum CRF. To accomplish this, we trained a suite of seventh-order polynomials on models with varying CRFs assumingthe two-layer iron-silicate models of Zeng & Sasselov (2013).We varied the CRF from 0 to 1 in 0.025 steps and performa linear least squares regression at each step. To illustratethis, a gradient of interpolations of the mass-radius relationfor the CRFs between 0 and 1 are shown in Figure 3.In order to have a general model, we are interested inparameterizing the CRF variable to understand the relationof the polynomials shown in Figure 3. To do this, we allowthe coefficients of the polynomials to be polynomials them-selves (but with respect to CRF rather than logarithmicmass), such that a i ( CRF ) = M i ∑ j = b i , j CRF j , (2)where M i is the polynomial order of the i th coefficientand b i , j is the j th “sub-coefficient” of the the i th coefficient.As an example, in the case of the a coefficient, wefirst graphed the coefficient for 40 steps as a function ofthe corresponding CRF (see top-left panel of Figure 4). Itwas immediately apparent that the points followed a smoothfunction that could be stably approximated by another poly-nomial. The polynomial-order was initially third-order andthen we stepped through until the polynomial appeared togo through almost all of the data, ranging from fifth to tenthorder. The same process was repeated for the remainingeight coefficients. The resulting functions are presented inFigure 4, with the coefficients tabulated in Table 1. MNRAS000
CRF max is far simpler and is briefly explained laterin Section 2.3. In what follows, we use the models of Zeng& Sasselov (2013), which are made available as a regulargrid of theoretical points. Whilst we could simply performa nearest neighbor look-up, this is unsatisfactory since ourprecision will be limited to the grid spacing and resulting posteriors would be rasterized to the same grid resolution.Instead, we seek a means to perform an interpolation of thegrid.The first successful literature interpolation of the Zeng& Sasselov (2013) models comes from Kipping et al. (2013),who found that for a specific fixed CRF, each of the varioustwo-layer models of Zeng & Sasselov (2013) are very well-approximated by a seventh-order polynomial of radius withrespect to the logarithm of mass, given by RR ⊕ = ∑ i = a i ( CRF ) × log (cid:16) MM ⊕ (cid:17) i (1)Temperature does not feature in this expression, as Zeng& Sasselov (2013) find its effects on the density profile aresecondary compared to pressure and can be safely ignored.Equation 1 is attractive since it is parametric, linear (andthus can be trained using linear least squares) and extremelyfast to execute as an interpolative model once the coefficientshave been assigned. Accordingly, we are motivated to pursuea similar strategy in what follows.Consider first the mass-radius relation for a 100% sili-cate planet. Plotting the radius against the logarithm of themass indeed reveals a series of points that are well describedby a seventh-order polynomial, as shown in Figure 2. Wefound that the range for which this interpolation works bestis M > . M ⊕ and thus we set this as a truncation pointduring training. Note that any planet that lies beneath thispolynomial curve must have an iron core.However, our goal is to not only determine whether ornot a planet has an iron core, but also quantify the mini-mum CRF. To accomplish this, we trained a suite of seventh-order polynomials on models with varying CRFs assumingthe two-layer iron-silicate models of Zeng & Sasselov (2013).We varied the CRF from 0 to 1 in 0.025 steps and performa linear least squares regression at each step. To illustratethis, a gradient of interpolations of the mass-radius relationfor the CRFs between 0 and 1 are shown in Figure 3.In order to have a general model, we are interested inparameterizing the CRF variable to understand the relationof the polynomials shown in Figure 3. To do this, we allowthe coefficients of the polynomials to be polynomials them-selves (but with respect to CRF rather than logarithmicmass), such that a i ( CRF ) = M i ∑ j = b i , j CRF j , (2)where M i is the polynomial order of the i th coefficientand b i , j is the j th “sub-coefficient” of the the i th coefficient.As an example, in the case of the a coefficient, wefirst graphed the coefficient for 40 steps as a function ofthe corresponding CRF (see top-left panel of Figure 4). Itwas immediately apparent that the points followed a smoothfunction that could be stably approximated by another poly-nomial. The polynomial-order was initially third-order andthen we stepped through until the polynomial appeared togo through almost all of the data, ranging from fifth to tenthorder. The same process was repeated for the remainingeight coefficients. The resulting functions are presented inFigure 4, with the coefficients tabulated in Table 1. MNRAS000 , 1–9 (2017)
Suissa et al.
Table 1.
Coefficients for the terms in Equation 2. These are implemented in our package
HARDCORE available at this URL.j b , j b , j b , j b , j b , j b , j b , j b , j
10 100010010.1
Figure 2.
Theoretical mass-radius grid points of a pure silicateplanet from Zeng & Sasselov (2013), shown in blue. For compar-ison, we show our th order polynomial fit. We only recommendthis interpolation for masses above 0.1 M ⊕ . To validate our model, we re-trained our model of allof the original training data Zeng & Sasselov (2013) butomitting a random single datum each time, serving as ahold-out validation point. We then computed the relativedifference between the prediction for that point using there-trained model, and the actual value. Repeating for random hold-out points, we find that the mean error of ourmodel is 0.045% and the maximum error is 0.24%.To assist the community, we make our model, which wedub HARDCORE , publicly available at this URL.
CRF max
Determining
CRF max is far more straight-forward than
CRF min . One may simply take the 100% iron mass-radiusmodels, in our case from Zeng & Sasselov (2013), and di-rectly compute the expected radius of a pure iron planetgiven an observed mass, R iron ( M obs ) .The maximum core radius fraction is then easily com-puted as CRF max = R iron ( M obs ) / R obs . In practice, this inver-sion is far simpler than CRF min since we need not interpolateacross intermediate mixtures of iron and silicate composi- core radius fraction0 1 [ M ⊕ ] r a d i u s [ R ⊕ ] Figure 3.
Interpolated theoretical mass-radius relations for asilicate-iron two-layer solid planet for various core radius fractions(CRFs), based off Zeng & Sasselov (2013). All interpolations forCRFs between 0 and 1 are seventh-order polynomials. We are thenmotivated to describe the dependence of the polynomials withrespect to the CRF, by making the these coefficients polynomialfunctions themselves. tions, but rather cam simply directly solve for
CRF max if ananalytic expression for R iron ( M ) is available. This could beestimated using our interpolative model but is more directlyaccessible by simply fitting Equation 1 to the Zeng & Sas-selov (2013) grid for the specific points corresponding to apure iron composition, which were previously presented inthe final column of Table 1 of Kipping et al. (2013). Thus far we have described a method to solve for
CRF max but not for
CRF min . The silicate-iron interpolative model is aforward model in which we begin with knowledge of both theplanet’s mass and minimum core radius fraction and com-pute the corresponding radius. In practice, however, we areinterested in the inverse model, where we wish to determine
CRF min from the mass and radius.The nested coefficient structure makes the problem non-linear with respect to
CRF min , yet it is one dimensional and
MNRAS , 1–9 (2017)
HARDCORE model for constraining an exoplanet’s core size a a Figure 4.
Interpolated functions for the coefficients of seventh-order polynomials. Given a specific CRF, one can use this parameterizationto solve for all eight coefficients, to then use for a seventh-order polynomial. This seventh order polynomial can then describe the mass-radius function corresponding to all CRFs, not just the ones tabulated by Zeng & Sasselov (2013). found to be unimodal in practice. For these reasons, it isamenable to a large number of possible optimization algo-rithms, but in what follows we adopted Newton’s method,since we are able to directly differentiate our functionsthanks to their parametric nature.Specifically, in our implementation we minimize the fol-lowing cost function, J , with respect to one degree of free-dom, the CRF : J = ( R Fe − Si ( CRF; M obs ) − R obs ) , (3)where M obs and R obs are the observed mass and radiusof the planet, and R Fe − Si ( CRF; M ) is the radius of a two-layer iron + silicate model with core radius fraction CRF andmass M . The latter function is determined using the smoothparametric interpolation model described in Section 2.2. Ourinversion algorithm, starts at an initial guess of CRF min = . and then iterates by computing CRF i + = CRF i − J ( CRF i )[ d J / dCRF ]( CRF i ) . (4)To improve speed and stability, we impose a check asto whether R obs is below that of a pure iron planet of mass M obs , in which case we fix CRF min = , or if the radius exceedsthat of an pure silicate planet of mass M obs (where again weuse the interpolative model of Kipping et al. 2013), in whichcase we fix CRF min = . Let us use the Earth itself as an example of our method. Wetook a 1 M ⊕ and R ⊕ planet and used the methods describedearlier to solve for CRF min and
CRF max , giving
CRF min = . and CRF max = . .In reality, the Earth is not perfectly described by the Zeng & Sasselov (2013) model and the core in particularis only ∼ % iron, with nickel and other heavy elementscomprising the rest. The mantle-core boundary occurs ata radius of 3480 km relative to the Earth’s mean radius of6371 km (Dziewonski & Anderson 1981), meaning that itsactual CRF = . . Accordingly, our CRF bounds correctlybracket the true solution, as expected.We go further by treating these limits as being boundson a prior distribution for CRF. Adopting the least infor-mative continuous distribution for a parameter constrainedby only two limits corresponds to a uniform distribution.Sampling from said distribution yields a marginalized CRFof CRF marg = . ± . , which is again fully compatiblewith the true value.To test our inversions in a probabilistic sense, wedecided to create a mock posterior distribution of anEarth-like planet where M ∼ N [ . M ⊕ , . M ⊕ ] and R ∼ N [ . R ⊕ , . R ⊕ ] (we also apply a truncation to the dis-tributions at zero to prevent negative masses/radii). This isclearly an optimistic assumption but a more detailed inves-tigation of sensitivity for different relative errors is tack-led later in Section 3.3. Generating samples, we in-verted each sample as described earlier to produce a pos-terior for CRF min and
CRF max . Our experiment returns near-Gaussian like distributions for both terms with a meanand standard deviation given by
CRF min = ( . ± . ) and CRF max = ( . ± . ) . This establishes that the inver-sions are stable against perturbations around physical solu-tions.As a brief aside, we argued earlier in Section 2.1 thatthe principle of exploiting the boundary condition of the-oretical models to infer CRF min does not explicitly requiresolid planets and works for mini-Neptunes too. To demon-strate this point with a specific example, let us return tothe earlier thought experiment of the Earth as a gaseous
MNRAS000
MNRAS000 , 1–9 (2017)
Suissa et al. planet consisting of a solid iron core surrounded by a lightH/He envelope. We consider that the mass and radius of theplanet remain the same as the real Earth, and that the en-velope significantly influences the radius but has negligiblemass. Using the 100% iron-model of Zeng & Sasselov (2013)and the corresponding th order polynomial interpolation ofKipping et al. (2013), we estimate that a 1 M ⊕ iron corewould have a radius of 0.77 R ⊕ and therefore the remain-ing 0.23 R ⊕ is given by the light, H/He envelope as depictedearlier in Figure 1. Accordingly, such a body would have acore radius fraction exceeding our inferred minimum valueof CRF (which was
CRF > . ), which is expected and self-consistent with our definition of CRF min .We highlight that the counter-example described aboveis highly unphysical though; R ⊕ are not expected to retainsignificant volatile envelopes in mature systems, both froma theoretical perspective (Lopez & Fortney 2014; Owen &Wu 2017) and an observational one (Rogers 2015; Chen &Kipping 2017; Fulton et al. 2017). A basic and important question to ask is what kind of pre-cisions on a planet’s mass and radius lead to meaningfulconstraints on
CRF min ? In other words, what is the corre-spondence we might expect between { ( ∆ M / M ) , ( ∆ R / R ) } and ( ∆ CRF min / CRF min ) ? This is key for designing future obser-vational surveys, where primary science objectives may cen-ter around inferring internal compositions. To investigatethis, we repeated the retrieval experiment described in Sec-tion 3.2, but varied the fractional error on mass and radiusaway from the fixed 1% value previously assumed.In total, we generated = experiments, where foreach one we generated a new mock posterior of mass-radius samples, which was then converted into a posteriorof CRF min and
CRF max . For each experiment, we record thestandard deviation of the resulting posteriors as ∆ CRF min and ∆ CRF max . The errors on the mass and radius were inde-pendently varied with a fractional error given by x , where x was varied across a regular grid from -4 to in 0.05 steps,giving 81 grid points in each dimension, and thus 6561 acrossboth. In all experiments, the underlying mass and radiusposteriors are generated assuming a mean of µ M = M ⊕ and µ R = R ⊕ .Figure 5 displays the results of this effort for each com-bination of mass and radius error. Given the shape of darkerareas of the color plot, one can see that radius is the dom-inant constraint, and that for the same fractional error onmass and radius, it is the radius term which mostly stronglyconstrains CRF min . For example, we find that in order toobtain a precision of 10% on
CRF min , we require a measure-ment on the mass better than 11% and a measurement onthe radius better than 3%.The ratio of these two numbers, close to three-to-one,led to us hypothesize that density was the underlying drivingterm. This can be seen by calculating error on density forindependent mass and radius via ∆ ρρ = ∆ MM + ∆ RR . (5)This is verified in the lower panels of Figure 5, where we find that although density doesn’t perfectly capture thedependency, it describes the vast majority of the variance.For precise densities ( (cid:46) %), the dependency is strictly lin-ear where we give the coefficients in the panels. This lineardependency breaks down as the errors grow, likely as a resultof the truncated normals used to generate the masses andradii becoming increasingly skewed and the finite support in-terval (zero to unity) of the CRF itself causing a saturationeffect.The marginalized CRF behaves quite different to theother two. The upper-central panel of Figure 5 alone looksfairly consistent with the previous, just with inflated errors.This is to be expected by the very act of marginalization.However, the bottom-central panel does not exhibit a simplelinear dependency, even at precise densities. In contrast, atprecise densities the marginalized CRF appears to saturateto ∼ %. This implies that no better than 10% precisioncan ever be obtained on the CRF using just mass and radiusalone.
Thus far, we have assumed a M = M ⊕ , R = R ⊕ planet. Inorder to generalize the scalings found, we decided to varythese inputs and repeat the entire process described above.We varied the mass from 1 to 10 Earth masses logarithmi-cally and the CRF from 0.2 to 0.8 uniformly, exploring over1000 different realizations. For the CRF max term, we find that ∆ CRF max (cid:39) α max (cid:32) ∆ ρρ (cid:33) (6)provides an excellent fit across all simulations, wherethe best-fitting value of the coefficient term ranged from . < α max < . . The relationship is sufficiently tightthat it is reasonable to simply adopt α max (cid:39) . as a generalrule of thumb. Repeating for the minimum limit on the CRFwe find that the function ∆ CRF min (cid:39) α min (cid:32) ∆ ρρ (cid:33) (7)again provides excellent fits, but now the coefficientterm varies dramatically from . < α min < across all simu-lations, or logarithmically a range of − . < log α min < . .The behaviour of α max appears to display a peculiar and pe-riodic dependency with respect to the dependent variables,CRF and mass. We were unable to identify a simple phys-ically motivated relation after substituting for terms suchas density and surface gravity, but were able to capture themost of the variance using an empirical approximate formulagiven by log α min (cid:39) − . + . β min , (8)where β min (cid:39)(cid:98) . ( CRF − log M − R [ . − M ] (cid:99)− . R [ ( log M − . )] − . + . M . (9) MNRAS , 1–9 (2017)
HARDCORE model for constraining an exoplanet’s core size -1.0 -1.5 -2.0 -2.5 -3.0 CRF min Δ CRF min ≃ α min ( Δρ / ρ ) Δ CRF min ≃ α min ( Δρ / ρ ) α min, ⊕ = α min, ⊕ = - - - - - - - - - - - - - - Δρ / ρ Δ C R F m i n -0.85 -0.90 -0.95 -1.00 CRF marg Δ CRF marg ≳ α marg, ⊕ Δ CRF marg ≳ α marg, ⊕ α marg, ⊕ = α marg, ⊕ = - - - - - - - - - - - - - Δρ / ρ Δ C R F m a r g -1.5 -2.0 -2.5 -3.0 -3.5 CRF max Δ CRF max ≃ α max ( Δρ / ρ ) Δ CRF max ≃ α max ( Δρ / ρ ) α max, ⊕ = α max, ⊕ = - - - - - - - - - - - - - - - Δρ / ρ Δ C R F m a x Figure 5.
Upper: Contour plots depicting the a-posteriori standard deviation on the minimum (left), marginalized (center) and maximum(right) CRF, as a function of the fractional errors on mass and radius. Lower: Re-parameterization of the above plots by combining themass and radius axes into a single density term, demonstrating a strong dependency in all cases. where R is a rounding function. We find the above func-tion has a robust (via the median absolute deviation) esti-mate of the RMS in the residuals of ∆ ( log α min ) = . .The marginal CRF plateau, which was 10% in the caseof the Earth, also varies in a non-trivial manner, rangingfrom 2.8% to 15.7% across our suite of simulations. Fortu-nately, the location of the plateau appears to be directlyrelated to α min , and thus we can use our earlier empiricalfunction to predict this term too with ∆ CRF marg (cid:39) α marg (10)where log α marg (cid:39) − . − . ( − . α min )+ . α marg . (11)We briefly comment that this saturation-behavior canbe understood as follows. With imprecise data, the poste-riors on the minimum and maximum limits will be broad,and so sampling a point between them will yield an evenbroader distribution. As the data become more precise, theposteriors distributions for CRF min and
CRF max converge to-wards sharp delta functions, but these two limits have noreason or expectation to be on top of one another. Accord-ingly, when we draw samples between them uniformly, wewill still get a broad distribution, albeit one less broad thanthat obtained when
CRF min and
CRF max were also broad.
To demonstrate the value of our prescription, we show herean example application to a real exoplanetary system. Asestablished in Section 3.3, precise measurement errors onboth planetary mass and radius are required in order to infer
CRF min at a meaningful level, certainly better than 10% onboth.One of the best laboratories for precise measurementsof planetary properties comes from dynamically interact-ing systems, particularly those in near mean motion reso-nances (MMRs) where transit timing variations (TTVs) maystrongly constrain planetary mass (Holman & Murray 2005;Agol et al. 2015). A good example comes from the Kepler-36 system, where two planets gravitationally perturb oneanother, enabling a measurement of both masses to betterthan 8% precision (Carter et al. 2012). Coupled with pre-cise planetary radii at precisions better than 2.5%, and thelow-radius, high-density nature of Kepler-36b in particular( R b = ( . ± . ) R ⊕ and ρ b = . + . − . g cm − ), Kepler-36 offers an excellent test case for our technique.To implement our code, we follow a similar procedureto that used in Section 3.2, except that we now use thereal mass-radius joint posterior distribution of Kepler-36bderived by Carter et al. (2012), which features indepen-dent samples. Kepler-36c was not used as the planet morelikely resembles a gas giant than a terrestrial body, as es-tablished using the MAH technique in Kipping et al. (2013), MNRAS000
CRF min at a meaningful level, certainly better than 10% onboth.One of the best laboratories for precise measurementsof planetary properties comes from dynamically interact-ing systems, particularly those in near mean motion reso-nances (MMRs) where transit timing variations (TTVs) maystrongly constrain planetary mass (Holman & Murray 2005;Agol et al. 2015). A good example comes from the Kepler-36 system, where two planets gravitationally perturb oneanother, enabling a measurement of both masses to betterthan 8% precision (Carter et al. 2012). Coupled with pre-cise planetary radii at precisions better than 2.5%, and thelow-radius, high-density nature of Kepler-36b in particular( R b = ( . ± . ) R ⊕ and ρ b = . + . − . g cm − ), Kepler-36 offers an excellent test case for our technique.To implement our code, we follow a similar procedureto that used in Section 3.2, except that we now use thereal mass-radius joint posterior distribution of Kepler-36bderived by Carter et al. (2012), which features indepen-dent samples. Kepler-36c was not used as the planet morelikely resembles a gas giant than a terrestrial body, as es-tablished using the MAH technique in Kipping et al. (2013), MNRAS000 , 1–9 (2017)
Suissa et al. who find that a volatile envelope comprises at least 36% ofKepler-36c’s bulk radius.The resulting posterior distributions on the
CRF min and
CRF max for Kepler-36b are shown in Figure 6. The
CRF min posterior indicates strong evidence for the presence of aniron core with
CRF min = . + . − . . We obtain a tighterconstraint on CRF max , as predicted by our earlier sensitiv-ity analysis, of
CRF max = . ± . , which is exactly whatwould be predicted from our scaling expression in Equation 6for a 10% error on the planet’s bulk density, which is indeedapproximately the reported value.By drawing random samples between CRF min and
CRF max on a sample by sample basis, we can constructthe
CRF marg posterior, revealing that
CRF marg = . ± . .Given that the fractional density error exceeds ∼ − . , thisprecision does not lie on the sensitivity plateau depictedin Figure 5, and therefore we predict that the precision onKepler-36b’s CRF marg could be improved. Since Carter et al.(2012) only considered ten quarters, it should be possible toconsiderably improve the uncertainties by re-analyzing theentire
Kepler time series.For comparison, the Earth’s
CRF is 0.55 but inverting
CRF marg for a synthetic Earth yields
CRF marg = . , show-ing a slight bias in the marginalization result (future workmay be able to correct for this by experimenting with differ-ent sampling schemes). Accordingly, our inference broadlyagrees with the conclusions of previous authors studyingKepler-36b’s interior- that the planet appears to be com-patible with having an Earth-like interior (Dorn et al. 2015;Unterborn et al. 2016; Owen & Morton 2016). It is worthemphasizing this is by no means guaranteed and one distin-guishing feature of our approach is that it infers an Earth-like composition rather than assuming it (e.g. see Lopez &Fortney 2013).Recall that our approach is actually relatively simplecompared to more sophisticated approaches, such as lever-aging stellar chemical proxies, strictly two-layer modeling orhierarchical Bayesian methods. All we’ve done is exploitedbasic boundary conditions in the problem. Despite this, wehave demonstrated its ability to reproduce the same infer-ence as other (often more complicated) models, giving usconfidence that the technique described here is a powerfuland effective tool for the community to constrain the interiorstructure of solid planets. In this work, we have presented a novel method for infer-ring the minimum and maximum iron core radius fraction(
CRF min & CRF max ) of an exoplanet using just it’s measuredmass and it’s radius. Building upon the earlier work of Kip-ping et al. (2013), we exploit two boundary conditions inthe theoretical models describing solid exoplanet interiors.Our method is valid under the assumptions of the specificunderlying theoretical model employed (we used the modelof Zeng & Sasselov 2013 for example) and the assumptionsthat the planet is differentiated, does not possess a high massvolatile envelope (although light envelopes are fine) and thatcores heavier than iron are not permitted.Although we used the theoretical models of Zeng & Sas-selov (2013) in this work, the method can be easily adapted
CRF = - + CRF = - + CRF max = - + CRF max = - + CRF min = - + CRF min = - + ( CRF ) p r o b . d e n s i t y Figure 6.
Posterior distribution of the minimum CRF (left),maximum CRF (right) and marginalized CRF (center) for Kepler-36b, based off the joint mass-radius posterior from Carter et al.(2012) and the model presented in this work. Posterior heightsnormalized to be equivalent. for whatever suite of theoretical models is preferred. Whilst
CRF max is reasonably straight-forward to compute (see Sec-tion 2.3),
CRF min has no simple solution. In order to facili-tate the inversion of the Zeng & Sasselov (2013) models, wehave developed a parametric interpolation of the silicate-iron two-layer model from that work, producing a functionable to predict a planet’s radius for any given mass and coreradius fraction. By iteratively inverting this expression withNewton’s method, we demonstrate that how the minimumCRF can be practically derived too (details are provided inSection 3).By drawing samples between the minimum and max-imum limits, a marginalized
CRF can be inferred. In thiswork, we have drawn samples assuming a uniform distri-bution although we have highlighted that this may not beoptimal, as it appears to impose a slight positive bias inthe results (see Section 4). As an applied example, we inferKepler-36b’s internal core size to be
CRF marg = . ± . ,compatible with, but slightly larger than, that of the Earth(see Figure 6).By inverting our model on a suite of planets with dif-fering measurement precisions, we have produced a detailedsensitivity map for CRF min , CRF max and
CRF marg (see Fig-ure 5). We find that both terms have standard deviationsproportional to the fractional density uncertainty, althoughcuriously
CRF marg saturates in precision once the density er-ror hits ∼ N > layer model will be degenerate ifinversion is attempted using just a mass and radius, and ourmethod circumvents this by exploiting boundary conditionsin the problem. However, it has been also proposed to ex-ploit stellar metallicity constraints from the parent star has MNRAS , 1–9 (2017)
HARDCORE model for constraining an exoplanet’s core size a third datum to break the degeneracy (Dorn et al. 2015),or simply adopt a two-layer model under certain reasonablecircumstances Zeng et al. (2016). We encourage the com-munity to view these approaches as complementary, eachoperating under different assumptions but ideally arrivingat consistent inferences.To aid the community it using our technique, we makeour code fully public ( HARDCORE available at this URL) andit is extremely fast to execute and perform inversions, pro-viding an efficient tool for others to use in their analyses.This work highlights the challenging nature of invert-ing exoplanet interiors, yet the potential to infer meaning-ful physical constraints. Precise masses in particular will becrucial in future work in this area, with surveys aiming todeliver sub m s − radial velocities and long-term transit tim-ing variations being particularly valuable resources for thisenterprise. ACKNOWLEDGMENTS
Thanks to members of the Cool Worlds Lab for helpful dis-cussions in preparing this manuscript.
REFERENCES
Agol, E., Steffen, J., Sari, R. & Clarkson, W. 2005, MNRAS, 359,567Baraffe, I., Chabrier, G., Fortney, J. & Sotin, C. 2014, in Proto-stars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ.Arizona Press), 763Batygin, K., Bodenheimer, P. & Laughlin, G. 2009, ApJ, 704, L49Batygin, K. & Stevenson, D. J. 2013, ApJ, 769, L9Burrows, A. S. 2014, Nature, 513, 345Carter, J. A. & Winn, J. N. 2010, ApJ, 709, 1219Carter, J. A., Agol, E., Chaplin, W. J., et al. 2012, Science, 337,556Chen, J. & Kipping, D. M. 2017, ApJ, 834, 17Dorn, C., Khan, A., Heng, K., Connolly, J. A. D., Alibert, Y.,Benz, W. & Tackley, P., 2015, A&A, 577, 83Dorn, C., Venturini, J., Khan, A., Heng, K., Alibert, Y., Helled,R., Rivol- dini, A. & Benz, W. 2017, A&A, 597, 37Dziewonski, A. M. & Anderson, D. L. 1981, Physics of the Earthand Planetary Interiors, 25, 297Fortney, J. J. & Nettelmann, N. 2010, Space Sci. Rev., 152, 423Fulton, B, J., Petigura, E. A., Howard, A. W., et al. 2017, AJ,154, 109Holman, M. J. & Murray, N. W. 2005, Science, 307, 1288Kaltenegger, L., Henning, W. G. & Sasselov, D. D. 2010, AJ, 140,1370Kaltenegger, L. 2017, Ann. Rev. of Astronomy & Astrophysics,55, 433Kipping, D. M., Spiegel, D. S. & Sasselov, D. D. 2013, MNRAS424, 1883Kipping, D. M. 2014, In search of exomoons. In Frank N.Bash Symposium 2013: New Horizons in Astronomy, Octo-ber 6ˆa ˘A¸S8, 2013, held at the University of Texas at Austin.arXiv:1405.1455.Kipping, D. M., Schmitt, A. R., Huang, X., Torres, G., Nesvorn´y,D., Buchhave, L. A., Hartman, J. & Bakos, G. ´A. 2015, ApJ,813, 14Kramm, U., Nettelmann, N., Fortney, J. J., Neuh¨auser, R. & Red-mer, R. 2012, A&A, 538, 146Lopez, E. D. & Fortney, J. J. 2013, ApJ, 776, 2Lopez, E. D. & Fortney, J. J. 2014, ApJ, 792, 1 McTier, M. & Kipping, D. M. 2017, MNRAS, submittedNoak, L., Godolt, M., von Paris, P., Plesa, A.-C., Stracke, B.,Breuer, D. & Rauer, H. 2014, Planetary and Space Science,98, 14Owen, J. E. & Morton, T. D. 2016, ApJ, 819, L10Owen, J. E. & We, Y. 2017, AAS Journals, submiited (arXiv e-print:1705.10810)Rogers, L. A. 2015, ApJ, 801, 41Seager, S. & Sasselov, D. D. 2000, ApJ, 537, 916Seager, S. & Hui, L. 2002, ApJ, 574, 1004Seager, S., Kuchner, M., Hier-Majumder, C. A. & Militzer, B.2007, ApJ, 669, 1279Spiegel, D. S., Fortney, J. J. & Sotin, C. 2014, PNAS, 111, 12622Unterborn, C. T., Dismukes, E. E. & Panero, W. R. 2016, ApJ,819, 32Valencia, D., O’Connell, R. J. & Sasselov, D. 2006, Icarus, 181,545Wolfgang, A., Rogers, L. A. & Ford, E. B. 2016. ApJ, 825, 19Zeng, L. & Sasselov, D. D. 2013. PASP, 125, 227Zeng, L., Sasselov, D. D. & Jacobstein, S. B. 2016, ApJ, 819, 127Zhu, W., Huang, C. X., Zhou, G. & Lin, D. N. C. 2014, ApJ, 796,67This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000