Tidally-Distorted, Iron-Enhanced Exoplanets Closely Orbiting Their Stars

TL;DR

This paper models the shape and composition of ultra-short-period exoplanets, revealing that tidally distorted planets like KOI 1843.03 can be highly elongated and iron-rich, with implications for planetary formation.

Contribution

First three-dimensional interior structure simulations of tidally distorted, iron-enhanced exoplanets, showing extreme elongation and specific composition requirements to avoid tidal disruption.

Findings

KOI 1843.03 may be elongated with an aspect ratio up to 1.79.

It must be as iron-rich as Mercury to avoid tidal disruption.

Nearly half of studied ultra-short-period planets are iron-enhanced.

Abstract

The transiting planet candidate KOI 1843.03 ($0.6 R_\oplus$ radius, 4.245 hour orbital period, $0.46 M_\odot$ host star) has the shortest orbital period of any planet yet discovered. Here we show, using the first three-dimensional interior structure simulations of ultra-short-period tidally distorted rocky exoplanets, that KOI 1843.03 may be shaped like an American football, elongated along the planet-star axis with an aspect ratio of up to 1.79. Furthermore, for KOI 1843.03 to have avoided tidal disruption (wherein the planet is pulled apart by the tidal gravity of its host star) on such a close-in orbit, KOI 1843.03 must be as iron-rich as Mercury (about 66% by mass iron compared to Mercury's 70% by mass iron, Hauck et al. 2013). Of the ultra-short-period ($P_\mathrm{orb} \lesssim 1$ day) planets with physically-meaningful constraints on their densities characterized to date, just under half (4 out of 9) are iron-enhanced. As more are discovered, we will better understand the diversity of rocky planet compositions and the variety of processes that lead to planetary iron enhancement.

Tables

Table 1. Table 1: Approximate values for the orbital velocities and fluxes at the orbital surface for the four planets we consider in this paper. We compare these values to the corresponding values for Mercury, the planet with the shortest orbital period in our own Solar System.

Planet	Orbital velocity (km s$-1$)	Flux (erg cm$-2$ s$-1$)
KOI 1843.03	320	$3\times {10}^9$
K2 137b	270	$1\times {10}^9$
K2 106b	270	$2\times {10}^{10}$
K2 229b	240	$3\times {10}^9$
Mercury	48	$9\times {10}^6$

Equations

$$a\simeq 2.44R_{\star }{\left(\frac{{\rho }_{\star }}{{\rho }_p}\right)}^{1/3}$$

$$P_{\mathrm{orb}}\simeq 12.6\text{leavevmode}\text{hr}{\left(\frac{{\rho }_p}{1\text{leavevmode}\text{g}\text{leavevmode}\text{cmtextsuperscript−3}}\right)}^{-1/2}.$$

$$H=\int {\rho }^{-1}\mathrm{d}P,$$

$${\hat{H}}_{i,j,k}=\hat{C}-{\hat{\Phi }}_{i,j,k}-{\hat{\Omega }}^2{\hat{\Psi }}_{i,j,k},$$

$$\left|\max \left({\hat{H}}_{n+1}-{\hat{H}}_n\right)/\max \left({\hat{H}}_{n+1}\right)\right|\le \delta ,$$

$$\left|\left({\hat{\Omega }}_{n+1}^2-{\hat{\Omega }}_n^2\right)/{\hat{\Omega }}_{n+1}^2\right|\le \delta ,$$

$$\left|\left({\hat{C}}_{n+1}-{\hat{C}}_n\right)/{\hat{C}}_{n+1}\right|\le \delta .$$

$${\Phi }_p\left(\mathbf{r}\right)=-G\underset{V}{\int }\frac{\rho \left({\mathbf{r}}^{\prime }\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}\mathrm{d}{V}^{\prime }.$$

$${\Phi }_p\left(r,\mu ,\varphi \right)=-4G\sum _{\ell =0}^{\infty }\sum _{\begin{array}{c}m=0\\ \ell +m\text{leavevmode}\mathrm{even}\end{array}}^{\ell }{ϵ}_m\frac{\left(\ell -m\right)!}{\left(\ell +m\right)!}{P}_{\ell }^m\left(\mu \right)\cos m\varphi \underset{0}{\overset{\infty }{\int }}\mathrm{d}{r}^{\prime }\text{leavevmode}{f}_{\ell }\left(r,r^{\prime }\right)\underset{0}{\overset{1}{\int }}\mathrm{d}{\mu }^{\prime }\text{leavevmode}{P}_{\ell }^m\left({\mu }^{\prime }\right)\underset{0}{\overset{\pi }{\int }}\mathrm{d}{\varphi }^{\prime }\text{leavevmode}\rho \left(r^{\prime },{\mu }^{\prime },{\varphi }^{\prime }\right)\cos m{\varphi }^{\prime }$$

$$f_{\ell }\left(r,r^{\prime }\right)=\{\begin{array}{ll}{r^{\prime }}^{\ell +2}/r^{\ell +1},& \mathrm{if}\text{leavevmode}{r}^{\prime }<r\\ r^{\ell }/{r^{\prime }}^{\ell -1},& \mathrm{if}\text{leavevmode}r<r^{\prime }\end{array},$$

$${ϵ}_m=\{\begin{array}{ll}1,& \mathrm{if}\text{leavevmode}m=0\\ 2,& \mathrm{if}\text{leavevmode}m\ne 0\end{array}.$$

$${\hat{\Phi }}_{p;i,j,k}=-4\sum _{\ell =0}^{{\ell }_{\text{max}}}\sum _{\begin{array}{c}m=0\\ \ell +m\text{leavevmode}\mathrm{even}\end{array}}^{\ell }{ϵ}_m\left(\frac{4\pi }{2\ell +1}\right){\tilde{P}}_{\ell }^m\left({\mu }_j\right)\cos m{\varphi }_k\times {\Lambda }_{i,\ell ,m}^{\hat{r}},$$

$${\Lambda }_{i,\ell ,m}^{\hat{r}}=h_{\hat{r}}\sum _{s=0}^{N-1}{S}_s^N{f}_{\ell }\left({\hat{r}}_i,{\hat{r}}_s\right){\Lambda }_{s,\ell ,m}^{\mu },$$

$${\Lambda }_{s,\ell ,m}^{\mu }=\sum _{t=0}^{P-1}{G}_t^P{\tilde{P}}_{\ell }^m\left({\mu }_t\right){\Lambda }_{s,t,m}^{\varphi },$$

$${\Lambda }_{s,t,m}^{\varphi }=\sum _{u=0}^{Q-1}{G}_u^Q{\hat{\rho }}_{s,t,u}\cos m{\varphi }_u.$$

$${\tilde{P}}_{\ell }^m\left(\mu \right)=\sqrt{\frac{2\ell +1}{4\pi }\frac{\left(\ell -m\right)!}{\left(\ell +m\right)!}}{P}_{\ell }^m\left(\mu \right)$$

$${\hat{\Psi }}_{i,j,k}=-\frac{1}{2}{\varpi }_{i,j,k}^2$$

$${\hat{\Psi }}_{i,j,k}=-\frac{1}{2}\left({\left({\hat{x}}_{i,j,k}-{\hat{x}}_{\text{cm}}\right)}^2+{\hat{y}}_{i,j,k}^2\right)$$

$${\hat{x}}_{\text{cm}}=\frac{{\hat{M}}_p{\hat{x}}_p+{\hat{M}}_{\star }{\hat{x}}_{\star }}{{\hat{M}}_p+{\hat{M}}_{\star }}$$

$${\hat{\Phi }}_{i,j,k}={\hat{\Phi }}_{p;i,j,k}+{\hat{\Phi }}_{\star ;i,j,k}$$

$${\hat{\Phi }}_{\star ;i,j,k}=\frac{-{\hat{M}}_{\star }}{\sqrt{{\left({\hat{x}}_{i,j,k}-{\hat{x}}_{\star }\right)}^2+{\hat{y}}_{i,j,k}^2+{\hat{z}}_{i,j,k}^2}}.$$

$${\hat{\Omega }}^2=-\left(\hat{\Phi }{|}_A-\hat{\Phi }{|}_B\right)/\left(\hat{\Psi }{|}_A-\hat{\Psi }{|}_B\right)$$

$$\hat{C}=\hat{\Phi }{|}_A+{\hat{\Omega }}^2\hat{\Psi }{|}_A.$$

$$\hat{H}\equiv H/G{r}_{{}_A}^2{\rho }_{\text{max}}.$$

$$\rho \left(P\right)=c{P}^n,$$

$$\begin{gathered} P_{\mathrm{en}}\left(\eta \right)=\frac{3}{2}{K}_{0,\mathrm{en}}\left({\eta }^{7/3}-{\eta }^{5/3}\right)\times \\ \left[1+\frac{3}{4}\left(K_{0,\mathrm{en}}^{\prime }-4\right)\left({\eta }^{2/3}-1\right)\right]. \end{gathered}$$

$${\tilde{H}}_{\mathrm{en}}\left(\rho \right)=\frac{3K_0{\rho }^{2/3}\left(9(K_0^{\prime }-4){\rho }^{4/3}+7\left(14-3K_0^{\prime }\right){\left(\rho {\rho }_0\right)}^{2/3}+5(3K_0^{\prime }-16){\rho }_0^{4/3}\right)}{16{\rho }_0^3}.$$

$$H_{\mathrm{en}}\left(\rho \right)={\tilde{H}}_{\mathrm{en}}\left(\rho \right)-{\tilde{H}}_{\mathrm{en}}\left({\rho }_{0,\mathrm{en}}\right).$$

$$H_{\mathrm{pv}}\left(P\right)=H_{\mathrm{en}}\left({\rho }_{\text{trans}}\right)+\left({\tilde{H}}_{\mathrm{pv}}\left(P\right)-{\tilde{H}}_{\mathrm{pv}}\left(P_{\text{trans}}\right)\right)$$

$$H_{\mathrm{fe}}\left(P\right)=H_{\mathrm{pv}}\left(P_{\text{cmb}}\right)+\left({\tilde{H}}_{\mathrm{fe}}\left(P\right)-{\tilde{H}}_{\mathrm{fe}}\left(P_{\text{cmb}}\right)\right).$$

Full text

methods: numerical — planets and satellites: composition — planets and satellites: interiors

††software: george (Ambikasaran et al., 2014), GNU Scientific Library (Galassi et al., 2009), isochrones (Morton, 2015), scipy (Jones et al., 2001–)

1 Introduction

The compositions of rocky planets reflect a combination of the compositions of their host star, the condensation sequence that concentrates elements heavier than hydrogen and helium into solids, and processing during planet formation and subsequent evolution. To leading order, the Earth is comprised of an iron-dominated core (32% by mass) and silicate mantle (68% by mass), with roughly the same relative elemental abundances as in the solar photosphere (Lineweaver & Robles, 2009). Most rocky exoplanets with measured masses and radii also follow this trend and are consistent with Earth’s composition with some scatter (Dressing et al., 2015). In contrast, Mercury, at 70% by mass iron (Hauck et al., 2013), is significantly enhanced in iron relative to solar abundances.

For planets on very short orbital periods ($\lesssim 1$ day), tides can be used to constrain the planets’ bulk densities and compositions. Planets in orbit around a star will experience a tidal force, as the day side of the planet feels a stronger attractive gravitational force than the night side. Planets, by the IAU definition111https://www.iau.org/static/resolutions/Resolution_GA26-5-6.pdf (Boss et al., 2007), are sufficiently massive for their self-gravity to overcome their rigid body forces and to achieve hydrostatic equilibrium shapes. As a result, a tidal field causes an orbiting planet to become elongated in the direction of the planet-star axis (e.g., de Pater & Lissauer, 2010). If the tidal forces are too strong (the planet is too close to its star), the planet may be tidally disrupted and pulled apart, thus becoming a ring around the host. The minimum distance at which a fluid planet can avoid tidal disruption is called the Roche limit. For an incompressible fluid, this limiting distance is given by

[TABLE]

(Roche, 1849), where $a$ is the orbital semi-major axis, $R_{\star}$ is the stellar radius, $\rho_{\star}$ is the stellar density, and $\rho_{p}$ is the planet density. Following Rappaport et al. (2013), we can rewrite Equation 1 using Kepler’s third law to express $a$ in terms of the orbital period $P_{\mathrm{orb}}$; then, the expression has no dependence on the stellar density, and we find

[TABLE]

The Roche limit is a familiar concept in the context of the rings and satellites of Saturn as well as interacting binary stars. The discovery of exoplanets that are very close (orbital period, $P_{\mathrm{orb}}<1$ day) to their host stars — found around $0.5\%$ of Sun-like stars (Sanchis-Ojeda et al., 2014) — open the opportunity to apply the Roche limit to Earth-mass-scale planets.

The transiting exoplanet candidate KOI 1843.03 has the shortest orbital period known to date. For KOI 1843.03 to have avoided tidal disruption on such a close-in orbit, previous estimates suggest that it must have a mean density of at least 7 g cm$-3$ (Rappaport et al., 2013). This density lower limit, however, relies upon interpolating the Roche limits of single-component polytrope models, wherein the pressure $P$ and density $\rho$ within the planet are related by a power-law $P\propto\rho^{\gamma}$ with constant $\gamma$. These do not accurately capture the density profiles of differentiated rocky bodies with sizes $>1000$ km. A more accurate calculation of the Roche limit is needed to constrain the composition of KOI 1843.03.

In this work, we develop the first self-consistent three-dimensional models for the interior structures of tidally-distorted rocky planets on ultra-short period (USP) orbits ($P_{\mathrm{orb}}<1$ day). We apply these models to refine calculations of the Roche limit for USP rocky planets and to explore the diversity of USP planet compositions. The paper is structured as follows: In Section 2 we describe the methods used. We outline the primary results in Section 3 and discuss in Section 4.

2 Methods

2.1 Modeling technique

Calculation of the Roche limit for generic equations of state (EOS) must rely on a numerical solution. Treating the planet as an extended body necessitates computing three forces: the gravitational force from the star, the gravitational force that the planet exerts on itself, and the centrifugal force in the planet’s rotating rest frame. The sum of all these forces influence the shape of the planet, which changes the mass distribution and, by extension, the forces on all points inside the planet. There is no simple, analytic way to capture the circular nature of this problem.

We use a relaxation method developed by Hachisu (1986a, b) to model the three-dimensional structure of USP rocky planets. Starting from an initial guess for the planet density distribution, the method iteratively adjusts the enthalpy and density distribution until a self-consistent solution is reached, within a tolerance of $1\times 10^{-5}$. We expand the Hachisu method to include the gravitational potential of a point source star at a fixed distance from the planet. We also modify the equation of state (which describes how the density of a material varies with pressure) to more accurately capture the behavior of silicate rocks and iron, which have nonzero density at zero pressure.

We model two-layer planets consisting of a silicate mantle (enstatite upper mantle and perovskite lower mantle) surrounding an iron core. We simulated more than $280000$ planet configurations over a grid of central pressures $P_{\textrm{max}}$, core-mantle boundary pressures $P_{\textrm{cmb}}$, and scaled star-planet orbital separations $a/R_{p}$. At each grid point, we begin by simulating a nearly spherical planet, self-consistently solving for the host star mass. We then simulate planets that are successively more distorted (elongated along the star-planet axis). Once material begins to fly off the planet, the Roche limit has been surpassed.

2.2 Coordinate system

To model a planet with an unknown shape, we define a three-dimensional coordinate system as shown in Figure 1, where any point may be specified by a radial coordinate $\hat{r}$, polar coordinate $\theta$ measured from the $z$-axis, and azimuthal coordinate $\phi$ in the $x$-$y$ plane.

Following Hachisu (1986a, b), we establish two points, $A$ and $B$, along the $x$- and $y$-axes, respectively, that lie on the surface of the planet. The simulation is conducted in scaled, dimensionless coordinates such that the distance from the origin to $A$ is $\hat{r}_{\!{}_{A}}\equiv 1$, and, similarly, the scaled distance from the origin to $B$ is $\hat{r}_{\!{}_{B}}$; we denote the corresponding dimensionful quantities as $r_{\!{}_{A}}$ and $r_{\!{}_{B}}$.

To approximate physical quantities that are continuous over all space, we define a grid of values of $\hat{r}$, $\mu\equiv\cos{\theta}$, and $\phi$, sampling each quantity at every grid point; for our simulations, we use $N,P,Q=\left(129,17,33\right)$ divisions in $\hat{r}$, $\mu$, and $\phi$, respectively. The symmetries inherent in this system allow us to limit the simulation space to $\mu\in\left[0,1\right]$ and $\phi\in\left[0,\pi\right]$; following Hachisu (1986a), we use $\hat{r}\in\left[0,16/15\right]$ to ensure that the planet does not exceed the simulation volume. For a dimensionless physical quantity $\hat{X}$, we use a notation such that $\hat{X}_{i,j,k}=\hat{X}\!\left(\hat{r}_{i},\mu_{j},\phi_{k}\right)$.

The star is treated as a point mass on the $x$-axis with mass $M_{\star}$ at a dimensionless distance $\hat{a}=a/r_{\!{}_{A}}$, where $a$ is the radius of the planet’s circular orbit, measured from the planet’s center of mass to the location of the star. We do not make any assumptions about the coordinate of the planet’s center of mass, so it does not necessarily coincide with the origin.

The planet rotates about this axis with Keplerian angular velocity given by $\Omega^{2}=G\left(M_{\star}+M_{p}\right)/a^{3}$. Assuming that the planet is tidally locked to the star, we may work in the rest frame of the planet, in which the star is stationary.

2.2.1 Relaxation method

A single iteration of the relaxation method begins with a proposal dimensionless density distribution $\hat{\rho}$. We convert the density distribution to a dimensionless enthalpy $\hat{H}$ (see Equation 24). Enthalpy in this context is defined as

[TABLE]

with $\rho$ the density and $P$ the pressure. Enthalpy must meet all of the boundary conditions — zero enthalpy at $A$ and $B$ and a dimensionless rotation rate $\hat{\Omega}$ consistent with Kepler’s third law — according to

[TABLE]

where $\hat{C}$ is a scalar constant. $\hat{\Phi}$ is the total dimensionless gravitational potential, including influence from both the star and the planet, and $\hat{\Omega}^{2}\hat{\Psi}$ is the dimensionless centrifugal potential. We then convert the new dimensionless enthalpy to a new dimensionless density distribution. Iterations of this procedure continue until a relative tolerance condition between consecutive iterations is reached. We use a relative tolerance value of $\delta=1\times 10^{-5}$ that must be satisfied for $\hat{H}$, $\hat{C}$, and $\hat{\Omega}^{2}$, such that, between iterations $n$ and $n+1$,

[TABLE]

[TABLE]

and

[TABLE]

These are the same metrics employed by Hachisu (1986a).

A two-layer planet may be uniquely specified by setting the values of $P_{\textrm{cmb}}$, $P_{\textrm{max}}/P_{\textrm{cmb}}$, $\hat{a}$, and $\hat{r}_{\!{}_{B}}$. For any set of values $\{P_{\textrm{cmb}},P_{\textrm{max}}/P_{\textrm{cmb}},\hat{a}\}$, we begin with the largest value of $\hat{r}_{\!{}_{B}}$ less than 1 and an ansatz dimensionless density distribution $\hat{\rho}$ that is a uniform density ellipsoid with radii $\left(\hat{r}_{\!{}_{A}},\hat{r}_{\!{}_{B}},\hat{r}_{\!{}_{B}}\right)$ in the $x$, $y$, and $z$ directions, respectively, ensuring that the ansatz satisfies the boundary conditions. After the system converges, we reduce the value of $\hat{r}_{\!{}_{B}}$ to the next grid point, effectively increasing the distortion each time, and use the previous converged result as the input ansatz to the next relaxation procedure.

2.2.2 Potential solver

By far, the most computationally difficult and expensive component of this method is finding the gravitational potential due to the extended planet itself at every point in space,

[TABLE]

Since the system is symmetric in $y$ and $z$, we may expand and simplify Equation 8 as

{widetext}

[TABLE]

where $P_{\ell}^{m}$ are the associated Legendre polynomials,

[TABLE]

and

[TABLE]

We employ Simpson’s rule, following Hachisu (1986b), in the $\hat{r}$ dimension but use Gauss-Legendre quadrature in the $\mu$ and $\phi$ dimensions. We develop the discrete, dimensionless equivalent of Equation 9, {widetext}

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Here,

[TABLE]

is the normalized associated Legendre polynomial222As computed by the function gsl_sf_legendre_sphPlm from the GNU Scientific Library (GSL) (Galassi et al., 2009). $h_{\hat{r}}$ is the interval in the $\hat{r}$ coordinate between successive grid points. The coefficients $S_{i}^{n}$ are the alternative composite Simpson’s rule coefficients (Press et al., 1988), and the coefficients $G_{i}^{n}$ are the fixed-order Gauss-Legendre quadrature weights333As computed by the function gsl_integration_glfixed_point from the GSL (Galassi et al., 2009), which depend on the integration interval. Since Gauss-Legendre quadrature is an “open” integration scheme, the endpoints of the integration interval in $\mu$ and $\phi$ are excluded; as a matter of computational convenience, we simply inject the endpoints with zero integration weight since they are needed to define the locations $A$ and $B$. We employ the same integration scheme to compute the dimensionless mass of the planet $\hat{M}_{p}$ and the coordinate of its center of mass $\hat{x}_{p}$.

Once the planet’s center of mass is computed, the coordinate of the star is easily determined by $\hat{x}_{\star}=\hat{x}_{p}+\hat{a}$. $\hat{\Psi}$ is given by

[TABLE]

where $\varpi_{i,j,k}$ is the distance from the point to the rotation axis; in our case,

[TABLE]

where

[TABLE]

is the coordinate of the center of mass of the entire system. The total gravitational potential $\hat{\Phi}$ is given by

[TABLE]

where

[TABLE]

However, at this point, the mass of the star is an unknown. We solve the system of equations that gives the mass consistent with the boundary conditions on the enthalpy and the dimensionless form of Kepler’s third law, $\hat{\Omega}^{2}=\left(\hat{M}_{\star}+\hat{M}_{p}\right)/\hat{a}^{3}$. Then,

[TABLE]

and

[TABLE]

Equation 4 then gives the enthalpy at every point.

We follow (Hachisu, 1986a) by defining a dimensionless enthalpy $\hat{H}$ and its maximum ${\hat{H}}_{\textrm{max}}$, where

[TABLE]

$\rho_{\textrm{max}}$ is fixed in the simulation because it can be determined directly from $P_{\textrm{max}}$ and the equation of state. To generate the non-analytic function that maps dimensionless enthalpy to dimensionless density, we finely sample the pressure $P$ from [math] to $P_{\textrm{max}}$ and calculate the density at each pressure with our equation of state. Equation 3 gives the enthalpy $H$ at each pressure, and $H_{\textrm{max}}$ is just the maximum of these values. Letting $\hat{\rho}\equiv\rho/\rho_{\textrm{max}}$, we linearly interpolate $\hat{\rho}$ as a function of $H/H_{\textrm{max}}$ and obtain a function which maps scaled, dimensionless enthalpy $\hat{H}/{\hat{H}}_{\textrm{max}}$ to dimensionless density $\hat{\rho}$.

2.2.3 Equation of state

Previous work modelling distorted stars by Hachisu (1986a) assumes a polytropic equation of state,

[TABLE]

where $\rho$ is density and $P$ is pressure. This equation of state is not appropriate for a rocky planet because it does not allow for nonzero density at zero pressure. In this work, we consider only planets with two layers: an iron core and a silicate mantle; our method may be extended to planets with different compositions and any number of layers, however.

At low pressures ($P\leq P_{\textrm{trans}}=23\times 10^{9}$ Pa), we apply the Birch-Murnagham equation of state (BME) for enstatite (Seager et al., 2007). For $\eta\equiv\rho/\rho_{0,\mathrm{en}}$, we have (Seager et al., 2007),

[TABLE]

Above $P_{\textrm{trans}}$, we use a tabulated equation of state for perovskite and iron or FeS (depending on the core composition assumed) (Seager et al., 2007). The transition between perovskite and the iron-dominated core is defined to occur at a core-mantle boundary pressure $P_{\textrm{cmb}}$. We have adopted room temperature (300 K) EOSs. Including thermal expansion, which we do not do in this work, will make the lower limits derived on the iron mass fraction of KOI 1843.03 even more severe. Figure 2 shows the full equation of state that we have adopted.

The relaxation method (described above) requires the conversion of enthalpy to density, which depends solely on the equation of state. Formally, enthalpy $H$ is given by Equation 3. The indefinite integral that corresponds to substituting Equation 26 in Equation 3 is {widetext}

[TABLE]

The enthalpy as a function of density $\rho$ is then given by

[TABLE]

The enthalpy as a function of density for the tabulated equations of state is computed by interpolating a cumulative trapezoidal integration approximating Equation 3, which we denote $\widetilde{H}_{\mathrm{pv}}$ and $\widetilde{H}_{\mathrm{fe}}$ for perovskite and iron, respectively. This gives us, for $P_{\mathrm{en}}\!\left(\rho_{\textrm{trans}}\right)=P_{\textrm{trans}}$,

[TABLE]

and

[TABLE]

2.3 Model validation

To validate our method, we reproduce the classical Roche limit for an incompressible fluid body. For this test case, we achieve $\left|\Delta P_{\mathrm{orb}}\right|/P_{\mathrm{orb}}\approx 0.5\%$, where $P_{\mathrm{orb}}$ is the analytic Roche limiting orbital period (Equation 2), and $\left|\Delta P_{\mathrm{orb}}\right|$ is the absolute difference between the analytic expectation and the Roche limit we numerically derive following the method described above. In our validation experiements, the relative error in $P_{\mathrm{orb}}$ is observed to be independent of the density of the fluid, as is expected because the simulation is run with a dimensionless, normalized density.

2.4 Model interpolation procedure

After all models have been computed, we distill meaningful results by smoothly interpolating within the model grid. When interpolating our model grids, we use the Gaussian process code george (Ambikasaran et al., 2014). Our chosen kernel is a constant kernel multiplied by a squared exponential kernel. We allow for “white noise,” which in this case is not observational but rather computational noise. We also use a convex hull algorithm as a safeguard against extrapolation. This reduces the extent to which our interpolation code can extrapolate outside our models’ bounds.

3 Results

3.1 KOI 1843.03

Our self-consistent 3D models show that KOI 1843.03 must be very iron-rich to avoid tidal disruption on an orbital period of 4.245 hours. Assuming the planet is composed a pure iron core surrounded by a magnesium-silicate mantle, we find that the $R_{p}=0.61^{+0.12}_{-0.08}\leavevmode\nobreak\ R_{\earth}$ radius constraints (Rappaport et al., 2013) imply that KOI 1843.03 must be at least $60_{-8}^{+7}\%$ iron by mass (Figure 3). Since rocky planets become compressed to higher densities with increasing size, larger values of the planet radius within the $1\sigma$ range translate into more relaxed lower bounds on the iron mass fraction of the planet. Based on our planet interior models, we expect KOI 1843.03’s mass to fall in the range $0.32$ – $1.06\leavevmode\nobreak\ M_{\earth}$ (Figure 4).

Figure 4 displays interpolated planet masses for a range of orbital periods and core mass fractions in systems consistent with KOI 1843.03’s host star mass and transit radius. The boundary of the colored contours in the lower left-hand corner corresponds to the Roche limit; as orbital period decreases, the core mass fraction is more tightly constrained. As anticipated, considering a fixed orbital period in Figure 4, increasing the planet’s iron mass fraction increases the planet’s mass. Less intuitively, at fixed core mass fraction, decreasing the planet’s orbital period also leads to an increase in the inferred planet mass. This is due to the tidal distortion of the planet; at shorter orbital periods, the volume of the planet exceeds $4/3\pi R_{\rm transit}^{3}$ by larger and larger factors.

As it orbits less than one stellar radius from the host star’s photosphere ($a/R_{\star}=1.9$), KOI 1843.03 will be significantly elongated in the direction of the star due to tidal distortion. Based on the parameters of KOI 1843.03, our models predict that it must be at least $31\%$ longer along the star-planet line than along the perpendicular axes (aspect ratio of about $1.3$), and our models support a value up to nearly twice as long along the star-planet line (aspect ratio of almost $1.8$); various possibilities are illustrated in Figure 5. For comparison, Saturn has an aspect ratio of about $1.1$.

The Earth’s core is not pure iron; it contains an unknown mixture of light elements. To determine what effect these light elements might have, we generated a second grid of models with an EOS appropriate to an FeS core (see Figure 6). KOI 1843.03 would need a core mass fraction of at least $80\%$ if it has a core comprised of FeS, but this measurement is only valid for a transit radius of $2\sigma$ greater than the mean. Both pure Fe and pure FeS are end-member core compositions; the true core density of KOI 1843.03 likely lies somewhere in between.

3.2 Iron-enhanced USP planets

KOI 1843.03 is one of a growing class of iron-enhanced, closely-orbiting planets discovered.

K2 137b is remarkably similar to KOI 1843.03 but has a slightly longer orbital period (by 4 minutes), a larger transit radius ($0.89\leavevmode\nobreak\ R_{\earth}$), and a more massive host star ($0.463\leavevmode\nobreak\ M_{\sun}$) (Smith et al., 2018). Our 3D models show that K2 137b must be at least $42\pm 5\%$ iron by mass to have avoided tidal disruption. Its mass lies somewhere between $1.01$ and $2.80\leavevmode\nobreak\ M_{\earth}$ (Figure 7), which is consistent with the radial velocity upper limit of $3\leavevmode\nobreak\ M_{\text{Jup}}$. We find that K2 137b’s aspect ratio is bounded between $1.21$ and $1.66$ (provided the true transit radius is within the 1$\sigma$ measurement, see Figure 8). The constraint on K2 137b’s aspect ratio is not as extreme as the constraint on KOI 1843.03’s.

Two additional transiting exoplanets — K2 229b ($M_{p}=2.59\leavevmode\nobreak\ M_{\earth}$, $R_{\mathrm{transit}}=1.164\leavevmode\nobreak\ R_{\earth}$, $P_{\mathrm{orb}}=14.0$ hr, Santerne et al. 2018) and K2 106b ($M_{p}=8.36\leavevmode\nobreak\ M_{\earth}$, $R_{\mathrm{transit}}=1.52\leavevmode\nobreak\ R_{\earth}$, $P_{\mathrm{orb}}=13.7$ hr, Guenther et al. 2017) — have been inferred to have iron-rich compositions based on their radial velocity measured masses. Based on our models, we infer iron mass fractions of $0.605^{+0.204}_{-0.181}$ and $0.691^{+0.206}_{-0.169}$, and aspect ratios of $1.02^{+0.002}_{-0.003}$ and $1.01^{+0.003}_{-0.003}$, for K2 229b and K2 106b respectively. (Note that the errorbars reported here are not one standard deviation, as the distributions tend to be non-Gaussian, but rather the 25% and 75% percentiles; the full distributions are shown in Figures 9 and 10.)

Of the roughly 7 rocky ultra-short period exoplanets ($R_{p}\leq 1.7\leavevmode\nobreak\ R_{\earth}$, $P_{\mathrm{orb}}\leq 1$ day) with masses and radii measured to within 20% precision to date (CoRoT-7b, Kepler-10b, Kepler-78b, K2 106b, K2 141b, K2 229b, HD-3167b), 2 are iron-enhanced. Including the planets with density upper limits from the Roche limit (KOI 1843.03 and K2 137b), we find that just under half (4 out of 9) of the ultra-short period exoplanets with physically-meaningful constraints on their densities characterized to date are iron-enhanced.

3.3 Analytic approximation to the Roche limit

Our numerical Roche limits can be approximated by modifying the power-law parameters of the well-known expression for the classical Roche limit, re-expressed in terms of orbital period using Kepler’s third law (Rappaport et al., 2013); this is given by Equation 2. Including additional terms up to quadratic order that encapsulate a dependence on the degree of central concentration of the planet (i.e., the ratio of the planet’s maximum/central density to its mean density), we obtain,

[TABLE]

where x\equiv\rho_{p}/\left(1\leavevmode\nobreak\ \text{\text{g}}\leavevmode\nobreak\ \text{\text{cm}\textsuperscript{-3}}\right), $y\equiv\rho_{\textrm{max}}/\rho_{p}$. Using sequential least squares programming444SLSQP, implemented in SciPy, (Jones et al., 2001–), using an algorithm originally developed by Dieter Kraft, we find the following best fit parameters to the 2D surface defining the Roche limit: $C=12.013$, $\alpha=-0.571$, $\beta=0.047$, $\gamma=0.108$, $\delta=0.562$, and $\epsilon=-0.527$.

Rappaport et al. (2013) report an interpolation formula without the quadratic terms, where $C=12.6$, $\alpha=-0.5$, and $\beta=-0.16$. If we set the quadratic terms to zero, our best-fitting revised Roche limit has $C=11.86$, $\alpha=-0.52$, and $\beta=0.086$.

While this approximation is a useful tool, interpolation within the model grid is our suggested approach for using the models.

4 Discussion

4.1 Formation scenarios

How did KOI 1843.03, K2 137b, K2 229b and K2 106b form and/or evolve to such short orbital periods and iron-rich compositions?

Several theories for the origin of Mercury’s high iron content involve impacts that remove the outer silicate layers of a differentiated planet. Scenarios involving a single giant impact (e.g., Benz et al., 2008), a hit-and-run collision (e.g., Asphaug & Reufer, 2014), or the cumulative effect of multiple collisions can all feasibly lead to a Mercury-like outcome, though a single-giant impact or hit-and-run impact require highly tuned collision geometries to reproduce Mercury’s mass and iron mass fraction (Chau et al., 2018). In the context of exoplanets, Marcus et al. (2010) used smoothed particle hydrodynamics (SPH) simulations of giant impacts to infer that iron mass fractions of up to than 80% can be achieved with a single impact for planets less than $2\leavevmode\nobreak\ M_{\earth}$. Marcus et al. (2010) neglected to track the dynamical evolution of the impact ejecta, however, and reaccretion of the ejected mantle is likely to dilute the iron-enhancement of a giant impact (Gladman & Coffey, 2009). USP planets would be susceptible to high-velocity erosive collisions due to the extreme orbital velocities along their orbits (for example, about 320 km s$-1$ for KOI 1843.03 compared to about 48 km s$-1$ for Mercury, Table 1). Their proximity to their stars, with shorter orbital timescales and stronger stellar irradiation environments (Table 1), would also affect the reaccretion of ejected silicates. Further work is needed to evaluate the effect of collisions on the compositions of USP planets.

Alternatively, these closely-orbiting iron-enhanced planets could have initially formed from iron-rich material. Both the condensation sequence (wherein iron condenses at a higher temperature than magnesium silicates) (Lewis, 1972) and photophoresis (which separates high-thermal-conductivity iron dust grains from lower-thermal-conductivity silicate grains) (Wurm et al., 2013) can lead to an enhancement of iron in the solid phase at the inner edge of the protoplanetary disk. These fractionation processes that operate primarily at the disk inner edge could imprint themselves as a statistical iron enhancement of the ultra-short-period planet population.

A third possibility is that KOI 1843.03 and K2 137b are right at their Roche limits and have been gradually losing their outer silicate layers to Roche lobe overflow as their orbits tidally decay (Jia & Spruit, 2017). If KOI 1843.03 started with a chondritic or Earth-like iron-to-silicate ratio and mass of about $0.7\leavevmode\nobreak\ M_{\earth}$ (an intermediate value between our estimated limits), the planet’s initial mass would have been about $1.4\leavevmode\nobreak\ M_{\earth}$. The orbital period precision achieved over the 4-year baseline of the Kepler mission is insufficient to resolve expected decay in KOI 1843.03’s orbit. This scenario does not explain the compositions of the longer-orbital-period K2 229b and K2 106b, however, since they are outside their Roche limits.

4.2 Thermal effects

We have not modeled the interior temperature profiles of these planets. Indeed, we have adopted room temperature (300 K) equations of state. The surface temperature of KOI 1843.03 and other USP planets can exceed 2000 K, with temperature increasing further toward the center. Thermal expansion would cause the planet of specified mass and composition to have a larger volume and lower mean density compared to the models presented here. Temperature may also affect the pressure of the phase transition between enstatite and perovskite, which we have fixed to 23 GPa, following Sotin et al. (2007). There are two common ways of incorporating temperature into the EOS (Jackson, 1998): One may either regard the typical EOS coefficients as being temperature-dependent, or one may add a “thermal” pressure at every point. The effect of temperature is more severe for lower mass bodies. Including thermal expansion will make the constraints on the iron mass fraction of KOI 1843.03 even more severe, strengthening our conclusions.

4.3 Effect of material strength

Our models provide the first self-consistent constraints on the hydrostatic equilibrium shapes and Roche limits of ultra-short period rocky planets. The effect of material strength in the planets’ shapes is not taken into account in these calculations.

Looking to the Solar System bodies for inspiration, we see that once bodies are roughly 200 km (for icy materials) to 300 km (for rocky materials) in radius their self gravity is sufficient to overcome their material strength and they achieve a rounded shape. Iapetus (mean radius $734.5\pm 2.8$ km, Roatsch et al. 2009) is the largest Solar System body measured to have significant deviations from a hydrostatic equilibrium shape (Thomas, 2010). Since KOI 1843.03 is about 4000 km in radius, it is safely in the regime where self-gravity dominates the material forces and hydrostatic equilibrium determines its leading-order shape, satisfying the minimum mass criterion in the IAU definition of a planet.

To leading order, Earth-mass scale planets (such as KOI 1843.03, K2 137b, K2 229b, and K2 106b) are in hydrostatic equilibrium, with rigidity representing a minor correction. In our models, the central pressure of KOI 1843.03 are on the order of $10^{11}$ – $10^{12}$ Pa which is orders of magnitude larger than the shear strength of iron (Clatterbuck et al., 2003) and peridotite (Handy et al., 1999).

The high instellations of USP planets can lead to molten surfaces (e.g., Léger et al., 2011; Kite et al., 2016), which further limit deviations from hydrostatic equilibrium shapes. In the extreme of no heat redistribution, the substellar point of KOI 1843.03 could exceed 2000 K, computed from $T_{ss}=T_{\star}\left(1-\alpha\right)^{1/4}\sqrt{R_{\star}/a}\approx 2500$ K  assuming a basalt-like planet surface albedo $\alpha=0.1$ (e.g., Kite et al., 2016). This temperature is sufficiently hot to melt metallic iron (melting point 1811 K) and is hotter than the liquidus of peridotite (Takahashi, 1986), the dominant rock in Earth’s upper mantle. Thus, the planet’s surface would have too little strength to sustain topography that could significantly influence the transit depth.

In using the Roche limit to constrain the bulk compositions of KOI 1843.03 and K2 137b, we have followed Rappaport et al. (2013) and Jia & Spruit (2017) and neglected the effect of material strength. It is unclear whether material strength or friction would help the planet to survive intact inside its Roche limit for gigayear timescales (e.g., Davidsson, 1999; Holsapple & Michel, 2006). As highlighted by Winn et al. (2018), further work is needed to model the destruction of USP planets that exceed their Roche limits.

4.4 Planet mass loss

The surface of KOI 1843.03 could be actively sublimating. Kite et al. (2016) models the exchange between atmospheric silicate, surface magma pools, and interior material for a hot, rocky exoplanet. If a rock vapor atmosphere is contributing to the transit depth in the Kepler bandpass, that only makes our constraints on the the iron fraction in KOI 1843.03 even more severe.

Though they do not themselves show evidence of evaporation in Kepler photometry, KOI 1843.03, K2 137b, K2 229b and K2 106b could be more massive cousins to the catastrophically evaporating rocky planet discovered orbiting KIC12557548. KIC12557548 shows asymmetric and variable transit shapes that have been interpreted as evidence of a dusty outflow of vaporized material driven by a thermal wind (Rappaport et al., 2012). Even for KOI 1843.03, the smallest among these close-orbiting iron-enhanced planets, with a mass in excess of $0.3\leavevmode\nobreak\ M_{\earth}$, the escape velocity from the surface is too high to drive a substantial hydrodynamic wind of sublimated silicates (Perez-Becker & Chiang, 2013). Detailed models of radiative hydrodynamic winds from evaporating rocky USP planets (Perez-Becker & Chiang, 2013) show that a $0.1\leavevmode\nobreak\ M_{\earth}$ rocky planet could survive at a surface temperature of $\sim 2200$ K with negligible mass loss for tens of gigayears.

4.5 Potential for follow-up observations

Due to the red colors of the early M host-star, the near infrared (NIR) and infrared (IR) offer the best opportunities for further observational characterization of KOI 1843.03. Given our constraints on its mass and bulk composition, the possible range for KOI 1843’s radial velocity semi-amplitude ($K_{1}=0.60$ – $1.98$ m s*-1*) spans the current state of the art precision of $1$ m s$-1$. The host star is too faint for precision radial velocity follow-up in the visible with any existing telescope or instrument but may be a feasible candidate for radial velocity follow-up in the NIR. Photometric follow-up in the infrared, for instance with the Spitzer Space Telescope or the James Webb Space Telescope could confirm the planetary nature of KOI 1843.03 (Désert et al., 2015).

Photometric follow-up could also provide a longer time baseline to reveal evidence of tidal evolution of the orbit. In principle, given sufficient time sampling, the detailed shape of the transit lightcurve (as the projected cross-section of the planet changes viewing angle during transit) may further constrain KOI 1843.03’s aspect ratio and bulk composition. Such an effect has been studied for distorted giant planets (Leconte et al., 2011).

To predict the IR transit signal-to-noise ratio (SNR), we use the isochrones software (Morton, 2015) modified for Spitzer bandpasses using data from Hora et al. (2008) and Indebetouw et al. (2005) to compute predicted apparent magnitudes for KOI 1843 in each of the four Spitzer bandpasses. In the 3.6 µm, 4.5 µm, 5.8 µm, and 8.0 µmbands, respectively, we predict magnitudes of $11.0\pm 0.05$, $11.0\pm 0.08$, $10.9\pm 0.08$, and $10.9\pm 0.07$. Scaling from the SNRs obtained by Désert et al. (2015) on stars with similar Spitzer magnitudes as KOI 1843, a signal-to-noise ratio for a single KOI 1843.03 transit could range from about $0.36$ to $0.93$. Observing multiple transits could improve the SNR.

KOI 1843.03 may not hold the records for the shortest orbital period and most distorted known exoplanet for long. The Transiting Exoplanet Survey Satellite (TESS, launched in 2018) should find several ultra-short-period transiting planets as it surveys the brightest stars over the entire sky; recent simulations by Barclay et al. (2018) predict detection for $52$ planets with orbital periods of $P_{\mathrm{orb}}<1$ day. Our models, which provide the first self-consistent constraints on the Roche limits of Earth-mass-scale rocky planets (Figure 3), will enable composition constraints on these future ultra-short period planet discoveries.

E.M.P. and L.A.R. thank Dr. Benjamin T. Montet for his analysis of the KOI 1843.03 Kepler light curve and and constraint on the rate of change in the planet’s orbital period. We also thank Drs. Andrew Vanderburg, Saul Rappaport, Josh Winn, and Darin Ragozzine for helpful discussion, and reviewer Dr. Erik Asphaug for valuable suggestions. This material is based upon work supported by a National Science Foundation Graduate Research Fellowship under Grant Nos. DGE1144152 and DGE1745303. L.A.R. acknowledges NSF grant AST-1615315. The computations in this work were carried out with resources provided by the University of Chicago Research Computing Center.