Giant Outer Transiting Exoplanet Mass (GOT 'EM) Survey. II. Discovery of a Failed Hot Jupiter on a 2.7 Yr, Highly Eccentric Orbit^*

The Kepler spacecraft observed Kepler-1704 at 30 minutes cadence in all 18 quarters of its primary mission. We accessed the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP; Jenkins et al. 2010; Smith et al. 2012; Stumpe et al. 2012) through the Mikulski Archive for Space Telescopes (MAST), stitching together the light curves from individual quarters into one time series with a common baseline flux using lightkurve (Lightkurve Collaboration et al. 2018). We further cleaned the photometry by removing all data points flagged for "bad quality" and dividing out the median background flux to produce a normalized light curve. We then measured a preliminary time of conjunction, duration, and period for the transiting planet using a box least-squares (BLS) transit search (Kovács et al. 2002), identifying only two transit events in quarters 2 and 13. The time separating these two transits was 989 days, the suspected orbital period of Kepler-1704 b. However, the data gap between quarters 7 and 8 occurred precisely in between these transits, introducing an ∼494 day orbital period alias.

We used the BLS results to mask out both transits and detrend any variability in the light curve without risk of obscuring the signal. Interpolating over the masked transit events, we fit a smoothed curve to systematics in the photometry using a Savitzky–Golay filter (Virtanen et al. 2020) and then subtracted out this additional structure to produce our final data product. Before unmasking the transit events, we clipped any remaining individual outliers with residuals that were greater than 3σ discrepant to the smoothed fit.

We present the binned, detrended Kepler transits of Kepler-1704 b in Figure 1. Under the assumption of a circular edge-on orbit, the mean transit duration of Kepler-1704 b and stellar properties reported by the NASA Exoplanet Archive suggest an orbital period of approximately 11 days. This scenario is thoroughly ruled out by the extensive Kepler data set. Instead, we explored the possibility that orbital eccentricity affected the duration of the transit.

Zoom In Zoom Out Reset image size

2.1.1. Photoeccentric Transit Modeling

The observed transits of Kepler-1704 b have a duration of ∼6 hr, which is nearly five times shorter than would be expected for a Jovian-size planet with a 989 day orbital period. The two plausible sources of this discrepancy are a high impact parameter (b) or high eccentricity (e), but a preliminary transit fit reveals that high b alone cannot account for the anomalously short transit duration. We instead developed a model to account for both of these properties through a photometric transit fit that takes into consideration the photoeccentric framework of Dawson & Johnson (2012), as shown in Figure 1.

We modeled the standard transit parameters, including orbital period (P), time of conjunction (T_C ), planet–star radius ratio (R_p /R_⋆), and b, along with the expected stellar density assuming a circular orbit, ρ_⋆,circ, to obtain a model that encodes information about the true orbital eccentricity of the planet according to Kipping et al. (2012b). We derived this dynamical information from our results by comparing our modeled ρ_⋆,circ to the true stellar density, ρ_⋆, represented by the median of our EXOFASTv2 ρ_⋆ posterior (Section 3). A value of ρ_⋆,circ greater than ρ_⋆ would imply that the planet transited faster than expected and vice versa, given an initial assumption of e = 0. Breaking from this assumption, however, we calculated which values of e and the argument of periastron (ω) were necessary to account for the unusually fast transit, subsequently bringing ρ_⋆,circ into agreement with ρ_⋆. For both parameters, we calculated posterior probability distributions using the log-likelihood function (Dawson & Johnson 2012),

$$\begin{eqnarray}&&\mathrm{log}P(e,\omega | {\rho }_{\star },{\rho }_{\star ,\mathrm{circ}})=-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{g{\left(e,\omega \right)}^{3}{\rho }_{\star }-{\rho }_{\star ,\mathrm{circ}}}{g{\left(e,\omega \right)}^{3}{\sigma }_{{\rho }_{\star }}}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&g(e,\omega )=\displaystyle \frac{1+e\sin \omega }{\sqrt{1-{e}^{2}}},\end{eqnarray} \tag{ 2 }$$

following the notation of Kipping (2010) and Kipping et al. (2012b).

Constraints on ω using this method tend to be broad, but they are sufficient to determine if a transit occurs closer to periastron (as is the case for Kepler-1704 b) or apastron. On the other hand, we were able to constrain the eccentricity of Kepler-1704 b here with high certainty. We found that the 68% credible interval for eccentricity is 0.901–0.970 (Figure 2).

Zoom In Zoom Out Reset image size

In a previous analysis of the photoeccentric effect in Kepler transit data, Dawson et al. (2015) found ${0}_{-0}^{+1}$ giant planets on highly eccentric orbits that are likely undergoing tidal circularization. This nondetection refuted the hypothesis of Socrates et al. (2012) that approximately four such planets should be detected if HEM is the dominant hot Jupiter migration mechanism. Dawson et al. (2015) only considered planet candidates with three or more transits, to more accurately account for the completeness of Kepler pipeline detections (e.g., Christiansen et al. 2020), so Kepler-1704 was not included in their analysis.

Assuming tidal decay at constant angular momentum, the highest allowed values of eccentricity from our photoeccentric modeling would produce a final orbital period below 10 days, the canonical threshold for hot Jupiters. Therefore, based on just this photoeccentric effect analysis, Kepler-1704 b is a candidate proto–hot Jupiter. However, additional orbital characterization via RV monitoring of the host is needed to refine the eccentricity and the nature of Kepler-1704 b.

We acquired 15 high-resolution spectra of Kepler-1704 with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on the Keck I telescope in support of our Doppler monitoring of the Kepler-1704 system. The baseline of these observations spans nearly a decade. For each observation, the starlight passed through a heated iodine cell before reaching the slit to enable the precise wavelength calibration of each RV measurement.

We did not acquire a high signal-to-noise ratio (S/N) template spectrum, as is typical for HIRES RV observations (e.g., Howard et al. 2010). Instead, we identified a preexisting "best-match" template spectrum in the HIRES spectral library following Dalba et al. (2020a). The best-match star was HD 203473, a brighter G6V star with similar spectroscopic properties to Kepler-1704 according to a SpecMatch–Emp ¹⁶ analysis (Yee et al. 2017). The use of a best-match template incurs extra uncertainty in addition to internal RV errors. Following conservative estimations by Dalba et al. (2020a, their Table 2), we added 6.2 m s⁻¹ to our internal RV errors in quadrature to account for this method. After swapping in the template of HD 203473, the RV extraction proceeded following the standard forwarding techniques employed by the California Planet Search (e.g., Howard et al. 2010; Howard & Fulton 2016).

We provide the full RV data set for Kepler-1704 in Table 1. The uncertainties listed include the additional uncertainty incurred by the matched-template method of RV extraction (Dalba et al. 2020a). We also include corresponding S_HK activity indicators derived from the Ca ii H and K spectral lines (Wright et al. 2004; Isaacson & Fischer 2010).

We note that the first RV measurement (from BJD = 2,455,669) is the least precise observation in the series. Its uncertainty is three standard deviations above the mean. This larger error is not surprising, as the exposure time for the spectrum used to measure that RV was substantially shorter than the others. The resulting best-fit velocity in each 2 Å chunk of spectrum, which typically only contains one stellar and one iodine line, was less precise, leading to the larger error in RV. When folded on the ephemeris of Kepler-1704 b, this data point occupies a noncritical phase in the orbit. However, this data point extends the baseline of the RV observations by 826 days and is critical to our consideration of acceleration in the Kepler-1704 system (Section 4.2). Although there is no obvious reason to exclude this data point from our analysis besides its larger uncertainty, we will treat this data point with skepticism moving forward.

