# Kepler Planet Masses and Eccentricities from TTV Analysis

## Abstract

We conduct a uniform analysis of the transit timing variations (TTVs) of 145 planets from 55 Kepler multiplanet systems to infer planet masses and eccentricities. Eighty of these planets do not have previously reported mass and eccentricity measurements. We employ two complementary methods to fit TTVs: Markov chain Monte Carlo simulations based on -body integration and an analytic fitting approach. Mass measurements of 49 planets, including 12 without previously reported masses, meet our criterion for classification as robust. Using mass and radius measurements, we infer the masses of planets’ gaseous envelopes for both our TTV sample as well as transiting planets with radial velocity observations. Insight from analytic TTV formulae allows us to partially circumvent degeneracies inherent to inferring eccentricities from TTV observations. We find that planet eccentricities are generally small, typically a few percent, but in many instances are non-zero.

###### Subject headings:

planets and satellites: detection## 1. Introduction

The Kepler mission’s census of exoplanetary systems has provided us with a statistical picture of the properties of planetary systems (Borucki et al., 2010). Small planets with radii in the range on short-period orbits, days, are among the most abundant, occurring around roughly half of Sun-like stars (e.g., Petigura et al., 2013; Fressin et al., 2013). The distributions of planet sizes and periods measured by Kepler provide important constraints for theories of planet formation and evolution. Mass measurements, which constrain planets’ compositions, and eccentricity measurements, which reveal the current dynamical states of planetary systems, are also essential clues for understanding formation and evolution processes.

The radial velocity (RV) method has provided the majority of exoplanet mass and eccentricity measurements to date (e.g., exoplanets.org; Wright et al., 2011). Mass measurements of transiting super-Earth and sub-Neptune planets are of particular interest given the ubiquity of such planets and their absence from our own solar system. Many sub-Jovian transiting planets have had their masses measured through radial velocity follow-up (e.g., Borucki et al., 2010; Bordé et al., 2010; Bakos et al., 2010; Hartman et al., 2011; Gautier et al., 2012; Gilliland et al., 2013; Pepe et al., 2013; Barros et al., 2014; Dumusque et al., 2014; Moutou et al., 2014; Alonso et al., 2014; Kostov et al., 2014; Marcy et al., 2014; Esteves et al., 2015; Dressing et al., 2015; Gettel et al., 2016; Sinukoff et al., 2016). These efforts have yielded important constraints on the mass-radius relationship of super-Earth and sub-Neptune planets: Weiss & Marcy (2014) infer a mass-radius relationship of planets smaller than using available RV mass determinations, supplemented with a handful of masses determined from transit timing. Rogers (2015) uses the sample of Kepler planets with Keck HIRES RV follow-up to infer that planets transition from mainly rocky to volatile-rich compositions above a size of . However, the RV method is of limited applicability to Kepler’s many sub-Neptune planets since the radial velocities induced by such planets often require intense follow-up efforts to detect. With a handful of exceptions, RV mass determinations of sub-Neptunes have been limited to planets with orbital periods shorter than 20 days.

The limitations of RV are even more acute for measuring sub-Neptunes’ eccentricities since high signal-to-noise ratios are required to accurately infer eccentricities with RV (Shen & Turner, 2008; Zakamska et al., 2011). Eccentricities of transiting planets can also be measured by modeling transit light curves. However, light curve modeling yields useful constraints only in special circumstances: giant planets with large eccentricities (Dawson & Johnson, 2012), planets with occultation detections (Shabram et al., 2016), or host stars with strong density constraints from asteroseismology (Kipping, 2014). Statistical analyses of the transit durations of planets around well-characterized host stars have been used to constrain planet samples’ overall eccentricity distributions (Van Eylen & Albrecht, 2015; Xie et al., 2016), but inferring the eccentricities of individual sub-Neptune planets from light curve modeling remains difficult.

