The Origin of Kepler-419b: A Path to Tidal Migration Via Four-body Secular Interactions

TL;DR

This paper investigates whether an undiscovered fourth body could explain the high-eccentricity orbit of Kepler-419b through secular interactions, using N-body simulations to explore parameter space and observational constraints.

Contribution

It demonstrates that a hidden fourth planet could induce high eccentricity in Kepler-419b without system destabilization, within specific mass and orbital parameters.

Findings

A fourth planet with certain parameters can excite Kepler-419b's eccentricity.

Such a planet could remain undetectable with current observational methods.

The proposed mechanism is consistent with system stability over 1 Gyr.

Abstract

We test the high-eccentricity tidal migration scenario for Kepler-419b, a member of the eccentric warm Jupiter class of planets whose origin is debated. Kepler-419 hosts two known planets (b,c). However, in its current configuration, planet c cannot excite the eccentricity of planet b enough to undergo high-eccentricity tidal migration. We investigate whether the presence of an undiscovered fourth body could explain the orbit of Kepler-419b. We explore the parameter space of this potential third giant planet using a suite of N-body simulations with a range of initial conditions. From the results of these simulations, coupled with observational constraints, we can rule out this mechanism for much of the parameter space of initial object d conditions. However, for a small range of parameters (masses between 0.5 and 7 $m_{\rm{Jup}}$, semi-major axes between 4 and 7.5 AU, eccentricities between 0.18 and 0.35, and mutual inclinations near 0$^{\circ}$) an undiscovered object d could periodically excite the eccentricity of Kepler-419b without destabilizing the system over 1 Gyr while producing currently undetectable radial velocity and transit timing variation signals.

Tables

Table 1. Table 1: Planet Parameters for Kepler-419b and Kepler-419c at Epoch BJD 2455809.4009671761741629 (Dawson et al., 2014 )

Parameter	Value
Stellar mass, $m_{\star }\left(m_{☉}\right)$	${1.39}_{-0.07}^{+0.08}$
Stellar radius, $R_{\star }\left(R_{☉}\right)$	${1.75}_{-0.07}^{+0.08}$
Mass b, $m_b\left(m_{\mathrm{Jup}}\right)$	$2.5\pm 0.3$
Period b, $P_b\left(\mathrm{days}\right)$	${69.7546}_{-0.0009}^{+0.0007}$
Semi-major axis b, $a_b\left(\mathrm{AU}\right)$	${0.370}_{-0.006}^{+0.007}$
Eccentricity b, $e_b$	$0.833\pm 0.013$
Argument of pericenter b, ${\omega }_b({}_{}{}^{\circ })$	${95.2}_{-1.2}^{+1.0}$
Mean anomaly b, $M_b({}_{}{}^{\circ })$	$68.69\pm 0.05$
Inclination b, $i_b({}_{}{}^{\circ })$	${88.95}_{-0.17}^{+0.14}$
Longitude of ascending node b, ${\Omega }_b({}_{}{}^{\circ })$	0 (fixed)
Longitude of pericenter b, ${\varpi }_b({}_{}{}^{\circ })$	${95.2}_{-1.2}^{+1.0}$
Mass c, $m_c\left(m_{\mathrm{Jup}}\right)$	$7.3\pm 0.4$
Period c, $P_c\left(\mathrm{days}\right)$	$675.47\pm 0.11$
Semi-major axis c, $a_c\left(\mathrm{AU}\right)$	$1.68\pm 0.03$
Eccentricity c, $e_c$	$0.184\pm 0.002$
Argument of pericenter c, ${\omega }_c({}_{}{}^{\circ })$	${275.3}_{-1.0}^{+1.2}$
Mean anomaly c, $M_c({}_{}{}^{\circ })$	$345.0\pm 0.3$
Inclination c, $i_c({}_{}{}^{\circ })$	${88}_{-2}^{+3}$
Longitude of ascending node c, ${\Omega }_c({}_{}{}^{\circ })$	$4\pm 12$
Longitude of pericenter c, ${\varpi }_c({}_{}{}^{\circ })$	$279\pm 12$
Mutual inclination, $i_{\mathrm{mut}}({}_{}{}^{\circ })$	$9_{-6}^{+8}$
${\omega }_b-{\omega }_c({}_{}{}^{\circ })$	$179.8\pm 0.6$
${\varpi }_b-{\varpi }_c({}_{}{}^{\circ })$	$176\pm 12$

Table 2. Table 2: Range of initial parameters explored for object d

Parameter	Range	Grid spacing
$m_d\left(m_{\mathrm{Jup}}\right)$	$0.5-115$	LOG10
$a_d\left(\mathrm{AU}\right)$	$1.0-10.0$	Linear
$e_d$	$0-0.45$	Linear
$i_{\mathrm{mut}}({}_{}{}^{\circ })$	$0-180$	Linear
${\omega }_d({}_{}{}^{\circ })$	$0-360$	Linear
${\Omega }_d({}_{}{}^{\circ })$	$0-360$	Linear

