The California-Kepler Survey. VIII. Eccentricities of Kepler Planets and Tentative Evidence of a High-metallicity Preference for Small Eccentric Planets

The exquisite photometric data from the Kepler mission has revolutionized our knowledge of exoplanet demographics. Studies revealed the occurrence rate of planets as a function of orbital period and planet size (Petigura et al. 2013; Zhu et al. 2018), the low mutual inclination of multiplanet Kepler systems (Fabrycky et al. 2014; Zhu et al. 2018), and the orbital period ratio distribution of the Kepler planets (Fabrycky et al. 2014; Steffen & Hwang 2015). Recently, large spectroscopic surveys enabled improved measurements of Kepler planet host star properties (e.g., De Cat et al. 2015; Johnson et al. 2017). These improved stellar properties revealed new details in the planet population, such as the bimodal distribution of planets between 1 and 4 R_⊕ (Fulton et al. 2017), and their eccentricity (e) distribution (Xie et al. 2016). Measuring the eccentricity distributions of different planet populations is important because they are relics of the processes which occurred during the epoch of planet formation and migration.

The eccentricity of planets detected with radial velocity (RV) can be measured via the shape of the Keplerian signal, and a large range of eccentricities have been observed (Winn & Fabrycky 2015). However, due to the difficulty in constraining planet eccentricity at low RV amplitudes (Shen & Turner 2008), generally only in giant planet systems with high signal-to-noise ratio (S/N) RV data can eccentricity be well constrained.

On the other hand, exoplanet eccentricity can be determined in systems where multiple planets transit via a dynamical analysis. The gravitational interaction among planets causes transit timing variations (TTVs), which are sensitive to the orbital period ratios, masses, and eccentricities of the planets (Agol et al. 2005; Lithwick et al. 2012). However, measurements of eccentricity by TTVs often suffer from two limitations: (1) it nearly always requires multitransiting systems and (2) the eccentricity is often degenerate with planet mass (Lithwick et al. 2012). Thus the Kepler data can uniquely determine eccentricities only in rare cases where this degeneracy is broken (Deck & Agol 2015). However, a planet's eccentricity also subtly affects the duration of a planet's transit regardless of any dynamical interactions.

The duration of a planetary transit is determined by the length of the transit chord across the face of its host star divided by the planet's orbital velocity. The transit chord length is given by the radius of the star (R_⋆) and the impact parameter (b) of the transit. The velocity of a planet on a circular orbit is uniquely determined by the planet's orbital period and the stellar mass, assuming M_p ≪ M_⋆. However, an eccentric planet's velocity depends additionally on the eccentricity and the phase of the planet during transit (since the orbital velocity is not constant throughout the orbit). This results in a degeneracy between b and e. Because the impact parameter also affects the transit shape (Winn 2010), careful transit modeling may uniquely determine b and therefore also e. Caution must be taken to account for TTVs, which may also alter the apparent shape of transits phased at constant period (Kipping 2014). This has been done for a subset of Kepler planets which have precisely characterized stellar hosts from asteroseismology (Van Eylen & Albrecht 2015; Van Eylen et al. 2019), resulting in a measurement of the low (≲0.06) eccentricities of ≲4R_⊕ planet pairs and confirming previous results of high mean eccentricity (e ∼ 0.2) for single transiting planet systems (Xie et al. 2016). An alternative approach to the precise characterization of individual planet transits is to use a statistical methodology which eliminates the need for individually measured impact parameters (Ford et al. 2008). This technique has previously been applied to the Kepler data with varying levels of stellar host property precision (Moorhead et al. 2011; Kane et al. 2012; Plavchan et al. 2014; Xie et al. 2016). The most recent result, Xie et al. (2016, X16 hereafter), reveals that single transiting planet systems from Kepler are drawn from a significantly broader distribution of eccentricities than the multiple transiting planet systems. A similar effect has been observed in systems of giant planets in single and multiple planet systems detected via RVs (Howard 2013; Limbach & Turner 2015).