Transit timing variations (TTVs) are a powerful tool for measuring masses and eccentricities in multi-transiting systems (Agol et al., 2005; Holman, 2005). The large TTV amplitudes induced in planets near mean motion resonances (MMRs) can probe the masses and eccentricities of small planets at relatively long orbital periods which would otherwise be difficult or impossible to measure via the RV method. However, inverting TTVs to infer planet properties poses a difficult parameter inference problem: it requires fitting a large number parameters, often with strong degeneracies, to noisy data. Statistical analyses of samples of TTV systems can overcome some of these difficulties (Wu & Lithwick, 2013; Hadden & Lithwick, 2014). Alternatively, the parameter inference challenge can be addressed with Markov chain Monte Carlo (MCMC) simulations when fitting TTVs of individual systems. MCMC is well suited for high-dimensional parameter inference problems and has been used frequently in TTV studies (e.g., Sanchis-Ojeda et al., 2012; Huber et al., 2013; Masuda et al., 2013; Schmitt et al., 2014; Jontof-Hutter et al., 2015, 2016; Hadden & Lithwick, 2016; Mills et al., 2016).

While MCMC can efficiently sample planet masses and orbits consistent with TTV observations, interpreting MCMC results is often complicated by strong parameter correlations and sensitivity to priors. Analytic TTV formulae identify degeneracies inherent to inverting TTVs and aid the interpretation of -body MCMC results. Components of the TTV signal responsible for mass and eccentricity constraints can be identified with analytic formulae and lend support to the robustness of -body results to different prior assumptions.

In this paper we compute MCMC fits to the TTVs of 55 Kepler multiplanet systems exhibiting significant TTVs, 33 of which do not have -body TTV fits reported previously in the literature. In addition, our work provides a uniform treatment of TTV systems that have previously been analyzed elsewhere. We complement our MCMC fits, which rely on -body integrations, with an analytic approach to TTV modeling. The paper is organized as follows: we review analytic TTV formulae in Section 2 and describe our fitting methods in Section 3. Results of our fits are described in Section 4 and discussed in Section 5: in Section 5.1 we use planet radii and masses derived from our TTV fits to infer the masses of planets’ gaseous envelopes. In Section 5.2 we briefly discuss implications of the eccentricities inferred from TTVs. We conclude in Section 6.

## 2. analytic TTV

### 2.1. The Analytic TTV Formula

A number of authors have derived analytic formulae to approximate TTVs using perturbative methods (e.g., Agol et al., 2005; Nesvornỳ & Morbidelli, 2008; Nesvornỳ, 2009; Lithwick et al., 2012; Deck & Agol, 2015, 2016; Agol & Deck, 2016; Hadden & Lithwick, 2016; Nesvornỳ & Vokrouhlický, 2016). Analytic formulae aid the interpretation of -body fitting results by elucidating degeneracies and identifying TTV features that constrain planet masses and eccentricities. In this paper we apply the analytic formulae derived in Hadden & Lithwick (2016, hereafter Paper I) as part of our analysis of each system’s TTVs. The main features of these formulae are summarized below.

For clarity, we focus our discussion on a planet near a :-1 first-order MMR with an exterior perturber. Other configurations are discussed in Section 2.3 The analytic formulae express a planet’s TTV as a sum of harmonic terms:

(1) |

where the three successive terms are called ‘fundamental’, ‘second-harmonic’, and ‘chopping’ TTV’s, and ’’ denotes complex conjugate. The frequencies of the harmonics are expressed in terms of the ‘super-period’, , and synodic period, . These depend only on the period of the planet and its perturber, and hence may be considered to be “known” for transiting planets. The dependence on masses and eccentricities is contained entirely in the amplitudes , , and . Since the first two of these are complex numbers, and the third real, there are in total five observables that can be used to infer masses and eccentricities. (The coefficients depend only on the period ratio of the two planets).

The fundamental TTV is typically the dominant component. Its complex amplitude depends on mass and eccentricity and is given by

