The Possible Tidal Demise of Kepler's First Planetary System

The Kepler spacecraft observed Kepler-1658 for 12 quarters at 30 minute cadence and three quarters at 1 minute cadence. We downloaded the Kepler light curve using lightkurve (Lightkurve Collaboration et al. 2018) and modeled this data set using exoplanet (Foreman-Mackey et al. 2021). We fit the light curve quarter by quarter to obtain midquarter transit timings (defined to be the first transit after the midpoint of the quarter). Except for the midquarter time and the limb darkening coefficients, we used the posteriors from the fit in Chontos et al. (2019) as priors for our analysis to inform the model for each quarter of the average transit shape. We simultaneously fit the rotational variability of the star using a Gaussian process (GP), defining a SHOTerm in celerite2 (Foreman-Mackey et al. 2017; Foreman-Mackey 2018) with a free amplitude scale σ, a fixed timescale of 5.66 days (the known rotation period of the star), and a fixed quality factor of 1.

We ran four chains for all of the fits in this paper. Each chain was tuned for 2000 steps before 1000 posterior draws were taken. For each fit, we verified that the Gelman–Rubin (Gelman & Rubin 1992) statistic was $\hat{R}\ll 1.01$ for all sampled parameters. The transit times are given in Table 1, and the light curves are plotted in Figure 1.

Zoom In Zoom Out Reset image size

Table 1. Transit Times for Kepler-1658b

Data Set	Transit Time (BJDTDB)
Kepler LC Quarter 0	${2454959.7314}_{-0.0015}^{+0.0014}$
Kepler LC Quarter 1	${2454982.82835}_{-0.00061}^{+0.00061}$
Kepler SC Quarter 2	${2455048.26751}_{-0.00022}^{+0.00021}$
Kepler LC Quarter 3	${2455140.65189}_{-0.00042}^{+0.00040}$
Kepler LC Quarter 4	${2455233.03736}_{-0.00035}^{+0.00033}$
Kepler LC Quarter 5	${2455325.42133}_{-0.00036}^{+0.00035}$
Kepler SC Quarter 7	${2455510.19192}_{-0.00023}^{+0.00023}$
Kepler SC Quarter 8	${2455602.57708}_{-0.00027}^{+0.00027}$
Kepler LC Quarter 9	${2455691.11211}_{-0.00031}^{+0.00033}$
Kepler LC Quarter 11	${2455883.58121}_{-0.00033}^{+0.00032}$
Kepler LC Quarter 12	${2455975.96583}_{-0.00036}^{+0.00034}$
Kepler LC Quarter 13	${2456064.50087}_{-0.00036}^{+0.00036}$
Kepler LC Quarter 15	${2456256.97026}_{-0.00035}^{+0.00037}$
Kepler LC Quarter 16	${2456349.35438}_{-0.00039}^{+0.00038}$
Kepler LC Quarter 17	${2456410.94385}_{-0.00064}^{+0.00065}$
Palomar/WIRC Visit 1	${2459097.8002}_{-0.0015}^{+0.0015}$
TESS Sector 41	${2459436.5407}_{-0.0023}^{+0.0023}$
Palomar/WIRC Visit 2	${2459790.6819}_{-0.0013}^{+0.0015}$
TESS Sector 54	${2459786.8359}_{-0.0030}^{+0.0028}$
TESS Sector 55	${2459813.7791}_{-0.0027}^{+0.0029}$

Note. For the Kepler data sets, LC and SC refer to long cadence (30 minute exposures) and short cadence (1 minute exposures), respectively.

Download table as:  ASCIITypeset image

