THE DYNAMICS OF THE MULTI-PLANET SYSTEM ORBITING KEPLER-56

Over the past few years, measurements of the sky-projected obliquity of exoplanets have found that large obliquities and even retrograde systems are common among hot Jupiters (e.g., Fabrycky & Winn 2009; Triaud et al. 2010; Morton & Johnson 2011; Moutou et al. 2011; Albrecht et al. 2012; Hébrard et al. 2013). Recently, Hirano et al. (2012), Sanchis-Ojeda et al. (2012), Albrecht et al. (2013), Chaplin et al. (2013), and Van Eylen et al. (2014) have measured the obliquity of six transiting multi-planet systems discovered by the NASA Kepler mission, and found they all have low obliquities. However, Huber et al. (2013), using asteroseismology, showed that large obliquities are not confined to Hot Jupiter systems. In fact Kepler-56 has two low mass inner planets whose orbit normal is tilted with respect to the stellar spin axis.

Several mechanisms have been suggested to explain the formation of misaligned planets. These theories can be divided into two categories. The first is based on tilting the orientation of an inner planet compared to the stellar spin axis. This category includes scattering and secular dynamical effects between a planet and a companion, or other planets in the system that can produce large obliquities (e.g., Fabrycky & Tremaine 2007; Chatterjee et al. 2008; Nagasawa et al. 2008; Naoz et al. 2011, 2012, 2013; Wu & Lithwick 2011; Li et al. 2014a, 2014b; Valsecchi & Rasio 2014a, 2014b). These mechanisms predict that an inner planet with a large obliquity has an outer perturber which is inclined with respect to the plane of the inner planet, the perturber can be either a stellar companion or a planet, or even multiple planets. In the second category, planets move inward from their birthplaces beyond the snow line by migrating inward through the protoplanetary disk (e.g., Lin & Papaloizou 1986; Masset & Papaloizou 2003). Large obliquities can then be produced either by tilting the stellar spin axis with respect to the orbital angular momentum (e.g., Winn et al. 2010; Lai et al. 2011; Rogers et al. 2012, 2013; Spalding & Batygin 2014), or by tilting the protoplanetary disk (Bate et al. 2010; Batygin 2012). This second category of models predicts that the various planets in a system should lie roughly in the same plane since they were confined to the same flattened disk.

Here we focus on the dynamical mechanism that produced the large obliquities in the Kepler-56 planetary system. Most of the theoretical studies investigating large obliquities focused on Hot Jupiters, mainly because these were observed to have large obliquities. The underlining physics of producing a misalignment in the presence of a perturber is very similar. Thus, such studies are relevant for investigating the Kepler-56 system (as we will show below).

Kepler multiple systems are typically packed, small sized (∼1–10 R_⊕; e.g., Lissauer et al. 2011; Swift et al. 2013), and close-in (∼1–100 days; e.g., Steffen & Farr 2013) systems. At face value these configurations may indicate that dynamical and secular processes are suppressed, since these systems better resemble the theoretical outcome of planet migration in the protoplanetary disk, given their low mutual inclinations (Lissauer et al. 2011; Fang & Margot 2012). Therefore, a large obliquity in a multi-planet system may be used as a laboratory to test the two categories of models summarized above. In other words, since it seems that these planets form in a disk, a tilt of the protoplanetary disk or of the star will cause the multiple planets to show the same obliquity.

Kepler-56 is an evolved star at the base of the red giant branch in the HR diagram with m_⋆ = 1.32 M_☉, R_⋆ = 4.23 R_☉, and an age of 3.5 Gyr (Huber et al. 2013). Furthermore, Huber et al. (2013) showed that the innermost planet (m_b = 0.07 M_J, R_b = 0.65 R_J, hereafter planet "b") has a period of 10.5 days, and a period of 21.4 days for the other planet (m_c = 0.57 M_J, R_c = 0.92 R_J, hereafter planet "c"). The mutual inclination between these two planets is measured to be <5°. Kepler-56 is an interesting system as it raises many questions regarding its formation and future evolution. Most importantly, Huber et al. (2013) measured the obliquity of the system using asteroseismology and placed a lower limit on the true obliquity of the two inner planets of ψ > 37°. The dynamical analysis of Huber et al. (2013) favors the scattering and later torquing scenario.

