The Discovery of the Long-Period, Eccentric Planet Kepler-88 d and System Characterization with Radial Velocities and Photodynamical Analysis

We obtained RVs of Kepler-88 on the HIRES spectrograph (Vogt et al. 1994) at the W. M. Keck Observatory between the years 2013 and 2020. We used the standard HIRES setup of the California Planet Search (see Howard et al. 2010 for details). Spectra were obtained using HIRES in the red-collimation mode with a warm molecular iodine gas cell in the light path for wavelength calibration. We used the C2 decker (, R = 60,000) to enable sky-subtraction for this relatively faint (V = 13.8) target. Since the target was faint, we only observed in good conditions (seeing <15, clear to thin clouds). For each spectrum, we achieved a signal-to-noise ratio of at least to ensure that our Doppler pipeline would deliver RVs with errors of (Howard & Fulton 2016).

We observed an iodine-free template spectrum bracketed by observations of rapidly rotating B-type stars to enable a deconvolution of the stellar spectrum from the spectrograph point-spread function (PSF). We then forward modeled our RV spectra with the deconvoled template stellar spectrum plus a night-specific model of our PSF convolved with an atlas iodine spectrum. We also used the blue HIRES chip to extract a Mt. Wilson value for each HIRES observation. Our Keck-HIRES RVs and values, plus the SOPHIE RVs from the literature, are presented in Table 1.

Table 1.  Kepler-88 RVs

Time	RV		SHK	Inst
(BJD — 2454900)	(m s−1)	(m s−1)
1575.047908	−35.7	2.8	0.142	HIRES
1575.40947	42.0	10.0		SOPHIE
1575.992695	−36.7	2.4	0.155	HIRES
1577.892731	−52.9	2.5	0.147	HIRES
1579.027079	−81.7	2.6	0.122	HIRES
⋯	⋯	⋯	⋯	⋯

Note. Times are in BJD — 2454900.0. SOPHIE RVs are from Barros et al. (2014); HIRES RVs are from this work. The SOPHIE RVs have had 20465.0 added, with respect to the values published in Barros et al. (2014), for easier zero-point calibration. The RV uncertainties do not include RV jitter.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The stellar properties of Kepler-88 were determined based on our high signal-to-noise template spectrum in combination with the Gaia parallax and Two Micron All Sky Survey (2MASS) photometry (Fulton & Petigura 2018). The stellar temperature is 5466 ± 60 K and with [Fe/H] = , the star is slightly metal-rich. The star has a similar mass but slightly smaller radius than the Sun ( = , = ).

Because the transit of planet b combined with dynamical information about the planet constrains the density of the star, we used a photodynmical fit to update the stellar characterization (e.g., Vanderburg et al. 2017, see Section 5.2). We used the best-fit values and uncertainties for the stellar mass and radius from Fulton & Petigura (2018) as priors in our photodynamical fit. After our photodynamical fit, the best-fit stellar mass and radius are = and = . The precision of the stellar radius determination was improved through the photodynamical fit, suggesting that the transits provide information about the stellar density and hence the stellar radius.¹¹

Of the stellar parameters reported here, only the stellar mass is dependent on isochrone fitting (see Fulton & Petigura 2018 for details). We caution that the formal error in the mass reported here does not account for systematic differences between the stellar isochrones formulated by different research groups, and so the reported error in the stellar mass (and hence density) might be underestimated.

In systems that are not very close to resonance and/or have large TTV uncertainties compared to the timing precision, only the low-frequency TTV super-period is detected. This low-frequency signal contains information about the mass ratio of the planets, but the absolute masses are degenerate with the eccentricities of the planets (Lithwick et al. 2012). However, in systems with high signal-to-noise TTV measurements like Kepler-88, it is possible to detect a higher-frequency signal: the synodic chopping signal. This signal abruptly changes direction after conjunctions between the transiting and perturbing planet (Agol & Fabrycky 2018). The chopping signal is therefore expected to occur at the synodic period, or