Equations

$$\frac{a_d(1-e_d)}{a_c(1+e_c)}>2.4[\max ({\mu }_c,{\mu }_d){]}^{1/3}{\left(\frac{a_d}{a_c}\right)}^{1/2}+1.5$$

$$a_d(1-e_d)$$

$$\left[1.7{\left(\frac{\max ({\mu }_d,0.007)}{0.01}\right)}^{1/3}{\left(\frac{a_d}{10\mathrm{AU}}\right)}^{1/2}+2.1\right]\mathrm{AU}$$

$$A_{bd}=\frac{1}{4}{n}_b\frac{m_d}{m_{\star }}\alpha \bar{\alpha }{b}_{3/2}^{(1)}(\alpha )⩾\frac{1}{10}{A}_{bc},$$

$$b_{3/2}^{(1)}(\alpha )=\frac{1}{\pi }{\int }_0^{2\pi }\frac{\cos \psi d\psi }{(1-2\alpha \cos \psi +{\alpha }^2{)}^{\frac{3}{2}}}.$$

$$\left(\frac{{\mu }_d}{0.01}\right){\left(\frac{a_d}{10\mathrm{AU}}\right)}^{-2}⩾{\left(\frac{b_{3/2}^{(1)}(\alpha )}{1.19}\right)}^{-1}.$$

$$\min [a_b(1-e_b^2)]<0.1$$

Full text

planets and satellites: dynamical evolution and stability — planets and satellites: individual (Kepler-419)

1 Introduction

Both hot and warm Jupiters are believed to form beyond the ice-line and migrate inward to their current semi-major axes (Bodenheimer et al. 2000; Rafikov 2006; see Dawson & Johnson 2018 for a review of hot Jupiter origins theories). Although disk migration (Goldreich & Tremaine, 1980; Ward, 1997; Baruteau et al., 2014) could deliver hot and warm Jupiters, it is difficult to reconcile disk migration with warm Jupiters’ eccentricity distribution (Figure 1) because planet-disk interactions tend to damp eccentricities (Bitsch et al., 2013; Dunhill et al., 2013). Planet-disk interactions can sometimes excite eccentricities, but they typically saturate at a random velocity equal to the sound speed ($e\lesssim 0.03$ for a 100 day orbit; Duffell & Chiang 2015). Eccentricity growth could occur after the disk migration stage via scattering, but such growth is limited by $v_{\rm{escape}}/v_{\rm{keplerian}}$, making it difficult for planets on close in orbits to attain large eccentricities (Goldreich et al., 2004; Ida et al., 2013; Petrovich et al., 2014). The disconnect between formation theories and observations makes warm Jupiters an interesting parameter space for testing migration mechanisms.

Eccentric warm Jupiters may have arrived through high-eccentricity tidal migration (e.g., Hut 1981; Wu & Murray 2003). Petrovich & Tremaine (2016) concluded that high-eccentricity tidal migration triggered by an outer planetary companion can account for $\sim 20\%$ of warm Jupiters and most warm Jupiters with $e\geq 0.4$. In this mechanism, planets are excited to large orbital eccentricities such that tidal friction with their host star works to circularize the planet’s orbit and reduce the semi-major axis. The planet’s angular momentum remains constant once it is decoupled from other perturbing bodies. Via high-eccentricity migration, many warm Jupiters would end their lives as hot Jupiters on circular orbits. Several theories have been proposed to explain the original excitation of eccentricity, including planet-planet scattering (Rasio & Ford, 1996), secular chaos (Wu & Lithwick, 2011), and stellar flybys (e.g., Kaib et al. 2013). These mechanisms are successful in producing a wide distribution of giant planet eccentricities from multi-planet systems with initially circular orbits (e.g., Jurić & Tremaine 2008; Chatterjee et al. 2008).

Here we focus on secular interactions between planets and their ability to induce high-eccentricity tidal migration. While some planets, such as HD 80606 b (Moutou et al., 2009) and HD 17156 b (Barbieri et al., 2009), have very large eccentricities and may currently be undergoing tidal migration en route to becoming a hot Jupiter, the intermediate eccentricity warm Jupiters would all need smaller angular momenta (or, equivalently, larger eccentricities at a given semi-major axis) than we observe to be undergoing tidal migration. One effect that could induce this migration is for a perturbing body to periodically excite the planet’s eccentricity through secular interactions (e.g., Dong et al. 2014) leading to migration by tidal friction during the high-eccentricity phase. We would be observing the planet during the low eccentricity portion of its eccentricity oscillation cycle. There is some evidence for this mechanism in the trend that eccentric warm Jupiters are more likely to have companions (Dong et al., 2014; Bryan et al., 2016) than circular warm Jupiters. The Kozai-Lidov mechanism (Wu & Murray, 2003; Naoz et al., 2011) is an extreme case of this type of secular eccentricity modulation, but typically requires large mutual inclinations (see Li et al. 2014 for the case of an eccentric, coplanar Kozai perturber).