Here we use Kepler-56's current observations to compute the probability distribution for its obliquity. (The Huber et al. 2013 reported observations already give enough information to calculate such a distribution.) This enables us to also put strong constraints on the probability distribution of the outer planet's inclination with respect to the innermost two. Furthermore, we estimate that the two inner planets will be engulfed in ∼129 Myr and ≲ 155 Myr, respectively. The engulfment of the inner planets is consistent with the deficit in short period planets around retired A stars (e.g., Johnson et al. 2007; Sato et al. 2008; Bowler et al. 2010; Schlaufman & Winn 2013).

The paper is structured as follows. We calculate the obliquity distribution function from observations, and show that the current observations give more information than just a lower limit (Section 2). We then discuss the current obliquity precession as a function of the system initial conditions (Section 3.2) and show that combining the physical understanding and the observed distribution, we can infer the outermost planet orbital inclination with respect to the innermost two as a function of the initial configuration (Section 3.3). We also calculate the orbit and obliquity future evolution as the star further ascends the giant branch (Section 4). We finally offer our discussion (Section 5).

Huber et al. (2013) analyzed the stellar oscillations observed in the Kepler photometry and used the splitting of the observed oscillation frequencies to measure the inclination between the stellar spin axis and the line of sight, finding $i_{\rm ls}^\star = 47^\circ \pm 6$. With the transit photometry, Huber et al. also measured the inclination of the inner planet's orbit with respect to the line of sight, finding $i_{\rm ls}^{b} =83{.\!\!^\circ}84$$^{+0.26}_{-0.25}$. Together, these angles place a lower limit on the three-dimensional angle between the stellar spin axis and planetary orbital plane of ψ > 37°.

The angle between the normal of the orbit and the stellar spin is not simply $i_{\rm ls}^{b}+i_{\rm ls}^\star$ since, for example, the angle $i_{\rm ls}^{b}$ can have different values on the sky plane (different values of α as shown in Figure 1).

Zoom In Zoom Out Reset image size

In this simple geometrical configuration (see Figure 1, left panel) and defining L_in and S as the angular momentum of the innermost orbit and stellar spin, respectively, the obliquity is defined by the scalar product between the three-dimensional spin axis unit vector ${{\bf n}}_s = {\bf S}/S=(\sin i_{\rm ls}^\star,0,\cos i_{\rm ls}^\star)$ and the three-dimensional normal to the innermost orbit ${\bf n}_{\rm in}={\bf L}_{\rm in} / L=(\sin i_{\rm ls}^{b},0,\cos i_{\rm ls}^{b})$, in random orientation with each other:

$$\begin{equation} \cos \psi ={{\bf n}}_s\cdot {\bf R}_{\rm ls} {{\bf n}}_{\rm in}. \end{equation} \tag{ 1 }$$

Here

$$\begin{equation} {\bf R}_{\rm ls} (\alpha) = \left(\begin{array}{@{}ccc@{}} \cos \alpha & -\sin \alpha & 0 \\ \sin \alpha & \cos \alpha & 0 \\ 0 & 0 & 1 \end{array}\right) \end{equation} \tag{ 2 }$$

is the rotation matrix in the azimuthal direction around the line of sight. We assume that α, the angle between the stellar spin and the orbital angular momentum in the azimuthal direction around the line of sight, is uniformly distributed. It is sufficient to multiply only once by the rotation matrix, with the random angle. Therefore, from Equation (1) we can estimate the cumulative distribution function of ψ. As shown in the right panel of Figure 1, the lower limit on ψ is of course the same one found by Huber et al. (2013), i.e., ψ > 37°, but an upper limit of 131° also exists and both these values have the same probability, which is larger than the probability of the angles in the range of 37° < ψ < 131°. We use ψ_obs to denote the observationally constrained value of ψ. Note that due to the degeneracy in the asteroseismology measurements, $i_{\rm ls}^\star$ could also be 133° ± 6. Setting $i_{\rm ls}^\star = 133^\circ$, ψ_obs is in the range of 49° < ψ < 143° (see the dashed line in Figure 1). Therefore, adding these two pieces together, the distribution of ψ_obs is symmetric over 37° < ψ < 143°. This distorts ψ_obs only slightly, because $i_{\rm ls}^{b}$ is almost 90°. Accordingly, we adopt $i_{\rm ls}^\star = 47^\circ \pm 6$, and have ψ_obs constrained in the range of 37° < ψ < 131° for the following discussion.