(2) |

where and are coefficients that depend only on the planets’ periods, is the perturber’s planet-star mass ratio, and

(3) | |||||

(4) |

where and are the free complex
eccentricities^{1}

### 2.2. Degeneracies

In attempting to extract planet parameters from the dominant (fundamental) harmonic of the TTV, there are at least two kinds of degeneracies (Lithwick et al., 2012). The first is between the two planets’ individual complex eccentricities — the eccentricity-eccentricity degeneracy. Since the fundamental TTV depends on the planets’ eccentricities only through the single linear combination , the and of each planet is not separately measurable. This degeneracy persists even with the inclusion of the chopping and second harmonic terms.

The second degeneracy is between and , which we call the mass-eccentricity degeneracy: a smaller can be compensated for by a larger (Eq. (2)). Usually in Equation (2) so that the fundamental TTV amplitude is sensitive to and the degeneracy spans a large range in .

The mass-eccentricity degeneracy can be broken if the chopping component can be resolved in the TTV signal (Nesvornỳ & Vokrouhlický, 2014; Deck & Agol, 2015). The chopping TTV is insensitive to eccentricities and provides a measurement of the perturbing planet’s mass. In fact, the chopping TTV amplitude can simply be equated to the perturber’s planet-star mass ratio,

(5) |

with appropriate choice of scaling for the coefficients in Equation (1). If a perturbing planet’s mass can be measured from the chopping TTV, then the fundamental TTV amplitude gives the planets’ combined eccentricity (Eq. (2)).

The mass-eccentricity degeneracy may alternatively be broken by detecting the second-harmonic component, the amplitude of which is given by

(6) |

with coefficients and depending only on planet periods.^{2}

### 2.3. Other configurations

With minor modifications, Equation (1) also describes the TTVs of planets subject to an interior perturber. Thus, breaking the mass-eccentricity degeneracy for one planet determines the pair’s and thereby automatically breaks the degeneracy for both planets.

Equation (1) without the fundamental TTV term also describes the TTVs of planets near a : second-order MMR. Proximity to a second-order resonance enhances the term in Equation (1), which in this context we refer to as the ‘second-order resonance’ component rather than second-harmonic component. TTVs near second-order MMRs exhibit essentially the same mass-eccentricity and eccentricity-eccentricity degeneracies as planets near first-order MMRs and the mass-eccentricity degeneracy can be broken by measuring a chopping component in either planet’s TTV.

We find that, in practice, measuring (only) fundamental TTVs of three or more planets generally does not resolve any of the degeneracies inherent to the two-planet case.

## 3. TTV Fit Methods

We analyze TTVs of 145 planets from 55 different Kepler systems. Figure 1 shows the TTV systems fit in this paper. Our analysis is based on transit times computed by Rowe et al. (2015) from long cadence Kepler data spanning Quarters 1-17. Each system is fit with both -body MCMC simulations and the analytic TTV formulae. Systems are discussed individually in Appendix A. Results for their masses and eccentricities are summarized in Tables 1 and 2.

Our MCMC simulations sample the most likely planet masses and orbital elements of a multi-transiting system given the planets’ transit times by using -body integrations. Specifically, for a system of planets our MCMC samples the posterior distribution of each planet’s planet-to-star mass ratio, , eccentricity vector components and , initial osculating period, , and time of first transit, where . The complete details of our MCMC implementation are given in Paper I.

As in Paper I, in order to assess how robustly planet properties are constrained by the TTVs we run two MCMC simulations for each system we fit with two different priors for masses and eccentricities. Our first (default) prior is logarithmic in planet masses () and uniform in eccentricities (). Our second (high mass) prior is uniform in planet masses () and logarithmic in eccentricity (). We will refer to the posterior distributions computed with the respective priors as the default and high mass posteriors. Inferred planet masses are classified as ‘robust’ if the 1 credible region of the default posterior excludes and includes the peak of the high mass posterior. Masses of 49 of the 145 planets are classified as robust. Twelve planets with robust masses do not have -body fits previously reported in the literature.