In Kepler-88, the expected chopping period is P_chop = 21.5 days.¹² We identified the chopping signal by fitting the Holczer et al. (2016) TTVs with a high-order polynomial¹³ (degree 18) and subtracted this polynomial fit from the TTVs (Figure 5, middle panel, blue points). The high-frequency variations in the TTVs show a characteristic chopping signal. Note that the strength of the chopping comes and goes at different phases of the TTV super-period. This episodic strength of the chopping signal is well reproduced by our best-fit model (orange crosses).

We used a Lomb–Scargle periodogram to identify the super-period of the TTVs at 611 days (Figure 6, top). After we removed the high-degree polynomial, the Lomb–Scargle periodogram of the TTV residuals had its highest peak at P = 20.9 or 23 days (which is the mirror reflection of 20.9 days about the Nyquist frequency), and a peak at the expected chopping signal at 21.5 days (bottom). Note that peak periods are reflected about the Nyquist frequency, day⁻¹.

Zoom In Zoom Out Reset image size

To better understand the origin of the peak at 20.9 or 23 days, we examined the behavior of the TTV chopping signal on longer timescales. We simulated TTVs for 50 yr using the best-fit parameters and the N-body integrator rebound (Rein & Tamayo 2016; Figure 7, top panel). The long-term TTVs have both the a super-period at 611 days, and a super–super period at 20.5 yr (7500 days). We used a fast Fourier transform to construct a low-pass filter and subtracted the filter from the model, thus obtaining the chopping signal (bottom panel). In Figure 8, we show the Lomb–Scargle periodogram of the chopping signal based on different numbers of consecutive transits. A periodogram of the 221 transits from the middle of the simulated chopping signal yields peaks at the chopping frequency (21.56 days) and 20.9/23 days. Including more transits increases the power at these distinct periods, and reveals splitting of the 20.9/23 days period. The peak at 20.91 days is at , i.e., the 2:3 resonance.¹⁴ The splitting appears consistent with the super–super period, as the peaks occur at 20.91 and days.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To investigate the origin of the super–super period, we plotted the evolution of and in a simulation with just the inner two planets (Figure 9), and confirmed that the super–super period (20.5 yr) is the timescale of precession of planet b due to planet c.¹⁵ There are approximately 11 retrograde epicycles during the 20.5 yr precession; each of these epicycles corresponds to a TTV super-period of ∼611 days.

Zoom In Zoom Out Reset image size

We downloaded the photometry of Kepler-88 (KOI-142, KIC 5446285) obtained during the Kepler prime mission from the MAST archive.¹⁶ Where available (quarters 4–17), we used short-cadence data; we used long-cadence data elsewhere. We detrended the photometry in the manner of Mills & Fabrycky (2017b). First, we segmented the lightcurve into chunks of approximately one day, masking any transits within each chunk. We then fit the photometry in each chunk with a cubic polynomial to model the continuum, including both systematic effects and stellar rotation. We divided the observed flux by our continuum model to obtain normalized photometry. We multiplied all the uncertainties by a scale factor such that the out-of-transit reduced χ² is 1.0, for both long cadence and short cadence independently.

Since no stellar companions are known, we fixed the transit dilution at zero for all the planets, and we fixed planet-to-star radius ratios of the two non-transiting planets, the longitude of ascending node for planet b of Ω_b = 0.0 (since this is an arbitrary angle on the sky plane), and the inclination and longitude of ascending node for planet d at i_d = 89°, Ω_d = 0.0 (the RVs only give information for this planet, and the TTVs do not help constrain its inclination). All other parameters were varied.

We arrived at our best estimate of the dynamical parameters in the following manner. First, we used the parameters published in N13 in conjunction with the software package TTVFast (Deck et al. 2014) to minimize our fit to the long-cadence determined TTVs reported in Holczer et al. (2016). When our fit to the long-cadence TTVs was optimized, we used our best fit as input orbital elements for Phodymm. We then ran 40 differential evolution MCMC (DE-MCMC) chains 10⁶ steps each to obtain improved values and formal uncertainties for each of the Phodymm variables. The chains were well mixed, based on both a visual inspection of the chain for each parameter and a maximum Gelman-Rubin statistic of 1.05 (Gelman & Rubin 1992). Our best-fit parameters and uncertainties for the photodynamical N-body model are in Table 5.