We observed transits of Kepler-1658b with Palomar/WIRC on UT 2020 September 5 and UT 2022 July 26. We collected data in the J band using a beam-shaping diffuser (Stefansson et al. 2017; Vissapragada et al. 2020; Greklek-McKeon et al. 2022). On the first night, we acquired 60 s exposures from 04:28 to 09:51 UT (airmass 1.0–2.2), and on the second night, we acquired 40 s exposures from 05:04 to 11:37 UT (airmass 1.1–1.7). Images from both nights were dark-corrected, flat-fielded, and background-subtracted per the methods in Vissapragada et al. (2020) and Greklek-McKeon et al. (2022). We then performed aperture photometry on Kepler-1658 along with 10 comparison stars on the first night and nine on the second night. We used uncontaminated annuli with 25 pixel inner radii and 50 pixel outer radii for local background subtraction. We tested aperture sizes from 5 to 25 pixels on both nights and chose the apertures that minimized the scatter in the final light curves, which was 18 pixels (45) for the first night and 12 pixels (30) for the second night.

We fit the Palomar/WIRC light curves using exoplanet, again using the results from Chontos et al. (2019) as priors for the fit for all parameters except the midtransit times and limb darkening coefficients. We modeled the systematics as a linear combination of comparison star light curves and the mean-subtracted times (Vissapragada et al. 2020; Greklek-McKeon et al. 2022). We also tried including combinations of additional decorrelation vectors into the systematics model, including the point-spread function (PSF) centroid offset, the background level of the image, and the airmass. We fit the light curve with all subsets of these three vectors and chose the systematics model that minimized the Bayesian information criterion (BIC). For the first night, the BIC was minimized when including only the background level, whereas on the second night the BIC was minimized when using only the airmass. The light curves are shown in Figure 2, and the transit times are reported in Table 1.

Zoom In Zoom Out Reset image size

TESS obtained photometry for Kepler-1658 (TOI-4480.01) in sectors 41, 54, and 55, all at a 2 minute cadence. We fit the TESS photometry sector by sector using the Chontos et al. (2019) priors for all values except the midsector transit time and the limb darkening coefficients. To handle the stellar variability in the TESS bandpass, we used the same GP parameters from the Kepler fit. The results are shown in Figure 2, and the transit times are reported in Table 1. To ensure our results were robust to choices in modeling methodology, three of us (F.D., M.G.-M., N.S.) reduced and fit the TESS data using independent pipelines and obtained consistent results.

The transit times are shown in Figure 3 with the ephemerides from Chontos et al. (2019) subtracted off. In agreement with their work, we found no compelling evidence for a changing orbital period in the Kepler data alone. However, the transits observed by Palomar/WIRC and TESS arrived early. In light of these new data, we sought to quantify the evidence for an evolving orbital period.

Zoom In Zoom Out Reset image size

We fit the transit times as a function of observing epoch t(N) using two models (Yee et al. 2020). The first model fit the data with a constant orbital period:

$$\begin{eqnarray}&&t(N)={t}_{0}+{NP},\end{eqnarray} \tag{ 1 }$$

where t₀ was the transit time at the reference epoch from Chontos et al. (2019). The second model fit the data including a constant period derivative dP/dN:

$$\begin{eqnarray}&&t(N)={t}_{0}+{NP}+\displaystyle \frac{1}{2}\displaystyle \frac{{dP}}{{dN}}{N}^{2}.\end{eqnarray} \tag{ 2 }$$

We used the nested sampling tool dynesty (Speagle 2020) to estimate the posteriors and Bayesian evidences ${ \mathcal Z }$ for these models. We performed each nested sampling run using single ellipsoid bounds, 1000 live points, and the random walk sampling method, terminating each run when the estimated log-evidence remaining was less than 0.01. The priors, posteriors, and evidences are given in Table 2. We then computed the Bayes factor $\mathrm{ln}B=\mathrm{ln}{{ \mathcal Z }}_{2}-\mathrm{ln}{{ \mathcal Z }}_{1}=17.5$ for this model comparison, representing decisive evidence (Trotta 2008) for the period derivative model over the constant period model. We conclude that the orbital period of the planet appears to be decreasing at a rate of ${131}_{-22}^{+20}$ ms yr⁻¹.

