Mass, Density, and Formation Constraints in the Compact, Sub-Earth Kepler-444 System including Two Mars-Mass Planets

TL;DR

This study uses transit-timing analysis of Kepler data to determine the masses of five sub-Earth planets in the ancient Kepler-444 system, revealing insights into its formation and orbital evolution.

Contribution

It provides the first mass measurements of two Mars-mass planets in a compact system, challenging smooth disk migration as the sole formation mechanism.

Findings

Two planets have confidently detected masses near Mars.

The mass ratio and tidal effects suggest complex formation history.

The system's architecture likely involved significant post-formation perturbations.

Abstract

Kepler-444 is a five planet system around a host-star approximately 11 billion years old. The five transiting planets all have sub-Earth radii and are in a compact configuration with orbital periods between 3 and 10 days. Here we present a transit-timing analysis of the system using the full Kepler data set in order to determine the masses of the planets. Two planets, Kepler-444 d ($M_\mathrm{d}=0.036^{+0.065}_{-0.020}M_\oplus$) and Kepler-444 e ($M_\mathrm{e}=0.034^{+0.059}_{-0.019}M_\oplus $), have confidently detected masses due to their proximity to resonance which creates transit timing variations. The mass ratio of these planets combined with the magnitude of possible star-planet tidal effects suggests that smooth disk migration over a significant distance is unlikely to have brought the system to its currently observed orbital architecture without significant post-formation perturbations.

Tables

Table 1. Table 1 Kepler-444 Planet Periods and Period Ratios

	Planet b			Planet c			Planet d			Planet e			Planet f
Period (d)	3.600105			4.545876			6.189437			7.743467			9.740501
Period Ratio			1.262707			1.361550			1.251078			1.257899
TTV Period (d)aafootnotemark:			89.5			73.1			1780.3			308.5

Table 2. Table 2 Photodynamic DEMCMC Posterior Median Values and 68.3% (1- σ 𝜎 \sigma equivalent) uncertainties.

Planet Parametersaafootnotemark:
	Planet b	Planet c	Planet d	Planet e	Planet f
$P$ (days)	${3.600105}_{-0.000037}^{+0.000031}$	${4.545876}_{-0.000031}^{+0.000030}$	${6.189437}_{-0.000037}^{+0.000053}$	${7.743467}_{-0.00010}^{+0.000060}$	${9.740501}_{-0.000026}^{+0.000078}$
$T_0$ (days)	${815.08383}_{-0.00055}^{+0.00052}$	${819.13903}_{-0.00044}^{+0.00042}$	${816.70059}_{-0.00072}^{+0.00072}$	${819.21772}_{-0.00083}^{+0.00087}$	${817.89759}_{-0.00032}^{+0.00038}$
$\sqrt{e}\cos \omega$	$-{0.03}_{-0.10}^{+0.14}$	${0.01}_{-0.13}^{+0.12}$	${0.098}_{-0.12}^{+0.065}$	$-{0.035}_{-0.090}^{+0.12}$	$-{0.059}_{-0.078}^{+0.12}$
$\sqrt{e}\sin \omega$	${0.048}_{-0.15}^{+0.099}$	$-{0.02}_{-0.11}^{+0.13}$	$-{0.014}_{-0.091}^{+0.10}$	${0.038}_{-0.11}^{+0.074}$	${0.052}_{-0.12}^{+0.075}$
$i$ (∘)	${92.00}_{-0.30}^{+0.26}$	${92.79}_{-0.11}^{+0.12}$	${91.95}_{-0.10}^{+0.11}$	${90.62}_{-0.35}^{+0.27}$	${92.087}_{-0.054}^{+0.058}$
$\Omega$ (∘)	0 (fixed)	0 (fixed)	0 (fixed)	0 (fixed)	0 (fixed)
$M_p/M_{\star }$ ($\times {10}^{-7}$)	${2.3}_{-1.6}^{+1.6}$	${4.5}_{-3.2}^{+3.5}$	${1.45}_{-0.81}^{+2.6}$	${1.34}_{-0.74}^{+2.35}$	${4.5}_{-3.5}^{+12}$
$R_p/R_{\star }$ ($\times {10}^{-3}$)	${4.967}_{-0.067}^{+0.070}$	${6.380}_{-0.087}^{+0.090}$	${6.613}_{-0.077}^{+0.079}$	${6.799}_{-0.076}^{+0.078}$	${9.39}_{-0.12}^{+0.13}$