The work in this paper combines the unprecedented population-wide precision and accuracy of the The California-Kepler Survey (CKS)-Gaia stellar sample (Fulton & Petigura 2018) with the most recent Kepler data release (DR25; Coughlin et al. 2017) to improve our knowledge of the eccentricities of the Kepler planets. We simulate populations of exoplanet systems with different eccentricity distributions and compare the resulting transit duration distribution to the distribution of durations observed with Kepler to determine the most likely eccentricity of this population. We use these results to search for trends in planet and host star properties as a function of eccentricity.

The population approach to measure eccentricity distributions from transit durations assumes that the 3D orientation of exoplanet systems is isotropic. Therefore the distribution of impact parameters will be uniform, subject to the observational bias that planets with high impact parameters have shorter and shallower transits and are thus less likely to be detected. Thus by assuming a given population of planets has randomly oriented invariant planes, the deviation from the expected distribution of observed durations reveals the eccentricity distribution of the population. A highly eccentric population will have more short-duration transits than a circular population due to the higher orbital velocities near pericenter combined with the increased transit probability due to the decreased planet-to-star separation (Burke 2008). Numerically integrating Equation (16) from Burke (2008) reveals that only ∼10%–20% of transit durations will be longer than expected for any eccentricity of ∼0.1–0.8 randomly observed from different orientations many times. This long-duration subset also favors durations just longer than the b = 0, e = 0 geometry, and can therefore be difficult to distinguish from a circular orbit when measurement errors are considered. If, however, a planet's transit is longer than expected for the b = 0, e = 0 case even when accounting for all uncertainties, a high eccentricity with transit near apocenter is a unique conclusion. On the other hand, the plentiful short-duration transits may each individually be explained by a high b rather than e. But by enforcing the assumption of random viewing orientations for a planet population, the population's high-e distribution may be determined by the surplus of these short-duration events over that expected from an e = 0 population. This method is also not biased by transit timing variations since time-shifts in the transits do not change the average transit duration.

To construct and compare our simulated populations to the observed Kepler population, we closely follow the methodology described in X16. We summarize the steps as follows.⁶

2.1. Single Transiting Planet Systems

We adopt a truncated Rayleigh (TR) distribution of eccentricities in our model since it is commonly used for eccentricities (Ford et al. 2008; Moorhead et al. 2011; Fabrycky et al. 2014; Xie et al. 2016) and thus readily comparable to results from the literature. This choice is physically motivated (Shabram et al. 2016), and has support over the allowed range of [0, 1]. The probability distribution function (PDF) of the TR is

$$\begin{eqnarray}{\mathrm{PDF}}_{\mathrm{TR}}({\sigma }_{e})=\left\{\begin{array}{ll}\tfrac{e}{{\sigma }_{e}^{2}}\exp \left(-\tfrac{{e}^{2}}{2{\sigma }_{e}^{2}}\right) & 0\lt =\,e\lt 1\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 1 }$$

For each system, we (1) draw an eccentricity from a the TR distribution with a specified width parameter, σ_e. (2) We draw an invariant plane and pericenter (ω) direction of the system at random. We exclude cases where the planet is guaranteed not to transit. (3) We compute the resulting transit duration. If the duration is unphysical due to, e.g., the pericenter of the planet hitting the star, we begin again from step 1. We also consider if the resulting transit would be detectable by computing the expected S/N of the simulated transit as $\sqrt{{d}_{\mathrm{sim}}/{d}_{\mathrm{obs}}}\ast {\mathrm{SN}}_{\mathrm{obs}}$, where d is the transit duration with subscripts for the simulated (sim) and observed (obs) cases, and SN_obs is the measured S/N of the Kepler object of interest (KOI). If the simulated transit S/N does not meet the threshold for inclusion in the Kepler catalog (>7.1), we start again from step 1. (4) We then compute the duration ratio (r) of the simulated transit to a transit with e = 0 and b = 0, as well as the observed duration ratio to a nominal e = 0, b = 0 transit. The likelihood of the observed transit duration ratio is computed as