In Section 3, we model the RVs and transits simultaneously, confirming that the orbital period of Kepler-1704 b is accurately represented by the time elapsed between the two Kepler transits (988.88 days) and not half of that value. Visual inspection of the RV data folded on an orbital period of 494.44 days suggests no Keplerian signal at this periodicity. Therefore, we hereafter do not consider the possibility that another transit occurred during the gap in observations between Kepler quarters 7 and 8.

The Kepler-1704 system has previously been observed in several imaging campaigns (for a summary, see Furlan et al. 2017). To explore the existence of bound or background stellar neighbors, we present three data sets acquired from the Exoplanet Follow-up Observing Program. ¹⁷

The first imaging data set comprises AO images from the PHARO instrument (Hayward et al. 2001) at the 200 inch telescope at Palomar Observatory, as published by Wang et al. (2015b). This work used a three-point dither pattern to obtain a set of images in the K_s band that were combined and searched for stellar companions (Figure 3, left panel). Wang et al. (2015b) claimed two detections: one source with ΔK_s = 3.3 with a separation and position angle (PA) of 547 and 1570, respectively, and another source with ΔK_s = 4.6 with a separation and PA of 319 and 3055, respectively. Both detections are visible in the left panel of Figure 3. The source with PA = 1570 (indicated by a green vertical arrow) is resolved by Gaia (Gaia Collaboration et al. 2016, 2021) and has the EDR3 source ID 2136191732305041920 (hereafter Gaia-213, for simplicity). The parallax and proper motion of Kepler-1704 and Gaia-213 as measured by Gaia definitively show that these two stars are not gravitationally associated. The other source identified by Wang et al. (2015b), as well as a brighter source near the upper edge of the image that was not identified by Wang et al. (2015b), are not present in the Gaia EDR3 catalog. In the left panel of Figure 3, these unidentified sources are shown with horizontal yellow arrows.

Zoom In Zoom Out Reset image size

The second imaging data set also comprises AO images from the PHARO instrument but in the Brγ filter, as published by Furlan et al. (2017). Surprisingly, only Kepler-1704 and Gaia-213 (at PA = 1570) are visible despite deeper magnitude limits near 319: ΔK_s = 4.9 versus ΔBrγ = 7.0 (Wang et al. 2015b; Furlan et al. 2017). The time elapsed between the epochs of imaging, roughly 1 month, is too short to explain the discrepancy.

The solution to this conundrum lies in the relative positioning of the two sources in question with respect to the positioning of Kepler-1704 and Gaia-213. The separation and PA between these two pairs are identical. Visual inspection suggests that the contrast between the stars in each pair is also similar. Thus, our conclusion is that the two sources identified by yellow horizontal arrows in the left panel of Figure 3 are spurious duplications of Kepler-1704 and Gaia-213 caused by an accidental image alignment error.

The third imaging data set comprises AO images from the NIRC2 instrument (Wizinowich et al. 2000) at the Keck II telescope, as published by Furlan et al. (2017). Observations were taken in the Brγ filter, and the field of view was too small to include any of the other sources (astrophysical or spurious) mentioned previously (Figure 4). The NIRC2 data yield a nondetection of a stellar neighbor within 2'' with delta-magnitude limits of 8.4 and 8.7 at separations of 05 and 1 0, respectively (Furlan et al. 2017). Since the NIRC2 observations of Kepler-1704 provide the strongest constraints on neighboring stars, we continue our analysis using only these data.

Zoom In Zoom Out Reset image size

We used the NIRC2 contrast curve (i.e., 5σ limiting delta magnitude as a function of separation) to derive the corresponding limiting mass for a bound companion. First, we downloaded a MESA Isochrones and Stellar Tracks (MIST) isochrone (Paxton et al. 2011, 2013, 2015; Choi et al. 2016; Dotter 2016) from the MIST web interpolator. ¹⁸ We provided values of initial stellar metallicity, extinction, and age based on the system modeling described in Section 3. This isochrone provide a numerical relationship between stellar mass and absolute K_s magnitude, which we treated interchangeably with Brγ. After converting absolute magnitude to apparent magnitude (using the distance from Section 3), we interpolated the ΔK_s values with those measured by NIRC2 to calculate an upper limit of companion mass as a function of projected separation (Figure 4). At wider separations, the delta-magnitude values exceeded those in the MIST isochrone. For those separations, we instead interpolated a 5 Gyr brown dwarf isochrone from Baraffe et al. (2003). Beyond a projected separation of ∼200 au, we find that any companion in the Kepler-1704 system must have a mass below ∼32 M_J.

The converged EXOFASTv2 fit yielded bimodal posterior probability distributions for the stellar mass (M_⋆) and age (Figure 7). The region of parameter space preferred by the MIST stellar evolution models, as influenced by all of the Kepler-1704 data, exists near the subgiant branch, as we suspected based on the SpecMatch radius estimation. EXOFASTv2 found that multiple stellar ages and surface gravity values ($\mathrm{log}g$) correspond to the T_eff prior, meaning that the bimodality is astrophysical and not due to inadequate posterior sampling. The bimodality propagates to the semimajor axis (a) of Kepler-1704 b and, to a lesser extent, its mass (M_p ; Figure 7).

Zoom In Zoom Out Reset image size

Since we could not distinguish between the two families of solutions with the data of the Kepler-1704 system in hand, we adopted the strategy of Ikwut-Ukwa et al. (2020) and divided the solutions at a fiducial M_⋆ value of 1.21 M_⊙, which corresponds to the trough between the posterior probability peaks in Figure 7. The lower stellar mass, older age solution contains 51.8% of the posterior samples, which we treated as a slight preference over the higher stellar mass, younger solution. Therefore, in Tables 2 and 3, we only publish the parameters for the preferred, lower stellar mass solution. The only exception is the planet mass, M_p , for which we determine the 68% confidence interval using the full posterior distribution. None of our interpretations of the nature or formation history of Kepler-1704 b are changed by considering the alternate solution.

A photometric dimming event with a depth corresponding to a giant planet transit can be created by substellar or stellar objects or various systematic signals (e.g., Brown 2003; Torres et al. 2005; Cameron 2012; Foreman-Mackey et al. 2016; Dalba et al. 2020a). False-positive signals can be harder to identify for longer (compared to shorter) orbital periods owing to the difficulty in quantifying the reliability of genuine transit events from similarly long-period exoplanets (e.g., Thompson et al. 2018). Indeed, Santerne et al. (2016) measured a 55% false-positive rate for Kepler giant planets within 400 days of the orbital period. For these reasons, long-period giant planet candidates like Kepler-1704 b must be vetted with Doppler spectroscopy before any weight is placed upon their standing as a genuine planet.

Our 10 yr baseline of RV measurements for Kepler-1704 confirmed the genuine planetary nature of Kepler-1704 b. It also confirmed the 988.88 day orbital period, placing Kepler-1704 b among the top five longest-period (noncontroversial) transiting exoplanets with precisely measured periods known to date. ²¹ With a semimajor axis of 2.03 au and an orbital eccentricity of 0.92, its elongated orbit brings it within 0.16 au of its host star and then slingshots it out to 3.9 au—the largest apastron distance of any transiting exoplanet with a known orbital period and eccentricity. Figure 8 is a diagram showing the orbit of Kepler-1704 b relative to those of Jupiter, the solar system terrestrial planets, and HD 80606 b. The RV data also contain a slight, although tentative, acceleration (${0.0031}_{-0.0027}^{+0.0029}$ m s⁻¹ day⁻¹) that possibly indicates the existence of an outer companion.

Zoom In Zoom Out Reset image size

The equilibrium temperature (T_eq) for Kepler-1704 b as shown in Table 3 is calculated following

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}={T}_{\mathrm{eff}}\sqrt{\displaystyle \frac{{R}_{\star }}{2a}},\end{eqnarray} \tag{ 3 }$$