Equations

$$K=\bigg{(}\frac{2\pi G}{P}(M_{\star}+M_{\mathrm{p}})\bigg{)}^{1/3}\frac{M_{\mathrm{p}}}{(M_{\star}+M_{\mathrm{p}})}\frac{\sin(i)}{\sqrt{1-e^{2}}},$$

$$\delta =\frac{n_j}{n_{j+1}}-\frac{(k+1)}{k},$$

$$t_{c,j}=\frac{4}{63}\frac{M_j{a}_j^{13/2}}{(G{M}_{\star }^3{)}^{1/2}{R}_j^5}\frac{3Q}{2k_2},$$

Full text

1. Introduction

Probing the mass-radius relationship for planets smaller than Earth is interesting to theorists as it may be used to constrain the formation and composition of these bodies, a topic of debate in the current literature (e.g Armitage, 2010; Chambers, 2010; Sinukoff et al., 2013; Chatterjee & Tan, 2014; Dupuy et al., 2016). A few planets in this size regime have been characterized (e.g. Sinukoff et al., 2013; Rappaport et al., 2013; Jontof-Hutter et al., 2015; Gillon et al., 2017); however, due to the small number of characterizable systems, little is yet known about the masses or compositions of the smallest ($\lesssim 1R_{\oplus}$, $\lesssim 1M_{\oplus}$) planets, despite them being among the most common in the galaxy (Malhotra, 2015).

Recent work has demonstrated the effectiveness of using photodynamic modeling to extract transit timing variations (TTVs) and planetary properties from systems with a low signal-to-noise ratio (SNR) (e.g., Carter et al., 2012; Barros et al., 2015; Mills et al., 2016). This technique takes advantage of the many transits of short-period planets observed in the *Kepler *data by fitting the entire light curve and all transits simultaneously. Here we apply this technique to Kepler-444.

Kepler-444’s planets (b, c, d, e, and f from inside to out) range in radii from 0.4 to 0.8 $R_{\oplus}$ and in orbital period from 3.6 to 9.8 days (Rowe et al., 2015; Campante et al., 2015). Their period ratios are near, but not exactly on, mean motion resonances (MMRs; see Table 1). Despite the compact architecture of the system, it is around a star $11.2\pm 1.0$ Gyr old (Campante et al., 2015) and therefore has likely been in a stable configuration for billions of years. A tight binary pair of M-dwarf stars also orbit together around Kepler-444 with a period of approximately 460 years and a distance of $\sim$ 60 AU (Campante et al., 2015). Such a configuration poses a puzzle regarding the early history of the Kepler-444 system, as planetary formation and migration in a truncated protoplanetary disk in the presence of a very nearby binary star pair is not well understood, with several effects newly proposed (e.g., Touma & Sridhar, 2015; Xu & Lai, 2016). Recent studies have attempted to understand the possible histories of the system and use it to place constraints on formation mechanisms (Dupuy et al., 2016; Papaloizou, 2016). However, such studies were unable to use the actual compositions or masses of the Kepler-444 planets since they were hitherto unknown. In this paper, we use photodynamics to put constraints on the masses of the planets in the Kepler-444 system and report mass detections for two of the planets: $M_{d}=0.036^{+0.065}_{-0.020}M_{\oplus}$ and $M_{e}=0.034^{+0.059}_{-0.019}M_{\oplus}$.

2. Methods

We initially identified potential transit timing variations in the Kepler-444 system by simultaneously fitting the raw *Kepler *light curve with a planet transit model (Mandel & Agol, 2002) and a 1-day wide polynomial to take into account systematic effects and stellar activity. We also compute the expected period of the TTV signal between each pair of planets analytically (Lithwick et al., 2012) (Table 1), noting that period of the expected signal for planets d and e matches the TTV observations well (Fig. 1). We find statistically significant TTV between planets d and e, but the signal for all other planets is undetectably low as theoretically expected. The same conclusion was reached independently by Hadden & Lithwick (2016), a survey of many *Kepler *systems showing TTVs.

