Eccentric Companions to Kepler-448b and Kepler-693b: Clues to the Formation of Warm Jupiters

TL;DR

This paper reports the discovery of non-transiting companions to two warm Jupiters, Kepler-448b and Kepler-693b, using transit timing and duration variations, providing insights into their formation and migration processes.

Contribution

It presents the first detection of massive, eccentric, non-transiting companions to warm Jupiters via dynamical modeling of TTVs and TDVs, suggesting possible in situ formation or high-eccentricity migration.

Findings

Companions have masses of 22 and 150 Jupiter masses.

Both companions are on eccentric orbits with high eccentricities.

Mutual inclinations suggest potential eccentricity oscillations affecting migration.

Abstract

I report the discovery of non-transiting close companions to two transiting warm Jupiters (WJs), Kepler-448/KOI-12b (orbital period $P=17.9\,\mathrm{days}$, radius $R_{\rm p}=1.23^{+0.06}_{-0.05}\,R_{\rm Jup}$) and Kepler-693/KOI-824b ($P=15.4\,\mathrm{days}$, $R_{\rm p}=0.91\pm0.05\,R_{\rm Jup}$), via dynamical modeling of their transit timing and duration variations (TTVs and TDVs). The companions have masses of $22^{+7}_{-5}\,M_{\rm Jup}$ (Kepler-448c) and $150^{+60}_{-40}\,M_{\rm Jup}$ (Kepler-693c), and both are on eccentric orbits ($e=0.65^{+0.13}_{-0.09}$ for Kepler-448c and $e=0.47^{+0.11}_{-0.06}$ for Kepler-693c) with periastron distances of $1.5\,\mathrm{au}$. Moderate eccentricities are detected for the inner orbits as well ($e=0.34^{+0.08}_{-0.07}$ for Kepler-448b and $e=0.2^{+0.2}_{-0.1}$ for Kepler-693b). In the Kepler-693 system, a large mutual inclination between the inner and outer orbits ($53^{+7}_{-9}\,\mathrm{deg}$ or $134^{+11}_{-10}\,\mathrm{deg}$) is also revealed by the TDVs. This is likely to induce a secular oscillation of the inner WJ's eccentricity that brings its periastron close enough to the host star for tidal star-planet interactions to be significant. In the Kepler-448 system, the mutual inclination is weakly constrained and such an eccentricity oscillation is possible for a fraction of the solutions. Thus these WJs may be undergoing tidal migration to become hot Jupiters (HJs), although the migration via this process from beyond the snow line is disfavored by the close-in and massive nature of the companions. This may indicate that WJs can be formed in situ and could even evolve into HJs via high-eccentricity migration inside the snow line.

Tables

Table 1. Table 1: Transit Times, Durations, and Radius Ratios of Kepler-448b

Transit number	Transit time (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)	Transit duration (day)	Planet-to-star radius ratio
$-1$	$128.7420\pm 0.0001$	$0.2808\pm 0.0003$	$0.0920\pm 0.0004$
$\cdots$

Table 2. Table 2: Transit Times, Durations, and Radius Ratios of Kepler-693b

Transit number	Transit time (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)	Transit duration (day)	Planet-to-star radius ratio
$0$	$173.609\pm 0.002$	$0.115\pm 0.007$	$0.115\pm 0.008$
$\cdots$

Table 3. Table 3: Transit Parameters of Kepler-448b Based on Short-cadence Data

Parameter	Value
(Parameters from the Mean Transit Light Curve)aaMedian and $68\%$ credible interval of the MCMC posteriors from the mean transit light curve. Here the circular orbit is assumed to relate mean stellar density to the semi-major axis over stellar radius, $a/R_{\star }$. The limb-darkening coefficients $u_1$ and $u_2$ are converted from $q_1$ and $q_2$.
$q_1={\left(u_1+u_2\right)}^2$	$0.199\left(7\right)$
$q_2=\frac{u_1}{2\left(u_1+u_2\right)}$	$0.39\left(2\right)$
Center of the mean transit (day)	$-0.00001\left(2\right)$
$R_p/R_{\star }$	$0.09044\left(7\right)$
Normalization	$0.999999\left(3\right)$
Transit impact parameter	$0.373\left(6\right)$
Mean stellar density ($g{\mathrm{cm}}^{-3}$)	$0.393\left(3\right)$
Photometric jitter	$0.000036\left(2\right)$
${\left(a/R_{\star }\right)}_{e=0}$	$18.77\left(4\right)$
$u_1$	$0.348\left(9\right)$
$u_2$	$0.10\left(2\right)$
(Mean Orbital Period and Transit Epoch)bbObtained by linearly fitting the series of transit times in Table 1. Errors are scaled by $\sqrt{{\chi }^2/d.o.f.}$.
$t_0$ (${\mathrm{BJD}}_{\mathrm{TDB}}$)	$2454979.5961\left(2\right)$
$P$ (day)	$17.855234\left(4\right)$

Table 4. Table 4: Transit Parameters of Kepler-693b Based on Long-cadence Data

Parameter	Value
(Parameters from the Mean Transit Light Curve)aaMedian and $68\%$ credible interval of the MCMC posteriors from the mean transit light curve. Here the circular orbit is assumed to relate mean stellar density to the semi-major axis over stellar radius, $a/R_{\star }$. The limb-darkening coefficients $u_1$ and $u_2$ are converted from $q_1$ and $q_2$.
$q_1={\left(u_1+u_2\right)}^2$	${0.6}_{-0.2}^{+0.3}$
$q_2=\frac{u_1}{2\left(u_1+u_2\right)}$	${0.4}_{-0.2}^{+0.3}$
Center of the mean transit (day)	$0.0000\left(3\right)$
$R_p/R_{\star }$	${0.117}_{\left(-6\right)}^{\left(+4\right)}$
Normalization	$1.00000\left(8\right)$
Transit impact parameter	${0.63}_{-0.15}^{+0.06}$
Mean stellar density ($g{\mathrm{cm}}^{-3}$)	${3.0}_{-0.5}^{+1.0}$
Photometric jitter	$0.00074\left(7\right)$
${\left(a/R_{\star }\right)}_{e=0}$	${34}_{-2}^{+4}$
$u_1$	${0.7}_{-0.3}^{+0.2}$
$u_2$	${0.1}_{-0.4}^{+0.3}$
(Mean Orbital Period and Transit Epoch)bbObtained by linearly fitting the series of transit times in Table 2. Errors are scaled by $\sqrt{{\chi }^2/d.o.f.}$.
$t_0$ (${\mathrm{BJD}}_{\mathrm{TDB}}$)	$2455006.613\left(1\right)$
$P$ (day)	$15.37566\left(3\right)$

Table 5. Table 5: Parameters of the Kepler-448 System from the Dynamical TTV and TDV Analysis

Parameter	Solution 1	Solution 2	Combined	Prior
Fitted Parameters
(Inner Orbit)
1. Time of inferior conjunction	${128.7418}_{\left(-3\right)}^{\left(+2\right)}$	${128.7413}_{\left(-2\right)}^{\left(+2\right)}$	${128.7415}_{\left(-3\right)}^{\left(+4\right)}$	$𝒰(128.73,128.75)$
$t_{\mathrm{ic},1}$ (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)
2. Orbital period $P_1$ (day)	${17.855183}_{\left(-6\right)}^{\left(+9\right)}$	${17.855181}_{\left(-5\right)}^{\left(+7\right)}$	${17.855182}_{\left(-6\right)}^{\left(+8\right)}$	${𝒰}_{\log }(17.8551,17.8553)$
3. Orbital eccentricity $e_1$	${0.34}_{-0.06}^{+0.07}$	${0.35}_{-0.08}^{+0.09}$	${0.34}_{-0.07}^{+0.08}$	$𝒰(0,0.7)$
4. Argument of periastron ${\omega }_1$ (deg)	$-{49}_{-21}^{+9}$	$-{131}_{-9}^{+18}$	$\cdots$	$𝒰(-180,180)$
5. Cosine of orbital inclination $\cos i_1$aaThere also exists a solution with negative $\cos i_1$. The solution is statistically indistinguishable from the one reported here, except that the signs of $\cos i_2$ and ${\Omega }_{21}$ are flipped (and thus $i_{21}$ remains the same; see Eqn. A15). In principle, the TTVs for the solutions with $\cos i_1>0$ and are not completely identical (see, e.g., A15 of Borkovits et al., 2015). However, the difference is in practice negligibly small for a transiting system because the effect is proportional to $\cot i_1$. I confirmed that this is indeed the case by performing the same numerical fit with the prior on $\cos i_1$ replaced by $𝒰(-0.02,0)$.	$0.013\left(1\right)$	$0.013\left(1\right)$	$0.013\left(1\right)$	$𝒰(0,0.02)$aaThere also exists a solution with negative $\cos i_1$. The solution is statistically indistinguishable from the one reported here, except that the signs of $\cos i_2$ and ${\Omega }_{21}$ are flipped (and thus $i_{21}$ remains the same; see Eqn. A15). In principle, the TTVs for the solutions with $\cos i_1>0$ and are not completely identical (see, e.g., A15 of Borkovits et al., 2015). However, the difference is in practice negligibly small for a transiting system because the effect is proportional to $\cot i_1$. I confirmed that this is indeed the case by performing the same numerical fit with the prior on $\cos i_1$ replaced by $𝒰(-0.02,0)$.
(Outer Orbit)
6. Time of the periastron passage	${1076}_{-11}^{+10}$	${1080}_{-11}^{+10}$	${1078}_{-11}^{+10}$	$𝒰(900,1200)$
${\tau }_2$ (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)
7. Periastron distance over	${9.5}_{-0.3}^{+0.4}$	${9.4}_{-0.3}^{+0.4}$	${9.4}_{-0.3}^{+0.4}$	${𝒰}_{\log }(5,20)$
inner semi-major axis $q_2/a_1$
8. Orbital eccentricity $e_2$	${0.61}_{-0.07}^{+0.11}$	${0.69}_{-0.10}^{+0.12}$	${0.65}_{-0.09}^{+0.13}$	$𝒰(0,0.95)$
9. Argument of periastron ${\omega }_2$ (deg)	$-{89}_{-7}^{+6}$	$-{89}_{-6}^{+5}$	$-{89}_{-6}^{+5}$	$𝒰(-180,180)$
10. Cosine of orbital inclination $\cos i_2$aaThere also exists a solution with negative $\cos i_1$. The solution is statistically indistinguishable from the one reported here, except that the signs of $\cos i_2$ and ${\Omega }_{21}$ are flipped (and thus $i_{21}$ remains the same; see Eqn. A15). In principle, the TTVs for the solutions with $\cos i_1>0$ and are not completely identical (see, e.g., A15 of Borkovits et al., 2015). However, the difference is in practice negligibly small for a transiting system because the effect is proportional to $\cot i_1$. I confirmed that this is indeed the case by performing the same numerical fit with the prior on $\cos i_1$ replaced by $𝒰(-0.02,0)$.	$-{0.3}_{-0.3}^{+0.3}$	$-{0.6}_{-0.2}^{+0.3}$	$-{0.5}_{-0.3}^{+0.4}$	$𝒰(-1,1)$
11. Relative longitude of	$3_{-4}^{+4}$	$-{178}_{-4}^{+5}$	$\cdots$	$𝒰(-180,180)$
ascending node ${\Omega }_{21}$ (deg)a,ba,bfootnotemark:
(Physical Properties)
12. Mass of Kepler-448 $m_{\star }$ ($M_{\odot }$)	$1.5\pm 0.1$	$1.5\pm 0.1$	$1.5\pm 0.1$	$𝒢(1.51,0.09,0.14)$
13. Mean density of	${0.79}_{-0.09}^{+0.09}$	${0.81}_{-0.09}^{+0.10}$	${0.80}_{-0.09}^{+0.10}$	$𝒢(0.76,0.11,0.12)$
Kepler-448 ${\rho }_{\star }$ ($g{\mathrm{cm}}^{-3}$)
14. Mass of Kepler-448b $m_b$ ($M_{\mathrm{Jup}}$)ccThis parameter is not constrained by the data; posterior is identical to the log-uniform prior.	${1.1}_{-0.8}^{+3.5}$	${1.1}_{-0.8}^{+3.5}$	${1.1}_{-0.8}^{+3.5}$	${𝒰}_{\log }(0.1,10)$
15. Mass of Kepler-448c $m_c$ ($M_{\mathrm{Jup}}$)	${21}_{-5}^{+6}$	${24}_{-5}^{+6}$	${22}_{-5}^{+7}$	${𝒰}_{\log }(0.001M_{\odot },0.1M_{\odot })$
(Jitters)
16. Transit time jitter ${\sigma }_{\mathrm{TTV}}$ (${10}^{-4}\mathrm{day}$)	$2.2\pm 0.3$	$2.2\pm 0.3$	$2.2\pm 0.3$	${𝒰}_{\log }(5\times {10}^{-2},5)$
17. Transit duration jitter ${\sigma }_{\mathrm{dur}}$ (${10}^{-4}\mathrm{day}$)	$3.8\pm 0.6$	$3.7\pm 0.6$	$3.7\pm 0.6$	${𝒰}_{\log }(0.1,10)$
Derived Parameters
Outer orbital period $P_2$ (day)	$\left({2.2}_{-0.5}^{+1.5}\right)\times {10}^3$	$\left({2.9}_{-1.0}^{+3.3}\right)\times {10}^3$	$\left({2.5}_{-0.7}^{+2.4}\right)\times {10}^3$	$\cdots$
Inner semi-major axis $a_1$ (au)	${0.154}_{\left(-3\right)}^{\left(+4\right)}$	${0.154}_{\left(-3\right)}^{\left(+4\right)}$	${0.154}_{\left(-3\right)}^{\left(+4\right)}$	$\cdots$
Outer semi-major axis $a_2$ (au)	${3.8}_{-0.6}^{+1.6}$	$5_{-1}^{+3}$	${4.2}_{-0.9}^{+2.4}$	$\cdots$
Periastron distance of	${1.47}_{-0.06}^{+0.07}$	${1.45}_{-0.06}^{+0.06}$	${1.45}_{-0.06}^{+0.07}$	$\cdots$
the outer orbit $a_2\left(1-e_2\right)$ (au)
Mutual orbital inclination $i_{21}$ (deg)	${20}_{-12}^{+17}$	${146}_{-16}^{+19}$	$\cdots$	$\cdots$
Physical radius of Kepler-448 ($R_{\odot }$)	$1.40\pm 0.06$	${1.39}_{-0.06}^{+0.07}$	${1.40}_{-0.06}^{+0.07}$	$\cdots$
Physical radius of Kepler-448b ($R_{\mathrm{Jup}}$)ddDerived from the posterior of $R_{\star }$ and mean and standard deviation of $R_p/R_{\star }$ from individual transits (Table 3).	${1.24}_{-0.05}^{+0.06}$	${1.23}_{-0.06}^{+0.06}$	${1.23}_{-0.05}^{+0.06}$	$\cdots$
Transit impact parameter of Kepler-448b	$0.372\left(6\right)$	$0.371\left(6\right)$	$0.371\left(6\right)$	$\cdots$
Occultation impact parameter of Kepler-448b	$0.22\pm 0.02$	$0.21\pm 0.02$	$0.22\pm 0.02$	$\cdots$
Log evidence $\ln 𝒵$ from Multinest	$-168.16\pm 0.09$	$-168.04\pm 0.09$	$-167.40\pm 0.09$