which assumes no albedo and perfect heat redistribution (Hansen & Barman 2007). However, including a factor of 1/($\sqrt{1\pm e}$) in this equation suggests that T_eq varies from ∼180 K at apastron to ∼900 K at periastron. This substantial ∼700 K swing in temperature likely affects the atmosphere on Kepler-1704 b.

In the following sections, we will investigate the possibility of companions, migration history, interior composition, and atmospheric characterization prospects for Kepler-1704 b. We take advantage of the fact that this planet's orbital period, eccentricity, and radius are known precisely, which is remarkable for an exoplanet with its orbital properties.

As described in Section 1 and extensively in the broader orbital dynamics literature (e.g., Naoz 2016), the presence of an outer planetary or stellar companion may have direct consequences on the migration history of a giant planet. For Kepler-1704, archival AO imaging data yield a nondetection of stellar companions beyond ∼100 au and upper mass limits on such a companion down to ∼50 au (see Figure 4 and Section 2.3). In the following sections, we exploit our long baseline of RV observations to improve upon these limits with an injection-recovery test (Section 4.2.1), an RV trend analysis (Section 4.2.2), and a chaos indicator analysis (Section 4.2.3).

4.2.1. RV Injection-recovery Test

We characterized the sensitivity of our RV data set to additional bound companions by running injection-recovery tests, in which we added synthetic signals to our RV data and converted the signal recovery rate into a map of search completeness. We used RVSearch (Rosenthal et al. 2021), an iterative periodogram search algorithm, to search for evidence of additional companions to Kepler-1704 b in the RV data and perform these tests. We initialized RVSearch with the best-fit Keplerian model for Kepler-1704 b and searched for additional companions with orbital periods spanning 2–10,000 days. We found no evidence for additional companions in this period range. Once the search was completed, RVSearch injected synthetic planets into the data and repeated the additional iteration to determine whether it recovered these synthetic planets. We ran 3000 injection tests for Kepler-1704. We drew the injected planet a and ${M}_{p}\sin i$ from log-uniform distributions and drew eccentricity from the beta distribution with shape parameters α = 0.867 and β = 3.03, which Kipping (2013) found represented the sample of RV-observed exoplanets. After RVSearch performed the injection-recovery tests, we measured the search completeness across a wide range of a and ${M}_{p}\sin i$ by determining the fraction of recovered synthetic signals in localized regions of a and ${M}_{p}\sin i$.

Figure 9 shows a pair of search completeness results, one of which includes the first low-S/N RV data point (left panel) and one of which excludes it (right panel). In both cases, our RV sensitivity to companions beyond the orbital separation of Kepler-1704 b is limited, dropping below 50% completeness at 4 M_J beyond 4 au. The sparsity and high rms of the RV data set drive the high lower limit on detectability in ${M}_{p}\sin i$, and the nearly 10 yr observational baseline sets the sharp change in completeness around 3 au.

Zoom In Zoom Out Reset image size

4.2.2. RV Trend Analysis

To build upon the injection-recovery test, we conducted a complementary analysis of acceleration (i.e., a long-term RV trend) in the Keck-HIRES RVs. This analysis focused specifically on partially sampled signals from giant planets, substellar objects, or stars that could be lurking undetected in the outer reaches of the Kepler-1704 system. When combined with a nondetection from the AO imaging, RV trends can greatly reduce the parameter space that a possible undetected companion could occupy (e.g., Crepp et al. 2012; Kane et al. 2019; Dalba et al. 2021a).

The EXOFASTv2 fit to the transit, RV, and SED (Section 3) included a parameter for "RV slope" ($\dot{\gamma }$), which quantified any acceleration measured from the RVs. As shown in Table 3, we made a low-significance detection of acceleration: $\dot{\gamma }\,={0.0031}_{-0.0027}^{+0.0029}$ m s⁻¹ day⁻¹. To refine the mass (M_c ) and orbital distance (a_c ) of the companion that could have caused this RV drift, we simulated RVs over a grid of scenarios broadly following the procedure of Montet et al. (2014).

First, we subtracted the maximum-likelihood EXOFASTv2 solution for Kepler-1704 b from the Keck-HIRES RV data but without including the acceleration (i.e., we set $\dot{\gamma }=0$). In doing so, we also inflated the RV uncertainties (${\sigma }_{{v}_{r}}(t)$) to account for the fitted RV jitter (Table 3). The resulting RV time series (v_r (t)) only contained the long-term trend.

Next, we defined a logarithmically spaced 30 × 30 grid in companion mass (1 M_J < M_c < 1 M_☉) and semimajor axis (4 < a/au < 200). The mass boundaries were chosen to complement the constraints from the RVsearch injection-recovery tests and the AO imaging (Section 2.3). The orbital distance boundaries were chosen to span the gap between the apastron distance of Kepler-1704 b and the stringent upper boundary from the AO imaging.

At each point along the M_c –a_c grid, we drew 500 sets of the orbital elements {ω, e, i}, which are the argument of periastron, eccentricity, and inclination, respectively. We drew ω randomly from a uniform distribution over the interval [0, 2π], and we drew i randomly from a uniform distribution in $\cos i$ over the interval [0, 1]. For e, we drew values from the beta distribution from Kipping (2013) mentioned previously (Section 4.2.1). These random draws were meant to account for the variety of orbital configurations a massive companion could have.

Then, for each of the individual orbits, we simulated 50 sets of RV time series (${\hat{v}}_{r}(t)$) with a cadence matching v_r . Each of the 50 sets started at a different orbital phase spaced evenly across the entire orbit. This accounted for the fact that the Keck-HIRES observations could have sampled any portion of the companion's orbit.

Finally, we used a least-squares regression routine to minimize the familiar statistic ${\chi }^{2}={\sum }_{t}{[{v}_{r}{(t)-{\hat{v}}_{r}(t)]}^{2}/{\sigma }_{{v}_{r}}(t)}^{2}$. This minimization was necessary because the Keck-HIRES RVs are relative, not absolute. Assuming uncorrelated errors, we converted the 50 χ² values for each individual orbit to relative probabilities following $P\propto \exp (-{\chi }^{2}/2)$, and we summed the probabilities to effectively marginalize over the portion of the orbit captured by the data. We also summed the probabilities of the 500 sets of orbits at each grid point to effectively marginalize over all orbital properties other than M_c and a_c . Lastly, we normalized the map of probabilities such that 22.5 million probability calculations summed to unity, ultimately yielding relative-likelihood values. Figure 10 (left panel) shows the resulting map.

Zoom In Zoom Out Reset image size

The slight acceleration detected in the full set of RVs prefers companions within roughly 30 au and less massive than a few hundred Jupiter masses, although some probable solutions are still present at wide separation and high mass. Incorporating the upper mass limit from the AO imaging (Section 2.3) trims a correlated region of parameter space at the highest masses and largest orbital separations. Also, assuming that any companion with at least a 50% RVsearch recovery rate should have been detected, the trend analysis further refines the likely parameter space of a companion.

We repeated this entire analysis but after removing the first Keck-HIRES RV data point, as its timing and quality may have inaccurately affected the measured RV trend (Section 2.2). The resulting map of relative likelihoods calculated without the first RV data point is shown in the right panel of Figure 10. For context supporting the second map, we also conducted a second EXOFASTv2 fit without the first Keck-HIRES data point that was otherwise identical to the fit described in Section 3. The only appreciable difference between the two EXOFASTv2 fits was the value of $\dot{\gamma }$, which decreased in significance to −0.0002 ± 0.0029 m s⁻¹ day⁻¹ in the latter case. This difference manifests in the relative-likelihood map as a preference toward larger orbits (a ≳ 30 au) rather than smaller ones. The map is again complemented by the AO imaging upper limit and the region with an over 50% RVsearch recovery rate.