$$\begin{eqnarray}&&{{ \mathcal L }}_{r}=\exp \left(-{\left(\displaystyle \frac{{r}_{\mathrm{sim}}-{r}_{\mathrm{obs}}}{{\sigma }_{\mathrm{obs}}}\right)}^{2}\right),\end{eqnarray} \tag{ 2 }$$

where σ_obs is the uncertainty of r_obs.

These steps are repeated 10,000 times for each KOI to probe the distribution of simulated durations for a given system and eccentricity distribution. The likelihoods are then multiplied to give a total likelihood for the e distribution given the observed transit duration. This procedure generates a Monte Carlo approximation of the integral

$$\begin{eqnarray}&&{{ \mathcal L }}_{{\sigma }_{e}}={\int }_{r}P({r}_{\mathrm{sim}}| {\sigma }_{e})P({r}_{\mathrm{obs}}){dr}.\end{eqnarray} \tag{ 3 }$$

We repeat this procedure 100 times and average the values of the trials at each σ_e for a large grid of σ_es to determine the overall ${ \mathcal L }({\sigma }_{e})$ for a given KOI.

Our catalog of KOIs includes all planets and planet candidates from Kepler Data Release 25 (Coughlin et al. 2017). To prevent spurious results from binary stars and other false positives, we remove any KOIs from our sample with R_p > 15 or that has either a false positive flag of 1 or a false positive probability (FPP) > 0.5 in the Morton et al. (2016) catalog. We find that the results are not sensitive to the exact FPP cutoff chosen. We also consider the effect of changing the S/N cutoff of 7.1 to reduce the number of false positives caused by noise. However, we find our results do not change significantly when higher S/N cutoffs are chosen (10 or 12), and therefore retain 7.1 as the nominal cutoff.

We restrict our sample to 0.5R_⊙ < R_⋆ < 2.0R_⊙, which includes the majority of Kepler targets and reduces the chance of calibration errors. We remove any systems whose uncertainty in R_p exceeds 0.02 R_⋆—this cutoff removes the handful of systems which may be unreliably measured in the Kepler data set due to limb-darkening degeneracies. Finally, we remove stars that Fulton & Petigura (2018) identify as having ≥5% contamination from nearby sources, as these stars may sometimes have biased or incorrect radii. The exact value of these cutoffs does not significantly affect our results. We are left with 439 single transiting KOIs whose ${{ \mathcal L }}_{{\sigma }_{e}}$ we sum to get the overall population likelihood as a function of σ_e.

Due to the possibility of unaccounted for systematics, we conservatively report 2σ equivalent uncertainties on e values.

2.2. Multiple Transiting Planet Systems

Our approach is fundamentally the same for the systems with two or more observed transiting planets as for single transiting planet KOIs, with the exception that we must also take into account the mutual inclination dispersion among the multiplanet systems which will also affect the transit durations. For instance, for extremely coplanar planetary systems, by geometric arguments alone, an interior planet must have a smaller impact parameter than an exterior planet. This is not true for a significantly mutually inclined system. Thus the mutual inclination of the systems affects the distribution of assumed underlying impact parameters, which is no longer uniform. To address this we modify step (2) above by drawing the mutual inclination and thus orbital plane of each planet in a system from a Gaussian distribution with width σ_i centered around the invariant plane. When we reach step (3), we check that every planet in the system meets the S/N threshold and restart the entire system from step (1) if any do not.

We then perform a grid search over both σ_e and σ_i to determine the likelihood surface in terms of both mutual inclination and eccentricity dispersion. To make the process more computationally tenable, for these systems we use 20 trials of the 10,000 point Monte Carlo integrals computing

$$\begin{eqnarray}&&{{ \mathcal L }}_{{\sigma }_{e},{\sigma }_{i}}={\int }_{r}P({r}_{\mathrm{sim}}| {\sigma }_{e},{\sigma }_{i})P({r}_{\mathrm{obs}}){dr}\end{eqnarray} \tag{ 4 }$$

for each system. The same cutoffs are then applied as for the singles (see Section 2.1) before computing the full population likelihood as a function of σ_e and σ_i. The multiplanet sample consists of 870 KOIs.