In order to perform a more robust, simultaneous fit for all planetary parameters, we first reprocessed the raw *Kepler *lightcurve data. We use short-cadence (58.8 second integration) data when it was available (*Kepler *observing quarters 4, 6, and 15-17) and long cadence data (29.4 minute integrations) otherwise. We first discarded points whose quality flag had a value equal to or greater than 16. We then detrended the light curves by masking out the expected transit times plus 20% of the transit duration to account for possible TTVs and then fit a cubic polynomial model with a 1000-minute width centered on photometric data points spaced by 30 minute intervals. We interpolated between these points to determine a baseline and divide the measured flux at each data point by these values. This detrending method produces two regions of extreme curvature in the lightcurve due to edge effects, so we discard the small regions with times BJD-2454900 = 1405.10 to 1405.18 days and 1490.88 to 1490.97 days. To account for certain Quarters showing higher noise levels than others despite all quarters having similar quoted uncertainties, we assign an uncertainty of 5.3030402e-05 to points in Quarter 12, 2.3470900e-04 in Quarter 16, and 6.5361999e-04 in Quarter 17, an increase over other regions by a factors of roughly 5, 4, and 11 respectively (the ratio of their out of transit standard deviation). Lastly, we increase the uncertainties in all points by a multiplicative factor of 1.38073 so that a fiducial fit to the light curve has a $\chi^{2}=1.000$. This approach substantially increases the uncertainties on the fitted parameters compared to using the values reported by *Kepler *, allowing unmodeled noise to propagate to the final uncertainties on our parameter posteriors. We used data from *Kepler *Data Release 21 (DR21) for this analysis because there is less scatter in the DR21 data compared to the DR25.

Our photodynamic model integrates Newtonian equations of motions for the star and five planets including the light travel time effect (which in this case is negligible). When any of the planets pass in front of the star along the line of sight, a synthetic light curve is generated (Pál, 2012), which can then be compared to the data. The parameters we include for each planet in the differential evolution Markov chain Monte Carlo (DEMCMC; Ter Braak, 2005) fit are $\{P,T_{0},e^{1/2}\cos(\omega),e^{1/2}\sin(\omega),i,\Omega,R_{p}/R_{\star},M_{p}/M_{\star}\}$, where $P$ is the period, $T_{0}$ is the mid-transit time, $e$ is eccentricity, $\omega$ is the argument of periapse, $i$ is inclination to the sky plane, $\Omega$ is nodal angle on that plane, and $R$ and $M$ are radius and mass respectively (with subscripts $p=b,c,d,e,f$ for the planets and $\star$ for the star). The star had five additional parameters: $\{M_{\star},R_{\star},c_{1},c_{2},dilute\}$, where $c_{i}$ are the two quadratic limb-darkening coefficients and $dilute$ is the amount of dilution from other nearby sources.

We put physically sensible, but permissive, minima ($\rho_{p}=0.0$) and maxima ($\rho_{p}=\rho_{\mathrm{Fe}}$) on the bulk planet densities, where $\rho_{p}$ is a planet’s bulk density and $\rho_{\mathrm{Fe}}$ is the density of iron for a body of planet p’s size. Taking values from Seager et al. (2007), the maximum densities for the 5 planets from b to f respectively are (9.5, 9.7, 10.3, 10.5, 12.2) g/cm3, differing due to the compressibility of iron. The prior on mass is otherwise flat between 0 and these values.

Since the mass and eccentricity implied by TTVs may be degenerate (Lithwick et al., 2012) and result in measured eccentricity values so high that the system go unstable on timescales much shorter than the age of the system Pu & Wu (2015), we use a Rayleigh prior on the eccentricity of all planets with width parameter $\sigma=0.02$. This is consistent with the values measured in other tightly packed planetary systems (Hadden & Lithwick, 2014; Fabrycky et al., 2014), and is consistent with long term stability because even moderate eccentricity has been shown to destabilize tightly packed systems (Pu & Wu, 2015).