We focus on an interesting system for testing theories of tidal migration, Kepler-419. Kepler-419 is a $1.39^{+0.08}_{-0.07}M_{\sun}$ star with two known planetary companions, Kepler-419b and Kepler-419c. Kepler-419b is a transiting warm Jupiter at semi-major axis $a=0.370^{+0.007}_{-0.006}$ AU with a large eccentricity of $e=0.833\pm 0.013$ measured via the “photoeccentric effect" (Dawson & Johnson, 2012) and confirmed via radial velocity (RV) measurements (Dawson et al., 2014; Almenara et al., 2018). Kepler-419c is a $m=7.3\pm 0.4M_{Jup}$ perturbing body with moderate eccentricity ($e=0.184\pm 0.002$) located at $a=1.68\pm 0.03$ AU.

Dawson et al. (2014) used RV measurements from Keck HIRES to measure the mass of planet b ($m=2.5\pm 0.3M_{Jup}$) and improve precision on the eccentricity value. They used transit timing variations (TTVs) of planet b to deduce precise physical and orbital parameters of planet c. The linear trend in the RV curve is consistent with the properties of planet c derived from the TTVs. Limits on transit duration variations provided no additional constraints on planet c’s properties. The full set of system parameters calculated by simultaneously fitting the RVs and TTVs are presented in Table 1, an adaptation of Table 4 in Dawson et al. (2014). Almenara et al. (2018) used RV measurements from SOPHIE (Perruchot et al., 2008; Bouchy et al., 2009) to independently confirm the presence of Kepler-419c. Saad-Olivera et al. (2018) uniquely constrained the masses of both planets using the TTV data alone.

At least one other system is consistent with the eccentricity modulated tidal migration scenario: Kepler-693, which hosts an eccentric warm Jupiter and, unlike Kepler-419c, a companion on a high mutual inclination orbit that induces large Kozai-Lidov eccentricity oscillations (Masuda, 2017). However, Dawson et al. (2015) found a paucity of high-eccentricity warm Jupiters compared with that expected if hot Jupiters obtain their low-$a$ orbits through tidal friction with their host star (Socrates et al., 2012b), suggesting high-eccentricity tidal migration might not be a common migration channel.

If high-eccentricity tidal migration can be ruled out for Kepler-419b, it could further support that this mechanism is not ubiquitous among the warm Jupiter population. Kepler-419b’s present day eccentricity is too low for tidal migration and Dawson et al. (2014) found that Kepler-419c cannot drive up the eccentricity high enough to induce tidal migration. Secular interactions with planet c do cause oscillations in the eccentricity of planet b, but it is currently in the high-$e$ phase of this cycle and never reaches $a(1-e^{2})<0.1$ (represented by the black dotted line in Figure 1 and the green shaded area in Figure 2, the time evolution of Kepler-419b in the two-planet solution). However, before ruling out tidal migration, we need to investigate the possibility that an undetected third planet or a brown dwarf in the system could allow for high-eccentricity migration by increasing the maximum value of Kepler-419b’s eccentricity.

We explore the parameter space for an unseen additional perturber in the Kepler-419 system that could excite the eccentricity of planet b enough for it to periodically reach the large values necessary to undergo tidal migration. We start with the two known planets in their present configuration and simulate the system forward in time. In the tidal migration scenario, the eccentricity and semi-major axis of planet b would decrease over time, while the parameters of planet c remain roughly constant. Any instabilities in the four-body system today would have been stronger in the past when the planets were more tightly packed, allowing us to assume that system configurations that go unstable in our simulations would not have been stable in the past. Similarly, if a set of system parameters fails to sufficiently excite the eccentricity of planet b in the present configuration, we can be confident that it did not do so in the past. Strong precession due to general relativity, tides, or rotational quadrupolar bulges would be needed to decouple planet b from secular interactions with the perturber (e.g., Wu & Murray 2003). With Kepler-419b’s current orbital configuration, none of these are strong enough to do so (Dawson et al., 2014), nor would they have been if Kepler-419b’s semi-major axis was larger in the past.

We test the feasibility of high-eccentricity tidal migration through a suite of $N$-body simulations of the Kepler-419 system with various initial conditions for a possible object d. We motivate our analysis and outline our approach in Section 2. We present the results of our simulations in Section 3. In Section 4, we assess how object d would affect the observables of the system. We summarize our conclusions, discuss alternative scenarios, and propose future efforts in Section 5.

2 Overview of Simulations

We seek to test the parameter space of a potential object d in this system for a set of initial conditions that could significantly contribute to the maximum eccentricity of planet b without destabilizing the system. We constrain our explored parameter space to the range of masses, semi-major axes, and eccentricities that could reasonably produce this effect.

To guide the construction of our parameter space, we use two approximate expressions based on the two requirements that guide our analysis: stability and secular perturbing strength. For stability, we require that any third planet or binary companion in the system must not disturb the orbits of the two known planets over the lifetime of the system via close encounters, collisions, or ejections. We also require that the secular interactions between planet b and a theoretical object d are strong enough to periodically boost the eccentricity of planet b such that it enters the tidal circularization regime. The two expressions we use have opposite demands and thus help to constrain our parameter bounds: stability requires a less massive and more distant object while secular perturbing strength requires a more massive and closer in object.