Table 5.  Phodymm MCMC Posteriors

Parameter	Units	Median ±1σ
Periodba	days	10.91647 ± 0.00014
T0,b	BJD−BJD0	55.08069 ± 0.00061
		−0.23578 ± 0.00031
		0.0044 ± 0.0027
ib	°	90.97 ± 0.12
	Jup	0.0300 ± 0.0036
Rb/Rs		0.03515 ± 0.00018
Periodc	days	22.26492 ± 0.00067
T0,c	BJD−BJD0	61.353 ± 0.025
		0.2392 ± 0.00095
		−0.0044 ± 0.0033
ic	°	93.15 ± 0.68
Ωc	°	-0.43 ± 0.19
	Jup	0.674 ± 0.016
Periodd	days	1403 ± 14
T0,d	BJD−BJD0	1335 ± 19
		0.63 ± 0.03
		0.08 ± 0.05
Msin	Jup	3.05 ± 0.16
Ms	solar	0.990 ± 0.023
Rs	solar	0.897 ± 0.016
c1		0.394 ± 0.062
c2		0.292 ± 0.096
Mb	Earth	9.5 ± 1.1
Mc	Earth	214.1 ± 5.2
Msin id	Earth	965 ± 44
Rb	Earth	3.438 ± 0.075
	g cm−3	1.29 ± 0.16
eb		0.05561 ± 0.00013
ec		0.05724 ± 0.00045
ed		0.41 ± 0.03
Ibc	°	2.23 ± 0.62

Notes. The MCMC parameters are above the line; derived parameters are below the line. All parameters are computed at epoch . Ω_b is an arbitrary reference angle and was fixed at 0.0. The photodynamical solution was not sensitive to the inclination of planet d, which we fixed at i_d = 89°.

^aNot the same as the linear ephemeris, which is days.

Download table as:  ASCIITypeset image

Figure 10 shows the photometry phase-folded to the individual transit times of Kepler-88 b. The sharp ingress and egress indicate that the individual transit times have been well determined. Furthermore, the distribution of the photometric residuals during transit are identical to the distribution of the photometric residuals outside of transit, and both are Gaussian, with a standard deviation of 550 parts per million per exposure. The residuals do not have strong correlated features, suggesting that the individual transit times, depths, shapes, and durations have been well modeled.¹⁷ No transits of planets c or d were detected: even a grazing transit of a 1 planet with an impact parameter of would produce a deeper transit than that of planet b, which would be easily identifiable in the photometry of Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One outcome of the photodynamical modeling is a model-dependent determination of the individual transit times. The transit midpoint times T_i, impact parameters b_i, and planet velocities during transit v_i and their uncertainties, are given in Table 6 for Kepler-88 b from the date of the first Kepler transit through 2022 November.

In both our RV-only and our photodynamical analysis, we confirm the existence of a giant planet at 22.26 days with a mass of 200 . A Lomb–Scargle periodogram of our RVs produces a very strong peak at 22.26 days, with no significant peaks at aliases or harmonics of this period, indicating that 22.26 days is in fact the period of the perturbing giant planet (Figure 12).

Zoom In Zoom Out Reset image size

In the Keck-HIRES RV data, we identify a third planet at with = . When we compute the Lomb–Scargle periodogram to the residual RVs of a two-planet fit (where the orbits are N-bodies), there is a significant peak at P = 1413 days, and there are no other peaks with comparable power (see Figure 12). We find evidence for the third planet in the significantly improved χ² statistic to the RVs, which we summarize in Table 4. Without the third planet, our best N-body fit to the TTVs + RVs has (see Figure 3). Including the third planet near P = 1413 days results in (see Figure 4). The inclusion of the third planet substantially improves the fit to the RVs while simultaneously reducing the HIRES RV jitter by a factor of ∼3. However, the goodness of fit to the TTVs does not change significantly between the two-planet and three-planet models, indicating that the TTVs provide essentially no evidence for the existence of planet d. This is unsurprising, since at , the outer planet causes negligible changes to the orbits of the inner planets on the timescale of the Kepler baseline. Therefore, the outer planet was only detected in RVs. The RV-led discovery highlights the importance of multi-method follow-up of the most architecturally interesting Kepler planetary systems.