The dilution is well-measured (Campante et al., 2015) so we fix $dilute=0.0394$ since it would otherwise be highly degenerate with $R_{p}/R_{\star}$. This implies that the $R_{p}/R_{\star}$ value uncertainties may be slightly underestimated, but since Campante et al. (2015) report $dilute=0.0394\pm 0.0001$, this will have only a very small effect on the reported posteriors. We also fix $\Omega_{p}=0$ for all planets since we expect very small mutual inclinations between the planets because we see five planets transit (see, e.g., Lissauer et al., 2011). Additionally, even modest values of $\Omega$ may greatly increase likelihood that the system becomes unstable over the system’s lifetime due exchange of eccentricity and inclination on secular time scales. Since the transit information gives only the stellar density and planet-to-star mass ratio (via TTVs), we model with a fixed $M_{\star}=0.758M_{\odot}$, which sets the overall scale of the system. We use generic flat priors in all other parameters.

3. Results

Median values and 68.3% confidence intervals from the of photodynamic model are reported in Table 2. The full data set of the DEMCMC chains can be downloaded from the online version of this article. We ran a 64-chain DEMCMC for 900,000 generations recording every 1,000 generation, conservatively throwing out the first 50,000 generations as a burn-in. The autocorrelation timescale for the slowest converging parameters was approximately 60,000 generations, thus we are left with $\gtrsim 850$ independent samples for each parameter. By numerically fitting the TTVs, the model produces mass constraints based on the *Kepler *data. Consistent with the measurement of individual transit times described in §2, planets b, c, and f do not induce significant TTVs on the other planets, which means their masses are not significantly detected. However, the TTVs in planets d and e are both significant enough to confidently place upper and lower bounds on the mass. The posteriors in mass are inconsistent with zero mass and fall off much more rapidly than the prior near $m=0$. Photodynamically measured TTVs are shown in Fig. 1, with the resulting mass constraints for planets d and e shown in Fig. 2 and reported for all planets in Table 2. Compared to the masses derived in Hadden & Lithwick (2016, $M_{d}=0.2^{+0.5}_{-0.1}M_{\mathrm{\oplus}}$ and $M_{e}=0.1^{+0.2}_{-0.1}M_{\mathrm{\oplus}}$), these new measurements ($M_{\mathrm{d}}=0.036^{+0.065}_{-0.020}M_{\oplus}$ and $M_{\mathrm{e}}=0.034^{+0.059}_{-0.019}M_{\oplus}$) are more precise. This is due in part to more a more restrictive and physical prior and in part to due to the photodynamic analysis method used in this study.

We compute the posterior of $Z_{j+1,j}$ for each neighboring planet pair by approximating the value as $|z_{j+1}-z_{j}|/\sqrt{2}$ (see, e.g., Hadden & Lithwick, 2016, Eq. 4), where $z_{j}=e_{j}e^{i\omega_{j}}$ for each planet $j$ and $i$ is the imaginary unit. We find the median and 68% confidence intervals or upper limits $Z_{c,b}=0.022^{+0.013}_{-0.011}$, $Z_{d,c}=0.021^{+0.013}_{-0.010}$, $Z_{e,d}\leq 0.023$, $Z_{f,e}\leq 0.020$. We note that the interior planets are consistent with the prior alone, but the planet pair with detected masses (d and e) has a smaller value preferring low free eccentricity (Fig. 2). The absence of measurable TTVs induced by planet f on planet e also constrains $Z_{f,e}$.

We also numerically integrate 100 draws from the DEMCMC posterior for 100 Myr to make sure we are exploring regions of parameter space stable for times comparable to a reasonable fraction of the system’s age. 95% of the samples remain stable. Since approximately equal numbers of systems are likely to go unstable in logarithmic bins of time (Pu & Wu, 2015), we expect $>80$% of our posterior to be stable for the measured system age of $\sim$11 Gyr.

Importantly, we note that at the 95% confidence level, both planets are inconsistent with being purely iron. Planet d requires a composition with a fraction of rock as least as great as Mercury ($\sim$30%), and, like Earth, planet e can be no more than 30% iron by mass. These measurements are plotted in Fig. 3 along with theoretical composition tracks taken from Seager et al. (2007).