The first expression we use (Equation 1) is an analytic approximation of stability, developed for two-planet systems in Petrovich (2015). Systems in which the inequality is satisfied are expected to evolve secularly without close encounters. Since our system hosts three satellites, this criterion is not directly applicable; however, if the inequality is not satisfied, we can expect the system to be unstable. If objects c and d would not be stable in a two-planet system, adding in planet b is unlikely to improve the stability.

[TABLE]

We can rewrite this criterion by plugging in the known values for Kepler-419c as

[TABLE]

where $\mu_{c}=m_{c}/m_{\star}$ and $\mu_{d}=m_{d}/m_{\star}$. Parameter sets that satisfy the inequality in Equation 2 are considered favorable for stability. The mass, semi-major axis, and eccentricity are limited by this constraint.

The second expression we use (Equation 3) is an approximation of the angular frequency of precession of a test particle induced by a body with some mass, assuming low eccentricities and low inclinations (Murray & Dermott, 1999). We use this expression to roughly approximate the secular effect of object d on planet b, despite our system breaking the test particle and low eccentricity assumptions and, in some iterations, the low inclination assumption. We can compare this to the empirical precession rate of planet b due to planet c measured from an N-body simulation of these two known planets (Figure 2). In order for the secular effect of object d to be competitive, we require that the frequency of planet b’s pericenter precession induced by object d must be greater than 1/10 that by planet c. The expression for this requirement is:

[TABLE]

where $A_{bd}$ is the angular frequency of precession of planet b due to planet d, $A_{bc}$ is the angular frequency of precession of planet b due to planet c, $n_{b}$ is the mean motion of planet b, $m_{\star}$ is the stellar mass, $\alpha=a_{b}/a_{d}$, $\bar{\alpha}=\alpha$ (for an exterior perturber), and $b^{(1)}_{3/2}(\alpha)$ is a Laplace coefficient of order unity,

[TABLE]

Substituting $A_{bc}=1.34$ kyr*-1* measured from the two-planet simulation (Fig. 2), as well as the mass and semi-major axis of planet b, into Equation 3 yields

[TABLE]

Equation 5 favors larger masses and smaller semi-major axes. Equation 2, however, favors larger semi-major axes and lower masses. Intermediate masses and semi-major axes are needed to satisfy both constraints. The range of properties for object d we include after considering these two constraints is provided in Table 2. Most notably, they are Jupiter,brown dwarf, and late M dwarf-mass objects with semi-major axes less than 10 AU and eccentricities less than 0.45. Due to their approximate nature, we do not use either of these expressions as hard constraints on the parameter space we explore; however, we show in Section 3 that they are useful guidelines because object d parameters that do not satisfy Equation 2 tend to lead to unstable systems and object d parameters that do not satisfy Equation 5 tend to not excite the eccentricity of planet b high enough.

Although approximate analytic expressions for 4-body secular evolution, such as Laplace-Lagrange excitation, exist, no existing approximation is appropriate here due to the combination of large eccentricities and small semi-major axis ratios. We could use these expressions as rough estimates of the long-term evolution of the system, but they are not reliable enough to determine which systems would be successful and which would not. Thus, we must numerically simulate the system to determine long-term stability and quantify the effect of object d on the eccentricity of planet b.

We run 3495 N-body simulations using the Mercury6 Bulirsch-Stoer integrator (Chambers, 1999) with the range of initial conditions for object d defined in Table 2. We search for stable systems with secular eccentricity oscillations in planet b that periodically result in a large enough eccentricity to tidally migrate.

3 Results

As discussed in the previous section, in order for a potential Kepler-419d to have contributed to the tidal migration of planet b, two conditions must be satisfied: (1) the amplitude of secular eccentricity oscillations in planet b induced by object d must be strong enough to significantly add to those induced by planet c so that it can get to a high enough eccentricity, and (2) the system must remain stable at least over the stellar lifetime. These two constraints compete with each other in that (1) favors a high mass, low semi-major axis object, while (2) favors a low mass planet widely separated from the other two planets in the system. If there is an area of parameter space that satisfies both constraints, we expect it to be in an intermediate mass and semi-major axis space where these constraints overlap.

We assess the results of our simulations by considering the satisfaction of these two conditions for each of the parameter sets for object d we test. In the following two subsections, we provide explanations of how we determine whether these conditions are satisfied.

3.1 Secular Excitation of Planet b

Following Socrates et al. (2012a) and Dong et al. (2014), we require that planet b must satisfy the inequality in Equation 6 in order for tidal migration to have occurred in this system. This inequality expresses the maximum orbital angular momentum that can result in a hot Jupiter ($e\approx 0$, $a\leqslant 0.1$ AU).

[TABLE]