Table 6. Table 6: Parameters of the Kepler-693 System from the Dynamical TTV and TDV Analysis

Parameter	Solution 1	Solution 2	Combined	Prior
Fitted Parameters
(Inner Orbit)
1. Time of inferior conjunction	${173.611}_{\left(-1\right)}^{\left(+1\right)}$	${173.610}_{\left(-1\right)}^{\left(+2\right)}$	${173.610}_{\left(-1\right)}^{\left(+1\right)}$	$𝒰(173.58,173.64)$
$t_{\mathrm{ic},1}$ (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)
2. Orbital period $P_1$ (day)	${15.37541}_{\left(-7\right)}^{\left(+10\right)}$	${15.37534}_{\left(-8\right)}^{\left(+10\right)}$	${15.37537}_{\left(-9\right)}^{\left(+10\right)}$	${𝒰}_{\log }(15.375,15.376)$
3. Orbital eccentricity $e_1$	${0.14}_{-0.04}^{+0.08}$	${0.3}_{-0.1}^{+0.2}$	${0.2}_{-0.1}^{+0.2}$	$𝒰(0,0.7)$
4. Argument of periastron ${\omega }_1$ (deg)	${41}_{-28}^{+55}$	${174}_{-39}^{+13}$	$\cdots$	$𝒰(-180,180)$
5. Cosine of orbital inclination $\cos i_1$aaThere also exists a solution with negative $\cos i_1$. The solution is statistically indistinguishable from the one reported here, except that the signs of $\cos i_2$ and ${\Omega }_{21}$ are flipped (and thus $i_{21}$ remains the same; see Eqn. A15). In principle, the TTVs for the solutions with $\cos i_1>0$ and are not completely identical (see, e.g., A15 of Borkovits et al., 2015). However, the difference is in practice negligibly small for a transiting system because the effect is proportional to $\cot i_1$. I confirmed that this is indeed the case by performing the same numerical fit with the prior on $\cos i_1$ replaced by $𝒰(-0.04,0)$.	$0.021\left(2\right)$	$0.022\left(3\right)$	$0.022\left(3\right)$	$𝒰(0,0.04)$aaThere also exists a solution with negative $\cos i_1$. The solution is statistically indistinguishable from the one reported here, except that the signs of $\cos i_2$ and ${\Omega }_{21}$ are flipped (and thus $i_{21}$ remains the same; see Eqn. A15). In principle, the TTVs for the solutions with $\cos i_1>0$ and are not completely identical (see, e.g., A15 of Borkovits et al., 2015). However, the difference is in practice negligibly small for a transiting system because the effect is proportional to $\cot i_1$. I confirmed that this is indeed the case by performing the same numerical fit with the prior on $\cos i_1$ replaced by $𝒰(-0.04,0)$.
(Outer Orbit)
6. Time of the periastron passage	${640}_{-17}^{+20}$	${640}_{-18}^{+22}$	${640}_{-17}^{+22}$	$𝒰(500,900)$
${\tau }_2$ (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)
7. Periastron distance over	${13}_{-1}^{+2}$	${13}_{-2}^{+1}$	${13}_{-2}^{+2}$	${𝒰}_{\log }(6,25)$
inner semi-major axis $q_2/a_1$
8. Orbital eccentricity $e_2$	${0.48}_{-0.06}^{+0.12}$	${0.46}_{-0.05}^{+0.10}$	${0.47}_{-0.06}^{+0.11}$	$𝒰(0,0.95)$
9. Argument of periastron ${\omega }_2$ (deg)	${41}_{-20}^{+20}$	${25}_{-21}^{+25}$	${30}_{-23}^{+24}$	$𝒰(-180,180)$
10. Cosine of orbital inclination $\cos i_2$aaThere also exists a solution with negative $\cos i_1$. The solution is statistically indistinguishable from the one reported here, except that the signs of $\cos i_2$ and ${\Omega }_{21}$ are flipped (and thus $i_{21}$ remains the same; see Eqn. A15). In principle, the TTVs for the solutions with $\cos i_1>0$ and are not completely identical (see, e.g., A15 of Borkovits et al., 2015). However, the difference is in practice negligibly small for a transiting system because the effect is proportional to $\cot i_1$. I confirmed that this is indeed the case by performing the same numerical fit with the prior on $\cos i_1$ replaced by $𝒰(-0.04,0)$.	$-{0.2}_{-0.2}^{+0.3}$	$-{0.3}_{-0.1}^{+0.2}$	$-{0.3}_{-0.2}^{+0.2}$	$𝒰(-1,1)$
11. Relative longitude of	${51}_{-10}^{+8}$	$-{138}_{-11}^{+12}$	$\cdots$	$𝒰(-180,180)$
ascending node ${\Omega }_{21}$ (deg)a,ba,bfootnotemark:
(Physical Properties)
12. Mass of Kepler-693 $m_{\star }$ ($M_{\odot }$)	${0.80}_{-0.03}^{+0.03}$	${0.80}_{-0.03}^{+0.04}$	${0.80}_{-0.03}^{+0.04}$	$𝒢(0.79,0.03,0.15)$
13. Mean density of	${2.2}_{-0.2}^{+0.3}$	${2.2}_{-0.3}^{+0.3}$	${2.2}_{-0.3}^{+0.3}$	$𝒢(1.93,0.18,0.53)$
Kepler-693 ${\rho }_{\star }$ ($g{\mathrm{cm}}^{-3}$)
14. Mass of Kepler-693b $m_b$ ($M_{\mathrm{Jup}}$)ccThis parameter is not constrained by the data; posterior is identical to the log-uniform prior.	${0.8}_{-0.6}^{+2.7}$	${1.0}_{-0.7}^{+3.2}$	${0.9}_{-0.7}^{+3.0}$	${𝒰}_{\log }(0.1,10)$
15. Mass of Kepler-693c $m_c$ ($M_{\mathrm{Jup}}$)	${167}_{-43}^{+59}$	${136}_{-34}^{+50}$	${145}_{-37}^{+58}$	${𝒰}_{\log }(0.001M_{\odot },0.3M_{\odot })$
(Jitters)
16. Transit time jitter ${\sigma }_{\mathrm{TTV}}$ (${10}^{-4}\mathrm{day}$)	$2_{-1}^{+4}$	$1_{-1}^{+4}$	$1_{-1}^{+4}$	${𝒰}_{\log }(5\times {10}^{-2},50)$
17. Transit duration jitter ${\sigma }_{\mathrm{dur}}$ (${10}^{-4}\mathrm{day}$)	$9_{-8}^{+14}$	$9_{-8}^{+13}$	$9_{-8}^{+14}$	${𝒰}_{\log }(0.1,{10}^2)$
Derived Parameters
Outer orbital period $P_2$ (day)	$\left({1.8}_{-0.4}^{+1.0}\right)\times {10}^3$	$\left({1.8}_{-0.3}^{+0.6}\right)\times {10}^3$	$\left({1.8}_{-0.3}^{+0.8}\right)\times {10}^3$	$\cdots$
Inner semi-major axis $a_1$ (au)	${0.112}_{\left(-1\right)}^{\left(+2\right)}$	${0.112}_{\left(-2\right)}^{\left(+2\right)}$	${0.112}_{\left(-1\right)}^{\left(+2\right)}$	$\cdots$
Outer semi-major axis $a_2$ (au)	${2.9}_{-0.5}^{+1.1}$	${2.8}_{-0.4}^{+0.7}$	${2.8}_{-0.4}^{+0.8}$	$\cdots$
Periastron distance of	$1.5\pm 0.2$	$1.5\pm 0.2$	$1.5\pm 0.2$	$\cdots$
the outer orbit $a_2\left(1-e_2\right)$ (au)
Mutual orbital inclination $i_{21}$ (deg)	${53}_{-9}^{+7}$	${134}_{-10}^{+11}$	$\cdots$	$\cdots$
Physical radius of Kepler-693 ($R_{\odot }$)	$0.80\pm 0.03$	$0.80\pm 0.04$	$0.80\pm 0.04$	$\cdots$
Physical radius of Kepler-693b ($R_{\mathrm{Jup}}$)ddDerived from the posterior of $R_{\star }$ and that of $R_p/R_{\star }$ from the mean transit light curve (Table 4).	$0.91\pm 0.05$	$0.91\pm 0.05$	$0.91\pm 0.05$	$\cdots$
Transit impact parameter of Kepler-693b	$0.58\pm 0.04$	$0.59\pm 0.05$	$0.59\pm 0.05$	$\cdots$
Occultation impact parameter of Kepler-693b	${0.68}_{-0.06}^{+0.06}$	${0.62}_{-0.08}^{+0.09}$	${0.64}_{-0.09}^{+0.08}$	$\cdots$
Log evidence $\ln 𝒵$ from Multinest	$-90.30\pm 0.08$	$-89.57\pm 0.08$	$-89.18\pm 0.08$

Table 7. Table 7: cos ⁡ n 1 , 2 subscript 𝑛 1 2 \cos n_{1,2} for sin ⁡ Ω 21 → 0 → subscript Ω 21 0 \sin\Omega_{21}\to 0

	$(\cos n_1,\cos n_2)$	$(n_1,n_2)$
$\cos {\Omega }_{21}\to 1$	$\mathrm{sgn}\left(\sin {\Omega }_{21}\right)(\frac{\sin \left(i_2-i_1\right)}{\sin i_{21}},\frac{\sin \left(i_2-i_1\right)}{\sin i_{21}})$	$(0,0)$ or $(\pi ,\pi )$ for $i_2-i_1\to 0$
$\cos {\Omega }_{21}\to -1$	$\mathrm{sgn}\left(\sin {\Omega }_{21}\right)(-\frac{\sin \left(i_2+i_1\right)}{\sin i_{21}},\frac{\sin \left(i_2+i_1\right)}{\sin i_{21}})$	$(0,\pi )$ or $(\pi ,0)$ for $i_2+i_1\to \pi$

Table 8. Table 8: Parameters of the Kepler-448 System from the Analytical and Numerical TTV Analyses

	Analytic		Numerical