Table 2. Model Selection for the Timing Data

Model	Parameter	Unit	Prior	Posterior
No Decay	t0	BJDTDB	${ \mathcal U }({t}_{{\rm{c}}}-1,{t}_{{\rm{c}}}+1)$	${2455005.92478}_{-0.00014}^{+0.00013}$
$\mathrm{log}({ \mathcal Z })$ = 84.3	P	days	${ \mathcal U }({P}_{{\rm{c}}}-0.001,{P}_{{\rm{c}}}+0.001)$	${3.84936720}_{-0.00000066}^{+0.00000060}$
Decay	t0	BJDTDB	${ \mathcal U }({t}_{{\rm{c}}}-1,{t}_{{\rm{c}}}+1)$	${2455005.92415}_{-0.00016}^{+0.00017}$
$\mathrm{log}({ \mathcal Z })$ = 101.8	P	days	${ \mathcal U }({P}_{{\rm{c}}}-0.001,{P}_{{\rm{c}}}+0.001)$	${3.8493733}_{-0.0000011}^{+0.0000012}$
	${\mathrm{log}}_{10}(-{dP}/{dN})$	log10(days orbit−1)	${ \mathcal U }(-10,-6)$	$-{7.796}_{-0.079}^{+0.061}$
Precession	t0	BJDTDB	${ \mathcal U }({t}_{{\rm{c}}}-1,{t}_{{\rm{c}}}+1)$	${2455005.8771}_{-0.0061}^{+0.0086}$
$\mathrm{log}({ \mathcal Z })=98.9$	Ps	days	${ \mathcal U }({P}_{{\rm{c}}}-0.001,{P}_{{\rm{c}}}+0.001)$	${3.8493936}_{-0.0000072}^{+0.0000085}$
	$e\cos {\omega }_{0}$		${ \mathcal N }$(-0.00840, 0.00080)	$-{0.00836}_{-0.00076}^{+0.00076}$
	$e\sin {\omega }_{0}$		${ \mathcal N }$(0.062, 0.019)	${0.059}_{-0.022}^{+0.021}$
	${\mathrm{log}}_{10}(d\omega /{dN})$	log10(rad orbit−1)	${ \mathcal U }(-8,2)$	$-{3.159}_{-0.071}^{+0.085}$

Note. ${ \mathcal U }(a,b)$ indicates a uniform prior with lower bound a and upper bound b. ${ \mathcal N }(a,b)$ indicates a normal prior with mean a and standard deviation b. P_c and t_c are shorthand for the best-fit ephemerides from Chontos et al. (2019).

Download table as:  ASCIITypeset image

There are a number of effects that can cause apparent decreases in the planetary orbital period on decade-long timescales (e.g., Patra et al. 2017; Bouma et al. 2020; Yee et al. 2020; Maciejewski et al. 2021; Ivshina & Winn 2022). Only when these effects are ruled out can we definitively attribute our observations to orbital decay.

3.2.1. Line-of-sight Acceleration

3.2.2. Apsidal Precession

Next, we consider apsidal precession of an eccentric orbit, which can mimic decay due to tides if the full precession cycle is not captured (e.g., Patra et al. 2017; Yee et al. 2020). Kepler-1658b's orbit is slightly eccentric, as its secondary eclipse arrives about half an hour early relative to an orbital phase of 0.5 in the Kepler data (Chontos et al. 2019), so we must consider this scenario carefully. We fit the data with a constant precession rate using the model from Giménez & Bastero (1995), used previously by Patra et al. (2017) and Yee et al. (2020):