Figure 1 shows a tidal circularization track (blacked dashed line) for a planet at the edge of this limit. After migration, a planet following this track would have eccentricity $e=0$ and period $P=10$ days. Observationally, hot Jupiter eccentricities tend to circularize around 0.06 AU (Socrates et al., 2012a). Since angular momentum is conserved by tidal circularization and tidal dissipation depends strongly on the planet-star separation, the $a(1-e^{2})<0.1$ cutoff shown here is a conservative limit (Dong et al., 2014). In reality, the cutoff for tidal migration on a timescale less than the age of the system could be much stricter; however, because tidal parameters for Jupiter-like planets such as Kepler-419b are not well known, we use the more conservative cutoff of requiring $a_{\rm{final}}=a(1-e^{2})$ for planet b to be in the hot Jupiter regime. Object d parameter sets that produce lower $\min[a_{b}(1-e_{b}^{2})]$ values are more likely to result in the migration of planet b within the age of the system

We use 1 Myr simulations to determine which properties of object d allow planet b to reach a sufficient eccentricity to satisfy Equation 6, given that the semi-major axis of planet b remains roughly constant. In order for planet b to undergo tidal migration in its current configuration, it must satisfy Equation 6 at some point within its secular cycle. Since, empirically, 1 Myr is enough time to encompass many secular timescales for all of our sets of object d initial conditions, we can rule out any solutions that do not satisfy our criterion in that time. Figure 3 shows a representative example of one simulation that satisfies Equation 6 over the first 1 Myr (bottom panel) and one that does not (middle panel). The top panel is discussed in Section 3.2.

For comparison of the effect on the eccentricity of planet b, Figure 2 shows $a_{b}(1-e_{b}^{2})$, simulated over 10 Gyr with only the known second planet included. Only the first $2*10^{5}$ years are plotted to show the timescale of the secular cycle. Although the planet’s eccentricity is oscillating, it does not satisfy Equation 6 in the two-planet solution. Our degree of confidence in that conclusion is very high because the orbital parameters for the two known planets are very precisely constrained (Table 1).

The first million years of our simulations partially constrain the allowed parameter space of object d, as demonstrated by the left column of Figure 4. The orange circles in this figure represent initial parameters of object d that do not sufficiently excite the eccentricity of the inner planet (Equation 6 is not satisfied), while blue circles represent successful initial parameter sets. Each panel in Figure 4 shows a projection of the many-dimensional space of object d parameters into two dimensions, which accounts for the overlapping points. We can rule out much of the low mass, large semi-major axis regime with the 1 Myr simulations since that region is dominated by orange circles and devoid of blue circles. Figure 4 also includes the results of 200 follow-up simulations that are discussed in Section 3.3, which accounts for the varying point density. The initial run of 3295 simulations were roughly uniformly sampled following Table 2.

3.2 System Stability

Because we are seeing Kepler-419b more than a Gyr after its formation (the age of the system is not well constrained beyond that, Dawson & Johnson 2012), the system configuration must be stable over the lifetime of the star. Thus, any simulations that go unstable before $\sim$1 Gyr can be ruled out. Unstable simulations are denoted by gray crosses in Figure 4 and an example is shown in the top panel of Figure 3.

In our simulations, a system configuration is deemed "unstable" if (a) a planet is ejected, (b) a planet collides with the central star, or (c) two planets experience a close encounter, defined as one object moving within a Hill radius of another object. Figure 5 provides a visual representation of how each unstable configuration went unstable in the parameter space of object d initial conditions. Green triangles represent systems where object d is ejected, purple stars represent systems where the inner planet collides with the star, and blue diamonds and red squares represent collisions between planets b/c and c/d, respectively. The different ways in which configurations go unstable are somewhat stratified in parameter space, with collisions between planets b and c tending to occur in systems with a high mass, large semi-major axis, small eccentricity object d; collisions between planets c and d tending to occur in systems where object d has a small semi-major axis; and collisions between planet b and the star populating the full parameter space except for the low mass, small semi-major axis object d regime. Planet ejections are rare in our simulations and mostly occur at intermediate semi-major axes and low masses for object d.

3.3 Simulation Results after 1 Gyr

After 1 Myr, the majority of object d parameters could be ruled out because they either went unstable (Section 3.2) or did not sufficiently excite the eccentricity of planet b (Section 3.1). The left column of Figure 4 shows the results of these simulations plotted in the parameter space of object d initial conditions. As stated, gray crosses represent unstable systems, small orange circles represent stable systems that do not satisfy Equation 6, and large blue circles represent systems that remain stable and satisfy Equation 6.