Parameter	Posterior	Prior	Posterior	Prior
Fitted Parameters
(Inner Orbit)
1. Time of inferior conjunction	${128.7414}_{\left(-3\right)}^{\left(+7\right)}$	$𝒰(128.7,128.8)$	${128.7418}_{\left(-2\right)}^{\left(+2\right)}$	$𝒰(128.73,128.75)$
$t_{\mathrm{ic},1}$ (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)
2. Orbital period $P_1$ (day)	${17.855219}_{\left(-13\right)}^{\left(+6\right)}$	${𝒰}_{\log }(17.85,17.86)$	${17.855179}_{\left(-6\right)}^{\left(+7\right)}$	${𝒰}_{\log }(17.8551,17.8553)$
3. Orbital eccentricity $e_1$	${0.5}_{-0.2}^{+0.3}$	$𝒰(0,0.95)$	${0.5}_{-0.2}^{+0.1}$	$𝒰(0,0.7)$
4. Argument of periastron ${\omega }_1$ (deg)	$-{43}_{-104}^{+158}$	$𝒰(-180,180)$	$-{15}_{-140}^{+166}$	$𝒰(-180,180)$
5. Cosine of orbital inclination $\cos i_1$	$0$ (fixed)	$\cdots$	$0$ (fixed)	$\cdots$
(Outer Orbit)
6. Time of the periastron passage	${1076}_{-10}^{+10}$	$𝒰(120,1600)$	${1072}_{-10}^{+12}$	$𝒰(900,1200)$
${\tau }_2$ (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)
7. Periastron distance over	${9.3}_{-0.3}^{+0.4}$	${𝒰}_{\log }(1,50)$	${9.4}_{-0.4}^{+0.4}$	${𝒰}_{\log }(5,20)$
inner semi-major axis $q_2/a_1$
8. Orbital eccentricity $e_2$	${0.54}_{-0.05}^{+0.10}$	$𝒰(0,0.95)$	${0.60}_{-0.07}^{+0.12}$	$𝒰(0,0.95)$
9. Argument of periastron ${\omega }_2$ (deg)	$-{91}_{-29}^{+52}$	$𝒰(-180,180)$	$-{116}_{-34}^{+58}$	$𝒰(-180,180)$
10. Cosine of orbital inclination $\cos i_2$	${0.1}_{-0.8}^{+0.7}$	$𝒰(-1,1)$	${0.1}_{-0.7}^{+0.6}$	$𝒰(-1,1)$
11. Relative longitude of	$-5_{-155}^{+162}$	$𝒰(-180,180)$	$1_{-131}^{+135}$	$𝒰(-180,180)$
ascending node ${\Omega }_{21}$ (deg)aaReferenced to the ascending node of the inner orbit, whose direction is arbitrary.
(Physical Properties)
12. Mass of Kepler-448 $m_{\star }$ ($M_{\odot }$)	${1.54}_{-0.09}^{+0.11}$	$𝒢(1.51,0.09,0.14)$	${1.54}_{-0.09}^{+0.11}$	$𝒢(1.51,0.09,0.14)$
13. Mass of Kepler-448b $m_b$ ($M_{\mathrm{Jup}}$)	$1$ (fixed)	$\cdots$	$1$ (fixed)	$\cdots$
14. Mass of Kepler-448c $m_c$ ($M_{\mathrm{Jup}}$)	${12}_{-4}^{+10}$	${𝒰}_{\log }({10}^{-4}{M}_{\odot },0.3M_{\odot })$	${15}_{-4}^{+7}$	${𝒰}_{\log }({10}^{-3}{M}_{\odot },0.1M_{\odot })$
(Jitters)
15. Transit time jitter ${\sigma }_{\mathrm{TTV}}$ (${10}^{-4}\mathrm{day}$)	$2.2\pm 0.3$	${𝒰}_{\log }(0.05,5)$	$2.2\pm 0.3$	${𝒰}_{\log }(0.05,5)$
Derived Parameters
Outer orbital period $P_2$ (day)	$\left({1.6}_{-0.2}^{+0.8}\right)\times {10}^3$	$\cdots$	$\left({2.0}_{-0.4}^{+1.4}\right)\times {10}^3$	$\cdots$
Inner semi-major axis $a_1$ (au)	$0.154\left(3\right)$	$\cdots$	$0.154\left(3\right)$	$\cdots$
Outer semi-major axis $a_2$ (au)	${3.1}_{-0.3}^{+0.9}$	$\cdots$	${3.6}_{-0.5}^{+1.5}$	$\cdots$
Periastron distance of	${1.44}_{-0.06}^{+0.07}$	$\cdots$	${1.45}_{-0.06}^{+0.07}$	$\cdots$
the outer orbit $a_2\left(1-e_2\right)$ (au)
Mutual orbital inclination $i_{21}$ (deg)	${118}_{-68}^{+31}$	$\cdots$	${54}_{-28}^{+81}$	$\cdots$

Table 9. Table 9: Parameters of the Kepler-693 System from the Analytical and Numerical TTV Analyses

	Analytic		Numerical
Parameter	Posterior	Prior	Posterior	Prior
Fitted Parameters
(Inner Orbit)
1. Time of inferior conjunction	${173.611}_{\left(-1\right)}^{\left(+2\right)}$	$𝒰(173.55,173.65)$	${173.609}_{\left(-2\right)}^{\left(+1\right)}$	$𝒰(173.58,173.64)$
$t_{\mathrm{ic},1}$ (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)
2. Orbital period $P_1$ (day)	${15.37565}_{\left(-7\right)}^{\left(+3\right)}$	${𝒰}_{\log }(15.37,15.38)$	${15.37549}_{\left(-7\right)}^{\left(+11\right)}$	${𝒰}_{\log }(15.375,15.376)$
3. Orbital eccentricity $e_1$	$0.6\pm 0.2$	$𝒰(0,0.95)$	$0.5\pm 0.1$	$𝒰(0,0.7)$
4. Argument of periastron ${\omega }_1$ (deg)	${64}_{-158}^{+72}$	$𝒰(-180,180)$	${84}_{-93}^{+62}$	$𝒰(-180,180)$
5. Cosine of orbital inclination $\cos i_1$	$0$ (fixed)	$\cdots$	$0$ (fixed)	$\cdots$
(Outer Orbit)
6. Time of the periastron passage	${636}_{-16}^{+28}$	$𝒰(170,1560)$	${660}_{-20}^{+23}$	$𝒰(500,900)$
${\tau }_2$ (${\mathrm{BJD}}_{\mathrm{TDB}}-2454833$)
7. Periastron distance over	${13}_{-2}^{+2}$	${𝒰}_{\log }(1,50)$	${15}_{-2}^{+2}$	${𝒰}_{\log }(6,25)$
inner semi-major axis $q_2/a_1$
8. Orbital eccentricity $e_2$	${0.50}_{-0.08}^{+0.18}$	$𝒰(0,0.95)$	${0.7}_{-0.1}^{+0.1}$	$𝒰(0,0.95)$
9. Argument of periastron ${\omega }_2$ (deg)	${60}_{-62}^{+58}$	$𝒰(-180,180)$	${86}_{-47}^{+42}$	$𝒰(-180,180)$
10. Cosine of orbital inclination $\cos i_2$	$0.0\pm 0.3$	$𝒰(-1,1)$	$-0.1\pm 0.4$	$𝒰(-1,1)$
11. Relative longitude of	$-6_{-157}^{+167}$	$𝒰(-180,180)$	$-{21}_{-139}^{+178}$	$𝒰(-180,180)$
ascending node ${\Omega }_{21}$ (deg)aaReferenced to the ascending node of the inner orbit, whose direction is arbitrary.
(Physical Properties)
12. Mass of Kepler-693 $m_{\star }$ ($M_{\odot }$)	$0.80\pm 0.03$	$𝒢(0.79,0.03,0.05)$	$0.80\pm 0.03$	$𝒢(0.79,0.03,0.05)$
13. Mass of Kepler-693b $m_b$ ($M_{\mathrm{Jup}}$)	$1$ (fixed)	$\cdots$	$1$ (fixed)	$\cdots$
14. Mass of Kepler-693c $m_c$ ($M_{\mathrm{Jup}}$)	${57}_{-21}^{+41}$	${𝒰}_{\log }({10}^{-4}{M}_{\odot },0.3M_{\odot })$	${94}_{-23}^{+37}$	${𝒰}_{\log }({10}^{-3}{M}_{\odot },0.3M_{\odot })$
(Jitters)
15. Transit time jitter ${\sigma }_{\mathrm{TTV}}$ (${10}^{-4}\mathrm{day}$)	$6_{-3}^{+4}$	${𝒰}_{\log }(0.05,50)$	$6_{-3}^{+4}$	${𝒰}_{\log }(0.05,50)$
Derived Parameters
Outer orbital period $P_2$ (day)	$\left({2.0}_{-0.6}^{+2.5}\right)\times {10}^3$	$\cdots$	$\left({5.3}_{-2.4}^{+8.0}\right)\times {10}^3$	$\cdots$
Inner semi-major axis $a_1$ (au)	${0.112}_{\left(-1\right)}^{\left(+2\right)}$	$\cdots$	${0.112}_{\left(-1\right)}^{\left(+2\right)}$	$\cdots$
Outer semi-major axis $a_2$ (au)	${3.0}_{-0.6}^{+2.1}$	$\cdots$	${5.8}_{-1.9}^{+4.9}$	$\cdots$
Periastron distance of	$1.5\pm 0.2$	$\cdots$	$1.6\pm 0.2$	$\cdots$
the outer orbit $a_2\left(1-e_2\right)$ (au)
Mutual orbital inclination $i_{21}$ (deg)	${133}_{-115}^{+30}$	$\cdots$	${133}_{-105}^{+30}$	$\cdots$

Equations

$$p(\mbox{\boldmath$\theta$}|\mbox{\boldmath$d$},I)=\frac{1}{\mathcal{Z}}\,p(\mbox{\boldmath$d$}|\mbox{\boldmath$\theta$},I)\,p(\mbox{\boldmath$\theta$}|I),$$

$$\mathcal{Z}\equiv\int p(\mbox{\boldmath$d$}|\mbox{\boldmath$\theta$},I)\,p(\mbox{\boldmath$\theta$}|I)\mathrm{d}\mbox{\boldmath$\theta$}$$

$${\mathcal{L}}_{\mathrm{TTV}}\equiv \prod _i\frac{1}{\sqrt{2\pi ({\sigma }_{t,i}^2+{\sigma }_{\mathrm{TTV}}^2)}}\exp \left[-\frac{(t_i-t_i^{\mathrm{model}}{)}^2}{2({\sigma }_{t,i}^2+{\sigma }_{\mathrm{TTV}}^2)}\right],$$

$$\mathcal{L}=$$

$$\times \prod _j\frac{1}{\sqrt{2\pi ({\sigma }_{T,j}^2+{\sigma }_{\mathrm{dur}}^2)}}\exp \left[-\frac{(T_j-T_j^{\mathrm{model}}{)}^2}{2({\sigma }_{T,j}^2+{\sigma }_{\mathrm{dur}}^2)}\right]$$

$$\times \frac{1}{\sqrt{2\pi {\sigma }_b^2}}\exp \left[-\frac{(b_{\mathrm{mean}}-b_{\mathrm{mean}}^{\mathrm{model}}{)}^2}{2{\sigma }_b^2}\right],$$

$$|\dot{T}|\le 3J_2{\left(\frac{a}{R_{\star }}\right)}^{-2}\frac{b}{\sqrt{1-b^2}}|\sin \lambda |.$$

$$\Delta T\lesssim 2\times 10^{-3}\mathrm{hr}\left(\frac{J_2}{10^{-4}}\right),$$

$$\dot{T}\le 4e{\left(\frac{a}{R_{\star }}\right)}^{-1}\sqrt{1-b^2}\times \frac{3}{1-e^2}{\left(\frac{na}{c}\right)}^2$$

$$\Delta T\lesssim 2\times 10^{-3}\mathrm{hr}\left(\frac{e}{1-e^2}\right)$$

$$\Delta t$$

$$\simeq \frac{P_2}{2}\left(1-\frac{4}{\pi }{e}_2\right)$$

$$\simeq \frac{P_1}{2}{\left[\frac{a_2(1-e_2{)}^{8/3\pi }}{a_1}\right]}^{3/2}\sim \frac{P_1}{2}{\left(\frac{q_2}{a_1}\right)}^{3/2},$$

$$\Delta A\simeq \frac{P_1}{2\pi }{A}_{L1}\simeq P_1\frac{15}{16\pi }{\left[\frac{q_2(1+e_2)}{a_1}\right]}^{-3/2}\frac{m_{\mathrm{c}}}{m_{\star }},$$

$${\Delta }_{\mathrm{LTTE}}=-\frac{a_2\sin i_2}{c}\frac{m_{\mathrm{c}}}{m_{\star }+m_{\mathrm{b}}+m_{\mathrm{c}}}\frac{(1-e_2^2)\sin u_2}{1+e_2\cos v_2},$$

$${\Delta }_1=\frac{P_1}{2\pi }{A}_{L1}(1-e_1^2{)}^{1/2}\left[{\delta }_{\mathrm{tidal}}+{\delta }_{ecc1}+{\delta }_{ecc2}+{\delta }_{\mathrm{noncopl}}\right].$$

$$A_{L1}=\frac{15}{8}\frac{m_{\mathrm{c}}}{m_{\star }+m_{\mathrm{b}}+m_{\mathrm{c}}}\frac{P_1}{P_2}(1-e_2^2{)}^{-3/2},$$

$${\delta }_{\mathrm{tidal}}=\left[\frac{8}{15}{f}_1+\frac{4}{5}{K}_1\right]\mathcal{M},$$

$${\delta }_{ecc1}=(1+\cos i_{21})\left\{K_{11}\mathcal{S}\left[2u_2-2(n_2-n_1)\right]-K_{12}\mathcal{C}\left[2u_2-2(n_2-n_1)\right]\right\},$$

$${\delta }_{ecc2}=(1-\cos i_{21})\left\{K_{11}\mathcal{S}\left[2u_2-2(n_2+n_1)\right]+K_{12}\mathcal{C}\left[2u_2-2(n_2+n_1)\right]\right\},$$

$${\delta }_{\mathrm{noncopl}}={\sin }^2{i}_{21}\left(K_{11}\cos 2n_1+K_{12}\sin 2n_1-\frac{2}{5}{f}_1-\frac{3}{5}{K}_1\right)\left[2\mathcal{M}-\mathcal{S}(2u_2-2n_2)\right],$$

$$\mathcal{M}=v_2-l_2+e_2\sin v_2,$$

$$\mathcal{S}(2u_2)=\sin 2u_2+e_2\left[\sin (u_2+{\omega }_2)+\frac{1}{3}\sin (3u_2-{\omega }_2)\right],$$

$$\mathcal{C}(2u_2)=\cos 2u_2+e_2\left[\cos (u_2+{\omega }_2)+\frac{1}{3}\cos (3u_2-{\omega }_2)\right],$$