Within the parameter space we are exploring, it is also helpful to consider which companions would be capable of overcoming precession caused by general relativity and exciting the eccentricity of Kepler-1704 b through Kozai–Lidov cycles. Dong et al. (2014) calculated an approximate strength criterion for warm Jupiters (their Equation (5)) that we apply to Kepler-1704 b. In the limiting case of an initially circular orbit that is much longer-period than that of a hot Jupiter, we identify the region of M_c –a parameter space with objects that are unable to have excited the eccentricity of Kepler-1704 b (Figure 10, yellow hatched region). By all of our other analyses, we cannot rule out the existence of a companion at the lowest masses and largest separation we consider. However, such a companion is also too low mass and orbits too far from Kepler-1704 b to overcome general relativity precession through Kozai–Lidov interactions.

Considering all of the companion analyses together yields three conclusions for possible outer companions in the Kepler-1704 system. First, for M_c ≳ 700 M_J, we should have either recovered the signal in the Keck-HIRES RVs (≥50% recovery rate) or detected the source directly in the AO imaging for nearly all values of a. Second, for 50 ≲ M_c /M_J ≲ 700, those with a ≲ 40 au should have been recovered by the RV data, and those with a ≳ 150 au should have been detected in the AO imaging. A companion with a separation between these values could go undetected. Third, we do not have sensitivity to companions with M_c ≲ 50 M_J within ∼150 au, so there could be substellar or planetary companions in this region. At these masses, our RV trend analysis reveals a preference for companions with a ≳ 30 au. However, some companions in this region of parameter space would be unable to excite the observed eccentricity of Kepler-1704 b.

4.2.3. MEGNO Simulations

To test whether additional constraints can be placed on the orbital configurations of the potential outer companion, we ran a dynamical simulation using the Mean Exponential Growth of Nearby Orbits (MEGNO) chaos indicator (Cincotta & Simó 2000). The MEGNO indicator demonstrates whether a specific system configuration would lead to chaos after a certain integration time by distinguishing between quasiperiodic and chaotic evolution of the bodies within the system (e.g., Hinse et al. 2010). The final MEGNO value returned for a specific orbital configuration is useful for determining the stochasticity of the configuration, where chaos is more likely to result in unstable orbits for planetary bodies. With a grid of orbital parameters, a MEGNO simulation can provide valuable information on the orbital configurations that are favored by dynamical simulations and reject configurations that return chaos results.

The MEGNO simulation to explore the dynamically viable locations for various outer companions was carried out within the N-body package REBOUND (Rein & Liu 2012) with the symplectic integrator WHFast (Rein & Tamayo 2015). We used the stellar and planetary parameters from Tables 2 and 3, respectively. We provided a linear-uniform grid in semimajor axis (20–60 au) and companion mass (1–100 M_J) that aligned with the higher-likelihood region in Figure 10 (right panel). The eccentricity of the outer companion was set to zero. The simulation was integrated for 20 million yr with a time step of 0.034 yr (∼12.4 days). This time step was chosen to be 1/80 of the orbital period of Kepler-1704 b, a fourth of the recommended value (Duncan et al. 1998), to increase the sampling near the periastron passage of this highly eccentric planet. The integration was set to stop and return chaos results if any of the planetary orbits started extending beyond 100 au.

Figure 11 shows the grid of results of the MEGNO simulation. Each grid point is color-coded according to the final MEGNO value for the orbital configuration of that outer companion. A MEGNO value around 2 (green) is considered nonchaotic (Hinse et al. 2010) and is thus a dynamically viable region where the outer companion could exist without making the system chaotic. Grid points in red indicate simulations that returned chaotic results, and those in white indicate irregular events, such as close encounters and collisions, all of which are unfavorable configurations for an outer companion.

Zoom In Zoom Out Reset image size

Only a few stripes of parameter space contain orbital configurations that lead to chaos. The stripes indicate systems where nonlinear eccentricity secular resonances overlap, leading to secular chaos (e.g., Lithwick & Wu 2011; Wu & Lithwick 2011). The stripes become less defined at wider orbits because our simulation has not been run long enough to capture the resonance evolution. Besides these narrow regions, this analysis fails to rule out any extra substantial area of parameter space where a massive companion could exist.

Continuing our discussion of results, we now shift the attention from the outer reaches of the Kepler-1704 system back to Kepler-1704 b itself.

With the measured mass and radius of Kepler-1704 b, along with other system properties, we retrieved the mass of its heavy elements or its bulk metals (M_z ) and calculated its bulk metallicity (Z_p ≡ M_z /M_p ) following Thorngren & Fortney (2019). Briefly, we modeled the thermal evolution of Kepler-1704 b using 1D structure models with a core composed of a rock/ice mixture at equal amounts, a homogeneous convective envelope made of an H/He rock/ice mixture, and a radiative atmosphere. The atmosphere models were interpolated from the grid of Fortney et al. (2007). Samples were drawn from the posterior probability distributions for planet mass, radius, and age (Section 3), and the heavy-element mass was adjusted in the structure models to recover the planet radius.

This analysis relied on two assumptions. First, we assumed that the planet radius is not inflated (e.g., Laughlin 2018) because of the average irradiation flux received by Kepler-1704 b (see Table 3), which is well below the canonical 2 × 10⁸ erg s⁻¹ cm⁻² empirical threshold for giant planets (Demory & Seager 2011; Miller & Fortney 2011; Sestovic et al. 2018). Second, we neglected any internal heating from circularization tides. We assumed that tides are an inefficient means of heating Kepler-1704 b, as evidenced by its tidal circularization timescale (τ_circ) given by Equation (3) from Adams & Laughlin (2006),

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{circ}} & = & 1.6\,\mathrm{Gyr}\,\left(\displaystyle \frac{Q}{{10}^{6}}\right)\left(\displaystyle \frac{{M}_{p}}{{M}_{J}}\right){\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-3/2}\\ & & \times \,{\left(\displaystyle \frac{{R}_{p}}{{R}_{J}}\right)}^{-5}{\left(\displaystyle \frac{a}{0.05\,\mathrm{au}}\right)}^{13/2},\end{array}\end{eqnarray} \tag{ 4 }$$

where Q is tidal quality factor that is assumed to be 10⁶ (similar to Jupiter). As listed in Table 3, τ_circ is 80,000 Gyr, much longer than the age of the universe, even considering the error introduced by our estimate value for Q.

The metallicity retrieval was complicated slightly by the bimodal probability distribution for age that we inferred from the comprehensive system modeling (Figure 7). Instead of using separate normal priors for stellar mass and age, we used a bivariate Gaussian kernel-density estimate. Then, we sampled the posterior with a Markov Chain Monte Carlo technique.

The results of the bulk metal mass retrieval are shown in Figure 12. Despite the bimodality in age, the marginalized posterior probability distribution for bulk metallicity is a near-normal distribution at Z_p = 0.12 ± 0.04, corresponding to M_z ≈ 150 M_⊕. To calculate the stellar metallicity (Z_⋆), we assumed that the iron abundance ([Fe/H]) scales with total heavy-metal content such that Z_⋆ ≡ 0.0142 × 10^[Fe/H] (Asplund et al. 2009; Miller & Fortney 2011), which yields Z_⋆ = 0.0229 ± 0.0031. Finally, we calculated the bulk metallicity enrichment relative to stellar for Kepler-1704 b as Z_p /Z_⋆ = 5.2 ± 1.9.

Zoom In Zoom Out Reset image size