We find that the single planet systems are best fit with ${\sigma }_{e}={0.167}_{-0.008}^{+0.013}$ at the 95% confidence level (see Figure 1). The maximum and uncertainties are found by interpolating likelihoods from the grid of σ_e with a degree three Savitzky–Golay filter (Savitzky & Golay 1964) using the seven nearest neighbors to each point. The best-fit value is driven by the balance between the majority of the population of planets whose transit durations are consistent with low eccentricities, and roughly a dozen systems which disfavor e = 0 to moderate to high significance (2 − 40σ). We list systems that individually strongly suggest a high eccentricity in Table 1. These systems have longer transit durations than possible with a circular orbit (see Section 2) and are unlikely to have arisen by chance given the uncertainties. We also hand inspected each of these light curves and found nothing anomalous about these candidates.

Zoom In Zoom Out Reset image size

The existence of only a few systems that strongly suggest high eccentricity motivated us to consider a two-population eccentricity model, as also investigated in X16 and Van Eylen et al. (2019). Our model is a simple combination of two Rayleigh distributions, σ_e,low and σ_e,high, and the fraction of systems in the low-e distribution (f). To fit for the σ_e values for each population, we chose a pair of σ_e values from our grid search, sort the planets by their best-fit σ_e values, and find the division of the sorted list such that the likelihood is maximized. We then iterate over all possible σ_e,low and σ_e,high pairs, and compare the highest likelihood values for all pairs. This results in a triangular 2D likelihood surface with a ridge peak indicating two populations with σ_e,low ≲ 0.05 and σ_e,high ≳ 0.3 best fit the data (Figure 2). In the best-fit solution, 69% of systems preferred the low-eccentricity distribution. We note that at high σ_e the eccentricity distribution changes very little due to the cutoff at e = 1. This causes the weak dependency on σ_e,high for σ_e,high ≳ 0.3 seen in Figure 2 and explains our decision to truncate our search at σ_e = 0.7 (σ_e can be arbitrarily large since the distribution is truncated at e = 1, but the distributions become unphysical at high σ_e).

Zoom In Zoom Out Reset image size

We compare the likelihood of the single-population and two-population model to determine if the two-population model is warranted. The introduction of the two additional free parameters (fraction of systems in the low-e population and σ_e for the low-e population) increased the natural log-likelihood by 108, a Δ AIC (Akaike Information Criteria; Akaike 1974) of 212, strongly favoring a two-population model (Burnham & Anderson 2003). We note that this is not trivial due to the small number of systems with a confidently detected high eccentricity, because for a single high-eccentricity population the expected number of systems with anomalously long eccentricities compared to a circular distribution would be small due to the decreased transit probability (Burke 2008). Our results are similar to the two-population results in X16, and for a more thorough analysis of the two-population models and eccentricity distributions see Van Eylen et al. 2019.

3.1. Eccentricity Subpopulation Comparisons

The presence of a few individual planets out of the 439 considered that determine the lower bound of the eccentricity distribution via their long durations makes subdividing the population into bins of stellar or planetary properties to identify trends difficult. If the bins are even moderately small, it is probable that some bins contain none of these long-duration planets, resulting in a wildly oscillating σ_e distribution between ∼0 and ∼0.2 as a function of the property. Therefore we adopt an approach which makes use of the previous division into high- and low-eccentricity planets and is agnostic to the functional form of the eccentricity dependence. We divide the systems into high- and low-eccentricity groups for the best-fit σ_e,low and σ_e,high described above and then compare the properties of the 5% of systems which most strongly favor high eccentricity with all of the systems in the low-eccentricity group. This method is adopted for two reasons. First, if we would compare every system in both groups, any differences between the two populations may be strongly diluted by systems which only slightly favor the high-eccentricity group due, for instance, to circular systems at the high end of the impact parameter distribution whose short transit durations very modestly (<1σ) favor a high-eccentricity distribution. By taking only the top 5%, we restrict our sample to systems which are very likely to be eccentric, while still having enough systems to make a meaningful comparison. Second, we include the entire low-eccentricity population to get a good sense of the underlying property distribution for comparison. We consider the properties M_⋆, R_⋆, T_eff, [Fe/H], P, and R_p.