$$f_1=1+\frac{25}{8}{e}_1^2+\frac{15}{8}{e}_1^4+\frac{95}{64}{e}_1^6+\mathcal{O}(e_1^8),$$

$$K_1(e_1,{\omega }_1)=$$

$$-\left(\frac{5}{16}{e}_1^4+\frac{3}{16}{e}_1^6\right)\cos 4{\omega }_1-\frac{3}{16}{e}_1^5\sin 5{\omega }_1+\frac{7}{64}{e}_1^6\cos 6{\omega }_1+\mathcal{O}(e_1^7),$$

$$K_{11}(e_1,{\omega }_1)=$$

$$-\frac{3}{16}{e}_1^3\sin 3{\omega }_1-\left(\frac{1}{16}{e}_1^4-\frac{1}{16}{e}_1^6\right)\cos 4{\omega }_1-\frac{1}{16}{e}_1^5\sin 5{\omega }_1+\frac{3}{64}{e}_1^6\cos 6{\omega }_1+\mathcal{O}(e_1^7),$$

$$K_{12}(e_1,{\omega }_1)=$$

$$-\frac{3}{16}{e}_1^3\cos 3{\omega }_1-\left(\frac{1}{16}{e}_1^4+\frac{5}{32}{e}_1^6\right)\sin 4{\omega }_1+\frac{1}{16}{e}_1^5\cos 5{\omega }_1+\frac{3}{64}{e}_1^6\sin 6{\omega }_1+\mathcal{O}(e_1^7).$$

$$\cos i_{21}=\cos i_1\cos i_2+\sin i_1\sin i_2\cos {\Omega }_{21}.$$

$$\cos i_{21}=\cos (i_2\mp i_1).$$

$$\frac{\sin i_{21}}{|\sin {\Omega }_{21}|}=\frac{\sin i_2}{\sin n_1}=\frac{\sin i_1}{\sin n_2}$$

Full text

planets and satellites: formation — planets and satellites: individual (KOI-12, Kepler-448, KIC 5812701) — planets and satellites: individual (KOI-824, Kepler-693, KIC 5164255) — techniques: photometric

††slugcomment: Accepted for publication in AJ on June 18

1 Introduction

Warm Jupiters (WJs), giant planets in moderately close-in orbits ($7\,\mathrm{days}<P<100\,\mathrm{days}$), pose a similar conundrum to that of hot Jupiters (HJs). A dozen WJs have been found to reside in circular orbits in multi-transiting systems (Huang et al., 2016), in which the orbital planes of the planets are likely well aligned. Such an architecture points to the formation via disk migration (Goldreich & Tremaine, 1980; Lin et al., 1996) as originally proposed for HJs, or in situ formation inside the snow line (Boley et al., 2016; Batygin et al., 2016). Alignments of the planetary orbits with their host stars’ equators, as confirmed for some of them (e.g., Sanchis-Ojeda et al., 2012; Hirano et al., 2012), may also support the absence of past dynamical eccentricity/inclination excitation via planet–planet scattering (Rasio & Ford, 1996) or secular chaos (Wu & Lithwick, 2011).

On the other hand, roughly half of WJs with measured masses from Doppler surveys have significant eccentricities that seem too large to result from disk migration or subsequent planet–planet scattering (Petrovich et al., 2014), but yet too small for their orbits being tidally circularized (Socrates et al., 2012b; Dawson et al., 2015). A possible explanation is that those moderately eccentric WJs are experiencing “slow Kozai–Lidov migration” (Petrovich, 2015): their eccentricities are still undergoing large oscillations driven by the secular perturbation from a close companion (Dong et al., 2014), without being suppressed by other short-range forces (Wu & Murray, 2003), and their orbits shrink only at the high-eccentricity phase. This scenario may indeed reproduce the observed eccentricity distribution of WJs with outer companions (Petrovich & Tremaine, 2016).

Observationally, long-period, massive companions to WJs are nearly as common as those of HJs (Knutson et al., 2014) and their orbital properties might be consistent with what is expected from this scenario (Bryan et al., 2016). Indeed, the apsidal misalignments of some of those eccentric WJs with outer companions provide statistical evidence for the oscillating eccentricity due to a large mutual orbital inclination (Dawson & Chiang, 2014). However, there has been no direct measurement of such a large orbital misalignment between WJs and their outer companions as to induce a large eccentricity oscillation, partly because they are mostly the systems detected with the radial velocity (RV) technique that yields no information on the orbital direction. Notable exceptions include the Kepler-419 system (Dawson et al., 2014) and the doubly-transiting giant-planet system Kepler-108 (Mills & Fabrycky, 2017), where the mutual inclinations were constrained via dynamical modeling of transit timing and duration variations (TTVs and TDVs), though their mutual inclinations are likely too small to drive secular eccentricity oscillations.

Transiting WJs without transiting companions provide a unique opportunity to search for close and mutually-inclined companions as direct evidence for the slow Kozai migration, because the full 3D architecture of the system can be dynamically constrained with a similar technique as used in the above systems. In addition, WJs on eccentric and intermediate orbits, unlike HJs, may still be interacting with the companion strongly, and their eccentricity can also help the TTV inversion for non-transiting objects by producing specific non-sinusoidal features. The TTV search for the outer companions on such “hierarchical” orbits is also complementary to the TTV studies so far, which have mainly focused on nearly sinusoidal signals typical for planets near mean-motion resonances (Hadden & Lithwick, 2016; Jontof-Hutter et al., 2016).

In this paper, I perform a search for non-transiting companions around transiting WJs using transit variations (Section 2) and report the discovery of outer companions to two transiting WJs, Kepler-448b (Bourrier et al., 2015) and Kepler-693b (Morton et al., 2016). Based on TTVs and TDVs of the WJs, I find that the companions are (sub-)stellar mass objects on highly eccentric, au-scale orbits (Sections 3–5). In particular, I confirm a large mutual orbital inclination between the inner WJ and the companion in the latter system, which can induce a large amplitude of eccentricity oscillation and the tidal shrinkage of the inner orbit (Section 6) — exactly as predicted in the “slow Kozai” scenario by a close companion. The companions’ properties, however, are not fully compatible with such migration starting from beyond the snow line. Thus I also assess the possibility of in situ formation (Section 6.3). I discuss possible follow-up observations in Section 7 and Section 8 concludes the paper.

2 Systematic TTV Search for Singly-Transiting Warm Jupiters

To identify the signature of outer companions, I analyzed TTVs of $23$ confirmed, singly-transiting WJs (Section 2.1) with the orbital period $7\,\mathrm{days}<P<100\,\mathrm{days}$ and radius $R_{\mathrm{p}}>8R_{\oplus}$ in the DR24 of the KOI catalog (Coughlin et al., 2016). Systems with multiple KOIs are all excluded, even though they consist of only one confirmed planet and false positives. I found clearly non-sinusoidal TTVs for Kepler-448/KOI-12b, Kepler-693/KOI-824b, and Kepler-419/KOI-1474b. The result is consistent with the TTV search by Holczer et al. (2016), who reported significant long-term TTVs for the same three KOIs in our sample.111Due to the update in the stellar radius in the DR25 catalog, the WJ sample defined above now includes Kepler-522/KOI-318b and Kepler-827/KOI-1355b, for which Holczer et al. (2016) reported long-term TTVs. Interpretation of their TTVs is less clear than for the two systems discussed in the present paper, because of the lack of clear non-sinusoidal features (Kepler-522b) and a low signal-to-noise ratio (Kepler-827b).

Among those planets, Kepler-419b’s TTVs have already been analyzed by Dawson et al. (2014) with RV data, and the companion planet Kepler-419c was found to be a super Jupiter on an eccentric and mutually-aligned orbit with the inner one. Therefore in this paper, I focus on Kepler-448b and Kepler-693b, which both show clear non-sinusoidal TTVs, and masses and orbits of the perturbers can be well constrained.

For the dynamical modeling taking into account the possible orbital misalignment, I reanalyze the transit light curves of Kepler-448b and Kepler-693b to derive both TTVs and TDVs (Section 2.2), consistently with the other transit parameters (Tables 1–4). Here I also fit the planet-to-star radius ratio, $R_{\mathrm{p}}/R_{\star}$, so that the possible transit depth modulation due to the different crowding depending on observing seasons does not mimic duration variations (cf. Mills & Fabrycky, 2017). As shown in Figure 3, I find significant TTVs consisting of a long-term modulation and a short-term, sharp feature (top panels). For Kepler-448b, the spike-like feature is more clear than that reported in Holczer et al. (2016), presumably owing to a more careful treatment of the local baseline modulation. No significant correlation is found between TTVs and the local light-curve slope or the fitted $R_{\rm p}/R_{\star}$, which supports the physical origin rather than due to star spots (Mazeh et al., 2015). I show in Sections 4 and 5 that this unusual feature is reproduced by a periastron passage of an outer non-transiting companion in an eccentric orbit. I also identify a significant ($\sim 5\sigma$) linear trend in the transit duration of Kepler-693b, which points to a mutual orbit misalignment. The duration change is also clear in Figure 6, where each detrended and normalized transit is overplotted around its center, along with the best-fit model.

2.1 Iterative Determination of Transit Times

In the systematic TTV search for WJs, I analyze Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) fluxes obtained with the long-cadence (LC) mode, retrieved from the NASA exoplanet archive. I follow the iterative method as described by Masuda (2015) to derive the consistent transit times and light-curve shapes. The method consists of (1) the determination of the central time for each transit and (2) the refinement of the transit parameters by fitting the mean transit light curve. The fitting in this subsection is performed by minimizing the standard $\chi^{2}$ (Markwardt, 2009) defined as the squared sum of the difference between the model and data divided by the PDCSAP-flux error.

In the first step, I fit the data segments of three times the total duration centered at each transit. I adopt the Mandel & Agol (2002) model for the quadratic limb-darkening law generated with pyTransit package (Parviainen, 2015), multiplied by a second-order polynomial function of time to take into account the local baseline modulation. The fitting is repeated iteratively removing $3\sigma$ outliers. Here a circular orbit is assumed and the values of central transit time $t_{\rm c}$ and three coefficients of the polynomial are fitted, while the other parameters (two limb-darkening coefficients $q_{1}$ and $q_{2}$ defined in Kipping (2013), mean stellar density $\rho_{\star}$, transit impact parameter $b$, planet-to-star radius ratio $R_{\mathrm{p}}/R_{\star}$, and orbital period $P$) are fixed.

In the second step, I shift each transit by the value of $t_{\rm c}$ to align all the transits around time zero. Then the data are averaged into bins of $0.05$ times the transit duration (to reduce the computational time), where the value and error of each bin is set to be the median and $1.4826$ times median absolute deviation of the flux values in the bin. The resulting “mean” light curve is again fitted with the same Mandel & Agol (2002) model for its central time $t_{0}$, normalization constant, $q_{1}$, $q_{2}$, $R_{\rm p}/R_{\star}$, $\rho_{\star}$, and $b$ assuming a circular orbit and fixing the period at the value refined in step (1); only in the first iteration, I adopt a quadratic function of time as a baseline and fit three coefficients rather than one normalization. The values of $q_{1}$, $q_{2}$, $R_{\rm p}/R_{\star}$, $\rho_{\star}$, and $b$ obtained from this fitting is used for the first step of the next iteration.

I perform five iterations for each KOI, starting with the second step based on the transit light curve phase-folded at the period given in the KOI catalog (Coughlin et al., 2016). I fit the resulting transit times with straight lines to search for any TTVs and identify the candidates mentioned above.

2.2 Transit Times and Durations of Kepler-448b and Kepler-693b

For Kepler-448b and Kepler-693b, I perform more intensive analyses taking into account possible variations of transit durations using the short-cadence (SC) data if available. For Kepler-448, I use the SC data for all the quarters, while for Kepler-693 I combine both the LC (Q$2$, Q$7$, Q$8$, Q$10$–Q$12$, Q$14$–$16$) and SC (Q$14$–Q$16$) data.

The method is the same as in Section 2.1 except for the following differences. I additionally fit $R_{\mathrm{p}}/R_{\star}$ and $b$ in the first step and repeated $10$ iterations. Here I fit $R_{\mathrm{p}}/R_{\star}$ so that the seasonal variation in the transit depths not be misinterpreted as the duration variation (e.g., Van Eylen et al., 2013; Masuda, 2015). The resulting $b$ and $R_{\rm p}/R_{\star}$ are combined with $\rho_{\star}$ and $P$ to yield the transit duration $T$ as the average of the total and full durations (see equations 14 and 15 of Winn, 2011). The errors in $t_{\rm c}$ and $T$ are scaled by the square root of the reduced chi-squared $\chi^{2}_{\nu}$ of the best-fit model to enforce $\chi^{2}_{\nu}=1$. When the mean light curve is fitted, the $\chi^{2}$ minimization is complemented by a short Markov Chain Monte Carlo (MCMC) chain (Foreman-Mackey et al., 2013), and the data are binned into one-minute bins for the SC data.

I also tuned the parameters specifically to each system as follows:

Kepler-448b — I fit the data of $1.6$ times duration for each transit. If more than $10\%$ of data is missing in the segment, the transit is excluded from the analysis. I use fourth-order polynomial because the light curve shows short-term wiggles likely due to the stellar rotation ($1.245\pm 0.124\,\mathrm{days}$; McQuillan et al., 2013).