$$\begin{eqnarray}\begin{array}{rcl}t(N) & = & {t}_{0}+{{NP}}_{{\rm{s}}}-\displaystyle \frac{{{eP}}_{{\rm{a}}}}{\pi }\cos \omega (N),\\ \omega (N) & = & {\omega }_{0}+\displaystyle \frac{d\omega }{{dN}}N\\ {P}_{{\rm{s}}} & = & {P}_{{\rm{a}}}\left(1-2\pi \displaystyle \frac{d\omega }{{dN}}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

where $d\omega /{dN}$ is the precession rate and P_s and P_a are the sidereal and anomalistic periods, respectively. The priors, posteriors, and evidence for this model are given in Table 2. We find that precession is capable of fitting the data just as well as the period decay model (Figure 3), with relatively insignificant evidence $\mathrm{ln}B=2.9$ in favor of the decay model. The required precession rate $d\omega /{dN}\approx 7\times {10}^{-4}$ rad orbit⁻¹ is about 4° yr⁻¹.

However, this rate is severely problematic from a theoretical standpoint. If the precession is due to the planetary tidal bulge, the rate of precession constrains the planet's Love number k_2p (Ragozzine & Wolf 2009; Patra et al. 2017):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\omega }{{dN}} & = & 15\pi {k}_{2{\rm{p}}}\left(\displaystyle \frac{{M}_{\star }}{{M}_{{\rm{p}}}}\right){\left(\displaystyle \frac{{R}_{{\rm{p}}}}{a}\right)}^{5}\\ & & \approx ({8.9}_{-2.1}^{+2.6}\times {10}^{-7}\,\mathrm{rad}\,{\mathrm{orbit}}^{-1}){k}_{2{\rm{p}}}.\end{array}\end{eqnarray} \tag{ 4 }$$

The Love number is related to internal structure and ranges from 0 and 1.5 (e.g., Russell 1928; Sterne 1939), but matching the required precession rate requires an unphysical k_2p ∼ 10³ for the planet. The precession rate from the planetary tidal bulge is not physically permitted to exceed values of order ∼10⁻⁶ rad orbit⁻¹.

For hot Jupiters orbiting evolved stars, the stellar tidal bulge can be more important (Ragozzine & Wolf 2009):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\omega }{{dN}} & = & 15\pi {k}_{2\star }\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}\right){\left(\displaystyle \frac{{R}_{\star }}{a}\right)}^{5}\\ & & \approx ({1.70}_{-0.35}^{+0.45}\times {10}^{-4}\,\mathrm{rad}\,{\mathrm{orbit}}^{-1}){k}_{2\star }.\end{array}\end{eqnarray} \tag{ 5 }$$

Though precession from the stellar bulge is much faster, an unphysical k_2⋆ ≈ 4 is still required. Precession induced by the rotational flattening of the star and planet are even smaller contributions than the tidal bulge terms, and the precession rate from general relatively is 2 orders of magnitude too small as well (Ragozzine & Wolf 2009).

Finally, we consider the potential impact of an unseen outer body in the system. Secular perturbations from the outer body with mass M_o, semimajor axis a_o and eccentricity e_o would drive precession of the inner planet; in the limit of an apsidally antialigned outer body with a large eccentricity, the precession rate can be written approximately as (Mardling 2007; Batygin et al. 2009)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\omega }{{dN}} & = & \displaystyle \frac{15}{8}\pi \left(\displaystyle \frac{{M}_{{\rm{o}}}}{{M}_{\star }}\right){\left(\displaystyle \frac{a}{{a}_{{\rm{o}}}}\right)}^{4}\displaystyle \frac{{e}_{0}}{e{\left(1-{e}_{{\rm{o}}}^{2}\right)}^{5/2}}\\ & & =({5.67}_{-0.36}^{+0.38}\times {10}^{-6}\,\mathrm{rad}\,{\mathrm{orbit}}^{-1})\left(\displaystyle \frac{{M}_{{\rm{o}}}}{10{{\rm{M}}}_{{\rm{J}}}}\right){\left(\displaystyle \frac{{a}_{{\rm{o}}}}{1\,\mathrm{au}}\right)}^{-4}\\ & & \times \,\displaystyle \frac{{e}_{{\rm{o}}}}{{\left(1-{e}_{{\rm{o}}}^{2}\right)}^{5/2}}.\end{array}\end{eqnarray} \tag{ 6 }$$