We have chosen our two MCMC priors to weight toward opposite extremes of the mass-eccentricity degeneracy in order to assess the significance of this degeneracy in each system. Additionally, both priors are chosen to be “uninformative” or broad and do not rely on any theoretical predictions for planet compositions or eccentricities. Aside from these considerations, our adopted priors are arbitrary and any number of other priors could be reasonably advocated. It is possible that a different choice of priors could affect the inferred masses and eccentricities significantly. Instead of exploring a wide variety of priors with -body MCMC, which is not computationally feasible, we complement MCMC simulations with fits using the analytic formulae. The analytic fits convert the amplitudes , , and measured from TTVs to constraints in the - plane. Constraints derived from our analytic fits are compared to the results of the -body MCMC simulations for each system in Appendix A. The analytic fits lend further confidence to the robustness of MCMC results by showing the degree to which various components of each planet’s TTV constrain planet parameters.

TTVs do not probe planet masses directly but rather planet-star mass ratios. Therefore we will often find it convenient to give planet masses in units of

(7) |

where is the ratio of host star mass to solar mass. Employing as the unit highlights mass uncertainties inherent to the TTVs separately from those arising from uncertain stellar properties.

Figure 2 compares planet masses inferred using our default versus high mass priors. Many mass inferences are seen to depend sensitively on the assumed prior, reflecting the degeneracies often inherent to inverting TTV observations.

## 4. Results: Masses and Eccentricities

Figure 3 shows inferred TTV planet masses on a mass-radius plot. A sample of transiting planets with RV-measured masses are shown as well for comparison. Most of the planets are less massive than and exhibit a wide diversity of radii. Many of the planets are less dense than a hypothetical pure water-composition planet, necessitating the presence of substantial gaseous envelopes to explain their low densities.

Compared to small transiting planets with RV mass measurements, our TTV sample typically finds lower density planets as illustrated by Figure 4. Although this trend has been noted before, (e.g., Weiss & Marcy, 2014; Dai et al., 2015), it is sometimes suggested that this discrepancy is due to mass errors in either the TTV or RV measurements. But Figure 4 demonstrates that much of the discrepancy can be explained if planets farther from the star tend to be less dense—perhaps because they are less affected by photoevaporation. One should note, however, that the distribution of planets in Figure 4 is also affected by a variety of selection effects. For example, the lack of dense planets at long orbital periods could reflect the decreased transit detectability of small planets with long periods (e.g., Gaidos & Mann, 2012). Also, differences between RV and TTV populations could reflect the different dependence of the two techniques on planet mass and radius (e.g., Steffen, 2016). Nonetheless, both TTVs and RVs find planets of similar density between periods of where there is significant overlap between the two samples. A simple two-sample Kolmogorov-Smirnov test comparing all densities measured with the two methods gives a probability of that they are drawn from the same underlying distribution, while the probability increases to if both samples are restricted to the period range .

We turn now to the planets’ eccentricities. As discussed in Section 2, although TTV observations do not strongly constrain individual planet eccentricities, they can constrain the combination . This combination may be considered as a surrogate for the individual planets’ complex eccentricities unless planets have , that is, comparable eccentricities and aligned orbits. We expect planets’ complex eccentricities have random relative orientations (see discussion in Paper I) so that, overall, values are a reliable surrogate for eccentricities. The inferred values of are summarized in Figure 5. The majority of eccentricities are inferred to be small: the median of all the posterior samples shown in Figure 5 is . While many planet pairs’ combined eccentricities are small, they are frequently inconsistent with zero.