Kepler-693b — I fit the data of twice the duration for each transit and use second-order polynomial because the light curve shows smaller variability than Kepler-448. I omit the transits with more than $10\%$ gaps for the SC data, while for the LC data I omitted those with more than $30\%$. The SC and LC data are analyzed independently and the resulting transit times and durations are averaged to give one measurement if both are available. For each data, I compute $\sqrt{\sigma_{\rm LC}^{2}+\sigma_{\rm SC}^{2}}/2$ and half of the difference between the LC and SC values, and assign the larger of the two as its error. The former was the case for most of the points, while the latter process helped mitigate the effect of a few outliers caused presumably by local features in the light curve.

The resulting transit times, durations, and radius ratios are summarized in Figure 3 and Tables 1 and 2. I also fit the mean transit light curve from the final iteration (binned version of those in Figure 6) with an MCMC including an additional Gaussian error in quadrature (denoted as “photometric jitter”) and derive the posteriors for the transit parameters summarized in Tables 3 and 4. For Kepler-448b, the result agreed well with those from the least-square fit with its standard errors scaled by the square root of $\chi^{2}_{\nu}$. For Kepler-693b, I obtain skewed posteriors and larger errors; the parameters are better determined by the LC data that include more transits.

3 Dynamical Modeling of TTVs and TDVs: Method

3.1 Model Assumptions

The TTVs observed in the two systems (Figure 3) are clearly different from the sinusoidal signal due to a companion near a mean-motion resonance (Lithwick et al., 2012). In particular, rapid timing variations on a short timescale suggest that the perturbing companions’ orbits are eccentric. In addition, such a feature is observed only once for each system, and so the companions must be far outside the WJs’ orbits. Thus I assume a hierarchical three-body system as the only viable configuration and model the observed TTVs and TDVs to derive the outer companions’ masses and orbits.

I only consider the Newtonian gravity between the three bodies regarded as point masses (see Section 3.5.1 for justification), as well as the finite light-travel time in computing timings. To better take into account the hierarchy of the system, I adopt the Jacobi coordinates to describe their orbits: the inner orbit (denoted by the subscript $1$) is defined by the relative motion of the inner planet around the host star, while the outer one (denoted by the subscript $2$) is the motion of the outer companion relative to the center of mass of the inner two bodies. The sky plane is chosen to be the reference plane, to which arguments of periastron and the line of nodes are referred. The direction of $+Z$-axis (which matters the definition of “ascending” node) is taken toward the observer.

3.2 Bayesian Framework

I derive the posterior probability density function (PDF) for the set of system parameters $\theta$ conditioned on the observed data $d$:

[TABLE]

where $I$ denotes the prior information. The normalization factor

[TABLE]

is called the global likelihood or evidence, which represents the plausibility of the model. The prior PDF p(\mbox{\boldmath\theta}|I) is given as a product of the prior PDFs for each parameter, which are assumed to be independent; they are presented in Tables 5 and 6. The likelihood \mathcal{L}\equiv p(\mbox{\boldmathd}|\mbox{\boldmath\theta},I) is defined and computed as described in Sections 3.4 and 3.5.

To invert the observed signals, I utilize the nested-sampling algorithm MultiNest (Feroz et al., 2009; Feroz & Hobson, 2008; Feroz et al., 2013) and its python interface PyMultiNest (Buchner et al., 2014), which allows us to sample the whole prior volume to identify multiple modes if any. I typically utilize $4800$ live points and default sampling efficiency of $0.8$, keep updating the live points until the evidence tolerance of $0.5$ is achieved, and allow for the detection of multiple posterior modes. When multiple modes are detected, MultiNest also computes evidence corresponding to each mode, which is referred to as “local” evidence.

3.3 Procedure for Finding Solutions

To check on the reliability of our numerical scheme and to find all the possible posterior modes, I adopt the following procedure. First, I use the analytic TTV formula for a hierarchical triple system (Borkovits et al., 2015, Appendix A) to fit only the TTV signal. This allows us to search a wide region of the parameter space and resulted in one global mode containing two peaks. The resulting solution was also found to be consistent with a rough analytical estimate based on the observed TTV features (see, e.g., Section 4.2). Then, I fit the same TTVs numerically adopting a slightly narrower prior range that well incorporates the global mode found from the analytic fit. The resulting posterior was found to be consistent with the one derived from the analytic fit. This agreement validates the numerical scheme I rely on. Finally, the same numerical method is used to model both TTVs and TDVs simultaneously to determine the physical and geometric properties of the outer companions. In Sections 4 and 5, I mainly report and discuss the results from the final fit, while the analytic and numerical posteriors from TTVs alone are found in Appendix B for comparison.

3.4 TTV Modeling

The likelihood for the TTV modeling is defined as follows:

[TABLE]

where $t_{i}$ and $\sigma_{t,i}$ are the transit times and their errors obtained in Section 2.2 (Tables 1 and 2). I also include $\sigma_{\rm TTV}$ as a model parameter that takes into account the additional scatter and marginalize over it to obtain more conservative constraint; it turns out that this parameter is not correlated with any other physical parameters (Figures 29 and 30) and hence does not affect the global shape of the joint PDF.

The model transit times $t^{\rm model}_{i}$ are computed both analytically and numerically as mentioned in Section 3.3; this is for cross-validation, as well as to obtain insights into how the parameters are determined. In both cases, the model includes $14$ parameters in addition to $\sigma_{\rm TTV}$ (see Tables 8 and 9): orbital period, orbital phase, eccentricity, and argument of periastron for both inner and outer orbits; cosine of the outer orbit’s inclination (inner one is fixed to be [math] for simplicity; this does not affect the result); difference in the longitudes of the ascending node; and masses of the host star and outer companion (here I fix the inner planet’s mass to be $1M_{\rm Jup}$ as it is not determined from TTVs). While I use the time of inferior conjunction $t_{\rm ic}$ and orbital period $P$ for the inner orbit, I choose the time of periastron passage $\tau_{2}$ and periastron distance relative to the inner semi-major axis $q_{2}/a_{1}$ to specify the outer orbital phase and period. This is because the latter two parameters are more directly related to the position and duration of the “spike” in the observed TTVs than $t_{\rm ic}$ and $P$, and thus better determined. I fit the mass scale of the whole system as well because I include the light-travel time effect; in practice, however, the observed TTVs are insensitive to this parameter and its value is solely constrained by the prior knowledge.

I first compute analytic transit times using the formula by Borkovits et al. (2011, 2015) developed for hierarchical planetary systems and triple-star systems (Appendix A). I include the light-travel time effect (LTTE) and $P_{2}$-timescale dynamical effect due to the quadrupole component of the perturbing potential; I did not find notable changes even with the octupole terms. The former is due to the motion of the inner binary (here the central star orbited by the inner planet), which changes finite time for the light emitted from the star and blocked by the planet to travel to us. The latter process involves the actual variation in the orbital period of the inner binary due to the tidal gravitational field exerted by the companion. Both of these effects are routinely observed in triple-star systems containing eclipsing binaries (e.g., Rappaport et al., 2013; Masuda et al., 2015).

Numerical TTVs are computed using TTVFast code (Deck et al., 2014) by integrating the orbits of the three bodies and finding the times at which sky-projected distance between the star and the inner planet, $d_{\rm sky}$, becomes minimum. I choose the time step to be $0.3\,\mathrm{days}$ for Kepler-448 and $0.1\,\mathrm{days}$ for Kepler-693, which are roughly $1/60$ and $1/150$ of the inner orbital periods, respectively. To compute transit times, I modify the default output of the TTVFast code to take into account the effect of light-travel time using the line-of-sight coordinate of the center of mass of the inner binary (i.e., central star and inner planet).

3.5 Joint TTV and TDV Modeling

The joint TTV and TDV modeling was performed using the following likelihood:

[TABLE]

where the second row is defined analogously to $\mathcal{L}_{\rm TTV}$ with the transit time $t$ replaced by the transit duration $T$. I also include the constraint on the transit impact parameter $b_{\rm mean}$ from the mean transit light curve.222Here we can use $b_{\rm mean}$ derived assuming a circular orbit because the parameter is based on the shape of the light curve alone. A possible bias due to the non-zero eccentricity is absorbed in the fitted stellar density, which is not used in any of the following analyses. See, e.g., Eqn. (28) and (29) of Winn (2011). The TDV modeling additionally requires the mean density of the host star $\rho_{\star}$ (equivalently stellar radius) and the “jitter” for durations $\sigma_{\rm dur}$, and $\theta$ consists of $17$ parameters listed in Tables 5 and 6 as fitted parameters. Here I fit inner orbit’s inclination and inner planet’s mass, although the latter is not constrained at all by the data.

This joint modeling is performed only numerically using TTVFast. In addition to the transit times above, I use another default output, $d_{\rm sky}$ at each transit center, to obtain the model transit impact parameter $b=d_{\rm sky}/R_{\star}$, where the stellar radius $R_{\star}$ is given by $(3m_{\star}/4\pi\rho_{\star})^{1/3}$; here I adopt $b$ at the transit around the center of the observing duration as $b^{\rm model}_{\rm mean}$. The transit durations are obtained from yet another default output of the code $v_{\rm sky}$, the sky-projected relative star–planet velocity at the transit center, via $T=2\sqrt{1-b^{2}}R_{\star}/v_{\rm sky}$.

3.5.1 Additional TDV Sources

While I only consider the Newtonian gravitational interaction between point masses, transit durations may also be affected by secular orbit variation due to stellar quadrupole moment and/or general relativistic effect. Indeed, we expect that Kepler-448 has a relatively large quadrupole moment $J_{2}$ due to its rapid rotation (McQuillan et al., 2013; Bourrier et al., 2015), and the inner orbits are moderately eccentric in both systems as will be shown later. Here we show that those effects on transit durations are negligible (the effects on transit times are even smaller; see Miralda-Escudé, 2002).

The duration drift due to the nodal precession is given by (Miralda-Escudé, 2002)

[TABLE]

where $\lambda$ is the sky-projected stellar obliquity. For Kepler-448b with $b\simeq 0.4$, $a/R_{\star}\simeq 20$, and $\lambda\simeq 10\arcdeg$, this results in the duration change over four years

[TABLE]

where the quoted value of $J_{2}$ is motivated by theoretical modeling and observational inference for a star with similar parameters (Szabó et al., 2012; Masuda, 2015). This $\Delta T$ is smaller than the measurement precision. The same is also true for Kepler-693b, for which $J_{2}$ is likely smaller and measurements of $T$ are less precise.

The drift caused by general relativistic apsidal precession is (Miralda-Escudé, 2002)

[TABLE]

or

[TABLE]

over four years; this is also negligible. Note that possible apsidal precession induced by the gravitational perturbation from the outer companion is already taken into account in our Newtonian model.

4 Results: Kepler-448/KOI-12

The resulting posterior PDF from the MultiNest analysis and the corresponding models are shown in Figures 7 (red solid lines) and Figure 29. The summary statistics (median and $68\%$ credible interval) of the marginal posteriors as well as the priors adopted for those parameters are given in Table 5. The solution consists of two separated posterior modes (denoted as Solution 1 and Solution 2) that reproduce the observed TTVs and TDVs equally well, without any significant difference in the local evidence computed for each mode.

In both solutions, the outer companion (we tentatively call it Kepler-448c) is likely to have a mass in the brown-dwarf regime ($22^{+7}_{-5}\,M_{\rm Jup}$) and reside in a highly eccentric orbit close to the inner WJ. In spite of the small outer orbit, its pericenter distance $1.46^{+0.07}_{-0.06}\,\mathrm{au}$ is more than nine times larger than the inner semi-major axis and well within the stable regime (Mardling & Aarseth, 2001). The system has one of the smallest binary separations compared to the planet apastron $0.21^{+0.05}_{-0.04}\,\mathrm{au}$ among the known planetary systems with (sub-)stellar companions; see figure 8 of Triaud et al. (2017).

I also find a modest eccentricity $0.34^{+0.08}_{-0.07}$ for the inner WJ from the combination of timing and duration information. While the TTVs alone favor an even larger value (Appendix B), transit duration combined with the spectroscopic prior on $\rho_{\star}$ lowers it. The inner WJ mass is also floated between $0.1\mathchar 45\relax 10\,M_{\rm Jup}$ but no constraint better than that from the RV data ($<10\,M_{\rm Jup}$ as the $3\sigma$ upper bound; Bourrier et al., 2015) is found. This is because the TTVs of the inner WJ does not depend on its own mass to the first order.

While the above constraints essentially come from TTVs, TDVs possibly allow us to constrain the mutual inclination $i_{\rm 21}$ between the inner and outer orbits. In fact, the two solutions are different in this respect: Solution 1 is “prograde”, i.e., the two orbits are in the same directions, while Solution 2 corresponds to a “retrograde” configuration. Non-zero mutual inclinations are slightly favored in both solutions with the posterior probability that $39\fdg 2<i_{21}<140\fdg 8$ being $27\%$. However, it is more or less due to the large prior volume corresponding to the misaligned solution and a wide range of the mutual inclination is allowed by the data (Figure 8). In particular, the data are totally consistent with the coplanar configuration, as visually illustrated in the right panels of Figure 7. I found the evidence for the coplanar model is not significantly different from the misaligned case. Only the configuration with $i_{\rm 21}\simeq 90\arcdeg$ is slightly disfavored due to the lack of a large TDV expected from the strong perturbation in the direction perpendicular to the orbital plane.