Parameters for object d that did not drive up planet b’s eccentricity (orange circles) were not continued beyond 1 Myr. Thus, in the right column of Figure 4 (snapshot at 1 Gyr), the only solutions that are plotted are ones that went unstable after 1 Myr but before 1 Gyr (gray crosses) and ones that satisfy both criteria through 1 Gyr (blue circles). Of the 3295 different sets of parameters for object d we test, 12 satisfy Equation 6 while remaining stable for 1 Gyr. These surviving systems are in low mutual inclination orbits with eccentricities between 0.175 and 0.35 and semi-major axes between 4 and 8 AU. Notably, more than half of the successful simulations are in violation of Equation 5 (upper cluster of blue circles in the bottom right panel of Figure 3). This approximation likely fails here due to the broken assumptions discussed in Section 2 (e.g., massless particle, low eccentricities). While Equation 5 was ultimately not used to strictly constrain our parameter space, it did provide a useful baseline when constructing the space.

Of the 12 sets of initial conditions of object d that survive beyond 1 Gyr and satisfy Equation 6, we can define two different populations: one low mass ($5\times 10^{-4}-7\times 10^{-3}m_{\sun}$) planetary companion group and one high mass ($5.5\times 10^{-2}-1.15\times 10^{-1}m_{\sun}$) stellar companion group with masses corresponding to brown dwarfs or late M dwarfs. We zoom in on these two areas of parameter space and run 100 simulations in each range with finer resolution in mass, semi-major axis, and eccentricity. These additional 200 simulations have already been included in Figure 4.

In the planetary-mass group, including both the original parameter sets and the 100 finer-resolution parameter sets, 16 simulations remain stable for 1 Gyr while satisfying Equation 6. In the stellar-mass group, 7 parameter sets are successful through 1 Gyr, bringing the total to 23 successful object d parameter sets.

Since there are areas of parameter space for a third planet in the Kepler-419 system that contribute strongly to the eccentricity oscillations of planet b without disturbing the stability of the system, we cannot rule out the high-eccentricity tidal migration mechanism for planet b’s origin based on these two criteria.

3.4 Chaotic System Evolution

We can divide our stable simulation results into two categories based on the nature of their secular evolution: (1) parameter sets that produce smooth secular oscillations, such as those in the bottom two panels of Figure 3, and (2) parameter sets that evolve chaotically. Figure 6 shows some representative examples of category (2). After 1 Myr of simulation runtime, 47% of our successful simulations (parameter sets that remain stable and satisfy Equation 6) qualitatively fall into the first category, while the remaining 53% fall into the second. However, after 1 Gyr, all of the surviving parameter sets are chaotic in their evolution. parameter sets that are not chaotic after 1 Myr either go unstable or are disrupted and become chaotic.

In all successful chaotic simulations, the orbits of planets b and c were not significantly disturbed from their measured values, suggesting the Kepler-419 system could have evolved in a chaotic manner and remain consistent with observations. In fact, some of the chaotically evolving systems reach a much lower minimum $a_{b}(1-e_{b}^{2})$ than the secularly evolving systems.

Chaotic systems that only passed below $a_{b}(1-e_{b}^{2})=0.1$ for a small amount of time are treated as systems that do not satisfy Equation 6 (orange circles in Figure 4) because an assumption for that criterion was that planet b would reach its maximum eccentricity periodically throughout the stellar lifetime. Thus we require the $a_{b}(1-e_{b}^{2})$ to periodically pass below 0.1 for at least half of the simulation runtime. This additional cut is already included in the number of successful simulations we quote in Section 3.3.

3.5 Longitude of Pericenter Constraints

The observed difference in longitude of pericenter $\varpi_{c}-\varpi_{b}$ ($\Delta\varpi$) between planet b and planet c is very close to $180^{\circ}$ (see Table 1). We are unlikely to have observed the orbits of these two planets in this configuration by chance, which suggests they are locked in a libration around $\Delta\varpi=180^{\circ}$ with a small amplitude. By simulating the two-planet case forward in time, we see that this angle does in fact librate about $180^{\circ}$ with an amplitude of $\sim 18^{\circ}$ in the best fit two-planet case and $\sim 16^{\circ}$ in the coplanar two-planet case. Any additional perturber must avoid disrupting this libration over the stellar lifetime for us to observe it today.

In Figure 7 we show the $\Delta\varpi$ for six representative three-planet configurations that survive longer than 1 Gyr. The systems in the left column have a planetary-mass object d, while the systems in the right column have a stellar-mass object d. Only the planetary-mass perturbers preserve the libration of $\Delta\varpi$ between b and c while satisfying the other conditions necessary for high-eccentricity tidal migration; for all stellar-mass perturbers, the angle circulates. The libration in two of the planetary-mass perturber systems gets disturbed before 1 Gyr, leaving 13 sets of object d parameters that preserve libration of $\Delta\varpi$ between planets b and c over the simulation runtime.

Almenara et al. (2018) suggest that any additional perturber in the system, would disrupt the libration of $\Delta\varpi$ between the two planets and dismiss four-body secular eccentricity excitation. However, we show that for planetary-mass objects, the $\Delta\varpi$ libration is preserved. Moreover, Since we only observe these planets in a snapshot in time, we cannot be certain that $\Delta\varpi$ is librating in the observed system. Thus, we cannot definitively rule out any areas of parameter space for a potential object d. However, because $\Delta\varpi$ circulates in all of the stellar-mass object d systems and we observe $\Delta\varpi$ very near $180^{\circ}$, we suggest that planetary-mass object d systems may provide a better explanation for the observations.