3.1. Overview of the System Architecture

3.2. Dynamics of Kepler-56

In the presence of a tilted outer orbit with inclination i_mut, the two inner planets will precess around the total angular momentum vector. Note that the precession of the orbit due to the oblateness of the star is negligible in this case. The torque felt by planet "b" due to stellar oblateness⁵ is more than two orders of magnitude smaller than the torque due to planet "c" (see Tremaine et al. 2009 and Tamayo et al. 2013). Therefore, the orbital evolution is not affected by the torque due to the stellar oblateness, and the system is in the "pure orbital regime" (Boué & Fabrycky 2014). We thus obtain the orbital evolution from an N-body simulation.

The obliquity angle is defined with respect to the innermost planet's orbital angular momentum, L_b. Thus, a natural coordinate choice for the spin is the Laplace–Runge–Lenz $({\hat{\bf q}}_b,{\hat{\bf h}}_b,{\hat{\bf e}}_b)$. Here, ${\hat{\bf e}}_b$ is the eccentricity vector (whose direction is toward the pericenter of planet "b" orbit), ${\hat{\bf h}}_b$ is the unit vector parallel to the orbital angular momentum of planet "b" (the vector h_b is the specific angular momentum vector, i.e., L_b = m_⋆m_b/(m_⋆ + m_b)h_b), and ${\hat{\bf q}}_b$ completes the right-hand triad of unit vectors. In this notation the precession of the stellar spin, S = (S_e, S_q, S_h), due to one planet is simply (Eggleton et al. 1998)

$$\begin{equation} \frac{d{\bf S}}{dt}_{\rm prec,a}={\bf S}\times {\bf K}_b + \frac{m_\star m_b}{m_\star +m_b} h_{b}/I_2 (-\tilde{Y}_b {\hat{\bf e}}_b+\tilde{X}_b{\hat{\bf q}}_b+\tilde{W}_b{\hat{\bf h}}_b), \end{equation} \tag{ 3 }$$

where $h_{b}=[G(m_\star +m_b)a_b(1+e_b^2)]^{1/2}$, G is the gravitational constant, and K_b = (X_b, Y_b, Z_b) represents the precession due to the orbital evolution:

$$\begin{eqnarray} X_b&=& \frac{di_b}{dt} \cos \omega _b +\frac{d\Omega _b}{dt}\sin \omega _b \sin i_b, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} Y_b&=& -\frac{di_b}{dt}\sin \omega _b + \frac{d\Omega _b}{dt}\cos \omega _b \sin i_b, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} Z_b&=& \frac{d \omega _b}{dt} + \frac{d\Omega _b}{dt} \cos i_b, \end{eqnarray} \tag{ 6 }$$

and $\tilde{X}_b$, $\tilde{Y}_b$, and $\tilde{W}_b$ represent the torque due to the stellar oblateness and the tidal dissipation:

$$\begin{eqnarray} \tilde{X}_b&=& -\frac{m_b k_\star R_\star ^5}{\mu \dot{l} a_b^5} \frac{S_h S_e}{(1-e_b^2)^2}\nonumber \\ &&-\frac{S_q}{2 \dot{l} t_{F_\star }} \frac{1+(9/2)e_b^2 + (5/8)e_b^4}{(1-e_b^2)^5}, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \tilde{Y}_b&=& -\frac{m_b k_\star R_\star ^5}{\mu \dot{l} a_b^5} \frac{S_h S_q}{(1-e_b^2)^2} \nonumber\\ &&- \frac{S_e}{2 \dot{l}t_{F_\star }}\frac{1+(9/2)e_b^2 + (5/8)e_b^4}{(1-e_b^2)^5}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \tilde{W}_b &=& \frac{1}{t_{F_\star }}\left[\frac{1+(15/2)e_b^2+(45/8)e_b^4+(5/16)e_b^6}{(1-e_b^2)^{13/2}}\right. \nonumber\\ &&\left.- \frac{S_h}{\dot{l}}\frac{1+3e_b^2+(3/8)e_b^4}{(1-e^2)^5}\right], \end{eqnarray} \tag{ 9 }$$