The photodynamical fit provides some advantages over the RVs alone. Although the RV alone identifies the period of planet c as P_c =  days, the photodynamical fit tightens the uncertainty by almost an order of magnitude, finding . Note that the precision on our determination of the orbit of planet c is a factor of three better than that of N13, which found days; the improvement in the precision of the orbital period must come from the additional Kepler photometry in our analysis, since the RVs alone did not determine the orbit of planet c as precisely as the TTV-based N13 work. Also, the RVs are only able to provide a mass upper limit for planet b, but the planet’s period, radius, and mass are determined with high confidence in the photodynamical analysis: (at epoch BJD = 2454954.62702), , and . The superior performance of the photodynamical model for the inner planets illustrates the complementary nature of transit photometry and RVs: together, these techniques reveal more about the 3D architecture of a planetary system than each of these techniques does alone. The superior mass determination of planets b and c in the photodynamical model can be traced to the chopping signal in the TTVs.

From the photodynamical analysis, we determined that the two inner planets, b and c, are apsidally anti-aligned. Our result agrees with N13. Apsidal anti-alignment is a predicted outcome of convergent or divergent Type-I migration in a viscous disk (Nelson 2018, and references therein). The combination of the near-resonant configuration for the inner planets and their apsidal anti-alignment could suggest a history of migration and resonant trapping, although the current anti-alignment is likely a short-lived coincidence (see Figure 9).

Planet d’s longitude of periastron passage is nearly aligned with that of planet b, and anti-aligned with that of planet c. The apsidal alignment of planet d is far less meaningful, since it is dynamically decoupled from the inner two planets. However, the eccentricity of planet d is large () compared to the inner planets (). The high eccentricity of planet d could be explained by planet–planet scattering or Kozai oscillations. By contrast, the modest eccentricities of planets b and c could possibly be explained by the equilibrium of disk and/or tidal circularization and N-body eccentricity pumping.

Because both the inclinations and relative longitude of ascending node are constrained, we can compute the mutual inclination between planets b and c:

We find that the mutual inclination of planets b and c is tightly constrained: . To test the extent to which the inclinations of the planets are constrained photodynamically, we initialized a DE-MCMC experiment with 40 walkers with our best fit, but flipped the inclination of planet c about the stellar meridian (i = 90°), allowing the same parameters to vary as in our control trial. In general, it is difficult for a DE-MCMC exploration to find this parameter space, because inclinations near 90° for planet c would produce deep transits, which are not observed. However, our inclination-flipped experiment performed substantially worse than our best fit, with ΔBIC = 33 (see Table 4). Therefore, our findings strongly disfavor the model in which planets b and c are on opposite sides of the star; rather, they seem to be on the same side of the star.

Constraints on the mutual inclinations of planets likely come from TDVs. The best-fit solution to our photodynamical model includes substantial TDVs for planet b (see Figure 13), though these are dominated by in-plane eccentricity precession (i.e., Nesvorný et al. 2013). We determined the TDVs based on the velocities and impact parameters in Table 6:

where is the stellar radius in units of astronomical units, b_i is the ith transit impact parameter in units of astronomical units, and is the velocity of the planet during the ith transit in units of astronomical units day⁻¹.¹⁸ In addition to long-term variation, the TDVs have a chopping signal.

Zoom In Zoom Out Reset image size