4. Followup Observations

4.1. Radial Velocities

The radial velocity (RV) signal induced on a host star by a planet is given by Cumming et al. (1999):

[TABLE]

where $K$ is the RV amplitude, $G$ the Newtonian gravitational constant, $P$ the planet’s period, $M_{\star}$ the stellar mass, and $M_{\mathrm{p}}$ the planet’s mass. Inserting values for for Kepler-444 planets, we see that the expected $K$ values range from $\sim$4-20 cm/s. This is below the current RV detection threshold (e.g., Plavchan et al., 2015).

4.2. PLATO

Because of the shallow transit depth, photometric follow-up is precluded for most existing instruments. However, the ESA’s Planetary Transits and Oscillations of Stars Mission (PLATO) has recently received approval with operational dates of 2024-2020111See PLATO SCIRD - http://sci.esa.int/plato/42730-scird-for-plato/.. The precision goal for PLATO is $3.4\times 10^{-5}$ in 1 hour for stars with $m_{\mathrm{V}}\leq 11$. Since Kepler-444 is 2 magnitudes brighter, we may expect a factor of $\sim$10 times more photons and thus a precision of $1\times 10^{-5}$ per hour. Each planned 50 second exposure should therefore have a precision of $\sqrt{3600/50}\times 10^{-5}\approx 8\times 10^{-5}$. Taking several solutions from the *Kepler *data posteriors, based on the planned observing strategy we produce 2-year sets of simulated PLATO transits beginning in 2025. We then add Gaussian noise to this data with $\sigma=8\times 10^{-5}$. Finally, we refit the combined actual *Kepler *data and simulated, noisy PLATO data to test how informative the PLATO measurements will be in further constraining the planet masses. We find that the mass constraints of planets d and e are improved to having $\sim$20% 1-$\sigma$ uncertainties. Such a measurement may allow tight constraints on the fraction of the planet which is iron, rocky, or volatile, potentially distinguishing a water-rich planet from an Earth-like composition. Additionally, we find that in some cases Planet b (the smallest radius planet $R_{b}=0.406\pm 0.013R_{\oplus}$) interacts with Planet c sufficiently to induce observable TTVs and a 99.7% confidence (3-$\sigma$ equivalents) non-zero mass detection of Planet b. Such a measurement would make it (as of right now) the smallest exoplanet with a detected mass orbiting a main sequence star. To conclude, we note that the results in this section are dependent on the true noise properties and observing strategy of PLATO, which are currently uncertain.

5. Implications for Formation and Tidal Evolution

Papaloizou (2016, hereafter P16) performs an in-depth analysis of the possible migration history of the Kepler-444 system, considering both migration and circularization effects due to planet-disk interactions. Since the planets are very low mass, P16 assumes they are in the Type I migration regime with migration timescale, $\tau_{\mathrm{mig}}\propto M_{\mathrm{p}}^{-1}$. If the planets migrate at different rates (due to mass and local disk density), then one would expect them to approach MMRs with other planets, at which point they would get trapped near those MMRs (Melita & Woolfson, 1996; Lee & Peale, 2002; Terquem & Papaloizou, 2007). Since the planets are up to 2% away from resonance, P16 speculate that significant relative contraction of the planets did not occur, although significant migration as a unit might have. In order to match the observed period ratios, P16 assumes that planet e is significantly (by a factor of $\sim$3) more massive than d. This allows e to easily migrate more quickly than, and thus contract and approach resonance with, planet d while the other planets remain relatively more distant from resonances. Our photodynamical fit finds that $M_{\mathrm{e}}/M_{\mathrm{d}}=0.93^{+0.14}_{-0.13}$, a significant departure from that assumption. This suggests that the present-day observed period ratios combined with smooth disk migration alone are generally insufficient for modeling specifics of the formation of the system. Many factors may have changed the migration of the planets while the disk was present, including local disk properties (Cossou et al., 2014) or turbulence in the disk (Oishi et al., 2007; Rein & Papaloizou, 2009). Alternately, the planets may have moved after the dispersal of the gas and dust disk, for instance via a combination of planetesimal crossings (Fernandez & Ip, 1984; Levison et al., 2007) or damping from tides raised by the star (Lee et al., 2013). Therefore, we caution against strict interpretations of observed exoplanet masses and architectures (or ensembles of these architectures) when it is likely that the systems have evolved substantially since their natal formation. We infer from the $M_{\mathrm{e}}/M_{\mathrm{d}}$ ratio that the system underwent significant orbital period changes after a migration formation, or formed in situ. We also note that very high, iron-like densities are disfavored, suggesting that large amounts of collisional stripping due to high velocity giant impacts likely did not occur (Marcus et al., 2010; Asphaug & Reufer, 2014).