This leaves only seven planets in the high-eccentricity group, however each is highly likely to be eccentric. Further, this a sufficient sample for comparison to the ∼300 in the low-eccentricity distribution via the Kolmogorov–Smirnov two-sample tests, which take into account the number of samples in each distribution. However, we caution against over-interpretation of these results due to the small number of objects involved. For a threshold of statistical significance, we initially select a nominal p-value of 0.01 (2.6σ equivalent detection). We also consider the problem of multiple hypothesis testing (see, e.g., Miller 1981), which can be alleviated via a Bonferroni Correction (Bonferroni 1936; Miller 1981). Dividing the threshold p-value by the number of tests performed, we compute a new threshold of 0.01/6 = 0.002. Correlations between stellar variables considered (e.g., [Fe/H]–M_⋆; Santos et al. 2003; Johnson et al. 2010) would require a less harsh reduction in the p-values (e.g., Sidak 1968), so this correction is conservative. The resulting 0.002 threshold is as conservative as generally recommended for evidence of new effects (even when assuming conservative prior odds of 1:10–0.005, Benjamin et al. 2018; when considering multiple tests and reproducibility—0.003, Berger & Sellke 1987; from a Bayesian testing perspective—0.001–0.005, Johnson 2013; and as commonly used as first evidence for a particle in physics—0.003). It corresponds to a Bayes Factor of >50 compared to the null hypothesis (a significant detection; Johnson 2013), and is equivalent to a >3σ individual detection.

We find that the high-eccentricity planets do not significantly differ in P or stellar T_eff from low-eccentricity planets via a Kolmogorov–Smirnov two-sample test.⁷ A trend for high-eccentricity planets to prefer larger radii than low-eccentricity planets is seen at the p = 0.045 level (see Figure 3). We caution against over-interpretation of this result because of the low significance (p-value outside our threshold) and since we are only taking the systems with the most significant eccentricities. It is possible this trend is due to the lower S/N of transits of very small planets, which may have greater uncertainties in transit duration and thus are less likely to strongly suggest a non-zero eccentricity. However, we do not find a significant correlation between planet detection S/N and eccentricity (Figure 3).

We find that high-eccentricity planets are preferentially found around high metallicity ([Fe/H] > 0; Figure 4). This may in principle be influenced by the marginally larger than expected radii of the high-eccentricity planets as discussed above, because a stellar metallicity-planet radius correlation has been shown to exist (Buchhave et al. 2014; Petigura et al. 2018). However, this sample includes only one planet with R > 4R_⊕, and the same trend is observed (KS test p value = 0.0016) when it is removed and the sample restricted to <4R_⊕ where metallicity is thought to be only weakly dependent on planetary radius (Buchhave et al. 2014). Thus we rule out any radius dependence as the cause of the observed metallicity–[Fe/H] correlation. There is also a slight preference for eccentric planets to have high-mass host stars (Figure 4), but it is not formally significant. We note that the correlation coefficient, ρ, between stellar mass and metallicity in our sample is 0.3. Thus the slight mass preference may be due to the larger number of metal-rich stars at high mass, but high-mass stars could not be the cause of the observed eccentricity–[Fe/H] correlation. Additionally, the preference for higher M_⋆ decreases in significance when considering only planets <4R_⊕ (p-value = 0.064; outside our threshold). Further work is needed to validate both trends, and the mass trend in particular is not confidently observed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Additional Validation

For an intuitive understanding of the comparison between the observed durations and the eccentricity, we show the distribution of impact parameters implied from the measured Kepler durations assuming a population with e = 0 in Figure 5. We then compare it with simulations of the distributions recovered by drawing from the stellar and orbital property uncertainties, and injecting a uniform distribution of impact parameters. This allows us to see how the addition of the uncertainties changes the nearly uniform distribution of impact parameters to a double-peaked distribution disfavoring b ≈ 0. In order to take into account longer-than-expected transit durations given the R_⋆ value and a circular planet's velocity, we analytically extend the formula for the e = 0 impact parameter (b_a) as