where $\dot{l} = \sqrt{G m_\star / a_b^3}$, and

$$\begin{equation} t_{F_\star } = \frac{t_{V_\star }}{9} \frac{m_\star ^2}{(m_\star +m_b)m_b} \left(\frac{a_b}{R_\star }\right)^8 \frac{1}{(1+2k_2)^2}. \end{equation} \tag{ 10 }$$

To calculate the orbital evolution due to the orbital precession (the K_b term), we take the time evolution of ω, Ω, and i of planets "b" directly from the N-body integration. This dominates the obliquity variation. The tidal effects are negligible until planet b is almost engulfed (see the discussion on the future evolution of Kepler-56 in Section 4). The timescale for the evolution of planet b's orbital separation due tidal dissipation in the star is defined in terms of the stellar viscous timescale $t_{V_\star }$. $t_{V_\star }$ is set to be 50 yr and kept constant, where $t_{V_\star }$ corresponds to Q ∼ 10⁶ for a 10 day orbit. The parameter k₂ is the apsidal precession constant, which is related to the Love parameter k_L via k₂ = 2 k_L (a similar equation exists for planet "b" and "c"). Note that the effects of tides in the planets are negligible. In fact, assuming a viscous timescale corresponding to Q = 12 and 10⁵ for planet "b" and "c" (Murray & Dermott 2000), respectively, the small planet radii yield much longer tidal timescales (see Equation (10)). In any case, the unconstrained nature of exoplanets makes it difficult to conclude how their tidal coefficients evolve.

Equations (3)–(10) imply that the time evolution of i_mut (and thus ψ) depends on the initial system's configuration. This can be constrained from the observed obliquity distribution. In Figure 2 we show the evolution of ψ assuming two possible initial configurations: S parallel to L_in and i_mut = 20° (left, hereafter "S || L_in" scenario), and S parallel to L_tot and i_mut = 40° (right, hereafter "S || L_tot" scenario), where L_in and $\bf L_{\rm tot}$ are the orbital angular momentum of the inner two planets and the total orbital angular momentum, respectively. We show below that these values for i_mut give a misalignment of at least 37° during the evolution (the minimum value constrained observationally).

Zoom In Zoom Out Reset image size

In the S || L_in scenario, ψ oscillates between well-aligned (ψ = 0°) and ∼2 × i_mut (∼382). In this case, we postulate that the system formed initially in a disk and planet "d" was perhaps scattered to large inclinations (e.g., Rasio & Ford 1996), causing the obliquity angle to precess between 0° and ∼2 × i_mut. Another possible case for this configuration is accretion of material onto the protoplanetary disk, which can tilt the outer parts of the disk and the total angular momentum (Bate et al. 2010; Tremaine 2011; Thies et al. 2011). Therefore, in the S || L_in scenario, ψ ∼ 40° can be produced by an initial inclination i_mut > 20. Note that a retrograde configuration with i_mut = 160° can also produce ψ ∼ 40°.

In the S || L_tot scenario, ψ remains close to the initial value. This configuration could have occurred if the inner parts of the disk were warped perhaps due to magnetic interactions with the inner disk edge (e.g., Lai et al. 2011). Therefore in the S || L_tot scenario, ψ ∼ 40° can be produced by an initial i_mut = 40°.

We show below that for the S || L_in scenario ψ is more likely to be detected in the maximum (at ∼382) where the derivative is closer to zero. For each possible obliquity value $\tilde{\psi }\in (0^\circ,180^\circ)$, we derived a cumulative distribution function of the mutual inclination, where CDF$(\tilde{\psi }|i_{\rm mut})=\Delta t(\psi <\tilde{\psi }|i_{\rm mut})/t$, where Δt is the time interval. This quantity will be used below to estimate the probability distribution of the system configuration for the actual observations. We run 35 N-body runs, for an array of initial inclinations i_mut between 5° and 175°, and calculate the cumulative probability for the two scenarios.

3.3. Inferring the Inclination Distribution Function from Observations