Table 2 lists combined eccentricities for individual planet pairs. In addition to , it is of interest to know which planets are consistent with ; such planets might have experienced significant damping by tides or other effects. Credible regions in cannot be used to address this question because must be non-negative. Therefore, following Zakamska et al. (2011), we define the signed quantity, , which is the projection of the s from the MCMC posterior onto the median of their distribution. More precisely, we define the median , by computing the median real and imaginary components of . Then, given , the value of is defined as

(8) |

where the ’’ indicates complex conjugate. Zakamska et al. (2011) show that an analogous quantity is useful for recovering solutions in the analysis of radial velocity data generated from circular orbits. Of the 90 adjacent planet pairs in our TTV sample, more than 60% have s inconsistent with 0 at 2 confidence for both the default and high mass posteriors.

## 5. Discussion

### 5.1. Gaseous Envelopes

Many of the planets shown in Figure 3 have such low densities that they must have massive gaseous envelopes. Such massive envelopes were likely accreted when the protoplanetary disk was still present; outgassing can only produce envelopes less massive than 2% of a planet’s total mass (Rogers et al., 2011). Here we convert observed masses and radii to envelope masses by interpolating the tables of Lopez & Fortney (2014). Those authors simulate envelopes on top of Earth-composition cores; as time proceeds, their envelopes cool, while being irradiated by the star. Thus the masses and radii of of their envelopes are essentially related by the fact that the cooling time should be comparable to the age. Although there are multiple unknowns that will affect the predicted envelope mass (especially atmospheric opacity and composition and, to a lesser extent, the core composition and heating/cooling) (see, e.g. Rogers et al., 2011; Howe et al., 2014; Lopez & Fortney, 2014), the results of their simulations are adequate for our purposes.

Inferred envelope mass fractions, , are plotted in Figure 6 for planets more massive than and smaller than with either RV or robust TTV mass measurements. Lopez & Fortney (2014) find that the radii of planets with envelope mass fractions 1% are quite insensitive to core mass. Therefore, the diversity of gas fractions in Figure 6 reflects the diversity of radii seen in Figure 3. On the other hand, inferred envelope mass fractions 1% are likely sensitive to our assumption of Earth-composition cores and some of the envelope mass fractions in Figure 6 may actually reflect a diversity of core compositions.

How do the inferred envelopes shown in Figure 6 compare to theoretical predictions of gas envelope accretion? A number of studies have examined envelope accretion by planetary cores embedded in a gaseous protoplanetary disk (e.g. Ikoma & Hori, 2012; Bodenheimer & Lissauer, 2014; Ginzburg et al., 2016; Lee & Chiang, 2016). These studies find that super-Earth cores readily accrete envelopes between a few to tens of percent of their total mass over the lifetime of a protoplanetary disk, similar to those shown in Figure 6. Planets in Figure 6 that are consistent with having no envelope generally receive higher fluxes than similar-mass planets covered by significant envelopes. This trend could result from photoevaporation removing the envelopes initially accreted by planets close to their stars (e.g. Lopez et al., 2012; Owen & Wu, 2013). Among planets that do have significant envelopes, there is no obvious trend in envelope fraction with planet mass. This suggests a wide diversity of factors that influence envelope accretion and retention. We plan to compare observed envelope masses with theoretical predictions of accretion and mass-loss in a future work.

### 5.2. Eccentricities

Overall, our results show that most planets have relatively small eccentricities of a few percent. This is consistent with the results of past TTV studies (Wu & Lithwick, 2013; Hadden & Lithwick, 2014) and transit duration analyses (Van Eylen & Albrecht, 2015; Xie et al., 2016) of multiplanet systems. These studies infer the overall eccentricity distribution of their samples by fitting Rayleigh distributions: Wu & Lithwick (2013) find a mean eccentricity for a sample of 44 planets, Hadden & Lithwick (2014) find for a sample of 139 planets, Van Eylen & Albrecht (2015) find for a sample of 74 planets, and Xie et al. (2016) find for a sample of 330 planets. In contrast to these population-level studies, we are able to measure the combined eccentricities of individual planet pairs, often with uncertainties on the order of a percent or better. Both Van Eylen & Albrecht (2015) and Xie et al. (2016) suggest that TTV systems’ proximity to resonance may imply eccentricities that are distinct from the larger population of multiplanet systems. While it is true that TTV systems are preferentially closer to resonance (see Figure 8 in Appendix A), they are not as radically distinct from the larger multiplanet population as is often suggested.