We place the bulk metal mass and metallicity enrichment in the context of other cool (T_eq ≲ 1000 K), weakly irradiated (〈F〉 < 2 × 10⁸ erg s⁻¹ cm⁻²) giant exoplanets from the Thorngren et al. (2016) sample ²² in Figure 13. By metal mass and enrichment, Kepler-1704 b is entirely consistent with the known trends. It contains more metal mass than its lower-mass counterparts, but it is broadly less enriched in metals relative to its host star. These findings are consistent with the theory of core accretion as its formation scenario, followed by a period of late-stage heavy-element accretion (e.g., Mousis et al. 2009; Mordasini et al. 2014). Planet Kepler-1704 b is similar to the other high-mass (M_p ≳ 2 M_J) giant planets in that it orbits a metal-rich star, something that has been predicted by population synthesis models (e.g., Mordasini et al. 2012) and likely relates to the correlation between host star metallicity and giant planet occurrence (e.g., Gonzalez 1997; Santos et al. 2004; Fischer & Valenti 2005).

Zoom In Zoom Out Reset image size

In Figure 13 (right panel), we also include a prediction from Ginzburg & Chiang (2020). They modeled concurrent gas accretion and mergers during giant planet formation as an alternate means of explaining the heavy-metal content of giant planets. The scatter in the data enclosed by the dotted black lines can be explained by the intrinsically chaotic nature of mergers, even if all systems evolve from nearly identical conditions, as quantified by a critical core mass of 10 M_⊕. We find that the mass and bulk metallicity enrichment of Kepler-1704 b are also consistent with the theory of concurrent gas accretion and mergers.

It is interesting to consider how trends in heavy-element mass, metal enrichment, and total planet mass relate to other orbital and stellar properties. In Figure 14, we show the relative residuals (calculated/best fit) of heavy-element mass and metallicity enrichment relative to stellar as a function of eccentricity for the Thorngren et al. (2016) sample of weakly irradiated giant exoplanets and Kepler-1704 b. As noted by Thorngren et al. (2016), there is no discernible trend in either quantity. However, given how sparsely populated the high-eccentricity region is, it is worthwhile to consider the (now) five systems with e > 0.6 individually. The residual heavy-element mass and metallicity enrichment of Kepler-1704 b and HD 80606 b are nearly identical, as are their orbital eccentricity and planet mass. However, HD 80606 b likely migrated via secular perturbations with HD 80607 (e.g., Wu & Murray 2003; Fabrycky & Tremaine 2007; Winn et al. 2009), whereas we are unsure if similar interactions with a planetary or stellar companion have influenced the migration history of Kepler-1704 b. If Kepler-1704 b and HD 80606 b followed different migration pathways, there is no evidence in their bulk metallicity to distinguish them. The residual heavy-element mass and the metallicity enrichment of these two planets are significantly different than those of HD 17156 b, which has e ≈ 0.67 (Fischer et al. 2007; Bonomo et al. 2017). Unlike HD 80606 b, HD 17156 b has no stellar companion, and its orbit is nearly aligned with its host star (Cochran et al. 2008; Narita et al. 2008; Barbieri et al. 2009). Also, HD 17156 b's orbital period is almost 2 orders of magnitude shorter than that of Kepler-1704 b. Therefore, it is perhaps not surprising that these planets experienced different formation histories that could account for the metallicity differences. The final two high-eccentricity planets in Figure 14 (KOI-1257 b and Kepler-419 b) have relatively imprecise residual heavy-element masses and metallicity enrichments. Planet KOI-1257 b is thought to be in a binary star system, possibly pointing to Kozai migration (Santerne et al. 2014). On the other hand, Kepler-419 b is joined by a massive outer giant planet that has a low mutual inclination, such that Kozai migration is likely not a viable migration theory (Dawson et al. 2014). The overall lack of a clear trend between heavy-element mass or metallicity enrichment and the presence of a companion and/or high stellar obliquity is likely in part a result of the small number of data points. However, it could also suggest that the heavy-element accretion occurs before or independently from the various channels of eccentricity excitement.

Zoom In Zoom Out Reset image size

The bulk heavy-element mass retrieval suggested that ∼150 M_⊕ of metals should exist within Kepler-1704 b. The distribution of those heavy elements within the planet can possibly affect the composition of its atmosphere. Specifically, Thorngren & Fortney (2019) showed that the bulk metallicity places an upper limit on atmospheric metallicity. For Kepler-1704 b, the core-free 2σ upper limit (Z_p = 0.2) for atmospheric metallicity is 35.7× solar. A natural next step in the exploration of Kepler-1704 b is to test this prediction via atmospheric characterization.

Considering only orbital period or semimajor axis, Kepler-1704 b is a rare opportunity for transmission spectroscopy (Seager & Sasselov 2000). However, long-period exoplanets pose specific challenges to this kind of technique. Not only are transits of such planets geometrically rare, but their timing is often uncertain. Since only two transits of Kepler-1704 b have been observed, the presence of extreme transit timing variations (TTVs; Wang et al. 2015a) cannot be ruled out (e.g., Dalba & Muirhead 2016; Dalba & Tamburo 2019; see Section 5.3). Furthermore, atmospheric temperature will (to first order) decrease with increasing orbital distance. As a result, atmospheres will be cooler, and scale heights and transmission spectrum features will be smaller. Surprisingly, this can be balanced by low surface gravity, as would be the case if Saturn was subject to transmission spectroscopy (Dalba et al. 2015). The transiting geometry of the long-period Kepler-1704 b also makes it a unique candidate for testing theories of atmospheric refraction (e.g., Sidis & Sari 2010; Dalba 2017; Alp & Demory 2018) that have not yet been observationally tested (Sheets et al. 2018).

However, considering the large radius of the subgiant Kepler-1704 and the high mass of Kepler-1704 b, this system is a challenging target for transmission spectroscopy. With a surface gravity of 86 m s⁻² and the average equilibrium temperature of 254 K from Table 3 (assuming no albedo), the atmospheric scale height is only ∼12 km, which corresponds to 1 part per million (ppm) in the transmission spectrum. The out-of-transit stellar mirage caused by refraction also scales with the atmospheric scale height, making such a detection similarly difficult (e.g., Dalba 2017).

On the other hand, we used Equation (3) to estimate that T_eq at periastron, which is within several days of transit, is ∼900 K. This suggests a 3.6× increase in the atmospheric scale height and transmission spectrum feature size. Although 4 ppm is still beyond the reach of current and future facilities, we caution that our intuition for predicting favorable transmission spectroscopy targets is largely based on our current understanding of hot, close-in exoplanet atmospheres. This possibly warrants skepticism. If Saturn were a transiting exoplanet, its warm stratosphere and active photochemistry would produce an ∼90 ppm absorption feature in its transmission spectrum at 3.4 μm (Dalba et al. 2015). Considering only Saturn's T_eq as defined in Equation (3) would underpredict its amenability to transmission spectroscopy. Other long-period giant exoplanets may prove surprising as well.

Even if transmission spectroscopy is not a viable atmospheric characterization technique, the 0.16 au periastron distance of Kepler-1704 b caused by its extreme eccentricity possible qualifies it for an IR phase-curve analysis.

4.4.1. IR Phase-curve Analysis

To predict the expected thermal signature of the planet during periastron passage, we calculated the IR phase curve for Kepler-1704 b during one complete orbital period. These calculations followed the methodology of Kane & Gelino (2011) using the stellar and planetary parameters provided in Tables 2 and 3, respectively. We assumed a passband of 4.5 μm and a Bond albedo of zero, and we calculated the flux ratio of planet to star using the "hot dayside" and "well-mixed" models. These models represent the extremes of heat redistributions, as they assume reradiated energy over 2π and 4π sr, respectively. The full IR phase curves for both models are shown in Figure 15, along with a zoomed panel that shows the location of the periastron passage.

Zoom In Zoom Out Reset image size