When the spin-orbit misalignment is due to the dynamical interaction between the planets, the obliquity distribution function derived from observations (Section 2, Figure 1) can be used to place strong constrains on the mutual inclination between the inner planets and planet "d," i.e., i_mut. We calculate the conditional probability distribution of i_mut given the observed distribution ψ_obs, i.e., p(i_mut|ψ_obs). This posterior probability can be written as

$$\begin{equation} p(i_{\rm mut}|\psi _{\rm obs})=\frac{p(\psi _{\rm obs}|i_{\rm mut}) p(i_{\rm mut})}{p(\psi _{\rm obs})}, \end{equation} \tag{ 11 }$$

where p(ψ_obs) is a normalization term, which we disregard because the shape of the distribution is of larger significance than the absolute probability here, and the absolute probability is out of the scope of this paper.

Furthermore, we use the distribution function of planet "d" line of sight inclination, $i_{\rm ls}^d$, to estimate the prior probability, p(i_mut). Note, that the actual value of $m_d \sin i_{\rm ls}^d$ only affects the normalization of the probability, but since we care about the shape of the probability we can ignore this. Note that if we assume the outer orbit to be isotropically distributed, the probability density function for $i_{\rm ls}^d$ takes the form of $\sin {i}_{\rm ls}^d$. This suggests that the most probable value for $i_{\rm ls}^d$ is 90°.

Following Ho & Turner (2011) we calculate the probability $p(i_{\rm ls}^d)$ assuming Cumming et al. (2008) mass function for m_d (see Figure 3). Note that the distribution in Cumming et al. (2008) is for msin i, not m. However, since the power law index is large, we use this power law for the mass distribution according to Ho & Turner (2011). The angle we are actually interested in is the angle between the normal to the outer orbit n_out and the normal to the inner orbit n_in. While $i_{\rm ls}^{b}$ has been measured to be 8384 (Huber et al. 2013) we have no information about the orientation of these two vectors on the plane of the sky. Consider first the case where the three-dimensional normal to the outer orbit, n_d, lies in the plane defined by n_in and the line of sight. This yields a simple relation between the different angles, i.e., $i_{\rm mut}=i_{\rm ls}^{b}-i_{\rm ls}^d$. Therefore, $p(i_{\rm mut})=p(i_{\rm ls}^{b}-i_{\rm ls}^d)$, where the latter is calculated from $p(\sin i_{\rm ls}^d)$ following Ho & Turner (2011). Their mass distribution function yields small m_d (compared to the measured $m_d \sin i_{\rm ls}^d$), thus $\sin i_{\rm ls}^d$ is more likely to be close to its maximum of 1. This suggests angles near 90 deg for $i_{\rm ls}^d$, which implies that $i_{\rm mut}=i_{\rm ls}^{b}-i_{\rm ls}^d$ is more likely to have a small value.

Zoom In Zoom Out Reset image size

However, another possible prior is that n_d has a random orientation (similar to the configuration depicted in the left side of Figure 1). Thus, as in Section 2, we multiply the normal to the orbit with the rotation matrix in Equation (2) assuming a random azimuthal angle α, i.e.,

$$\begin{equation} \cos i_{\rm mut} = {\bf n}_{\rm out}\left(i_{\rm ls}^d\right) \cdot {\bf R}_{\rm ls}(\alpha) {\bf n}_{\rm in}\left(i_{\rm ls}^{b}\right), \end{equation} \tag{ 12 }$$

where ${\bf n}_{\rm out}(i_{\rm ls}^d)$ is chosen with $p(\sin i_{\rm ls}^d)$ distribution, which gives p(i_mut). This prior also gives a high probability for large values of i_mut, as this case covers large parts of the parameter space. Below we consider these two cases.

The probability of ψ_obs for a given i_mut, i.e., p(ψ_obs|i_mut) can be calculated from

$$\begin{equation} p(\psi _{\rm obs}|i_{\rm mut}) = \int _0^\infty p(\tilde{\psi }|i_{\rm mut}) p_{\rm obs}(\tilde{\psi }) d\tilde{\psi }, \end{equation} \tag{ 13 }$$

where p_obs(ψ) was computed in Section 2, Figure 1. The probability p(ψ|i_mut) is calculated from theory for the two different cases, i.e., S || L_tot and S || L_in. In the discrete description we calculate the probability distribution of i_mut for each ψ. This can be easily derived from the cumulative distribution function calculated in Section 2, Figure 1.