We explored the extent to which we could constrain the inclination of planet d from the TTVs. By keeping constant but varying i_d and M_d, we found that the best-fit solution to the TTVs did not significantly degrade. We tried this experiment in forward modeling the TTVs with TTVFast and also with Phodymm. In both cases, a wide range of mutual inclinations between planet d and the inner planetary system are supported by the data. For example, initializing the inclination of planet d at i_d = 30° from the sky plane, i.e., about 60° from the inner planets, and thus requiring its mass to be , only increased (worsened) the BIC of our photodynamical model by six, which suggests only modestly better performance of a coplanar model. We consider the long-term orbital stability of such solutions in the next subsection. Simulations in similar planetary systems have found that a long-period giant planet is likely to be coplanar with the inner planets, as this configuration is usually stable for longer periods of time than highly mutually inclined geometries (Becker & Adams 2017). In an analysis of the Kepler-88 system, Denham et al. (2019) found that a planet at semimajor axis 2.4 au with would be stable so long as its mass was . At , Kepler-88 d could take on a wide range of inclinations without exceeding 20 and thereby disrupting the inner system.

The architecture of a long-period giant planet accompanying two closer-in planets reminded us of the Kepler-56 system. Kepler-56 is a red giant star hosting two coplanar, transiting planets whose orbits are misaligned with respect to the stellar rotation axis, which is determined from the asteroseismic modes of the star (Huber et al. 2013). RV monitoring of the system revealed a long-period, non-transiting giant planet with moderate eccentricity (Otor et al. 2016). Follow-up theoretical work (Gratia & Fabrycky 2017) suggested that the outer planet can become eccentric due to planet–planet scattering. In such cases, additional outer planet(s) are likely ejected, leaving the surviving outer planet on an eccentric and inclined orbit. The perturbations ripple to the inner system, not necessarily disrupting it, but possibly causing precession of the orbital plane that periodically misaligns those two planets from the host star’s equatorial plane.¹⁹

Here we run a long-term N-body simulation for Kepler-88 to observe the evolution. For the simulation we entered a fit of the data into the Mercury package (Chambers 1999) and used the Burlisch–Stoer integrator for 0.1 Myr to record the secular-timescale effects. We have assumed the outer planet is inclined 30° from our line of sight and ∼27° from the inner planets, but that is not constrained from the data. In Figure 14 we show that no substantial eccentricity is transferred from the outer giant to the inner planets on these timescales, but a long-term precession effect can excite the inner planets to a large inclination from its original plane, and hence relative to the star. The orbital planes of the two inner planets remain closely aligned to each other.

Zoom In Zoom Out Reset image size

We follow Boué & Fabrycky (2014) to evaluate the secular timescales. The frequency at which the inner planets would precess due to the outer planet, if the outer planet were not to back-react, is ν₃ = 2.7 × 10⁻¹² rad s⁻¹ (74 kyr period). The frequency that the outer planet would precess due to the inner planets, if they were not to back-react, is ν₄ = 1.5 × 10⁻¹³ rad s⁻¹ (1.33 Myr period). Together, the frequency of precession should be rad s⁻¹ (81 kyr period), where I_cd is the mutual inclination between planets c and d. The frequencies related to stellar spin precession are several orders of magnitude smaller. Therefore the inner planets are effectively precessing due to the outer planet, without much back-reaction (as is evident also in Figure 14, with a precessional frequency of 77 kyr, near the analytic value), a conclusion that requires only that the outer planet’s angular momentum dominates the system, which is true for any possible inclination of the outer planet. Also, the star’s precession cannot keep up with that relatively quick motion. So we expect that a spin orientation measurement of the star, with respect to the transiting planet’s orbit, would help diagnose whether the outer planet has a significant inclination with respect to the inner planets. For instance, if the spin–orbit misalignment is 20° currently, it is likely cycling between 0° and at least 20°, meaning the inclination of planet d with respect to the inner planets is at least 10°. There may, of course, be more to the dynamics than this simple picture, such as a non-trivial scattering history (Gratia & Fabrycky 2017).

For completeness, we ran two other simulations for 0.1 Myr, one in which planet d is nearly coplanar to the inner system, and one corresponding to the a tilt of 60° and a doubling of the planetary mass as described above. Both simulations remained stable, with quasi-periodic oscillations similar to those in Figure 14. We conclude that planet d does not threaten the stability of the system for a wide range of inclinations.