4 Observational Constraints

We can rule out most of the parameter space for object d by applying the two criteria described in Section 3 (Equation 6 and stability for 1 Gyr) to our simulations. However, we can also compare the solutions that did satisfy those criteria with current RV (Section 4.1) and TTV (Section 4.2) measurements of the system to constrain the parameter space even further. As additional observations are published, the allowed parameter space we present here will likely shrink as RV and TTV signals due to a long period planet are detected or ruled out.

4.1 Radial Velocities

We model the RV signal of the host star due to object d for each of the 23 sets of initial conditions that satisfy Equation 6 and survive for 1 Gyr. In Figure 8, we plot a subset of these RV curves as well as the Keck HIRES (Dawson et al., 2014) and SOPHIE (Almenara et al., 2018) observations of the host star with the best fit two-planet signal subtracted out. The error on these observations includes 40 ms*-1* of stellar jitter (Dawson et al., 2014). The residuals to the two-planet solution are flat, with a root mean square (RMS) spread of $\sim 47$ms*-1*. Since the eccentricity perturbations in planet b due to object d are secular, the mean anomaly $M_{d}$ is not constrained by our simulations. We know the shape and orientation of the orbit, but not the planet’s current position along that orbit. Thus, we shift $M_{d}$ such that the observations fall on the flattest part of the RV curve, minimizing the RV change ($\Delta v$). We plot The minimum-amplitude stellar signal (blue) as well as the minimum-amplitude (red) and maximum-amplitude (orange) planetary signals as dashed lines on Figure 8.

For the planetary-mass group of object d parameter sets, the RV signals due to object d are too small to be detectable with the current data. A small fraction of this group has a large enough RV semi-amplitude, $K\gtrsim 47$ms*-1*, to potentially be detected through continued RV monitoring with, for example, the SOPHIE spectrograph. A detailed quantitative analysis of the detectability of these potential planets is beyond the scope of this paper, but Figure 8 allows us to qualitatively compare the signal to the RMS scatter in the data.

The stellar-mass group produces much larger RV signals. In fact, each of the 7 stellar-mass parameter sets produces a signal with amplitude significantly larger than the RMS scatter of the observations. Thus, a stellar-mass object d that sufficiently excites the eccentricity of planet b cannot be present in the Kepler-419 system.

4.2 Transit Timing Variations

We independently model the variation in the time of transit for planet b due to object d for each of the 23 sets of initial conditions that satisfy Equation 6 and survive for 1 Gyr. In Figure 9, we plot this TTV signal for several representative examples of object d that encompass the spread of solutions in mass for both planetary and stellar objects. Here we show the TTV signal using the mean anomaly $M_{d}$ that minimizes TTV amplitude (the same approach used to minimize the RV amplitude in Section 4.1). We qualitatively interpret an amplitude larger than the error bars to be a detectable signal.

In the left column of Figure 9, we show that all of the planetary-mass solutions fail to produce detectable TTVs. In the right column, we show that the TTV amplitude for each of the stellar-mass solutions is larger than the spread of the error bars, independently ruling out stellar-mass objects.

We note that the amplitude of the TTV signal is strongly dependent on the mean anomaly $M_{d}$. Had we used the $M_{d}$ that maximized TTV amplitude, nearly all of the planetary-mass signals could be ruled out. To demonstrate this, the blue triangles in Figure 9 show one solution with an arbitrary $M_{d}$. Thus, while we cannot use our TTV analysis alone to rule out any additional solutions, some fine-tuning was required to reach this conclusion. In fact, for the more massive planetary perturbers, about half of the possible values for $M_{d}$ produce detectable signals, allowing us to rule out configurations with mean anomalies $M_{d}<150^{\circ}$ and $M_{d}>310^{\circ}$. For the less massive perturbers ($m_{d}<\sim 3m_{\rm{Jup}}$), the TTVs would not be detectable for any value of $M_{d}$

5 Summary and Discussion

High eccentricity tidal migration has been invoked to explain the origins and properties of hot and warm Jupiters, giant planets orbiting close to their host stars. Although Kepler-419b does not currently have a large enough eccentricity for tidal friction to circularize its orbit, if a perturber were to periodically drive up its eccentricity, it could periodically reach a small enough separation from the star for tidal forces to become important, leading to its current small semi-major axis orbit. Kepler-419c, a coplanar cold Jupiter with a moderate eccentricity, is perturbing the eccentricity of planet b, but not enough to cause high-eccentricity tidal migration on a timescale comparable to the age of the system. However, it is possible for an additional perturbing planet (or small star) to periodically boost the eccentricity of the warm Jupiter into the tidal migration regime.