The joint TTV/TDV analysis also allows for constraining the impact parameter during a possible occultation, $b_{\rm occ}$, based on that during the transit as well as the eccentricity and argument of periastron of the inner orbit. I find $b_{\rm occ}=0.22\pm 0.02$, which means that the secondary eclipse should have been detected if Kepler-448b was a star. The fact strengthens the planetary interpretation by Bourrier et al. (2015) who confirmed $m_{\mathrm{b}}<10\,M_{\rm Jup}$ at the $3\sigma$ level.

The following subsections detail the prior information and interpretation of TTVs and TDVs.

4.1 Adopted Parameters

I adopt $m_{\star}=1.51^{+0.14}_{-0.09}\,M_{\odot}$ and $\rho_{\star}=0.758^{+0.12}_{-0.11}\,\mathrm{g\,cm^{-3}}$ (converted from the radius $r_{\star}=1.41\pm 0.06\,R_{\odot}$) as the priors based on the spectroscopic values by Bourrier et al. (2015). The values agree with, and are more precise than, the latest KIC values (Mathur et al., 2016): $m_{\star}=1.386^{+0.093}_{-0.076}\,M_{\odot}$, $r_{\star}=1.367^{+0.267}_{-0.067}\,R_{\odot}$, and $\rho_{\star}=0.7647^{+0.09607}_{-0.2874}\,\mathrm{g\,cm^{-3}}$. Note that I do not use the values from the joint analysis in Bourrier et al. (2015) because they are derived assuming a circular orbit for the inner transiting planet, which turns out unlikely from our dynamical modeling. I adopt $b_{\rm mean}=0.373$ and $\sigma_{b}=0.006$ based on the posterior from the mean transit light curve (Table 3).

4.2 Constraints from TTVs

The mass and some elements of the outer orbit ($\tau_{2}$, $q_{2}/a_{1}$, $e_{2}$) are well determined from TTVs in spite of the non-transiting nature of the companion. The relationship between these parameters and observed TTV features can roughly be understood as follows, with the help of the analytic formula in Appendix A.

Most of the information comes from the spike-like feature around $\mathrm{BJD}\sim 2454833+1050$, which is caused by the close encounter of the outer body and thus defines the time of its periastron passage $\tau_{2}$ (the effect is represented by $\mathcal{S}$ and $\mathcal{C}$ functions in Equations A9 and A10; see also Figure 28). In addition, its duration $\Delta t$ is sensitive to both $a_{2}$ and $e_{2}$: the former determines the overall orbital time scale and the latter changes the fraction of time spent around the periastron. The $\Delta t$ may be roughly estimated using Kepler’s second law as follows:

[TABLE]

where its order of magnitude is not so sensitive to the rather arbitrary choice of the interval of integration ($-\pi/2$ to $\pi/2$). Since $\Delta t\sim 300\,\mathrm{days}$ and $P_{1}\sim 20\,\mathrm{days}$, this estimate gives $q_{2}/a_{1}\sim 10$. Finally, its amplitude $\Delta A$ is given by (Equations A2 and A3)

[TABLE]

assuming $m_{\mathrm{b}},m_{\mathrm{c}}\ll m_{\star}$. Combining $\Delta A\sim 2\,\mathrm{min}$ with the above estimate $q_{2}/a_{1}\sim 10$, we find $m_{\mathrm{c}}/m_{\star}\sim 10^{-2}$. These numbers are in reasonable agreement with those from the full dynamical modeling.

In addition to the spike, a long-term modulation due to “tidal delay” (Borkovits et al., 2003; Agol et al., 2005, represented by $\mathcal{M}$ function in Equation A8) is also visible (blue curves in Figure 28): tidal force from the outer companion delays the inner binary period, depending on the distance between the two. This effect, combined with the short-term spike, allows for further constraints on $\tau_{2}$, $\omega_{2}$, $e_{2}$, and $P_{2}$. In principle, this additional constraint enables the separate determination of $a_{2}$ and $e_{2}$ rather than only $q_{2}$, although $a_{2}$ is not well constrained because it is not clear whether the whole cycle of the outer orbit is covered given the only one periastron passage observed.

The analytic expressions A2–A7 show that the short-term modulation represented by $\mathcal{S}$ or $\mathcal{C}$ function is produced only when $e_{1}\neq 0$ or $i_{\rm 21}\neq 0$. In fact, the observed TTVs alone are found to be fitted well either by the non-zero $e_{1}$ ($\delta_{\rm ecc}$ in Equation A2) or non-zero $i_{\rm 21}$ ($\delta_{\rm noncopl}$) effects. This causes the anti-correlated degeneracy between $e_{1}$ and $i_{\rm 21}$. Since $e_{1}$ also affects the spike amplitude, we have two classes of solutions: high-$e_{1}$–low-$i_{21}$–low-$m_{\mathrm{c}}$ solution and low-$e_{1}$–high-$i_{21}$–high-$m_{\mathrm{c}}$ solution. The joint TTV/TDV model slightly favors the latter.

The directions of the orbits, $\omega_{1}$, $\omega_{2}$, $\cos i_{2}$, $\Omega_{21}$ are not well constrained from the TTVs alone. The complicated degeneracy between $\omega_{2}$ and $\Omega_{21}$ comes from the fact that TTVs are rather sensitive to such “dynamical” angles referred to the invariant plane as $n_{1}$ and $n_{2}$ in figure 1 of Borkovits et al. (2015). Nevertheless, I use $\cos i_{2}$ and $\Omega_{21}$ because of the simplicity in implementing the isotropic prior. The relationship between the two sets of angles are summarized in Appendix A.2.

Finally, the LTTE effect (Equation A1) is about $0.1\,\mathrm{min}(a_{2}/\mathrm{au})$ for $m_{\mathrm{c}}/m_{\star}\sim 10^{-2}$. This term therefore does not play a major role, as also seen in Figure 28.

4.3 Constraints from TTVs and TDVs

In the joint analysis, I find a smaller $e_{1}$ than that derived from the TTVs alone, because the combination of $T$, $b$, and $\rho_{\star}$ favors a small $e_{1}\sin\omega_{1}$ (see equations 16, 18, and 19 of Winn, 2011). This constraint is combined with those from TTVs to better determine $e_{1}$ and $\omega_{1}$. The distribution of $\omega_{1}$ is still bimodal because no constraint is available for $e_{1}\cos\omega_{1}$ from TDVs.

These additional constraints on $e_{1}$ and $\omega_{1}$ partly solve the degeneracies with other angles mentioned above. In particular, the value of $m_{\mathrm{c}}$ from the joint TTV/TDV analysis is larger than from the TTVs alone, because the low-$e_{1}$–high-$i_{21}$–high-$m_{\mathrm{c}}$ solution is favored.

Absence of significant TDVs only weakly constrains the mutual inclination. As mentioned, the solution with $i_{21}\sim 90\arcdeg$ is disfavored; if this were the case, a large tangential perturbation should have produced significant TDVs around the periastron passage, which are not present in the data. The solutions with $\Omega_{21}\sim 0\arcdeg$ and $\Omega_{21}\sim 180\arcdeg$ are basically indistinguishable and the log-likelihood is not sensitive to $\cos i_{2}$ either. In the TDV models (bottom left panel of Figure 7), both the positive and negative bumps are seen around the outer periastron because the inner inclination is perturbed in the opposite ways depending on the sign of $\cos i_{2}$. The current data do not favor either of the cases significantly.

5 Results: Kepler-693/KOI-824

The resulting posterior PDF from the MultiNest analysis and the corresponding models are shown in Figures 9 (red solid lines) and Figure 30. The summary statistics (median and $68\%$ credible region) of the marginal posteriors as well as the priors adopted for those parameters are given in Table 6. Again I found two almost identical solutions with prograde and retrograde orbits.

Similarly to Kepler-448c, the outer companion (we tentatively call it Kepler-693c) orbits close to the inner WJ in a highly eccentric orbit, except that it has a larger mass $1.5^{+0.6}_{-0.4}\times 10^{2}\,M_{\rm Jup}$ and can be a low-mass star. The outer pericenter distance of $1.5\pm 0.2\,\mathrm{au}$ is also similar to Kepler-448c and satisfies the stability condition (Mardling & Aarseth, 2001). Again non-zero eccentricity is favored for the inner orbit ($e_{1}=0.2^{+0.2}_{-0.1}$), and its mass is not constrained at all from the data.

A notable difference from the previous case is the clear TDV signal, which tightly constrains the mutual inclination to be $|i_{21}-90^{\circ}|\simeq 40\arcdeg$ (Figure 8, bottom). I checked that the data can never be explained by the aligned configuration (right column of Figure 9): I find the Bayesian evidence for the mutually-inclined model ($\ln\mathcal{Z}=-89.18\pm 0.08$) is larger than that of the coplanar one ($\ln\mathcal{Z}_{\rm copl}=-103.1\pm 0.07$) by a factor of $10^{6}$. The observed mutual inclination is the largest among those dynamically measured for planetary systems, and likely above the critical angle for the Kozai oscillation (posterior probability that $39\fdg 2<i_{21}<140\fdg 8$ is $80\%$).333Santerne et al. (2014) inferred the orbital inclination of the binary companion of the transiting WJ host KOI-1257/Kepler-420 to be $i_{2}=18\fdg 2^{+18\fdg 0}_{-5\fdg 4}$ ($68.3\%$ interval), based on the trend in the SOPHIE RVs and bisector and full-width half maximum of the stellar lines, combined with the Kepler light curve and stellar spectral energy distribution. If confirmed, the value translates into a large $i_{21}\gtrsim 90\arcdeg-i_{2}$ ($\Omega_{21}$ is not constrained at all) and a similar eccentricity evolution as discussed in the present paper is expected, although they caution that the outer binary inclination is poorly constrained, with the $99\%$ limit on $i_{2}$ being $[8\fdg 2,85\fdg 2]$. In Section 6, I will show that the perturbation from the inclined companion does cause a large eccentricity oscillation of the inner WJ.

Kepler-693b’s mass is not measured, and the planet was statistically validated by Morton et al. (2016). A possible concern is that the absence of secondary eclipse, which the statistical validation partly relies on (in addition to other factors including the transit shape and non-detection of the companion via the AO imaging by Baranec et al., 2016), may lose its meaning if $b_{\rm occ}$ is larger than unity due to the non-zero inner eccentricity. However, our dynamical modeling finds $b_{\rm occ}=0.64^{+0.08}_{-0.09}$, which excludes the possibility. In addition, the derived mean stellar density is also compatible with the KIC value. Therefore, the low false positive probability (less than $1/3000$) derived by Morton et al. (2016) is still valid in the light of our new constraints.

5.1 Adopted Parameters

I adopt $m_{\star}=0.793^{+0.054}_{-0.029}\,M_{\odot}$ and $\rho_{\star}=1.931^{+0.526}_{-0.1773}\,\mathrm{g\,cm^{-3}}$ (corresponding to $r_{\star}=0.833^{+0.033}_{-0.062}\,R_{\odot}$) from KIC as the priors. I adopt $b_{\rm mean}=0.6$ and $\sigma_{b}=0.1$ that roughly incorporates the $68\%$ credible region of the posterior from the mean transit light curve (Table 4).

5.2 Constraints from TTVs

The situation is basically similar to the Kepler-448 case, except that the TTV amplitude is larger by almost an order of magnitude and so is the companion’s mass. In addition, the whole shape of the short-term feature is not entirely observed due to the data gap, as shown in the large scatter of the models (Figure 9). This explains a weaker constraint on $q_{2}/a_{1}$, which is mainly determined by the duration of the feature (Section 4.2).

5.3 Constraints from TTVs and TDVs

It is also the case in this system that the duration data point to a lower inner eccentricity than favored by the TTVs alone and increase the estimated mass. A striking difference is the clear trend in the duration that points to non-zero $i_{21}$. The trend comes from the nodal precession of the inner orbit (Miralda-Escudé, 2002) and has also been observed in the Kepler-108 system recently (Mills & Fabrycky, 2017). While a similar duration drift can also be caused by the apsidal precession of the eccentric inner orbit, this is unlikely to be the case given the failure of the coplanar model. In fact, both precession frequencies are basically comparable, but the effect of the apsidal precession on the duration is smaller than that of the nodal precession by a factor of $R_{\star}/a_{1}$ (Miralda-Escudé, 2002). This difference can be significant for WJs with $a/R_{\star}>10$.

5.3.1 Mutual Inclination from the Likelihood Profile

The constraint on the mutual inclination in Figure 8 and Table 6 is based on the one-dimensional Bayesian posterior marginalized over the other parameters. This depends not only on the likelihood and prior function but also on volume of the parameter space with high posterior probabilities. Thus, the resulting credible interval is generally not the same as the likelihood-based confidence interval even for a uniform prior, and a low value in the marginal posterior probability does not necessarily mean a poor fit to the data (Feroz et al., 2011).

To see what value of the mutual inclination is “excluded” by the data alone in the frequentist sense, I also estimate the confidence interval for $i_{21}$ based on the likelihood profile. Specifically, I derive the maximum likelihood $\hat{\mathcal{L}}(i_{21})$ for each fixed value of $i_{21}$ optimizing the other parameters, and examine its form as a function of $i_{21}$. This has been done by searching for minimum $\chi^{2}$ solutions for a grid of $\cos i_{2}\in[-1,1]$ and $\Omega_{21}\in[-90\arcdeg,90\arcdeg]$ (prograde solutions). Here the constraints on $m_{\star}$ and $\rho_{\star}$ are also incorporated in $\chi^{2}$, and $\sigma_{\rm TTV}$ and $\sigma_{\rm dur}$ were fixed to be zero.