Using Equations (11)–(13) we can find the mutual inclination probability function given the observed obliquity distribution. This is depicted in Figure 4, bottom panels. We consider the two initial configurations scenarios, i.e., S || L_tot and S || L_in, and the two $p(i_{\rm ls}^d)$ cases, i.e., n_out random along line of sight (right panels) and n_out in the same plane as n_in (left panels). Since the obliquity distribution function derived from observations has two high probability peaks (ψ = 37° and ψ = 131°), the possible i_mut values that can produce such distribution function also have two peaks. In the case of S || L_in, the double peak distribution is also probable since the precession of a retrograde orbit can as well produce this configuration. Note that if we also consider the case when $i_{\rm ls}^\star = 133^\circ \pm 6$ (due to the degeneracy in the asteroseismology measurements), ψ_obs is symmetric, and p(i_mut|ψ_obs) would also be symmetric.

Zoom In Zoom Out Reset image size

Interestingly, better observations may help constraining i_mut but will not disentangle the degeneracy between the S || L_tot and S || L_in cases. We show this in the top panels of Figure 4, where we consider an example of ψ = 40 ± 3°. In the S || L_tot scenario, the symmetry is broken, since there is a direct link between the obliquity and i_mut in this case, as seen from the right panel of Figure 2. Note that the two different p(i_mut) cases produce slight differences in the probability peak. Assuming that n_out and n_in are coplanar produces a decreasing probability toward i_mut ∼ 45°, as in this case near polar configurations are less likely. On the other hand, assuming a random orientation for n_out produces an increasing probability toward the larger i_mut values. In fact, as mentioned above, this case yields a larger parameter space for near polar configurations. Having a precise observation also improves the i_mut estimation for the S || L_in but the double peak probability remains, because the same obliquity can be reached in a prograde and retrograde configurations. The degeneracy can be broken only for the case where S || L_tot, a more precise measurement of ψ will be available. This can be seen in the top panels of Figure 4, where ψ_bin represents 37° < ψ < 43°.

So far, we have assumed two possible priors for p(i_mut). These represent two extreme possibilities, one that favors low mutual inclinations and one that favors large values. The truth may lay in between. Thus, we have tested the possibility that n_d is randomly oriented within a small interval as a prior (see the left side of Figure 1, where α is now confined to a certain interval). In this case, different from what Figure 3 shows, we assumed an initial tilt of 37° between the stellar spin axis and the angular momentum of the inner orbit. This way, we consider the possibility that the source of the obliquity is not dynamical. We find that for α ≳ 10° Equation (11) and the observed obliquity distribution favors large mutual inclinations. In other words, the three planets will be aligned, and the observations will be consistent with tilt of the star or the disk in the migration scenario, if the random angle α < 10°.

Here we focus on the fate of the innermost planet and the future evolution of the obliquity as a result of tidal dissipation in the star and stellar evolution. We compute a detailed model of the host star with the publicly available stellar evolution code MESA (version 4798 Paxton et al. 2011, 2013). Specifically, we follow Huber et al. (2013) and consider a star with an initial mass and metallicity of 1.32 M_☉ and Z = 0.032, respectively. We evolve the stellar model with the same physical assumptions adopted in Valsecchi & Rasio (2014b). Briefly, we account for stellar wind mass loss following the test suite example provided with MESA for the evolution of a 1 M_☉ star, and we set the mixing length α_MLT parameter to 1.92, following the MESA star Standard Solar Model (Paxton et al. 2011, Table 10). Note that the mass loss is negligible in this case, since the star is only slightly evolved (as shown in Figure 5). This negligible mass loss explains why the planets' orbits are significantly expanding, different from the case of the Earth when the Sun evolves into a red giant. The model agrees with the observationally inferred stellar mass, radius, and effective temperature (within 1 σ) at 4.418 Gyr. The latter is consistent with the age quoted by Huber et al. (2013) within 1σ (3.5 ± 1.3 Gyr).

Zoom In Zoom Out Reset image size