$$\begin{eqnarray}\begin{array}{rcl}{b}_{1} & = & \left(1+\displaystyle \frac{{R}_{p}}{{R}_{\star }}\right)-\left(\displaystyle \frac{a}{{R}_{\star }}\sin \left(\displaystyle \frac{\pi d}{P}\right)\right)\\ {b}_{a} & = & \mathrm{sgn}({b}_{1})\sqrt{| {b}_{1}| }.\end{array}\end{eqnarray} \tag{ 5 }$$

Thus transit durations that are longer than theoretically possible for an (R_⋆, v_p) pair due to either eccentricity or measurement error are treated as negative impact parameters and the distribution is thus real for all possible input values. We note that although b is isomorphic to duration, the transformation is nonlinear and therefore the introduction of noise in the duration (and other parameter) measurements will not simply broaden an initial uniform distribution of impact parameters but also distort it. For b in [0, 1 + R_p/R_⋆), $| \tfrac{\partial b}{\partial d}|$ is strictly increasing and thus causes a bias toward high b. In other words, the density of durations as a function of b grows as b increases away from 0 and thus random noise preferentially biases bs to higher values. This effect increases with the level of noise. Additionally, noise may cause negative best-fit b values (which would create a pileup at b = 0 if negative b was not allowed). We emphasize this here as it has caused some confusion in the past (see Dong & Zhu 2013). Figure 5 shows that the single planet systems have both outliers at large negative values and a surplus of apparent high impact parameters—both indicators of eccentricity due to transit near apocenter and pericenter, respectively (Burke 2008). The routine described in Section 2.1 can be restated as a method of measuring how well the observed and theoretical b distributions match as a function of population eccentricity, while taking into account the uncertainties and observational bias which causes a decline in detection at very high b. The multiple planet case is more complicated due to the nonuniform underlying b distribution from mutual inclination constraints, so we do not show a similar simplified figure for such systems.

Zoom In Zoom Out Reset image size

The median radius uncertainties on the CKS-Gaia stars are ≈3%, and have been shown to be accurate compared to asteroseismology samples (Fulton & Petigura 2018). Nevertheless we conduct a test to determine how much unaccounted for systematic biases in the stellar properties would change our results. We refit the eccentricity distributions of the singles using a new set of stellar properties, where the stellar radii are all changed by 3% higher and lower with the reported uncertainties unchanged. The results do change the greatest likelihood σ_e and uncertainties, but still strongly rule out a σ_e near zero (see Figure 6). Since the transit duration is more weakly dependent on stellar mass, systematic biases are even less important there. Any biases for other stellar parameters would not affect the observed trends since the relative distributions of, for instance, metallicity between the two populations are measured rather than any absolute values being used. Thus we confirm that even if the stellar properties are systematically wrong by several σ in either direction, the qualitative results of our study will stand.

Zoom In Zoom Out Reset image size

We also consider systematics in the Kepler data. We compare transit durations reported in Data Release 24 (DR24; Coughlin et al. 2016) and Data Release 25 (DR25; Coughlin et al. 2017). The majority agree to within 1σ and an approximately normal distribution is generated by examining $({T}_{\mathrm{dur},\mathrm{DR}25}-{T}_{\mathrm{dur},\mathrm{DR}24}))/{\sigma }_{\mathrm{dur},\mathrm{DR}25}$ with σ = 1.05, excluding a few outliers.

The results of the transit duration simulations accounting for mutual inclination among multiple transiting planetary systems outlined in Section 2.2 indicate that the population of planets in multiple transiting systems have low eccentricities (σ_e = 0.0355 ± 0.012 at the 2σ level). We show a contour plot of the data with a cubic spline smoothing in Figure 7. The mutual inclination distribution is found to have ${\sigma }_{i}={2.45}_{-0.53}^{+0.65}$ degrees at the 2σ level. Similar to the single planet case, we consider individual systems whose transit durations rule out low eccentricities. Since the final multiplanet sample consists of 870 KOIs, we look for systems which prefer non-zero σ_e at greater than 3.7σ, as this leads to the expectation of fewer than 1 false positive. No systems reach this cutoff, suggesting a uniformly low-eccentricity distribution.