Figure 10 shows the resulting profile of $-2\ln(\hat{\mathcal{L}}(i_{21})/\hat{\hat{\mathcal{L}}})$, where $\hat{\hat{\mathcal{L}}}$ denotes the maximum likelihood found by optimizing all the model parameters including $i_{21}$. This is equivalent to the chi-squared difference from its minimum value, $\Delta\chi^{2}$, in our current setting. Here the $\Delta\chi^{2}$ value is scaled so that the minimum $\chi^{2}$ solution has $\chi^{2}_{\nu}=1$. The resulting $1\sigma$ “confidence” interval is $i_{21}=47\arcdeg^{+17\arcdeg}_{-15\arcdeg}$, and the $2\sigma$ lower bound is found to be $i_{21}=18\arcdeg$. Thus I conclude that a high mutual inclination is indeed favored by the data and the conclusion is insensitive to the prior information on $i_{21}$.

6 Long-term Orbital Evolution

A high mutual inclination confirmed for the Kepler-693 system may result in a large eccentricity oscillation of the inner orbit as long as the perturbation is strong enough to overcome other short-range forces. If this is indeed the case, the system may serve as a direct piece of evidence that some WJs are undergoing eccentricity oscillations. Even in the Kepler-448 system where highly-inclined solutions are not necessarily favored, significant eccentricities of both the inner and outer orbits may still lead to the eccentricity excitation due to the octupole-level interaction.

Given the full set of parameters constrained from the dynamical analysis, we can assess the future of the systems rather reliably by extrapolating the dynamical solutions. In this section, we explore the effect of the secular perturbation due to the outer companions on the inner planet’s orbit and its tidal evolution.

6.1 Oscillation of the Inner Eccentricity

I compute secular evolution of both the inner and outer orbits along with the spins of the star and inner planet. I use the code developed and utilized in Xue & Suto (2016) and Xue et al. (2017), which takes into account (i) gravitational interaction up to the octupole order and (ii) precessions due to general relativity as well as tidal and rotational deformation of the bodies. Here I neglect magnetic braking of the star and tidal dissipation inside the bodies, and assume zero stellar and planetary obliquities for simplicity. The rotation periods are set to be $1.25\,\mathrm{days}$ for Kepler-448, $10\,\mathrm{days}$ for Kepler-693, and $1\,\mathrm{day}$ for the inner planets, although the spin evolution does not affect the result significantly. I adopt standard values for the dimensionless moments of inertia (0.059 and 0.25 for the star and inner planet) and the tidal Love numbers (0.028 and 0.5). The orbits are integrated for $10\,\mathrm{Myr}$ (sufficiently longer than the oscillation timescale; see below), starting from $1000$ random sets of parameters sampled from the posterior distribution from the dynamical analyses (Sections 4 and 5).

Figure 13 shows the initial (black circles) and minimum (red diamonds) values of $a_{1}(1-e_{1}^{2})$ over the course of evolution, the latter of which correspond to the maxima of $e_{1}$. In both systems, significant eccentricity oscillations occur at least for some of the solutions. If we adopt $a(1-e^{2})<0.1\,{\rm au}$ as a conventional threshold for the migrating WJs (Socrates et al., 2012b; Dong et al., 2014; Dawson et al., 2015), the periastrons become close enough to drive the tidal migration for $12\%$ of the solutions for Kepler-448b and $96\%$ for Kepler-693b, excluding the tidally-disrupted cases shown with transparent colors. If we choose $a(1-e^{2})<0.05\,{\rm au}$ instead as a threshold (e.g., Anderson et al., 2016), the fractions become $5\%$ and $33\%$ for Kepler-448b and Kepler-693b, respectively. Large eccentricity oscillations (and hence small minimum periastrons) are observed mainly for $40\arcdeg\lesssim i_{21}\lesssim 140\arcdeg$; this explains why a larger fraction of solutions have sufficiently small $a(1-e^{2})$ for the Kepler-693 system.

The current eccentricities of the two WJs ($1\sigma$ region shown with vertical dotted lines in Figure 16) also turn out to be the most likely values to be observed over the course of oscillation if we are random observers in time; they are around the peaks of the inner eccentricity distribution during the $10\,\mathrm{Myr}$ evolution averaged over the dynamical solutions (gray histograms in Figure 16).

For highly inclined solutions, both libration and circulation of the argument of periastron (Kozai, 1962), modified by the octupole effects, are observed.444On the other hand, the libration of the difference in the apsidal longitudes $\Delta\varpi_{\rm inv}$, as discussed in Dawson & Chiang (2014), was not observed in the current simulations. In either of the two systems, the Kozai-Lidov timescale $\tau_{\rm KL}$ (Kiseleva et al., 1998) and the octupole one $\tau_{\rm Oct}$ are both found to be much shorter than that of general relativistic apsidal precession $\tau_{\rm GR}$ for any solution found from TTVs and TDVs: for reference, I find typical $\tau_{\rm KL}\sim 10^{-2}\,\mathrm{Myr}$, $\tau_{\rm Oct}\sim 10^{-1}\,\mathrm{Myr}$, and $\tau_{\rm GR}\sim 10^{0.5}\,\mathrm{Myr}$ for Kepler-448, and $\tau_{\rm KL}\sim 10^{-3}\,\mathrm{Myr}$, $\tau_{\rm Oct}\sim 10^{-1.5}\,\mathrm{Myr}$, and $\tau_{\rm GR}\sim 10^{0.5}\,\mathrm{Myr}$ for Kepler-693. Thus the perturbation from the “close friends” is indeed strong enough to overcome general relativity (Dong et al., 2014). Note that this is the case for the octupole effect as well, which explains why some low-$i_{21}$ solutions in the Kepler-448 system lead to a large eccentricity oscillation.

6.2 Migration Timescales

Are the two WJs currently migrating into HJs due to the eccentricity oscillation? If this is the case, the current migration timescale $\tau_{\rm mig}$ needs to be comparable to the system age $\tau_{\rm age}$. If $\tau_{\rm mig}\gg\tau_{\rm age}$ they are unlikely to be migrating, while if $\tau_{\rm mig}\ll\tau_{\rm age}$ the WJs should have evolved into HJs rapidly and we are unlikely to observe the system in the current state. Here we perform an order-of-magnitude comparison between the two timescales, given their large observational and theoretical uncertainties.

Figure 19 shows the distribution of tidal migration timescale given by equation 2 of Petrovich & Tremaine (2016) at the minimum periastron distance computed in Section 6.1. Strictly speaking, the relevant timescale is expected to be longer because oscillation of $e_{1}$ slows down the migration (Petrovich, 2015), but in our case the maximum eccentricity is not so close to unity (since $a_{1}$ is already small) that the modification is unlikely large (cf. figure 2 of Petrovich, 2015); thus we simply neglect the correction. The timescales are computed for three different values of the viscous time of the planet, which gives a characteristic timescale for dissipation: $t_{\rm V}=0.015\,\mathrm{yr}$, $t_{\rm V}=1.5\,\mathrm{yr}$, and $t_{\rm V}=150\,\mathrm{yr}$, while the dissipation inside the star is neglected (see Socrates et al. (2012a) for comparison with different parameterizations). The three values roughly correspond to (i) very efficient tidal dissipation required for some high-eccentricity migration scenarios to explain the observed HJs (Petrovich, 2015; Hamers et al., 2017), (ii) dissipation required to circularize the orbits of HJs with $P\leq 5\,\mathrm{days}$ within $10\,\mathrm{Gyr}$ (Socrates et al., 2012a), and (iii) values calibrated based on a sample of eccentric planetary systems (Hansen, 2010; Quinn et al., 2014), while the limit from the Jupiter–Io system ($t_{\rm V}\gtrsim 15\,\mathrm{yr}$) lies in between the latter two. I also indicate (rough) estimates for the ages of the two systems with vertical dashed lines: $1.5\,\mathrm{Gyr}$ for Kepler-448 based on spectroscopy (Bourrier et al., 2015) and $5\,\mathrm{Gyr}$ for Kepler-693 as a tentative value given that the host star has dimensions of a K dwarf.

The comparison between the histogram and the dashed line shows that the migrating solutions with $\tau_{\rm mig}\sim\tau_{\rm age}$ exist for a wide range of $t_{\rm V}$ for the Kepler-693 system. In particular, the eccentricity oscillation plays a crucial role for $t_{\rm V}\gtrsim 1.5\,\mathrm{yr}$ so that such solutions realize. The migrating solutions also exist for the Kepler-448 system, though they seem plausible only for a small fraction of significantly misaligned solutions that lead to the large eccentricity oscillation, or require efficient tidal dissipation with $t_{\rm V}\lesssim 1.5\,\mathrm{yr}$.

6.2.1 Possible Fates of the Inner Planets

The timescale arguments above indicate the inner WJs may evolve into HJs within the lifetime of the system. If this is the case, HJ systems with close substellar companions as found by a long-term RV monitoring (Triaud et al., 2017) may have been WJs like ours in the past.

As a proof of concept, I compute the evolution of the two systems for $1\,\mathrm{Gyr}$, including the tidal dissipation with the planetary viscous timescale $t_{\rm V}=1.5\,\mathrm{yr}$ for an illustration. I fix $t_{\rm V}=50\,\mathrm{yr}$ for the star. I stop the calculation when the inner orbit is circularized (both $a_{1}<0.1\,\mathrm{au}$ and $e_{1}<0.01$ are achieved) before $1\,\mathrm{Gyr}$. If we change $t_{\rm V}$, we expect things just happen on a different timescale corresponding to the change.

Figure 22 compares the initial and final semi-major axes from those simulations. We see that some solutions do survive the tidal disruption and evolve into HJs within $1\,\mathrm{Gyr}$. Such an outcome is rarer in the Kepler-693 system than in the Kepler-448 system. This is consistent with the expectation that the former system is likely older than the latter at least by a factor of a few. In fact, this kind of path may be even rarer than it seems in the right histograms of Figure 22, because some of the survived HJs (shown with bluer colors) have experienced the circularization on a much shorter timescale compared to the system age due to the rapid eccentricity surge: if we are random observers in time, it is a priori unlikely to observe a system with such a short remaining lifetime compared to the current age (Gott, 1993). However, I do not attempt to correct for this effect here, given that the outcome will be sensitive to the uncertain tidal parameter in any case. If correctly taken into account, this kind of argument will potentially allow for better constraints on the system parameters.

6.3 Implications for the In Situ Formation Scenario

While the observed properties of the two inner planets are consistent with those of migrating WJs, presence of the outer (sub)stellar companions as close as $1.5\,\mathrm{au}$ challenges the high-eccentricity migration scenario from beyond the snow line. Planets in S-type orbits around tight binaries with periastron distances less than $10\,\mathrm{au}$ have also been reported around KOI-1257 (Santerne et al., 2014), Kepler-444 (Dupuy et al., 2016), HD 59686 (Ortiz et al., 2016), and possibly $\nu$ Octantis (Ramm et al., 2016). If confirmed to be a low-mass star, Kepler-693c has the smallest peristron among such stellar companions.

A similar issue has also been discussed for WJs with outer planetary-mass companions (Antonini et al., 2016): the outer orbits in these systems, if primordial, are in most cases too small for the inner WJs to have migrated from $\gtrsim 1\,\mathrm{au}$. In addition, population synthesis simulations of high-eccentricity migration from $\gtrsim 1\,\mathrm{au}$, either via the companion on a wide orbit (Anderson et al., 2016; Petrovich & Tremaine, 2016) or secular chaos in multiple systems (Hamers et al., 2017), have difficulty in producing a sufficient number of WJs relative to HJs. These may also argue for the WJ formation via disk migration or in situ. Considering the prevalence of compact super-Earth systems revealed by Kepler, the latter can well be possible theoretically (Lee et al., 2014) and may also have observational supports (Huang et al., 2016).

The companions discovered in the Kepler-448 and Kepler-693 systems may further argue for the in situ origin. Such low-mass stellar or brown-dwarf companions on au-scale orbits may be formed via fragmentation at a larger separation followed by the orbital decay due to dissipative dynamical interactions involving gas accretion and disks, which proceed in $\lesssim 1\,\mathrm{Myr}$ (Bate et al., 2002; Stamatellos & Whitworth, 2009; Bate, 2012). This implies that giant-planet formation and migration must have completed very quickly if they preceded those of the outer companion. Alternatively, the companions may have arrived at the current orbit well after disk migration and disk dispersal via chaotic dissolution of an initial triple-star system or an impulsive encounter with a passing star (Marzari & Barbieri, 2007; Martí & Beaugé, 2012). While these scenarios are compatible with the eccentric and inclined outer orbit, they may suffer from the fine-tuning problem. In the former scenario, for example, a binary orbit typically shrinks only by a factor of a few, limited by energy conservation (Marzari & Barbieri, 2007). The outcomes are likely more diverse in the latter, but only a small fraction of them usually constitutes a suitable solution (Martí & Beaugé, 2012), and such close encounters as to alter the binary orbit significantly are likely rare when the planet formation is completed, even in a cluster with $\sim 10^{3}$ stars (Adams et al., 2006). Also note that tidal friction associated with the close encounter with the primary (e.g., Kiseleva et al., 1998) is unavailable to shrink the binary orbit in the presence of the inner planet. Considering these possible difficulties of the alternative scenarios, in situ formation in a tight binary seems to be an attractive possibility that provides simple solutions both for our two systems and for other theoretical and observational issues, although the disk migration followed by rearrangement of the outer orbit cannot be excluded.