To match the necessary precession rate, an outer body with 10M_J at 1 au would need an eccentricity of e ≈ 0.9. A body with this mass, separation, and eccentricity would have been readily observed in the RVs, with semiamplitude (540 m s⁻¹)$\sin i$, but the data do not indicate its presence. Unless the inclination of the outer body is very close to i = 0° (requiring some fine-tuning), this mechanism also seems unlikely.

We conclude that, while apsidal precession model can match the data, it requires an implausibly large precession rate and can thus be ruled out on physical grounds. Still, it would be helpful to strengthen this conclusion further with observations. For WASP-12b, the decisive evidence against apsidal precession came from the timing of secondary eclipses, which should arrive late for a precessing orbit (Patra et al. 2017; Yee et al. 2020). However, because the eclipse is quite shallow for Kepler-1658b (62 ppm in the Kepler bandpass, and ostensibly 200–300 ppm in the near-infrared), it is difficult to make the same measurement from the ground or with TESS. Now that we are a decade removed from the Kepler era, a secondary eclipse observation with the Hubble Space Telescope (HST) or JWST would be helpful for testing the precession model more definitively.

The best explanation for the early transit times observed by Palomar/WIRC and TESS is tidal decay of Kepler-1658b's orbit. Given our observed $\dot{P}$, we can calculate the rate at which orbital energy is being lost (e.g., Millholland & Laughlin 2018):

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}=\displaystyle \frac{{\left({{GM}}_{\star }\right)}^{3/2}{M}_{p}}{6\pi }{a}^{-5/2}\dot{P}=-{1.08}_{-0.21}^{+0.20}\times {10}^{31}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 7 }$$

To estimate the implied tidal quality factors, we assume the constant time lag model of Leconte et al. (2010) with zero stellar and planetary obliquity (though both could feasibly increase the dissipation if maintained over long timescales; Leconte et al. 2010; Millholland & Laughlin 2018; Su & Lai 2022). As the eccentricity is small, it is acceptable to truncate the model at second order in e (Equation (22) in Leconte et al. 2010)). Then, for inspiral dominated by dissipation in the star, the implied stellar tidal quality factor is

$$\begin{eqnarray}\begin{array}{rcl}{Q}_{\star }^{{\prime} } & = & \left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}\right){\left(\displaystyle \frac{{R}_{\star }}{a}\right)}^{5}\displaystyle \frac{27\pi }{\dot{P}}\left(\displaystyle \frac{{\omega }_{\star }}{n}-1\right)\\ & & ={2.50}_{-0.62}^{+0.85}\times {10}^{4}.\end{array}\end{eqnarray} \tag{ 8 }$$

Dynamical tides are required to make the star so dissipative. For WASP-12 (with ${Q}_{\star }^{{\prime} }\sim {10}^{5};$ Yee et al. 2020; Turner et al. 2021; Wong et al. 2022), it has been suggested that g modes deposit their energy efficiently via wave-breaking in the radiative core (Weinberg et al. 2017; Bailey & Goodman 2019; Barker 2020). This explanation is inadequate for Kepler-1658, for which the Barker (2020) models predict ${Q}_{\star }^{{\prime} }\sim {10}^{8}$ with this mechanism. Dissipation of inertial waves in the convective zone (Ogilvie & Lin 2007) are more effective for this rapidly rotating star. For their model closest to Kepler-1658, Barker (2020) find that inertial wave dissipation should result in ${Q}_{\star }^{{\prime} }\sim {10}^{4}$, in good agreement with our observations. If stellar dissipation drives the inspiral of Kepler-1658b, inertial waves are likely responsible.

If instead the inspiral is dominated by dissipation in the planet, then assuming the planet is tidally locked, the implied planetary tidal quality factor is