Since both planets c and d are gas giants, they must have formed early in the disk lifetime, when gas was abundant. The presence of multiple giant planets in this system is unsurprising since [Fe/H] = , and the occurrence of giant planets increases with stellar metallicity (Fischer & Valenti 2005). Perhaps additional giant planets were present earlier, or are still present. Planets c and d likely underwent viscous (Type I) migration in the protoplanetary disk. As the gas disk dissipated, planet–planet scattering would likely have increased, and low- and high-eccentricity migration likely became important at this time. The high eccentricity of planet d probably arose due to a significant exchange of angular momentum with another gas giant planet.

The formation of planet b could have been contemporaneous with the giant planets if the planet were somehow gas-starved, resulting in only a low-mass volatile envelope. Or perhaps planet b formed when gas was less abundant and was caught in mean-motion resonance with planet c during an epoch of inward migration of planet c.

Giant planets are present around a large number of the Kepler systems that host small, transiting planets (Marcy et al. 2014; Mills et al. 2019), and perhaps at greater frequency than giant planets occur around field stars (Zhu & Wu 2018; Bryan et al. 2019). Kepler-88 joins their ranks. Furthermore, Kepler-88 has two giant planets. Other systems with multiple giant planets in addition to small transiting planets include WASP-47 (Neveu-VanMalle et al. 2016) and 55 Cnc (McArthur et al. 2004). The stellar and planetary properties of these systems are summarized in Table 7.

Table 7.  Systems with Multiple Giant Planets and Small Transiting Planets

Parameter	Units	Kepler-88	55 Cnc	WASP-47
Stellar Parameters
	K	5466 ± 60	5196 ± 24	5552 ± 75
		0.99 ± 0.024	0.905 ± 0.015	1.040 ± 0.031
		0.897 ± 0.016	0.943 ± 0.010	1.137 ± 0.013
[Fe/H]	dex	0.27 ± 0.06	0.31 ± 0.04	0.38 ± 0.05
Age	Gyr	1.9 ± 1.6	10.2 ± 2.5	6.5 ± 2
Innermost Transiting Planet
Letter		b	e	e
Period	days	10.91649 ± 0.00014	0.736539 ± 0.000007	0.789592 ± 0.000012
		3.44 ± 0.08	1.91 ± 0.08	1.810 ± 0.027
		9.5 ± 1.1	8.08 ± 0.31	6.83 ± 0.66
		1.29 ± 0.16	6.4 ± 0.8	6.35 ± 0.64
Ecc.		0.05561 ± 0.00013	0.040 ± 0.027	0.03 ± 0.02
Innermost Giant Planet
Letter		c	b	b
Period	days	22.2649 ± 0.0007	14.65152 ± 0.00015	4.1591289 ± 0.0000042
		214.0 ± 5.23	264.0 ± 1.0	363.1 ± 7.3
Ecc.			0.0034 ± 0.0032	<0.002
Outermost Known Giant Planet
Letter		d	d	c
Period	days	1403 ± 14	4825 ± 39	588.5 ± 2.4
		965 ± 44	1232 ± 22	398.2 ± 9.3
Ecc.			0.019 ± 0.013	0.296 ± 0.017

Note. For Kepler-88: planetary parameters and are from this work, and other stellar parameters are from Fulton & Petigura (2018). For 55 Cnc: stellar parameters are from von Braun et al. (2011), and planetary parameters are from Demory et al. (2016, planet e) and Baluev (2015, other planets). For WASP-47: stellar and planetary parameters are from Vanderburg et al. (2017).

Download table as:  ASCIITypeset image