There are several caveats to this calculation. We assumed an instantaneous response of the planetary absorption and IR emission, whereas the radiative and advective timescales will determine the nature of the phase lags in the thermal emission profiles (Langton & Laughlin 2008; Cowan & Agol 2011a). This, combined with the blackbody emission and zero albedo assumptions, means that the calculations presented in Figure 15 may be considered as an upper limit on the expected IR emission. Furthermore, the variation in temperature would also alter the atmospheric composition. Some of the energy would be converted into latent heat to dissociate larger molecules or particulates. There would also be an interconversion between CO and CH₄ (e.g., Visscher 2012). The timescale of this reaction, and also the vertical mixing timescale, should be considered to produce a more accurate model of the phase curve. We leave these considerations for a future work and instead derive a first-order upper limit on the phase-curve emission.

Despite the various assumptions that apply to this phase-curve modeling, the order of magnitude (${ \mathcal O }$(10²) ppm) of the thermal flux increase is likely accurate. Several instruments on board JWST will have sensitivity in the near- to thermal-IR and, based on preliminary noise floor expectations (Greene et al. 2016), should be capable of detecting the Kepler-1704 b phase variation. Borrowing from solar system intuition, 4–6 μm is likely a promising wavelength for such an observation. Jupiter's and Saturn's atmospheres have a low opacity in this wavelength region that exists between bands of methane and phosphine where radiation from 5 to 8 bars can escape (Irwin et al. 2014). Jupiter's radiance near 5 μm even exceeds that at mid-IR wavelengths (e.g., Irwin et al. 1998, 2014; Fletcher et al. 2009a).

The periastron passage of Kepler-1704 b occurs over ∼10 days and includes the transit. Low-cadence time-series observations from the F444W filter of NIRCam, for example, could detect the peak flux ratio and the width of the feature, assuming that the visit-to-visit photometric variability does not overwhelm the astrophysical signal. Including at least one high-cadence, longer visit at or following periastron would also be valuable because the phase curve may exhibit "ringing" as the hot spot from periastron rotates in and out of view (e.g., Cowan & Agol 2011b). This effect is not featured in our simulation of Kepler-1704 b, which assumed pseudosynchronous rotation (Hut 1981). However, given the inefficiency of tides at the periastron distance of Kepler-1704 b, this assumption may be an oversimplification. The detection of a ringing oscillation in the IR phase curve would test this assumption and possibly directly yield the effective planetary rotation rate.

The detectability of the thermal phase curve for Kepler-1704 b should be explored in more detail with an atmospheric structure code (e.g., Mayorga et al. 2021). Such an effort is beyond the scope of this paper and should likely wait until JWST is launched and commissioned. In addition to a broadband detection of phase variability, the prospects for spectroscopic detection should also be investigated. Transmission spectroscopy may not be an effective tool to measure atmospheric composition (e.g., metallicity), so any other possible method would be extremely useful. Atmospheric metal enrichment (relative to stellar) is a specifically valuable property to measure because it can be compared to the planet's bulk metallicity enhancement (Section 4.3). One prediction would be an atmospheric metallicity less than the bulk metallicity if some heavy elements comprise a planetary core or there is otherwise an increasing gradient in metals with depth. However, for Jupiter and Saturn, recent high-precision gravity data and well-established atmospheric composition results suggest a more complicated picture (Niemann et al. 1998; Wong et al. 2004; Fletcher et al. 2009b; Wahl et al. 2017; Guillot et al. 2018; Iess et al. 2019; Müller et al. 2020b). More elaborate theories including inverse compositional gradients (e.g., Debras & Chabrier 2019) are needed to explain Jupiter and Saturn and could possibly be refined through atmospheric characterization of exoplanets like Kepler-1704 b.

Based on the optimistic prospect of JWST observations, we determined the timing of transits and periastron passages of Kepler-1704 b occurring in the next 10 yr (Table 4). For each event, we checked for visibility from JWST using the General Target Visibility Tool. ²³ This tool only predicts visibility through the end of 2023, but we assumed the same visibility of Kepler-1704 in later years. The 2023 transit of Kepler-1704 b will not be visible to JWST, but the 2028 transit will be visible. The 2025 transit will occur within 24 hr after the visibility window closes, and the 2031 transit will occur roughly 6 days after the visibility window opens. If the solar avoidance restrictions change after launch, these transits may or may not be visible to JWST. The periastron passage of Kepler-1704 b occurs several hours before transit, so its visibility is similar. However, as shown in Figure 15, the peak of the thermal flux ratio spans ∼10 days. Even if the exact moment of periastron is (or is not) visible, some portion of the event is expected to be visible to JWST.

Table 4. Future Transit and Periastron Timing Predictions

	Conjunction (Transit) Time		Periastron Time		JWST
Epoch a	BJDTDB	UTC	BJDTDB	UTC	Visibility b
5	2,460,016.0902 ± 0.0046	2023-03-12 14:10	2,460,015.78 ± 0.20	2023-03-12 06:37	None
6	2,461,004.9714 ± 0.0055	2025-11-25 11:19	2,461,004.66 ± 0.20	2025-11-25 03:46	Partial
7	2,461,993.8525 ± 0.0065	2028-08-10 08:28	2,461,993.54 ± 0.20	2028-08-10 00:55	Full
8	2,462,982.7336 ± 0.0074	2031-04-26 05:36	2,462,982.42 ± 0.20	2031-04-25 22:03	Partial

Notes. The times listed here do not account for possible uncertainty owing to yet-undiscovered TTVs (see Section 5.3).

^a Epoch = 0 is defined as the first transit observed by the Kepler spacecraft. ^b The JWST visibility after 2023 December 31 is based on previous years' visibility. Epochs for which the full periastron passage of Kepler-1704 b partially falls outside of the predicted visibility windows are labeled as "Partial" (see the text).

Download table as:  ASCIITypeset image

4.4.2. Radio Emission

Unlike for transmission spectroscopy, the relatively high mass of Kepler-1704 b is beneficial to attempts to measure planetary radio emission. Lazio et al. (2010) searched for radio emission from HD 80606 b during a periastron passage but measured only an upper limit. That experiment was based on the expectation that the variation in planet–star distance over an eccentric orbit would lead to a dramatic increase in magnetospheric emission. Assuming that luminosity scales with the planet–star distance as L ∝ d^−1.6 (e.g., Farrell et al. 1999), the factor of 24.3 change in distance for Kepler-1704 b would produce a 165× increase in luminosity. While this is slightly smaller than the 200× increase expected for HD 80606 b, a future radio search may be aided by the fact that Kepler-1704 b can possibly emit at higher frequencies. We estimate that the upper-limit emission frequency as determined by the local plasma frequency in the emission region for Kepler-1704 b is

$$\begin{eqnarray}&&\nu =24\,\mathrm{MHz}\left(\displaystyle \frac{{\rm{\Omega }}}{{{\rm{\Omega }}}_{{\rm{J}}}}\right){\left(\displaystyle \frac{{M}_{p}}{{M}_{{\rm{J}}}}\right)}^{5/3}{\left(\displaystyle \frac{{R}_{p}}{{R}_{{\rm{J}}}}\right)}^{3},\end{eqnarray} \tag{ 5 }$$

where Ω is the angular rotation rate (Farrell et al. 1999; Lazio et al. 2004, 2010). In this equation, all values are scaled to those of Jupiter. For HD 80606b, tidal forces are expected to force the planet into pseudosynchronous rotation with a period of 39.9 hr (Hut 1981; Lazio et al. 2010). It is unlikely that this would also apply to Kepler-1704 b, for which the larger periastron distance renders tides inefficient. Therefore, the assumption of a Jupiter-like rotation period (∼9.9 hr) is reasonable. In that case, evaluating Equation (5) gives 310 MHz. Lazio et al. (2010) argued that this equation may actually underpredict the cutoff frequency of exoplanets, as it does for Jupiter, and suggested that the upper limit may be 60% larger. In that case, the cutoff frequency for Kepler-1704 b would be 497 MHz, which is more accessible to existing radio observatories than HD 80606b's 55–90 MHz.