The advanced evolutionary stage of the host star (which is off its main sequence) poses the interesting possibility that, if Kepler-56 is similar to other Kepler multi-planet systems, it may have had planets that were engulfed as the star expanded (such possibilities have been investigated in the literature, Bear & Soker 2011, 2012). If this is the case, the observed small stellar rotation rate suggests that the host star in Kepler-56 did not engulf a large planet. In fact, to increase the stellar spin by more than 10%, the engulfed planet should have had a mass larger than 0.6 M_J (neglecting the possibility of core-envelope decoupling, see, e.g., Teitler & Königl 2014). However, we note that magnetic braking, stellar winds, and the expansion of the star as a result of natural stellar evolution might all contribute to spin down after engulfment. Nevertheless, it also seems unlikely (but not impossible) that a very massive planet could have migrated to the innermost configuration, with two lighter planets outside (planets "b" and "c") which also supports the notion that no inner planet was engulfed.

The inner planets' orbital evolution is affected by tides, whose efficiency changes as the star evolves. In the left panels of Figure 5 we show the forward evolution of the stellar mass, radius, and Love parameter. The latter was computed following Valsecchi et al. (2012).

Both the specific angular momentum (h_b) and the eccentricity undergo tidal dissipation, which leads to circularization and orbital shrinking. Following Eggleton et al. (1998) we have

$$\begin{eqnarray} \frac{de_b}{dt}&=&{-}\tilde{V}_b e_b, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} \frac{dh_b}{dt} &=& {-}\tilde{W}_b h_b. \end{eqnarray} \tag{ 15 }$$

The parameters $\tilde{W}_b$ and $\tilde{V}_b$ are the dissipation coefficients (see Eggleton et al. 1998), where $\tilde{W}_b$ is given by Equation (9), and $\tilde{V}_b$ is defined as:

$$\begin{eqnarray} \tilde{V}_b &=& \frac{9}{t_{F_\star }}\left[\frac{1+(15/4)e_b^2+(15/8)e_b^4+(5/64)e_b^6}{(1-e_b^2)^{13/2}}\right. \nonumber\\ &&\left.-\frac{11S_h}{18\dot{l}}\frac{1+(3/2)e_b^2+(1/8)e_b^4}{(1-e^2)^5}\right]. \end{eqnarray} \tag{ 16 }$$

We compute the evolution of the orbital separation, eccentricity, and inclination, using the extrapolated orbital parameters from the initial direct N-body integration, together with the equations mentioned above. We stop the integration when the innermost planet is engulfed (a_b = R_⋆) and neglect possible mass transfer events between the planet and the star (e.g., Trilling et al. 1998), for simplicity. The evolution is shown in the right two panels of Figure 5. During the first ∼0.1 Gyr of evolution, the star loses about 0.1% of its mass and its radius expands by about 40%. After this stage tidal effects become increasingly important and planet "b" is quickly engulfed. We note that the tidal treatment adopted here does not fully account for how the evolution of the star affects the efficiency of tides (i.e., the stellar viscous timescale $t_{V_\star }$ is kept fixed). However, a more consistent orbital evolution calculation with the method adopted in Valsecchi & Rasio (2014a, 2014b), but only accounting for the evolution of the innermost planet, yields similar results. Note also that the precession due to the stellar oblateness affects the final stages of the evolution (very close to the final engulfment of planet "b"). This occurs because tidal dissipation dominates the dynamics only toward the end of the evolution right before the engulfment, and thus, it does not change the overall orbital dynamics. Moreover, we find that the final semimajor axis of the planet is ∼0.03 AU during the engulfment, which is within twice the Roche limit (the Roche limit is 0.016 AU according to the prescription of Paczyński 1971). This suggests that the planet may be tidally distorted during the engulfment and that the accumulated heat due to tidal dissipation as the planet orbits the star multiple times may increase the chance of tidal disruption (Li & Loeb 2013). Past studies have investigated the engulfment of planets by their host stars (Nordhaus et al. 2010; Bear & Soker 2011; Kaib et al. 2011; Kratter & Perets 2012; Veras et al. 2013; Lillo-Box et al. 2014). Figure 5 shows that the innermost planet will be engulfed in ∼129 Myr. Similarly, the second planet (Kepler 56c) will be engulfed in less than ∼155 Myr.