This motivates us to examine whether the moderate eccentricities of our WJs can be explained in the in-situ framework, in which a near-circular orbit is normally expected. Here we consider a specific form of question, in the same spirit as Anderson & Lai (2017):555See Mustill et al. (2017) for other possible pathways of eccentric WJ formation. suppose that the inner WJs were produced into circular orbits with the current semi-major axes, can their observed non-zero eccentricities be explained by the perturbation from the detected companions? To see this, I perform a similar set of simulations as in Section 6.1 setting $e_{1}=\omega_{1}=0$ initially, while sampling the other parameters from the posterior. Note that this experiment is applicable to the disk migration case as well, whose outcome is also a short-period planet on a circular orbit.

Figure 25 summarizes the maximum inner eccentricities achieved during the $10\,\mathrm{Myr}$ evolution against the current mutual inclination. The maximum value exceeds the current best-fit value (horizontal dashed line) for $14\%$ and $76\%$ of the solutions that did not lead to tidal disruption in the Kepler-448 and Kepler-693 systems, respectively. Thus the current architecture can indeed be compatible with the initially circular orbits. The sequences of the maximum $e_{1}$ ($e_{\rm 1,max}$) and initial $i_{21}$ ($i_{\rm 21,init}$) roughly follows the relation $\sqrt{1-e_{\rm 1,max}^{2}}\sqrt{3/5}=\cos i_{\rm 21,init}$.

Considering the arguments in Section 6.2, it is also conceivable that the inner WJs were formed into circular orbits at larger semi-major axes than observed now (but still inside the snow line) and have migrated to the current orbits via the tidal migration. In this case, excitation of eccentricity should have been easier because the gravitational interaction with the companion was initially stronger. This process might serve as yet another path of HJ formation: some HJs may have been isolated WJs formed in situ or via disk migration, whose orbits were shrunk due to the tidal high-eccentricity migration driven by a close companion.

7 Discussions

7.1 Follow-up Observations

While the future observations of transit times and durations will surely improve the constraint on the outer period and eccentricity, I did not find any systematic dependence of the future TTV and TDV evolutions on the current mutual orbital inclination at least within about 10 years. Moreover, it is challenging to observe a full transit from the ground, due to the long transit duration ($\simeq 7\,\mathrm{hr}$) for Kepler-448b and faintness ($V=16.8$) of the host star for Kepler-693b.

How about RV observations? Based on the constraints from TTVs and TDVs, the RV semi-amplitude due to the outer companion is expected to be $2.7^{+0.7}_{-0.6}\times 10^{2}\,\mathrm{m\,s^{-1}}$ for Kepler-448 and $2.7^{+0.9}_{-0.6}\,\mathrm{km\,s^{-1}}$ for Kepler-693. The variation may be observable for Kepler-693 with a large telescope, while the detection is implausible for Kepler-448 with $v\sin i_{\star}\simeq 60\,\mathrm{km\,s^{-1}}$.

Instead, Gaia astrometry (Perryman et al., 2001) is promising to detect the outer binary motion and better determine the mutual inclination. Since the misalignment in the sky plane $\Omega_{21}$ is dynamically well constrained in both systems (and the inner planets are transiting), inclinations of the outer orbits $i_{2}$, if measured, significantly improve the constraint on $i_{21}$. Based on the dynamical modeling, the expected astrometric displacements due to the outer companions are $1.4^{+1.2}_{-0.5}\times 10^{2}\,\mu\mathrm{as}$ for Kepler-448 with $V=11.4$ (assuming the distance $d=426\pm 40\,\mathrm{pc}$ from Bourrier et al., 2015) and $4.4^{+3.4}_{-1.5}\times 10^{2}\,\mu\mathrm{as}$ for Kepler-693 with $V=16.8$ (assuming $d=1110\,\mathrm{pc}$ from the isochrone fit). They are both well above the astrometric precision expected after the nominal five-year mission (Perryman et al., 2014).

7.2 Frequency of Close and Massive Companions to WJs

It is admittedly difficult to evaluate the completeness of our TTV search due to the complex dependence of the signal on the system parameters. Nevertheless, our detection of the close and massive companions in two systems, among the sample of $23$ WJs, suggests such an architecture is not extremely uncommon. We also need to take into account that the two systems would not have been detected if the periastron passage did not occur during the Kepler mission. The ratios of the outer orbital periods ($P_{2}\sim 2500\,\mathrm{days}$ for Kepler-448c and $P_{2}\sim 1800\,\mathrm{days}$ for Kepler-693c) to the Kepler observing duration ($\sim 1400\,\mathrm{days}$) suggest that there may be one or two more WJs with similar companions in the Kepler sample that have evaded our search. This crude estimate seems compatible with the theoretical argument by Petrovich & Tremaine (2016) that roughly $20\%$ of WJs can be accounted for by high-eccentricity migration, although the observed architectures of our systems may not support the migration from $\gtrsim 1\,\mathrm{au}$ via this process as assumed in Petrovich & Tremaine (2016).

8 Summary and Conclusion

This paper reported the discovery of close companions to two transiting WJs via their TTVs and TDVs. The companions have masses comparable to a brown dwarf or a low-mass star ($22^{+7}_{-5}\,M_{\rm Jup}$ and $150^{+60}_{-40}\,M_{\rm Jup}$), and they are on highly eccentric orbits ($e\gtrsim 0.5$) with small periastron distances ($\simeq 1.5\,\mathrm{au}$). For the companion of Kepler-693b, a large mutual orbital inclination ($\simeq 50\arcdeg$) with respect to the inner planetary orbit is indicated by TDVs, while the constraint on the mutual inclination is weak for the Kepler-448 system. They are among the few systems with constraints on mutual inclinations, and that inferred for the Kepler-693 system is the largest ever determined dynamically for planetary systems. The value is indeed large enough for the eccentricity oscillation via the Kozai mechanism to occur: more than $90\%$ of the solutions inferred for Kepler-693b (and some $10\%$ for Kepler-448b) imply that the inner WJs’ eccentricities exhibit a large enough oscillations for tidal dissipation to affect the inner orbits significantly, by bringing $a(1-e^{2})$ less than $0.1\,\mathrm{au}$. The corresponding migration timescales can be compatible with the hypothesis that the inner WJs are tidally migrating to evolve into HJs, for a wide range of viscous timescales.

The architectures of the two systems support the scenario that eccentric WJs are currently undergoing eccentricity oscillation induced by a close companion and experiencing the slow tidal migration where the orbit shrinks only at the high-eccentricity phase. On the other hand, the origin of the current highly eccentric/inclined configuration is still unclear. Specifically, they may not fit into the classical picture that close-in gas giants have migrated from beyond the snow line, given the close and (sub-)stellar nature of the outer companions. Formation of gas giants within the snow line onto circular orbits, followed by eccentricity excitation by the companion, therefore seems another viable option to be pursued, although the companion may instead have arrived at the current orbit after disk migration of the inner WJ. Regardless of the origin of the current configuration, the long-term evolution simulation demonstrates a new pathway of HJ formation via “high-eccentricity” migration of a WJ formed in situ or via disk migration.

The author is grateful to the entire Kepler team for making the revolutionary data available. The author thanks Hajime Kawahara, Yasushi Suto, Yuxin Xue, Saul Rappaport, Josh Winn, Sho Uehara, Sarah Ballard, Gongjie Li, Tim Morton, Dong Lai, Kaloyan Penev, Adrian Price-Whelan, Chelsea Huang, George Zhou, Bekki Dawson, Cristobal Petrovich, Amaury Triaud, Kengo Tomida, Kaitlin Kratter, Eve Ostriker, and Ji-Ming Shi for helpful conversations and insightful discussions; and Tamás Borkovits for answering questions on the analytic ETV formula. Comments on the manuscript from Dong Lai, Josh Winn, and an anonymous referee were also very helpful in improving the paper and enriching its content. K.M. acknowledges the support by JSPS Research Fellowships for Young Scientists (No. 26-7182) and also by the Leading Graduate Course for Frontiers of Mathematical Sciences and Physics. This work was performed in part under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

Appendix A Analytic TTV Formula for Hierarchical Triple Systems

In a part of the TTV analysis, we adopt the analytic timing formula for hierarchical three-body systems by Borkovits et al. (2015), which takes into account the eccentricities of both inner and outer orbits and arbitrary mutual inclination between them (see also Borkovits et al., 2011, 2015, 2016, for its applications). Among various effects that could possibly affect the observed transit times, we include two of the most important effects: light-travel time effect (LTTE) and $P_{2}$-timescale dynamical effects up to the quadrupole order. The other effects including the octupole-level dynamical effects and short-term perturbations are at least an order-of-magnitude smaller than the quadrupole terms and thus neglected (see Borkovits et al., 2015, for a complete discussion of these effects). Note that, in this appendix, $+Z$-axis is taken away from the observer’s direction; this definition is opposite to the one in the main text, and so arguments of periastron in the formulae below differ by $\pi$ from the ones in the main text.

A.1 Formula

We model the timing variations from the linear ephemeris due the LTTE effect and the $P_{2}$-timescale dynamical effect. The LTTE term is given by

[TABLE]

where $v_{2}$ is the true anomaly and $u_{2}\equiv\omega_{2}+v_{2}$ is the true longitude. The $P_{2}$-timescale dynamical effect is modeled up to the quadrupole order as follows (Borkovits et al., 2015, equation 5):

[TABLE]

Here the overall amplitude is fixed by the factor

[TABLE]

and each TTV component is given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The angles $n_{1}$ and $n_{2}$ are the directions of $Z>0$ part of the line of intersection of the inner and outer orbits measured from the ascending nodes of the two orbits, defined as in figure 1 of Borkovits et al. (2015) between $[0,\pi]$, and $l_{2}$ is the mean anomaly. Note that $K_{11}(e_{1},\pi-\omega_{1})=K_{11}(e_{1},\omega_{1})$ and $K_{12}(e_{1},\pi-\omega_{1})=-K_{12}(e_{1},\omega_{1})$.

A.2 Conversion of the Angles

In the main body of the paper, we did not use the physically most natural parametrization of the angles because it complicates the assignment of the prior in the MultiNest fitting. Here we summarize how to relate our set of fitted angles ($i_{1}$, $i_{2}$, $\Omega_{21}$) to that of physical angles ($i_{21}$, $n_{1}$, $n_{2}$) in the analytic formula by Borkovits et al. (2015), since this case is not covered in their Appendix D.

The mutual inclination $i_{21}$ can be computed as usual:

[TABLE]

Note that $\sin i_{21}\neq 0$ for $\Omega_{21}\neq 0$ since $(i_{1},i_{2})\neq(0,0)$ nor $(\pi,\pi)$ for transiting systems as discussed in this paper. For $\cos\Omega_{21}=\pm 1$, the above equation reduces to

[TABLE]

Let us first consider the case when $\sin i_{21}\neq 0$. If $\sin\Omega_{21}\neq 0$, the sine and cosine rules of the spherical trigonometry yield

[TABLE]

and

[TABLE]

[TABLE]

for $\sin\Omega_{21}<0$ case. For $\sin\Omega_{21}>0$, the cosine rules are replaced by

[TABLE]

[TABLE]

In fact, the two cases can be written in a single form as

[TABLE]

[TABLE]

Thus, for a non-zero mutual inclination we have

[TABLE]

and

[TABLE]

[TABLE]

If $\sin i_{21}=0$, we may just define $n_{2}-n_{1}=0$ for $i_{21}=0$ and $n_{2}+n_{1}=\pi$ for $i_{21}=\pi$ to correctly compute the non-vanishing terms of TTVs. Other than $i_{1}=i_{2}=0$ or $\pi$, this includes the two cases: (i) $i_{1}-i_{2}=0$ and $\Omega_{21}=0$ and (ii) $i_{1}+i_{2}=\pi$ and $\Omega_{21}=\pi$. As shown in Table 7, we have $n_{1}-n_{2}=0$ for either $i_{2}-i_{1}\to\pm 0$ in case (i). Similarly in case (ii), we have $n_{1}+n_{2}=\pi$ for either $i_{2}+i_{1}-\pi\to\pm 0$. Although the other combination is ambiguous, it appears only in the vanishing terms of the TTV formula.

Appendix B Results of Analytic and Numerical TTV Analyses

B.1 Comparison of the Two Results

For TTVs, I performed both an analytic fit with a wider prior range and a numerical fit with a narrower prior range. As shown in Tables 8 and 9, I found consistent posteriors from the two analyses; this agreement validates our numerical procedure. The best-fit models are basically indistinguishable from those in Figures 7 and 9.

B.2 Decomposition of the TTV Solutions Using the Analytic Formula

The analytic formula allows us to understand how each physical effect in Equation A2 contributes to the observed TTVs. Figure 28 shows the decomposed signals for each of the (i) “LTTE” $\Delta_{\rm LTTE}$, (ii) “tidal” $\delta_{\rm tidal}$, (iii) “eccentric” $\delta_{\rm ecc1}+\delta_{\rm ecc2}$, and (iv) “non-coplanar” $\delta_{\rm noncopl}$ terms for $10$ solutions randomly sampled from the posterior obtained in the previous section. The plot shows that the $\delta_{\rm ecc}$ terms play a crucial role in producing the short-term feature, especially for Kepler-448b.

Appendix C Corner Plots for the Posteriors from Dynamical Analyses

Figures 29 and 30 show the corner plots of the posterior distributions obtained from the dynamical TTV and TDV analyses in Sections 4 and 5. The figures are generated using corner.py by Foreman-Mackey (2016).