A full simulation of the potential for radio emission from Kepler-1704 b is beyond the scope of this paper, and the ability to make such a detection, at least relative to previous attempts for HD 80606 b, will likely be hindered somewhat by the greater distance to the Kepler-1704 system. However, even our approximate calculation suggests that Kepler-1704 b is one of the best systems to investigate the magnetospheric response to a rapidly changing planet–star distance. Such an observation stands to extend the study of giant exoplanet magnetic fields beyond the innermost hot Jupiters (e.g., Cauley et al. 2019) and explore magnetic field generation in planets akin to Jupiter and Saturn.

Based solely on the measured orbital eccentricity, we discard disk migration as the explanation for the orbital properties of Kepler-1704 b. Papaloizou et al. (2001) showed that eccentricities up to ∼0.25 could be achieved through disk interactions for a variety of planet masses. However, eccentricity is generally damped by the disk for giant planets with M_p < 5 M_J (Bitsch et al. 2013). Recent work revisiting disk cavity migration argued for eccentricities up to 0.4 for giant planets (Debras et al. 2021), which is possibly a viable theory for other outer giant planets like Kepler-1514 b ($e={0.401}_{-0.014}^{+0.013}$), which also harbors an inner Earth-sized companion (Dalba et al. 2021b). Explaining the current orbit of Kepler-1704 b, however, requires excitation to high eccentricity by another body.

Through multiple analyses, we rule out stellar companions with mass greater than ∼700 M_J at most orbital separations in the Kepler-1704 system. As described in Section 4.2, undetected, less massive companions may still be present at a variety of separations and could have driven Kepler-1704 b to its high eccentricity through secular Kozai–Lidov perturbations (e.g., Wu & Murray 2003; Naoz et al. 2011). Also, star–planet Kozai migration from a stellar companion that was present when Kepler-1704 b formed but subsequently lost due three-body interactions also remains a possible explanation. However, motivated by our nondetection of a companion and only a tentative detection of acceleration in ∼10 yr of RV measurements, we discard secular perturbations as being the most likely explanation for the high eccentricity of Kepler-1704 b. We recommend that future dynamical simulations explore the areas of parameter space that we have not ruled out to see if a hidden companion could theoretically explain the properties of Kepler-1704 b (see Jackson et al. 2019).

This brings us to eccentricity excitation theories involving close, fast dynamical interactions. Specifically, could planet–planet scattering (e.g., Rasio & Ford 1996) provide an explanation for the eccentricity of Kepler-1704 b? Many aspects of the observed eccentricity distribution of giant exoplanets can be explained by planet–planet scattering (e.g., Moorhead & Adams 2005; Chatterjee et al. 2008; Raymond et al. 2010; Bitsch et al. 2020), including planets with eccentricities above 0.99 (Carrera et al. 2019). For Kepler-1704 b, we find that planet–planet scattering is consistent with its orbital properties, host star, and bulk interior properties. Object Kepler-1704 is metal-rich ([Fe/H] = 0.196 ± 0.057 from Table 2). Dawson & Murray-Clay (2013) demonstrated that metal-rich stars tend to host high-eccentricity hot Jupiters, which they interpreted as evidence supporting HEM by planet–planet scattering owing to the well-known correlation between stellar metallicity and giant planet occurrence (e.g., Santos et al. 2004; Fischer & Valenti 2005). Even though Kepler-1704 b is not a hot Jupiter, it is reasonable that it could have formed alongside other giant planets that were subsequently scattered. After many close encounters, possibly even tens of thousands (Carrera et al. 2019), Kepler-1704 b could have been driven to its current eccentricity. However, its final periastron distance was too far for tides to efficiently circularize the orbit, leaving Kepler-1704 b as the failed hot Jupiter that we have characterized here.

Now, we consider this migration pathway in the context of the bulk heavy-metal mass, which we found to be ∼150 M_⊕. This enrichment could have been acquired through late-stage accretion of planetesimals or core mergers with concurrent gas accretion (Section 4.3). For the former explanation, Shibata et al. (2020) found that migrating giant planets can capture tens of Earth masses worth of planetesimals that may otherwise not be available in situ (e.g., Shibata & Ikoma 2019). The amount of heavy elements accreted scales with increasing migration distance and decreasing migration timescale, both of which are expected for eccentricity excitement through planet–planet scattering. However, Shibata et al. (2020) also found that the most enriched giant exoplanets, containing more than ∼100 M_⊕ of heavy elements, likely require an additional source of enrichment.

The latter explanation mentioned above for how Kepler-1704 b acquired its heavy-metal mass involves the merger of cores during gas accretion (Ginzburg & Chiang 2020). The interaction of cores during this stage of giant planet formation, while gas is still present in the disk, is broadly consistent with the high eccentricity of Kepler-1704 b. Bitsch et al. (2020) found that such scattering events at this stage are common, and systems that start with more planetary embryos create giant planets with higher eccentricities as long as the damping rates for inclination and eccentricity are slow. Indeed, slow rates are required to reproduce the eccentricity distribution of the known giant planets (Bitsch et al. 2020). It is therefore possible that the same processes that led to the accretion of heavy elements for Kepler-1704 b also contributed to exciting its eccentricity. In reality, owing to the fact that planet mergers (or collisions) are less efficient at producing high eccentricities than scattering events (e.g., Ford & Rasio 2008; Jurić & Tremaine 2008; Anderson et al. 2020), some combination of the aforementioned theories along with planet–planet scattering after the dispersal of the gas disk likely produced the Kepler-1704 system as seen today.

Moving forward, it would be useful to compare this proposed formation history to other well-characterized long-period transiting giant planets. It will be particularly interesting to compare the bulk interior properties of giant planets in systems with and without outer companions that could have induced Kozai migration. If bulk heavy-element composition and migration mechanisms are linked, as may be the case for Kepler-1704 b, we might expect to find a correlation between interior properties, orbital properties, and the existence of companions.

A critical missing piece in our discussion of the migration of Kepler-1704 b is the stellar obliquity. A substantially misaligned orbit of Kepler-1704 b would warrant a reexamination of Kozai eccentricity oscillations, although planet–planet scattering can also cause misaligned orbits (e.g., Chatterjee et al. 2008; Naoz et al. 2012). Moreover, the effective temperature of Kepler-1704 makes this system a perfect laboratory for testing the theory that hot Jupiters preferentially realign cool stars (e.g., Schlaufman 2010; Winn et al. 2010). The effective temperature of Kepler-1704 is ${5745}_{-89}^{+88}$ K, which is well below the ∼6200 K Kraft break (Kraft 1967) that has been implicated by hot Jupiter obliquity observations. Since tidal forces are inefficient for failed hot Jupiters like Kepler-1704 b, we would expect that these planets would show a variety of obliquities and would not be preferentially aligned like hot Jupiters orbiting similarly cool stars.