The tidal evolution of the inner planets affects the stellar spin evolution (Equation (3)). The same equation holds for planet "c" (substituting subscript "b" with "c"). The stellar spin evolves due to the precession of planets "b" and "c," and their tidal torque.

We extrapolate the evolution of their precession directly from the N-body calculation. The evolution of the stellar spin direction and magnitude is shown in Figure 6. In particular, we show the evolution of the obliquity ψ and the angle between the stellar spin and the total angular momentum (ϕ). The magnitude of the spin decreases due to the mass lost and the expansion of the stellar radius (irrespective of the scenario considered). This exercise reveals that the obliquity behavior for the two cases does not vary much as the star evolves. In the S || L_in scenario with an initial i_mut = 20°, the obliquity oscillates between zero and ∼40°, the amplitude slightly decreases, and additional modulations due to tides appears. In the S || L_tot scenario with an initial i_mut = 40°, the obliquity monotonically decreases.

Zoom In Zoom Out Reset image size

We studied the configuration and obliquity of Kepler-56, a multi planet system with two coplanar inner planets that are misaligned with respect to their host star. Two main scenarios were proposed in the literature to explain the large obliquities observed for close-in exoplanets. The first model involves dynamical evolution between multi planetary members or stellar companion (e.g., Winn et al. 2010; Fabrycky & Tremaine 2007; Chatterjee et al. 2008; Nagasawa et al. 2008; Naoz et al. 2011, 2012, 2013; Wu & Lithwick 2011; Li et al. 2014b). The second model proposes disk migration as the main mechanism that controls the planetary configuration, while the star spin axis is tilted with respect to the planets by other mechanisms (e.g., Winn et al. 2010; Lai et al. 2011; Rogers et al. 2012, 2013; Spalding & Batygin 2014). The two scenarios lead to different configurations for the configuration of the planets with respect to each other and the star. The dynamical scenario predicts that large obliquities are associated with an inclined perturber, while the disk-migration scenario predicts aligned planetary systems. We showed that the large obliquity observed in Kepler-56 is consistent with a dynamical nature.

We showed that we can improve the Huber et al. (2013) lower limit on the obliquity (ψ > 37°). Specifically, using the Huber et al. (2013) current observations, we found the probability distribution of the observed true obliquity (see Figure 1). This probability has a large range with two main peaks at ψ = 37° and ψ = 131°. Furthermore, using this probability distribution we gave the probability distribution of the inclination of the third planet with respect to the inner two (i_mut). This is highly dependent on the system's initial conditions. For this reason, we explored two scenarios: S || L_in and S || L_tot. In the former, the initial spin axis of the star was set along the orbital angular momentum of the inner two planets. A possible origin for this configuration is that the system formed initially in a disk and the third Jupiter-like planet was perhaps scattered to a large inclination. Instead, in the S || L_tot scenario, the initial stellar spin was set parallel to the total orbital angular momentum. This initial condition may be ad hoc, and possibly caused by, e.g., magnetic interactions Lai et al. (e.g., 2011) warping the inner parts of the disk. For these two scenarios, we found the mutual inclination probability function for the observed obliquity distribution (see Figure 4, bottom panels). Both configurations have a double peak distributions, with zero probability of having aligned configuration between the two orbits. The degeneracy between the two probability peaks may be broken only for the S || L_tot case, with a more precise measurement of ψ. However, a precise measurement of ψ would not disentangle between the S || L_tot and S || L_in cases, as shown in the top panels of Figure 4.

We finally considered the effect of the stellar evolution on the system's parameters and, specifically, the obliquity. We evolved the host star using MESA (Paxton et al. 2011, 2013) and extrapolated the planet orbital evolution calculated with direct N-body integration (since the latter is rather regular and periodic). We have also included the spin precession and tidal evolution. This exercise revealed that the obliquity behavior for the two cases does not vary significantly as the star evolves. It also shows that planet "b" will be engulfed in ∼129 Myr.

S.N. is supported by NASA through an Einstein Post-doctoral Fellowship awarded by the Chandra X-Ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract PF2-130096. This work was also supported by NASA grant NNX12AI86G at Northwestern University. J.A.J. gratefully acknowledges funding from the Alfred P. Sloan and David & Lucile Packard foundations.