As noted above, a majority of our TTV sample (roughly ) have eccentricities that are non-zero at 2 confidence. This seemingly contradicts the inference that the planets were born in a disk of gas, which would presumably damp eccentricities. The non-zero eccentricities suggest that planets may have accreted their envelopes from a depleted disk in which gas damping was ineffective (Lee & Chiang, 2016).

Some of the planets in Table 1 orbit their host star at close enough separations that their eccentricities may be affected by tidal dissipation. Tidal dissipation depends steeply on planet period and the shortest period planets are expected to have experienced significant tidal damping. Figure 7 plots (Equation (8)) versus the periods of planet pairs’ inner planets. All but one of the seven shortest period planets within days have combined eccentricities consistent with 0, suggesting circularization by tides. The shortest period planet, Kepler-80 d, has a small but non-zero . Beyond periods of 5 days there is a mixture of combined eccentricities that are consistent with 0 or very nearly so, as well as combined eccentricities that are definitively non-zero.

## 6. Summary and Conclusion

We have presented a uniform analysis of the TTVs of 55 multiplanet systems in order to infer planet masses and eccentricities. We employ both -body and analytic fitting in a complementary approach in order identify the degeneracies inherent to TTV inversion and understand when and how they are broken by components of the TTV signal. We use two sets of MCMC simulations for each systems, each of which uses a prior weighted to opposite extremes of the mass-eccentricity degeneracy predicted by the analytic formulae, in order to classify inferred planet masses as robust or not.

Low mass () planets exhibit a wide diversity of sizes, with many of these planets less dense than a hypothetical pure ice composition planet, indicating the presence of significant gaseous envelopes. The wide diversity of planet sizes can be attributed to a diversity of planet envelope mass fractions. We use our TTV fits together with the evolutionary models of Lopez & Fortney (2014), to convert planet masses and radii to envelope mass fractions. We plan to use our sample of TTV-characterized planets in a future work examining planets’ accretion and retention of gas envelopes during and after dispersal of the protoplanetary disk.

With guidance from the analytic TTV model, we have focused on planets’ combined eccentricities, , rather than individual eccentricities, which are largely unconstrained by TTVs. We find that planets typically have eccentricities of a few percent or less, in agreement with past statistical studies of multiplanet systems (Wu & Lithwick, 2013; Hadden & Lithwick, 2014; Van Eylen & Albrecht, 2015; Xie et al., 2016). The shortest period planets have eccentricities consistent with 0 and thus may have experienced significant tidal eccentricity damping. Detailed modeling of individual TTV systems can potentially shed light on how dynamical processes such as migration and eccentricity damping may have shaped the systems we observe today.

We measure a number of planet masses for planets in the super-Earth/mini-Neptune size range where there is great theoretical interest in understanding the mass-radius relationship. The observational biases of our TTV sample are difficult to assess given the TTV signal’s complicated dependence on periods and eccentricities. We see incorporating TTV non-detections and quantifying selection effects as important directions for future work in understanding the mass-radius relationship.

Planet-star radius ratio from Borucki et al. (2011). b Planet-star radius ratio from Fabrycky et al. (2012). c Planet-star radius ratio from Masuda (2014). d Planet-star radius ratio from Steffen et al. (2013). e Planet-star radius ratio from Mills et al. (2016). f Planet-star radius ratio from Rowe et al. (2014). g Planet-star radius ratio from Kepler KOI Q1-17 data release DR24, hosted on the Exoplanet Archive.