In theory, the obliquity between Kepler-1704 b and its host star could be measured through the Rossiter–McLaughlin (RM) effect (McLaughlin 1924; Rossiter 1924). If successful, it would stand as the longest-period planet, by far, to have an obliquity measurement. In practice, an RM experiment will be challenging. Our SpecMatch analysis (Section 2.2) inferred a low stellar rotational velocity of 2.74 ± 1.0 km s⁻¹. By Equation (40) of Winn (2010), the maximum expected amplitude of the RM effect is only 11 m s⁻¹. Assuming 30 minute exposure times (as used for the current RV data), this would only allow 12 data points across the entire transit. With the ∼7 m s⁻¹ RV precision achieved using the best-match template (see Table 1 and Section 2.2), any detection of obliquity would likely be marginal. We recommend that any future effort to observe the spectroscopic transit of Kepler-1704 b should first acquire a high-S/N spectral template of Kepler-1704 to reduce the internal RV precision by several meters per second. Owing to the extreme eccentricity and the argument of periastron, the transit duration is short enough that a fortunately timed transit could be observed from a single site. For the Keck I telescope, only the second half of the 2023 transit (Table 4) will be visible. Around 14:30 UTC on 2023 March 12, Kepler-1704 will rise above the Nasmyth deck at Keck I at a favorable airmass of ∼1.5. Again assuming 30 minute exposure times, that would place roughly six data points across the second half of the RM signal. Even with the actual template and improved internal precision, a detection of obliquity would likely be moderate at best. Not until 2028 will the Keck I telescope have the optimal position for an RM detection. The mid-transit time of the 2028 August transit is almost perfectly timed with Kepler-1704 crossing the meridian, and the full transit (plus post-transit baseline) is observable. However, in the coming years, new precise RV facilities with the capability of achieving a few meters per second precision on faint (V = 13.4) stars, such as MAROON-X (Seifahrt et al. 2018) or the Keck Planet Finder (Gibson et al. 2016), should consider conducting RM measurements of long-period Kepler planets like Kepler-1704 b.

With only two transit events detected by Kepler, we cannot rule out large TTVs that could possibly preclude future transit observations (e.g., Wang et al. 2015a). The TTVs of this nature would require a massive perturber with an orbit that is sufficiently close to Kepler-1704 b to allow for gravitational interaction. The RV observations presented here largely rule out such a companion on orbits interior to Kepler-1704 b (Figure 9). However, giant planet companions on wider orbits could be present. One avenue of future work would be to apply the companion limits and stability results described here to a dynamical analysis of the Kepler-1704 system to set limits on TTV magnitudes.

Another avenue toward addressing this issue would be to observe a third transit. As has been done with Kepler-421 b (Dalba & Muirhead 2016), Kepler-167 e (Dalba & Tamburo 2019), and HIP 41378 f (Bryant et al. 2021), Kepler-1704 can be observed for a window of time surrounding the expected third transit according to a linear ephemeris. For this experiment, a missed transit places a lower limit on transit-to-TTVs. The ground-based observability of future Kepler-1704 b transits is not as restrictive as that described for JWST in Section 4.4.1. For example, in 2023 March, Kepler-1704 will rise to reasonable elevations in the last few hours of the night in the northern hemisphere. Any single site will likely struggle to detect the full transit. However, detections of ingress and egress from multiple sites spread out in longitude would refine the ephemeris of Kepler-1704 b and identify any TTVs.

Object Kepler-1704 represents an interesting comparison for the Kepler-167 system, in which an early K dwarf star hosts three inner super-Earth-sized planets and an outer transiting Jupiter analog on a P = 1071 day orbit (Kipping et al. 2016). Although the mass of Kepler-167 e—the outer giant planet—has not been measured, its orbital eccentricity has been constrained to ∼0.06 by the transit shape and duration. This low eccentricity combined with the presence of multiple inner super-Earth planets suggests that the migration mechanism for Kepler-167 e was likely gentle and driven by interactions with the disk. Kepler-167 is of solar metallicity, if not slightly metal-poor, so it is possible that Kepler-167 e was the only giant planet formed in the outer disk, so scattering events never occurred. Dalba & Tamburo (2019) ruled out the existence of TTVs in the ephemeris of Kepler-167 e, which further implied a lack of an outer massive companion. A mass and bulk metallicity measurement for Kepler-167 e would provide an interesting comparison with Kepler-1704 b, which likely experienced dynamical interactions with other bodies during and/or after its formation.

Giant transiting exoplanets with multiyear orbital periods are possibly exciting targets for dedicated exomoon searches (e.g., Kipping et al. 2012a; Heller et al. 2014; Teachey & Kipping 2018). Now that we have measured the mass and orbital properties of Kepler-1704 b, the plausibility of this planet hosting a system of exomoons should be investigated in more detail. Given the suspected active dynamical formation history of Kepler-1704 b, its ability to have maintained a system of exomoons is perhaps questionable. Indeed, the investigation of exomoon stability under tidal forces (e.g., Barnes & O'Brien 2002; Adams & Bloch 2016; Sucerquia et al. 2020), planet–planet scattering (e.g., Nesvorný et al. 2007; Gong et al. 2013; Hong et al. 2018), disk torques (e.g., Namouni 2010; Spalding et al. 2016), and secular migration owing to a stellar companion (e.g., Martinez et al. 2019; Trani et al. 2020) are active areas of theoretical research. Although any such study is beyond the scope of this work, we can approximate the Hill radius of Kepler-1704 b at periastron (where it is smallest):

$$\begin{eqnarray}&&{r}_{{\rm{H}},\mathrm{peri}}\approx a(1-e){\left(\displaystyle \frac{{M}_{p}}{3{M}_{\star }}\right)}^{1/3}.\end{eqnarray} \tag{ 6 }$$

For Kepler-1704 b, we find that r_H,peri ≈ 2.6 × 10⁶ km. Any exomoon on a circular prograde orbit around Kepler-1704 b would need a semimajor axis less than (roughly) half of this value to survive the close periastron passages (Hamilton & Burns 1991). For perspective, the semimajor axis of Callisto, Jupiter's most distant Galilean moon, is roughly 1.9 × 10⁶ km. Of course, this calculation neglects all of the processes that led to Kepler-1704 b reaching its current orbital configuration. We offer Kepler-1704 b as a potentially interesting case study for more detailed investigations of exomoon formation and stability in the future.

The fact that Kepler-1704 b swings through its host star's habitable zone on its eccentric orbit is also potentially interesting from an exomoon standpoint (e.g., Heller 2012; Heller & Barnes 2013; Hill et al. 2018). However, the plausibility of life developing on an exomoon that experiences such intense variation in stellar irradiation should be thoroughly scrutinized.

The vast majority of known giant planets on au-scale orbits have unknown radii because they either do not transit or are not known to transit. Without a radius, subsequent investigations of atmospheres and interiors are uncertain, if not altogether impossible. Measuring the masses of the modest sample of known transiting giant planets with au-scale orbits and discovering more such planets will be important to advancing our understanding of giant planet formation and migration. These discoveries will also drive new theoretical advances in giant planet interiors, which are needed given that changing model assumptions can substantially alter our conclusions about the interior structures of giant exoplanets (e.g., Müller et al. 2020a).

Only a handful of outer giant exoplanets like Kepler-1704 b exist within the Kepler sample, and they all orbit relatively faint stars. This creates two problems. First, their limited number means that unfortunate transit timing (see Table 4 and also Dalba & Tamburo 2019) can drastically slow progress to obtain new observations and advance our theoretical understanding. Second, their faintness must be overcome (if at all possible) by larger investments of highly competitive telescope time.

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), which is actively searching for transits of bright stars around the entire sky, presents solutions to both problems. The only drawback is the tendency of TESS's observing strategy to yield single transit events for most planets with orbital periods greater than a couple dozen days (e.g., Dalba et al. 2020b; Díaz et al. 2020; Gill et al. 2020; Lendl et al. 2020). If the Kepler mission had adopted the TESS mission's observing strategy, not only would Kepler-1704 b have been identified through a single transit, but its 6 hr transit duration could have easily been misconstrued as corresponding to a relatively short orbital period. This suggests (and more quantitative efforts have shown; Cooke et al. 2018; Villanueva et al. 2019) that given enough time and targets, TESS will identify transits from a unprecedented sample of long-period giant planets. Yet the advancement of giant planet theory and understanding will rely on continued challenging follow-up efforts to characterize these planet's masses, orbits, interiors, and atmospheres.