$$\begin{eqnarray}\begin{array}{rcl}{Q}_{{\rm{p}}}^{{\prime} } & = & -\displaystyle \frac{3}{2}\left(\displaystyle \frac{{M}_{\star }}{{M}_{{\rm{p}}}}\right){\left(\displaystyle \frac{{R}_{{\rm{p}}}}{a}\right)}^{5}\displaystyle \frac{171\pi {e}^{2}}{\dot{P}}\\ & & ={13.5}_{-3.5}^{+4.9}{\left(\displaystyle \frac{e}{0.06}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 9 }$$

A similar expression is given by Efroimsky & Makarov (2022), who found planetary eccentricity tides to be capable of explaining the orbital decay of WASP-12b. For Kepler-1658b, we require ${Q}_{{\rm{p}}}^{{\prime} }\sim 10$ if all the energy is dissipated in the planet. The quality factors for Jupiter and Saturn are much larger, ${Q}_{{\rm{p}}}^{{\prime} }\sim {10}^{5}$ (Goldreich & Soter 1966), but dynamical tides can lead to small effective quality factors via f mode diffusive growth and/or gravity wave dissipation in the radiative upper envelope of the planet (Lubow et al. 1997; Ogilvie & Lin 2004; Wu 2018). The former mechanism cannot work as it requires a pericenter distance of less than four tidal radii (about 5R_⊙ for this system), but it seems plausible that the outer part of the planet is radiative and could support the resonant excitation of g modes.

Given these estimates for the tidal quality factors, the majority of the energy is probably lost in the star and not the planet. But even a small amount of planetary dissipation would have observable consequences for the planetary energy budget as we describe in the next section. Finally, it is worth noting that the stellar spin period (5.66 days) and planetary orbital period (3.85 days) are close to a near-integer ratio, with ω_⋆/n ≈ 2/3. If this is not coincidental, the near commensurability may encode the system's history of tidal angular momentum exchange.

A secondary eclipse has been detected for Kepler-1658b in the optical, with depth 62 ± 4 ppm (Chontos et al. 2019). If this eclipse is due solely to reflected light, the geometric albedo would be A_g = 0.72 ± 0.09—by far the largest ever measured for a hot Jupiter. This is comparable to the geometric albedos of icy satellites in the solar system (Madden & Kaltenegger 2018) and is unexpected for hot Jupiters (e.g., Adams et al. 2022). We find it more likely that the planet is overluminous in the optical due to its own thermal emission, but the maximum dayside temperature of this planet is T_maximum = 2796 ± 73 K (Cowan & Agol 2011). At this blackbody temperature, the thermal contribution in the Kepler bandpass is negligible at ≲15 ppm; the dayside must be T_observed ≈ 3450 K to match the eclipse depth via thermal emission. The dayside temperature of the planet cannot be due to stellar irradiation alone.

The extra energy required to superheat the planet beyond its maximum dayside temperature could be provided by dissipation in the planet. This requires the luminosity from eccentricity tides to be a substantial fraction of the incident stellar irradiation (e.g., Bodenheimer et al. 2001; Jackson et al. 2008). To superheat the dayside of the planet, an additional luminosity is required:

$$\begin{eqnarray}&&{\rm{\Delta }}L={\sigma }_{\mathrm{SB}}\pi {R}_{{\rm{p}}}^{2}({T}_{\mathrm{observed}}^{4}-{T}_{\mathrm{maximum}}^{4})\approx 8\times {10}^{29}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 10 }$$

This is 1 order of magnitude smaller that the total rate at which orbital energy is being dissipated (Equation (7)); it is plausible that ≲10% of the energy from the shrinking orbit is being dissipated in the planet. Planetary dissipation is especially interesting considering that many gas giants orbiting evolved stars appear to be (re)inflated, which requires the deposition of additional energy beyond the incident stellar radiation (Grunblatt et al. 2016, 2017; Lopez & Fortney 2016). If a similar mechanism operates for these systems, tidal heating could be a natural explanation for the inflation of planets orbiting evolved stars.