Like Kepler-88, WASP-47 has a nearly circular hot Jupiter and a slightly eccentric longer-period giant planet (Sinukoff et al. 2017; Vanderburg et al. 2017; Weiss et al. 2017). Similarly, 55 Cnc has a close-in, nearly circular warm Jupiter at P = 14.7 days, and three other known giant planets at 44.4, 261, and 4800 days (Marcy et al. 2002; McArthur et al. 2004; Naef et al. 2004; Fischer et al. 2008; Wright et al. 2009; Dawson & Fabrycky 2010; Endl et al. 2012; Nelson et al. 2014; Baluev 2015). In general, systems with hot Jupiters do not tend to have companions within ∼1 au (Steffen et al. 2012b), although many such systems have companions from 5 to 20 au (Bryan et al. 2016). Perhaps systems with hot/warm Jupiters in proximity to small exoplanets and/or with metal-rich stars are an exception to these patterns in the broader population. Long-baseline RV studies of the Kepler and Transiting Exoplanet Survey Satellite (TESS) systems with hot and/or warm Jupiters in addition to small planets will reveal whether these systems also have distant giant planets.

Kepler-88 differs from 55 Cnc and WASP-47 in that the innermost known planet has an orbital period of 11 days, rather than <1 day. Both 55 Cnc e and WASP-47 e are examples of ultra-short period (USP) planets, which are defined as having P < 1 day (Dawson & Fabrycky 2010; Becker et al. 2015). USPs are generally small (; Sanchis-Ojeda et al. 2014); this is likely because their high equilibrium temperatures do not support a volatile envelope of hydrogen and helium, although they could support thin envelopes of heavier mean-molecular weight species like water or silicates (Lopez 2016). Kepler-88 b is nearly the size of Neptune, and, at a density of 1.1 , has abundant hydrogen and helium. However, its mass of is very similar to the masses of 55 Cnc e and WASP-47 e. We speculate that perhaps the nearby giant planet may eventually perturb Kepler-88 b into an orbit with P < 1 day, where photoevaporation would remove the H/He envelope. One mechanism that could accomplish this rearrangement is that the outer gas giant could perturb the 2:1 resonance to ever-widening libration amplitude, whereupon a scattering interaction between Kepler-88 b and c sends Kepler-88 b on an eccentric orbit which is tidally circularized at a much smaller orbital period. Kepler-88 is significantly younger than the other systems (age = 1.8 ± 1.6 Gyr, see Table 7). Thus, Kepler-88 might represent an early prototype of the 55 Cnc and WASP-47 systems, in particular their inclusion of a hot super-Earth.

The measurement of planetary masses with RV and transit timing may always be confused by the presence of additional planets perturbing the star (RV) or the transiting planet(s) (TTVs). In addition, both techniques can suffer from stellar systematics such as stellar jitter (RV) or stellar photometric inhomogeneities (TTVs), as well as instrumental systematics. This motivates a comparison of these two techniques when both are available. Unfortunately the number of systems for which this is possible is very small due to the small probability of transit of RV-detected systems, and the poor RV precision of most TTV systems found with Kepler. Mills & Mazeh (2017) found only nine planets which had both RV and TTV mass measurements, of which eight agree to better than 2σ, while one (Kepler-89d) may be influenced by the presence of additional undiagnosed planets (Mayo et al. 2017). Note that each technique has a slightly different dependence upon the stellar mass, so that precise stellar parameters are needed to carry out a comparison.

The Kepler-88 system adds another planet for which both RV and TTV measurements are available: Kepler-88c. In our RV analysis, we found a minimum mass of this planet of , while in an analysis of the TTVs only using a three-planet model we found a best-fit mass of . The photodynamical analysis indicates an inclination for planet c of 9314, so these two determinations agree to within 1σ. Thus, Kepler-88c adds to the number of cases in which consistent RV and TTV masses are obtained, building confidence in both techniques.

The presence of two planets in near-resonant orbits can sometimes be confused for a single planet with moderate eccentricity (Wittenmyer et al. 2019). Our RV baseline is too short, and our sampling of the periastron too sparse, to fully explore the possibility of a fourth planet in the system (see Figure 12). Our most recent RV just after the 2020 periastron passage of planet d shows no sign of a fourth planet.

The TESS spacecraft observed Kepler-88 during its northern hemisphere campaign this summer. The photometric precision of TESS should be adequate to detect Kepler-88 b, if the planet is still transiting (Christ et al. 2019).