Since the planets orbit very close to their host star, we consider the effects of tidal dissipation on the observed orbital period ratios. It is possible that tides on planets in or near a MMR causes their proximity to orbital resonance to change (generally spreading planets apart away from resonance) over Gyr timescales (Papaloizou, 2011; Lee et al., 2013). Following Papaloizou (2011, henceforth P11), we define $\delta_{j}$ as the distance from orbital resonance by

[TABLE]

where $n_{j}$ is the $j^{th}$ planet’s mean motion and $k$ is the degree of the near first order resonance between planets $j$ and $j+1$. P11 equation (40) gives the relation between the change in $\delta_{j}$ as a function of time and orbital parameters of the system. To determine analytically the amount tides would move planets away from exact resonance as a function of time (equation (42)), P11 integrates equation (40) from $t^{\prime}=0$ to $t^{\prime}=t$ and assumes $\delta_{j,t=0}=0$, i.e., the system begins in exact MMR. If, however, we integrate from $t^{\prime}=0$ to $t^{\prime}=11$ Gyr (the age of Kepler-444), and we know $\delta_{j,t=11Gyr}$ based on the observed system, we may solve for $\delta_{j,t=0}$ as a function of $Q/k_{2}$, the ratio of the tidal $Q$ factor and the love number. This factor enters via the tidal circularization time

[TABLE]

for the $j^{th}$ planet (Goldreich & Soter, 1966; Rasio et al., 1996). We solve for the total change in distance from resonance since the planets’ formation $\Delta_{j}=\delta_{j,t=11Gyr}-\delta_{j,t=0}$. For the inner pair of planets (b and c, $k=4$), we find that $\Delta_{1}\approx 7\times 10^{-4}-7\times 10^{-7}$ for values of $Q/k_{2}$ ranging from 1-1000, using the approximation that $(Q/k_{2})_{b}\approx(Q/k_{2})_{c}$ which is reasonable given their similar size and proximity in the system. In the solar system, the rocky planets and large, rocky moons have $10\lesssim Q/k_{2}\lesssim 500$ (Goldreich & Soter, 1966). Since the observed $\delta_{1,t=11Gyr}=1.27\times 10^{-2}$, we see that tidal dissipation was insufficient to have moved the innermost pair a significant distance from its current period ratio and rules out tidal dissipation breaking a natal MMR. These findings are confirmed by long-term numerical N-body integrations, following MacDonald et al. (2016).

The other pairs of planets have longer periods, and in the case of c and d, are further from resonance. They are thus generally less affected by tides. However, the period ratio of planets d and e are very close to resonance (Table 1) so even a small amount of dissipation may significantly impact their $\delta_{3}$. Following Lee et al. (2013), we can set a limit on the tidal $Q/k_{2}$ factor for the innermost planet by using their equation (18) with the observed system age and planet parameters. We find $(Q/k_{2})_{d}\gtrsim 12$. This limit is very near solar system values for rocky bodies, and possibly hints that the pair started in an exact MMR and was driven apart via this mechanism. This suggests that disk migration may have driven this pair of planets together, but the lack of tidally-broken commensurabilities among the other planets suggests the migration was not smooth or there were significant external perturbations after the disk dissipated.

We thank an anonymous referee for insightful comments which considerably added to the quality of this manuscript. This research is supported by Grant NNX14AB87G issued through NASA’s *Kepler *Participating Science Program.