Zoom In Zoom Out Reset image size

We may still divide the KOIs into groups that prefer high- and low-eccentricity values as we did with the single transiting planets, even if systems do not individually strongly prefer one or the other. In addition to the parameters considered for the singles, we also consider the period ratio of the nearest neighboring planet (P_rat), the size of the largest planet in the system in Earth radii (R_max), and the number of planets in the system (N_pl). However, we find no statistically significant differences between the distributions of any of these parameters using the same methodology as for the singles at the p = 0.05 level. This is perhaps unsurprising as it is far more difficult to separate high- and low-eccentricity planets in this sample.

The most important insight this eccentricity population study provides is that high-e planets of 1.4–4R_⊕ size are preferentially found around metal-rich stars. Taken together with previous work on larger planets (Dawson & Murray-Clay 2013; Buchhave et al. 2018), these results suggest that eccentric planets of all sizes are preferentially found around super-solar metallicity stars. We also confirm previous results that the single transiting planet systems in Kepler are drawn from a broader range of eccentricities than the multiple transiting systems. The Rayleigh σ_e values we find are slightly lower than reported for the singles ($\bar{e}\approx 0.3$) in X16, but agree with being inconsistent with e ≈ 0. We note that we also performed a validation study using only the CKS spectral data, which resulted in a slightly higher e distribution for both the single and multiple transiting systems. This suggests that the improved accuracy provided by the Gaia calibrations on the stellar properties (Fulton & Petigura 2018) may explain the slight discrepancy between the results. On the other hand, our multiplanet σ_e agrees well with the results of X16 as well as those found by Van Eylen & Albrecht (2015), suggesting low eccentricities around the majority of multiple transiting Kepler systems.

A well-known correlation exists between the existence of giant planets and a star's metallicity (Fischer & Valenti 2005) and mass (Johnson et al. 2010). It is notable that these are the two characteristics which are most strongly correlated with eccentricity in the single transiting planet case. High stellar mass is associated with more massive disks (Andrews et al. 2013), which may encourage giant planet formation (Kennedy & Kenyon 2008), but also decreased lifetimes which may inhibit it (Ribas et al. 2015). However, high-metallicity environments promote both more robust planet formation and longer disk lifetimes increasing both the frequency of giant planets (Ercolano & Clarke 2010; Yasui et al. 2010) and perhaps Earth to sub-Saturn planets (Wang & Fischer 2015; Petigura et al. 2018, but see Buchhave & Latham 2015). These giant planets may interact with compact inner planet systems to increase eccentricity while decreasing multiplicity or apparent multiplicity due to greater mutual inclinations (Hansen 2017; Huang et al. 2017; Pu & Lai 2018), whereas typical multiplanet systems may not reach high eccentricities and mutual inclinations via self-excitation (Becker & Adams 2016). If this hypothesis is correct, a search for giant planets around stars hosting an apparent single eccentric planet should recover companions at a high rate.

More work and greater statistical certainty is required to disentangle the different influences on formation with the distributions of planetary eccentricities. We therefore look forward to future transit surveys such as TESS (Ricker et al. 2015) and PLATO (Rauer et al. 2014) and their follow-up campaigns, which will provide many more transiting planet durations to validate the planetary eccentricity–stellar metallicity relationship.

We thank the Kepler and Gaia teams for years of work making these precious data sets possible. We are grateful to the time assignment committees of the University of Hawaii, the University of California, the California Institute of Technology, and NASA for their generous allocations of Keck observing time that enabled this large project. This research has made use of NASA's Astrophysics Data System, the Exoplanet Orbit Database, and the Exoplanet Data Explorer at exoplanets.org. A.W.H. acknowledges NASA grant NNX12AJ23G. L.M.W. acknowledges support from the Beatrice Watson Parrent Fellowship. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has long had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.