Using a suite of $\sim 3500$ N-body integrations, we explored the parameter space for a potential fourth body in the system that could be exciting the eccentricity of planet b high enough to induce tidal circularization of its orbit without destabilizing the system over the stellar lifetime. We ran these integrations for 1 Gyr and narrowed down the allowed parameter space for this additional perturber to a small region. We also produced RV and TTV curves for all sets of object d parameters that survived for 1 Gyr and compared them to the observations. From the results and analysis of our simulations, we conclude the following:

Successful sets of initial object d parameters (cases where the system remains stable for 1 Gyr and Kepler-419b reached $a(1-e^{2})<0.1$) are constrained to semi-major axes between 4 and 8 AU, eccentricities between 0.175 and 0.35, and mutual inclinations below a few degrees. 2. 2.

The simulated systems with object d initial conditions that survive for 1 Gyr while meeting the eccentricity criterion can be divided into two groups by mass: a planetary-mass group covering $\sim 0.5-7$ $m_{\rm{Jup}}$ and a stellar-mass group covering $\sim 50-115$ $m_{\rm{Jup}}$. 3. 3.

Many of the simulations that survived to 100 Myrs and all of the simulations that survived to 1 Gyr produced $a_{b}(1-e_{b}^{2})<0.1$ values that evolved chaotically (not in smooth, consistent secular cycles). The orbits of planets b and c were not significantly disturbed in these cases, so they remain consistent with the observations. 4. 4.

The measured difference in longitude of pericenter between Kepler-419b and Kepler-419c, $\varpi_{c}-\varpi_{b}$ ($\Delta\varpi$), is very near 180*∘, suggesting the angle is likely librating about this value with a small amplitude (Dawson et al., 2014). In the two-planet solution, this angle does librate. All of the stellar-mass group of successful parameter sets for object d failed to preserve this libration for the length of the simulation, along with 2 of the 16 parameter sets from the planetary-mass group. There is a small chance our measurement of $\Delta\varpi$ near 180∘* (Dawson et al., 2014) is just a coincidence; however, our results suggest that planetary mass perturbers may provide a better explanation for the observations. 5. 5.

The RV and TTV signals for all of the planetary-mass group proved to be consistent with the observations. However, the stellar-mass group could be ruled out because an object in that group would have been detected in both the RV and TTV observations. With the current RV precision and stellar jitter, a subset of the planetary-mass group could potentially be ruled out by extending the RV baseline to the present day if significant linear trends in the RVs are not seen. For the more massive planetary perturbers, TTVs constrain the mean anomaly to the range $150^{\circ}<M_{d}<310^{\circ}$. 6. 6.

Pending additional RV measurements, we are left with a potential planet between 0.5 and 7 Jupiter masses on a 4-7.5 AU orbit with an eccentricity between 0.18 and 0.35. If a planet in this range of parameter space were to exist in this system, it could produce eccentricity oscillations in Kepler-419b strong enough to periodically pass very close to the host star, allowing tidal friction and, therefore, inward migration.

As additional RV measurements of Kepler-419 are taken, we will be able to further constrain the parameter space allowed for a third planet that could perturb the eccentricity of Kepler-419b enough to be inducing tidal circularization of its orbit. The RV semi-amplitudes of possible solutions for object d are between 5 and 85 ms*-1*. Kepler-419 is an active F-star, making precise RV measurements difficult for this system; however, a longer time baseline combined with detailed stellar noise modeling might put additional constraints on the mass and semi-major axis of any additional hidden planets. Similarly, it may be possible to constrain the allowed parameter space for an additional perturbing planet using future TTV measurements of the system, although they would have to significantly improve upon Kepler’s precision.

Ruling out a four-body tidal migration scenario for Kepler-419b with future observations could lend support to new recent alternative hypotheses for the origin of its eccentric, close-in orbit. Petrovich et al. (2019) posited that Kepler-419b could have acquired its large eccentricity through adiabatic transport during the disk clearing stage. This scenario requires a massive disk and in situ formation of the warm Jupiter, but could feasibly explain the observed properties of the system. Almenara et al. (2018) suggest two alternative scenarios that could increase the eccentricity of planet b: (1) spin-orbit coupling (Correia et al., 2012) and (2) collision with another planet in a mean-motion resonance. These two scenarios could occur in situ or following disk migration. Although we present stable solutions for a four-body system in this paper, planets b and c are very near a stability boundary. Thus, planet b could only have a slightly larger semi-major axis in the past and could not have migrated far while remaining in a stable configuration. Antonini et al. (2016) claim that most observed Jupiter pairs are dynamically unstable if the inner planets are placed on orbits with semi-major axes larger than 1 AU. This result could suggest that Kepler-419b was and other warm Jupiters with companions were formed in situ.

We gratefully acknowledge support from grant NNX16AB50G awarded by the NASA Exoplanets Research Program and the Alfred P. Sloan Foundation’s Sloan Research Fellowship. This research has also made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. We thank Andrew Shannon for helpful discussions and Cristobal Petrovich, Eric Ford, Jason Wright, Caryl Gronwall, and Jim Kasting for helpful comments. We also thank the anonymous referee for their helpful suggestions.