DIFFICULTY IN THE FORMATION OF COUNTER-ORBITING HOT JUPITERS FROM NEAR-COPLANAR HIERARCHICAL TRIPLE SYSTEMS: A SUB-STELLAR PERTURBER

Since the first discovery of an exoplanet, 51 Peg b, more than 100 hot Jupiters (HJs) with semimajor axis <0.1 au have been detected around main-sequence stars. Nevertheless, their origin remains one of the important unsolved puzzles in this field. It is generally believed that such gas giants are unlikely to be formed in situ and, instead, that they first formed at a large distance from the central star, most likely beyond the ice line, and then migrated significantly inward to their current orbits (but see, e.g., Boley et al. 2015 and Batygin et al. 2015 for different ideas).

The migration mechanisms have not yet been established, but possible scenarios include (1) disk-planet interaction (e.g., Lin et al. 1996; Alibert et al. 2005), (2) planet–planet scattering (e.g., Rasio & Ford 1996; Nagasawa et al. 2008; Nagasawa & Ida 2011; Beauge & Nesvorny 2012), (3) Lidov–Kozai migration (e.g., Kozai 1962; Lidov 1962; Wu & Murray 2003; Fabrycky & Tremaine 2007; Naoz et al. 2011; Petrovich 2015a; Anderson et al. 2016), and (4) secular migration (Wu & Lithwick 2011).

In reality, each of these different migration mechanisms may have contributed to the observed HJ population to some degree. Each mechanism often predicts a different statistical distribution and correlations for the resulting orbital parameters of the planetary systems, and relevant observations may provide a potential clue to distinguish between the different mechanisms. For example, disk–planet interaction would imply that gas giants smoothly migrate inward in a gaseous disk, and thus the angle, ψ, between the stellar spin and planetary orbital axes would not significantly change from its initial value (most likely very close to zero, but it is possible that the spin axis of the central star is moderately misaligned with the normal vector of the primordial disk (e.g., Bate et al. 2010; Foucart & Lai 2011; Batygin 2013; Lai 2014). In contrast, the other migration mechanisms mentioned above rely on a dynamical process after the depletion of the gas disk, which can induce strong spin–orbit misalignment. For this reason, measurement of ψ can be a useful probe in understanding the origins of HJs.

Indeed, the Rossiter–McLaughlin effect has been very successful in measuring the sky-projected spin–orbit angle, λ, for transiting planetary systems (McLaughlin 1924; Rossiter 1924; Queloz et al. 2000; Ohta et al. 2005; Winn et al. 2005). Approximately one-third of the measured systems exhibit significant misalignment of λ > π/4, and a dozen systems even turn out to be in a retrograde orbit (λ > π/2); see Figure 7 of Xue et al. (2014) for example. Such unexpected and counter-intuitive discoveries imply that those HJs should have experienced violent dynamical processes.

While all three of the above dynamical migration mechanisms could produce retrograde HJs, it has proven difficult to produce counter-orbiting HJs (e.g., Fabrycky & Tremaine 2007; Naoz et al. 2011; Liu et al. 2015; Petrovich 2015b). For clarity, we refer to counter-orbiting planets as those with 160° < ψ < 180°, and retrograde planets are simply used to indicate ψ > π/2 throughout the present paper, even if the distinction may not be conventional.

In this context, we should note that the observed λ differs from the true spin–orbit angle ψ; they are related in terms of the orbital inclination i_orb and the obliquity of the stellar spin axis i_⋆ as

$$\begin{eqnarray}&&\mathrm{cos}\psi =\mathrm{cos}{i}_{\star }\mathrm{cos}{i}_{{\rm{orb}}}+\mathrm{sin}{i}_{\star }\mathrm{sin}{i}_{{\rm{orb}}}\mathrm{cos}\lambda \approx \mathrm{sin}{i}_{\star }\mathrm{cos}\lambda .\end{eqnarray} \tag{ 1 }$$

The above approximation holds for transiting systems with ${i}_{{\rm{obs}}}\approx \pi /2$. Since the stellar axis is usually defined so that 0 < i_⋆ < π/2, Equation (1) implies that $\psi \geqslant \lambda$ if 0 < λ < π/2, while ψ ≤ λ if $\pi /2\lt \lambda \lt \pi$.

The true spin–orbit angle ψ is not as easy to obtain, but can be estimated by combining the measurement of i_⋆ via asteroseismology (Unno et al. 1989; Gizon & Solanki 2003; Huber et al. 2013; Campante 2014; Christensen-Dalsgaard 2014). Benomar et al. (2014) performed the first quantitative determination of ψ for transiting planetary systems around main-sequence stars. For HAT-P-7, their asteroseismology analysis yields ${i}_{\star }\approx 30^\circ$, and they obtain ψ ≈ 120° from the joint analysis of the Rossiter–McLaughlin measurement of $\lambda \approx 180^\circ$. For Kepler-25c, they obtain ${i}_{\star }=65\_\_AMP\_\_fdg;{4}_{-6\_\_AMP\_\_fdg;\quad 4}^{+10\_\_AMP\_\_fdg;\quad 6}$ and $\psi =26\_\_AMP\_\_fdg;{9}_{-9\_\_AMP\_\_fdg;\quad 2}^{+7\_\_AMP\_\_fdg;\quad 0}$, which should be compared with λ = 94 ± 71. Indeed, these results demonstrate the importance of the projection effect mentioned above. More importantly, planetary systems with λ ≈ 180° may not necessarily be counter-orbiting, but just retrograde. This may also be the case for HAT-P-6b with λ = 165° ± 6° (Albrecht et al. 2012) and HAT-P-14b with λ = 1891 ± 51 (Winn et al. 2011).

Therefore, the existence of counter-orbiting planets has not yet been observationally established. Nevertheless, it is tempting to consider a dynamical model that can theoretically explain counter-orbiting HJs, if they exist at all. One promising mechanism was recently proposed by Li et al. (2014). They considered a near-coplanar hierarchical triple system and derived a flip condition that the inner planet becomes counter-orbiting under secular perturbation up to the octupole order of the gravitational potential of the outer object in a very eccentric orbit.

To be more specific, their flip condition is written as

$$\begin{eqnarray}&&\epsilon \gt {\epsilon }_{\mathrm{crit},i}=\displaystyle \frac{8}{5}\displaystyle \frac{1-{e}_{1,i}^{2}}{7-{e}_{1,i}(4+3{e}_{1,i}^{2})\mathrm{cos}({\omega }_{1,i}+{{\rm{\Omega }}}_{1,i})},\end{eqnarray} \tag{ 2 }$$

in terms of

$$\begin{eqnarray}&&\epsilon \equiv \displaystyle \frac{{a}_{1}}{{a}_{2}}\displaystyle \frac{{e}_{2}}{1-{e}_{2}^{2}}\end{eqnarray} \tag{ 3 }$$

which characterizes the ratio of the orbit-averaged octupole to quadrupole terms in the massless limit (${m}_{1}\ll {m}_{0},{m}_{1}\ll {m}_{2}$). In the above expressions, e, a, ω, Ω, and m denote the eccentricity, semimajor axis, argument of periastron, and longitude of the ascending node with the subscripts 1 and 2 indicating the inner and outer body, respectively. In the massless limit, a₁, a₂, and e₂ are conserved, and thus as defined by Equation (3) is also a constant of motion. The other orbital elements are time-dependent, and we use the subscript i in Equation (2) to indicate their initial values.

Petrovich (2015b) presented a more general form of the flip condition (2) on the basis of the conservation of energy (the orbit-averaged quadrupole and octupole potential terms) for the coplanar hierarchical triple configuration. His result, Equation (11) of Petrovich (2015b), can be written as

$$\begin{eqnarray}&&{\epsilon }_{{\rm{L}}}\lt {\epsilon }_{{\rm{oct}}}\equiv \displaystyle \frac{{m}_{0}-{m}_{1}}{{m}_{0}+{m}_{1}}\displaystyle \frac{{a}_{1}}{{a}_{2}}\displaystyle \frac{{e}_{2}}{1-{e}_{2}^{2}}\lt {\epsilon }_{{\rm{U}}},\end{eqnarray} \tag{ 4 }$$

which reduces to Equation (2) in the massless (or test-particle) limit of the inner planet. The lower and upper limits, _L and _U, defining the boundary of the flip region are determined by the value of the final angle ${\varpi }_{f}\equiv {\mathrm{cos}}^{-1}{\hat{{\boldsymbol{e}}}}_{1,f}\cdot {\hat{{\boldsymbol{e}}}}_{2,f}$ between the inner and outer orbital unit Lenz vectors. Specifically, _L and _U correspond to ϖ_f = 0 and π, respectively. The upper limit can be practically neglected for sub-stellar perturbers as considered in the present study, but is very important for planetary perturbers (Y. Xue et al. 2016, in preparation).

Li et al. (2014) numerically computed the evolution of such coplanar triple systems in the massless limit, and confirmed that the flip condition is very well described by Equation (2). Also, in the large inclination regime, that analytical flip criterion agrees well with the numerical results, even up to m₂/m₁ > 7 (Teyssandier et al. 2013). They found that e₁ increases monotonically and the mutual orbital inclination between the two bodies, i₁₂, remains low just before the flip, and then the orbital flip of the inner planet proceeds over a very short timescale as e₁ becomes very close to unity, $1-{e}_{1}$ ∼ 10⁻³–10⁻⁴. In that case, the angular momentum of the inner planet is roughly given as ${m}_{1}\sqrt{{{Gm}}_{0}(1-{e}_{1}^{2})}\approx {m}_{1}\sqrt{2{{Gm}}_{0}(1-{e}_{1})}$, and even a small perturbative torque may easily change the angular momentum of the inner planet, and thus flip its orbit if the value of $1-{e}_{1}$ is sufficiently small. Li et al. (2014) proposed a coplanar-flip mechanism for the formation of counter-orbiting HJs in which the inner planet flips by ∼180° before the tidal evolution dominates, and then its extremely eccentric orbit is quickly circularized due to the strong tidal interaction with the central star.

However, Liu et al. (2015) found that short-range forces, general relativity (GR), planetary tide (PT; non-dissipative), and rotational distortion, suppress the extreme value of e₁ that otherwise could be achieved due to the octupole term in hierarchical triple systems with large mutual orbital inclination (i.e., not coplanar configuration) in the Lidov–Kozai oscillation. These additional forces induce precession of the Lenz vector of the inner planet and impose a strict upper limit on the maximum achievable value of e₁; as the short-range forces become stronger, the orbital flips are more confined to the region where the mutual orbital inclination i₁₂ is close to 90°. This result strongly implies that one needs to incorporate those short-range forces in order to describe properly the dynamics of near-coplanar hierarchical triple systems, which is not taken into account in Li et al. (2014).

Petrovich (2015b) performed a series of such simulations for planetary perturbers including short-range force effects that induce pericenter precession of the inner orbit, such as GR, PT, and rotational distortion. All of the resulting HJs in his simulations turn out to be in a prograde and low-obliquity orbit. This is mainly because most of his initial conditions do not satisfy _oct < _U in the flip condition (4), even when they satisfy _oct > _L. Therefore, his set of simulations does not cover the relevant parameter space for the formation of counter-orbiting HJs, even though the simulations are for planetary perturbers, unlike for sub-stellar perturbers as we consider below.

These interesting previous results motivated us to systematically explore the fate of the inner planet in near-coplanar hierarchical triple systems including the quadrupole and octupole terms of the gravitational potential of the outer perturber and short-range forces. Our simulation is based on the orbit-averaged secular dynamics following the formulation of Correia et al. (2011), in which the stellar and planetary spin effects are also incorporated and the octupole-order effect is included following Liu et al. (2015). Because we are primarily interested in the origin of counter-orbiting HJs, we consider only those systems that initially satisfy the analytical flip condition (2). The present paper focuses on the stellar perturber, and the parameter space relevant to the planetary perturber will be discussed in a subsequent paper. In this case, we find that most of the systems are tidally disrupted and a fraction of surviving planets remain mainly as prograde HJs; the formation of counter-orbiting HJs is possible only in a very restricted parameter range.

The rest of the paper is organized as follows. Section 2 describes the basic configuration of the hierarchical coplanar triple systems that we simulate. The simulation results are presented and discussed in Section 3 where we consider the parameter dependence in detail. Section 4 presents the spin–orbit angle distribution in this scenario. Finally, Section 5 is devoted to a summary and the implications of the present paper. The set of equations that we employ is based on Correia et al. (2011) and Liu et al. (2015), but is summarized explicitly in Appendix A for convenience and clarity. The analytical expression for the short-range force effects are summarized in Appendix B.

A schematic configuration of near-coplanar hierarchical triple systems for our numerical simulations is illustrated in Figure 1. A central star of mass m₀ and radius R₀ is located at the origin of the coordinate. The normal vector of the invariable plane is the total orbital angular vector ${{\boldsymbol{G}}}_{{\rm{tot}}}={{\boldsymbol{G}}}_{1}+{{\boldsymbol{G}}}_{2}$ of the inner and outer bodies. Thus, the mutual orbital inclination angle of the two orbits is given by ${i}_{12}={i}_{1}+{i}_{2}$, where i₁ and i₂ are the inclinations of each orbit with respect to the invariable plane. Throughout the present paper, we adopt m₀ = 1 M_⊙, R₀ = 1 R_⊙, m₁ = 1M_J, and R₁ = 1R_J for definiteness.

Zoom In Zoom Out Reset image size

The equations of motions that we adopt are based on the work of Correia et al. (2011), in which short-range forces, GR, spin rotation, and tidal effects for both the star and inner planet are included in addition to the quadrupole term of the orbit-average gravitation potential of the outer body. We modify their equations so as to incorporate the octupole secular terms following Liu et al. (2015). The full equations of motion are described in Appendix A.

2.1. Model Parameters

2.2. Fiducial Case

Figure 2 plots the final states of the inner planet in our fiducial model for coplanar hierarchical triple systems. In this particular example, we first explore the range of 0.005 < _i < 0.15 to ensure the validity of the analytical flip conditions (Equation (2) by Li et al. 2014 and Equation (4) by Petrovich 2015b). Figure 2 clearly shows that the region below those flip conditions agrees with that of the non-migrating planets in our runs. Therefore, their conditions are accurate for distinguishing the significant migration and non-migration boundary, even if they do not necessarily lead to RHJs; see the discussion below. Its most important conclusion is that retrograde HJs are very difficult to form, despite the fact that the plotted region of ${e}_{1,i}\mbox{--}{\epsilon }_{i}$ is chosen so as to satisfy the flip condition (2) in the massless limit; ∼90% of the systems are tidally disrupted and ∼10% survive as prograde HJs. The fraction of retrograde HJs turns out to be less than 1%. Since this may be a rather unexpected result, we plot the dynamical evolution of typical systems for e_1,i = 0.9 in six panels of Figure 3. While we adopt 10 days as the spin rotation period of the inner planet throughout the current analysis, it may be more relevant to use 10 hr as in the case of Jupiter. In reality, however, the result turns out to be fairly insensitive to the value, as shown in Appendix C below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The lower limit of the analytical flip condition (2) by Li et al. (2014), _crit,i, is a very good approximation for the necessary condition, but obviously is not a sufficient condition because it is derived on the basis of orbital dynamics without short-range force effects. Our simulation shows that _crit,i becomes slightly larger, especially for large e_1,i (small a_1,i). The details of the short-range force effects are described in Section 2.3. One example of the behavior in the region between the _crit,i that we adopted and the real flip boundary including the short-range forces effects for ${\epsilon }_{i}=0.025({a}_{1,i}=13.21\;{\rm{au}})$ is illustrated in Figure 3(a). The system exhibits an oscillation both in $1-{e}_{1}$ and i₁₂, but the resulting pericenter distance q₁ is not small enough for the tidal effect to operate. Thus, the semimajor axis a₁ stays constant and no significant migration occurs for 10¹⁰ years. Therefore, all of the other systems with _i < _crit,i not simulated in the present paper show the same behavior.

If _i is slightly larger than _L, the amplitude of oscillation in $1-{e}_{1}$ becomes larger, as plotted in Figure 3(b) for _i = 0.028 (a_1,i = 14.81 au). In this case, the maximum eccentricity reaches 0.998 rather than 0.990 in the non-migrating example. This large-amplitude oscillation allows the inner planet to reach a minimum pericenter distance of ${a}_{1}(1-{e}_{1})\sim 0.03\;{\rm{au}}$ where tidal dissipation efficiently extracts orbital energy, which results in gradual damping of a₁ at each maximum eccentricity (minimum pericenter distance). Thus, PHJ systems form via multiple close approaches within a typical timescale of several 10⁹ years. The example of Figure 3(b) results in a HJ at a ∼ 0.065 au with i₁₂ ∼ 4°. Indeed, this slow coplanar migration has been systematically studied by Petrovich (2015b), who proposes this as a potential path to PHJs, and our results are in agreement with his proposal.

As _i increases further, the octupole potential starts to dominate and drives e₁ very close to unity. At the same time, the orbit flip happens if the dissipative tide is neglected. Along the line of e₁ = 0.9, we observe two continuous regions where PHJs form (9.0%) via the coplanar-flip mechanism. One example in this region is shown in Figure 3(c) for ${\epsilon }_{i}=0.034({a}_{1,i}=18.01\;{\rm{au}})$. It suggests that this path to PHJs happens over a much smaller timescale than that of Figure 3(b); note the different scales of time in each panel. In this case, $1-{e}_{1}$ monotonically decreases and comes close to ∼10⁻³ where the tidal effect becomes important. Therefore, in the middle of increasing i₁₂, the system starts to be circularized and becomes a PHJ with a_1.f ∼ 0.035 au within <10⁷ years. The mutual orbital inclination oscillates with gradually increasingly amplitude in the range i₁₂ ∼ 0°–30°, and then damps from ∼22° to ∼9° during the circularization stage. Since the eccentricity increases until the end of the orbit flip if no short-range force effects are taken into account, such HJs have relatively low i₁₂. In total, the resulting PHJs (PHJ 9.0%) are preferentially located in the low region. Most of them are formed through the coplanar-flip mechanism within a very short timescale (∼10⁷ years), while a few result from secular tidal damping via eccentricity-inclination oscillation.

Beyond that value of _i, the orbit of the inner planet is indeed flipped, but the fate is very sensitive to change due to the subtle competition between the flipped condition and the tidal disruption as illustrated in Figures 3(d)–(f). As a result, the system behavior looks chaotic and there seems to be no systematic parameter region for the formation of RHJ (see Figure 2).

The evolution necessary for the formation of RHJs similar to Figures 3(d) and (e) happens only in a very narrow parameter range; _i = 0.070 (a_1,i = 37.21 au) and ${\epsilon }_{i}=0.103({a}_{1,i}=54.81\;{\rm{au}})$, respectively. The former is circularized at the second closest point of $1-{e}_{1}$. The orbit suffers from tidal circularization during an orbit flip process within a timescale of a few 10⁷ years. Since the eccentricity of the inner orbit increases in the orbit flip stage, the system suffers from tidal circularization in the beginning of the orbit flip stage in order to not be tidally disruped. Thus, this system ends with i_12,f = 162°, which is only slightly smaller than the highest i₁₂ ever reached, 177°, while the latter is circularized at the first closest point due to the stronger perturbation of the outer body. The tidal circularization starts when the orbit flip process is completed. Since the tidal circularization does not modify i₁₂ significantly, i_12,f remains almost unchanged in the counter-orbiting regime, 172° with ±1° oscillation. Such a high value of i₁₂ suggests that the counter-orbiting HJ can be formed via the coplanar-flip mechanism, which supports the conclusion of Li et al. (2014).

Figure 3(f) presents an example of a tidally disrupted inner planet for _i = 0.113 (a_1,i = 60.10 au). Its pericenter falls into the Roche limit at the second extreme eccentricity approach when $1-{e}_{1}$ reaches ∼2 × 10⁻⁴. Such a state is preferentially found in systems in which the inner planet has a relatively large semimajor axis, since the gravitational interaction between two orbits is stronger when the inner orbit reaches the extreme eccentricity. The comparison among the panels (d), (e), and (f), as well as Figure 2, strongly indicates that the fate of the systems is very sensitive to the parameters. Nevertheless, the conclusion that most of the systems satisfying the flip condition (2) are tidally disrupted, instead of forming counter-orbiting HJs, is quite general.

Since most of the systems become disrupted via orbital flip in our simulation, the condition of forming retrograde or counter-orbiting HJs is fairly fine-tuned. Considering the two successful examples of RHJs as shown above, a subtle change in the initial conditions may significantly modify the evolution and tidally disrupt the system, as shown in Figure 2. Therefore, we may need to fine-tune the parameter sets in order to successfully create RHJs, which seems to be unlikely. Based on the low ratio (RHJ 0.4%) and such an uncertainty, it is difficult to form retrograde or counter-orbiting HJs via the coplanar-flip mechanism.

Before moving to the next subsection, we would like to note that there is an interesting pattern in Figure 2; there are a few branching structures in prograde HJs. These are more significant in Figures 10 and 11 below. Although we have not yet successfully explained this behavior, we suspect that they are related to some timescales in the orbital evolution. We hope to come back to this issue in our next paper.

2.3. Effect of Short-range Forces

Liu et al. (2015) showed that the pericenter precessions due to short-range force effects suppress the growth of the eccentricity of the inner planet and reduce the flip region of i₁₂ for systems under the Lidov–Kozai oscillation. In this subsection, we show that similar suppression also works for near-coplanar triple systems.

A small area around the bottom right region of Figure 2 corresponds to non-migrating planets, despite the fact that they initially satisfy the flip criterion Equation (2). Indeed, this comes from short-range force effects. In order to see their effects separately, we consider the NM (non-migrating) planet example of Figure 3(a) (_i = 0.025 (a_1,i = 13.21 au), e_1,i = 0.9).

The left and right panels of Figure 4 plot the evolution of the mutual orbital inclination, i₁₂, and the pericenter distance of the inner planet in units of its initial semimajor axis. Since this example corresponds to the NM case, the latter is almost equivalent to 1 – e₁. We show the results for the secular orbital perturbation effect alone, orbital and GR correction, orbital and planetary rotational distortion (PRD), orbital and PT,³ and orbital and all short-range force effects, from top to bottom.

Zoom In Zoom Out Reset image size

As expected, the case without the short-ranges forces (top panels) flips the orbital inclination each time $1-{e}_{1}$ becomes less than $\sim {10}^{-3}$. The flip repeats periodically since no other dissipational effects are included. If only the PRD is included, the system still shows the orbital flip, but the maximum value of e₁ is slightly suppressed relative to the purely orbital case.

Precession due to PT could effectively limit the orbital flip with a maximum eccentricity of less than 0.999. On the other hand, GR is very effective in suppressing the eccentricity; the maximum value of e₁ under GR correction barely reaches ∼0.99. Thus, the system stays outside the tidal circularization region for 10¹⁰ years, and the PT never becomes important in reality, as shown in the bottom panels of Figure 4.

The above behavior can be understood by comparing the precession timescales of the Lenz vector ${\hat{{\boldsymbol{e}}}}_{1}$ for those short-range forces, which we plot in Figure 5 on the basis of the expressions in Appendix B. Clearly, GR plays a dominant role for e₁ < 0.995, while PT becomes dominant for ${e}_{1}\gt 0.995;$ PRD is sub-dominant in either case. This is in good agreement with our simulation result shown in Figure 4, and therefore the precession induced by the short-range forces, in particular GR, prevents the orbital flip. Since the short-range forces become stronger for the smaller semimajor axis, the NM planets are located around the high-e_1,i and low-_i region.

Zoom In Zoom Out Reset image size

The previous section has presented the result for our fiducial model, and discussed the dynamical behavior for several examples. Next we consider the dependence of parameters employed in the fiducial model separately in each subsection below. The full list of different models is summarized in Table 1, and we plot the two models in each subsection as examples. Since we already confirmed that planets with ${\epsilon }_{i}\lt {\epsilon }_{{\rm{crit}}}$ do not migrate in practice, we run the models for _crit < _i < 0.15 in what follows.

3.1. Mass of the Outer Perturber

We adopt m₂ = 0.03 M_⊙ as our fiducial value, but one might wonder if the larger mass would be more relevant as (sub-)stellar perturbers. While this sounds reasonable, the larger m₂ significantly increases the tidal disruption ratio, and in practice there is no chance to form retrograde planets.

This is clearly shown in the left and right panels of Figure 6 for m₂ = 1 M_⊙ and 0.1 M_⊙, respectively. Since we focus on the parameter space satisfying the analytic flip condition (2), gravitational perturbation due to the outer body is sufficiently strong to produce potentially the orbital flip. Under such circumstances, the larger m₂ results in the larger e₁ (extremely closer to unity), which leads to the stronger tidal effect. Therefore, in order to survive the tidal disruption, the inner planet should have smaller a_1,i for larger m₂. This is why the fractions of both PHJ and RHJ decrease as m₂ increases.

Zoom In Zoom Out Reset image size

Thus, it is very difficult to form RHJs via the near-coplanar-flip mechanism if the outer perturber has a stellar mass m₂ > 0.1 M_⊙. This is why we adopt m₂ = 0.03 M_⊙ as our fiducial value.

3.2. Semimajor Axis of the Outer Perturber

Consider next the dependence on a₂. Again the proper choice of this parameter is not easy. If a₂ is larger, then the near-coplanar configuration is unlikely. On the other hand, the sub-stellar perturber closer to the central star may also be difficult to form. As a compromise, we select a₂ = 500, 200, 100, and 50 au in Table 1 with a₂ = 500 au being the fiducial value. Figure 7 presents the results for a₂ = 200 and 50 au. There are two important messages from Figure 7.

Zoom In Zoom Out Reset image size

First, RHJs do not form for a₂ ≤ 200 au. In order to become a RHJ, the inner planet needs to experience orbital flip before tidal circularization. This shows a preference for larger a_1,i because the inner planet suffers from less tidal dissipation before reaching the extreme eccentricity for the orbital flip. Even larger a_1,i, however, results in stronger gravitational perturbation from the outer body, and thus the inner planet is tidally disrupted. Due to that subtle competition, RHJs in our fiducial model are confined to the narrow region 0.07 < _i < 0.11. As a₂ decreases, the entire system becomes more compact for the same value of _i. Thus, the stronger gravitational perturbation of the outer body brings the inner planet to the orbit within the Roche limit more easily because the pericenter distance at the same maximum eccentricity is smaller. This is why RHJs disappear for the smaller a₂ models. For the same reason, PHJs are limited for the lower _i region.

Second, the NM planet fraction drops as a₂ decreases: 1.8%, 1.8%, 1.1%, and 0.0% for a₂ = 500, 200, 100, and 50 au, respectively. In the fiducial model, short-range forces suppress the maximum eccentricity and the inner planet does not flip or become tidally circularized around the high-e_1,i and very low-_i region. The same value of e_i, however, corresponds to the smaller a₁ for the smaller a₂ models. Thus, the pericenter distance for those systems becomes smaller, which enhances the tidal dissipation and thus circularizes the orbit. As a result, systems gradually migrate and finally become PHJs via secular eccentricity-inclination oscillation, as illustrated in Figure 3(b). In any case, the formation of RHJs is more difficult for the smaller a₂ than the fiducial model.

3.3. Eccentricity of the Outer Perturber

Sub-stellar perturbers may exhibit a broad range of eccentricity, and we run six simulation sets with e_2,i = 0.3, 0.4, 0.5, 0.6 (fiducial), 0.7, and 0.8, and two examples from these models are plotted in Figure 8. We find that the fraction of PHJs monotonically increases for the more eccentric outer perturber. PHJs tend to form preferentially in low a_1,i where tidal dissipation becomes effective. Since we consider the same range of _i for all of the models, the corresponding value of a_1,i for the same _i becomes smaller as e_2,i increases. Thus, the dependence of the fraction on e₂ is mainly due to scaling. While RHJs are very rare, their fraction also increases slightly as e₂ increases, but it would be mainly due to the scaling of a_1,i with respect to _i.

Zoom In Zoom Out Reset image size

3.4. Mutual Orbital Inclination of the Inner and Outer Orbits

The initial orbits of the inner and outer bodies are naturally expected to be inclined to some extent. While our fiducial model adopts i_12,i = 6°, we examine more inclined cases of i_12,i = 15° and 30° as well as an idealized coplanar case (i_12,i = 0°). The Lidov–Kozai mechanism starts to work for more inclined cases, but we do not consider it here because the orbital flip does not happen in those cases, as mentioned in the Introduction.

The left and right panels of Figure 9 present the results for i₁₂ = 0 and i₁₂ = 15°. In the exact coplanar case, the net force normal to the orbital plane always vanishes and the orbits cannot flip. Thus, RHJs cannot form but PHJs can.

Zoom In Zoom Out Reset image size

As i_12,i increases, the fraction of PHJs decreases monotonically and they are confined around the narrow region with high-e_1,i and low _i. In rare cases, RHJs form in a scatter manner over the ${e}_{1,i}\mbox{--}{\epsilon }_{i}$ plane, probably due to the chaotic nature of the system.

3.5. Viscous Timescale of the Inner Planet

Unfortunately, it is well known that the viscous timescale of planets, t_v,p (equivalently, the tidal delay time and tidal quality factor) is the most uncertain parameter in the equilibrium tidal theory. The observational data for the Jupiter–Io system put an empirical lower limit on that of Jupiter as t_v,J > 15 yr. On the other hand, Socrates et al. (2012) stated that t_v,p < 1.5 years for an initially highly eccentric planetary orbit with a semimajor axis of ∼5 au to be circularized into <0.06 au within 10 Gyr. They argued that the discrepancy between their upper limit and the empirical lower limit for Jupiter should not be taken seriously given various theoretical uncertainties concerning the tidal dissipation model and the diversity of physical properties of the exoplanets.

For instance, more recent work by Storch & Lai (2014) examined the possibility of tidal dissipation in the solid cores of giant planets, and claimed that tidal dissipation in the core can reconcile the Jupiter–Io tidal constraint and very efficient high-eccentricity migration simultaneously.

Given the somewhat confusing situation, we decided to adopt t_v,p = 0.03 as our fiducial value, simply following Li et al. (2014). Our purpose in the present paper is not to find a suitable value for t_v,p, but to understand the role of t_v,p in the orbit flip of near-coplanar triple systems. Thus, we examine the other three cases with t_v,p = 0.3, 0.003, and 0.0003 yr as well.

The results are plotted in Figure 10. As expected, the fate of the inner planet is very sensitive to the very uncertain value of t_v,p. When t_v,p is smaller, the tide on the planet becomes stronger and the planet suffers from very efficient circularization, even at a larger pericenter distance. Thus, the majority of the tidally disrupted planets for t_v,p = 0.3 yr survive as PHJs and RHJs for t_v,p = 0.003 yr. The lower right region of Figure 10 corresponds to planets at a relatively larger pericenter distance, which are thus insensitive to the value of t_v,p.

Zoom In Zoom Out Reset image size

Of course, the value of t_v,p = 0.003 yr is very extreme and unrealistic; even the paucity of observed RHJs is inconsistent with the choice. Nevertheless, Figure 10 clearly illustrates that the uncertainty of the tidal dissipation model is the key to understanding the formation and dynamical evolution of HJs in general.

3.6. The Proportional Constant for the Roche Limit

Finally, we consider the criterion of the tidal disruption itself. As discussed in Section 2.1, the proportional factor f of the Roche limit in Equation (5) is not precisely determined. While we adopt f = 2.7 following Guillochon et al. (2011) from hydrodynamical simulations, f = 2.16 is reported by Faber et al. (2005) and f = 1.66 is adopted in simulations by Naoz et al. (2012).

As shown in the previous subsection, the efficiency of the tidal disruption is the most important factor in determining the fate of the inner planet. Thus, we plot the cases of f = 2.16 and f = 1.66 in the left and right panels of Figure 11, respectively.

Zoom In Zoom Out Reset image size

Similar to Figure 10, short-range forces are effective and suppress the growth of the eccentricity of the inner planet in the lower right region of Figure 11. Thus, the pericenter distance of the inner planets around the region is larger than R_roche in any case, and the fudge factor f hardly changes the evolution of those planets.

On the other hand, tidally disrupted planets in our fiducial model are sensitive to the value of f. As is clear from Figure 11, those planets turn out to survive as PHJs and RHJs for the smaller value of f, and there are no tidally disrupted planets for f = 1.66. Indeed, the result with f = 1.66 is already virtually indistinguishable from the case where tidal disruption happens only when the inner planet falls into the central star.

So far, we have classified the surviving HJs into prograde or retrograde according to the mutual orbital inclination angle i₁₂ of the inner and outer orbits, i.e., i₁₂ < 90° or >90°, respectively. In reality, however, i₁₂ cannot be measured directly since the possible outer perturbers of the observed HJs are hardly identified. Thus, observationally, the distinction between prograde and retrograde HJs is made based on the the value of λ, the sky-projected angle of i_s1, obtained from the Rossiter–McLaughlin effect. Since our current simulation runs solve the evolution of the stellar spin axis as well, we can address the validity of a somewhat conventional assumption of i₁₂ = i_s1. The result is plotted in Figure 12, which basically confirms that i₁₂ can be used as a proxy for i_s1 as long as the stellar spin vector is initially completely aligned with the orbital angular momentum vector of the inner planet (i_s1,i = 0), as we adopted in the present runs.

Zoom In Zoom Out Reset image size

Now we show the distribution of i_s1 in Figure 13. These plots indicate that PHJs and RHJs in our simulations correspond almost exclusively to well-aligned (i₁₂ < 20°) and counter-orbiting ($180^\circ -{i}_{12}\lt 20^\circ$) planets. This is not the case, however, for models with very strong tidal interaction (t00030 and f216), which exhibit a very broad distribution of i₁₂, and thus of i_s1.

Zoom In Zoom Out Reset image size

In the coplanar-flip mechanism, the planetary orbit suffers from tidal circularization after the orbit flip. Thus, the system ends up with PHJs if the tidal circularization happens before the orbit flip, and RHJs if the orbit flip occurs before circularization. On the other hand, the tidally disrupted planets have a very broad distribution of i_s1 that we define at the epoch when the pericenter distance of the inner planet reaches the Roche limit. This result implies that those planets fall into the Roche limit over a very short timescale that is less than that of the orbit flip.

Observation of the Rossiter–McLaughlin effect has revealed a dozen possible retrograde planets, which has already challenged the conventional theory of planet formation. Although there exists no reliable candidate (yet), the presence of counter-orbiting planets would have an even stronger impact on the formation theory; somewhat conventional planetary migration scenarios including disk–planet interaction, planet–planet scattering, and the Lidov–Kozai migration are successful in producing retrograde planets, but fail to explain the counter-orbiting planets in general.

An interesting and attractive possibility is based on the extreme eccentricity evolution expected for the near-coplanar hierarchical triple system. Indeed, Li et al. (2014) and Petrovich (2015b) derived an analytical condition for the orbital flip of the inner planet which holds for the massless limit of the inner planet under the quadrupole and octupole gravitational potentials of the outer perturber while neglecting short-range forces (GR, star and inner PT, and rotational distortion).

In the present paper, we have performed a series of systematic simulations for the sub-stellar outer perturber case, including short-range forces, and examined in detail the conditions for orbital flip in a more realistic situation.

Our main findings are summarized as follows.

(1) Most of the near-coplanar hierarchical triple systems that satisfy the analytical flip condition do not produce counter-orbiting planets. Instead, the inner planets in those systems are tidally disrupted. A small fraction of the systems end up with prograde HJs and very few retrograde HJs are produced. Systems that do not satisfy the analytical flip condition do not exhibit any significant migration of the inner planet.

(2) The break-down of the the analytical flip condition is due to the short-range forces, which suppresses the extreme eccentricity evolution of the inner planet that is required for the orbital flip.

(3) The results are almost independent of the model parameters, and thus are fairly generic unless unrealistically strong tidal effects are assumed.

(4) The mutual orbital inclination angle between the inner planet and outer perturber, and the spin–orbit angle between the central star and the inner planet are almost the same. Their distribution of the surviving HJs is bimodal; ∼0°–20° for prograde, and ∼160°–180° for retrograde planets, and virtually nothing in-between.

Our simulation runs span the parameter space that satisfies the analytical flip condition, and more importantly uniformly sample the _1,i–e_1,i plane without assuming any prior distribution for their realistic values. Therefore, the predicted statistics for the fate of the inner planet under such configurations are significantly biased. Having emphasized such warnings, however, it might be instructive to present some statistics simply illustrating the difficulty of forming counter-orbiting planets in a near-coplanar hierarchical triple system.

Figure 14 plots the fraction of four different final outcomes of the inner planet; NM (non-migrating planet) in black, PHJ (prograde HJ) in red, RHJ (retrograde HJ) in blue, and TD (tidally disrupted planet) in green. The left panel corresponds to the number fraction of each fate simply based on the numbers out of 1800 runs for each model summarized in Table 1. The right panel is computed from their sub-sample with $10\;{\rm{au}}\lt {a}_{1,i}\lt 30\;{\rm{au}}$ so as to sample the a_1,i–e_1,i plane assuming eccentric inner gas giant planets orbiting at reasonable distances from the central star for comparison purposes.

Zoom In Zoom Out Reset image size

In any case, our basic conclusion remains the same even if the statistics shown here just for example may be highly biased; it is very difficult to produce the retrograde planet in the present scenario, while some fraction of prograde HJs might have formed through this channel. This implies that the formation of counter-orbiting planets imposes an even more serious challenge for the theory. Instead, it could simply be that counter-orbiting planet candidates with the projected spin–orbit angle λ ≈ 180° are mildly misaligned, with their true spin–orbit angles ψ being much less than 180° as suggested for HAT-P-7b (Benomar et al. 2014). In this respect, future observational searches for counter-orbiting planets combined with the Rossiter–McLaughlin effect and asteroseismology continue to be important, and hopefully will provide an exciting puzzle for planet formation.

Finally, we note that the presence of numerous tidally disrupted planets is not specific to near-coplanar hierarchical triple systems, but is a fairly generic outcome in planetary migration models and in spin–orbit realignment models (Lai 2012; Rogers & Lin 2013; Xue et al. 2014; Li & Winn 2015). Thus, it is of vital importance to look for possible signatures of such tidal disruption events observationally. Indeed, recent studies on the determination of the orbital decay rate (Jiang et al. 2015) and on the unsual photometric signals in KIC 8462852 (Bodman & Quillen 2015; Boyajian et al. 2015), for instance, are closely related to this important direction.

We thank Shoya Kamiaka and Kento Masuda for useful discussions. We are also grateful to an anonymous referee for several important suggestions that improved our earlier manuscript. This research is supported by the Grant-in Aid for Scientific Research by Japan Society of Promotion of Science No. 24340035.

To be self-contained, we write the secular equations of motion used in the present paper for a hierarchical triple system. We consider gravitational interaction up to the octupole expansion of the outer body as described by Liu et al. (2015). In addition, we include the general relativistic correction, the spin effect of the central star and the inner planet, and tidal effect following Correia et al. (2011). In addition, we incorporate the damping of the stellar spin due to magnetic braking following Barker & Ogilvie (2009).

The subscripts 0, 1, and 2 distinguish the quantities for the central star, the inner planet, and the outer perturber, respectively. The mass and radius of those objects are denoted by m and R. The spin rate ω_i and gravity coefficients ${J}_{{2}_{i}}$ for the star (i = 0) and inner planet (i = 1) are written as

$$\begin{eqnarray}&&{J}_{{2}_{i}}={k}_{{2}_{i}}\displaystyle \frac{{\omega }_{i}^{2}{R}_{i}^{3}}{3{{Gm}}_{i}},\end{eqnarray} \tag{ 6 }$$

where ${k}_{{2}_{i}}$ is the second Love number that characterizes the deformation property of each body.

All of the equations are written in Jacobi coordinates with ${{\boldsymbol{r}}}_{1}$ being the relative to the position from m₀ to m₁, and a and e are the semimajor axis and eccentricity, respectively. Then, the evolution of the spin and orbit can be tracked in the octupole approximation using three parameters. Spin angular momentum:

$$\begin{eqnarray}&&{{\boldsymbol{L}}}_{i}={C}_{i}{\omega }_{i}{\hat{{\boldsymbol{s}}}}_{i},\end{eqnarray} \tag{ 7 }$$

where ${\hat{{\boldsymbol{s}}}}_{i}$ is the unit vector of ${\hat{{\boldsymbol{L}}}}_{i}$ and C_i is the principal moment of inertia; orbital angular momentum:

$$\begin{eqnarray}&&{{\boldsymbol{G}}}_{i}={\beta }_{i}\sqrt{{\mu }_{i}{a}_{i}(1-{e}_{i}^{2})}{\hat{{\boldsymbol{k}}}}_{i},\end{eqnarray} \tag{ 8 }$$

where ${\hat{{\boldsymbol{k}}}}_{i}$ is the unit vector of ${\hat{{\boldsymbol{G}}}}_{i}$ with ${\beta }_{1}={m}_{0}{m}_{1}/({m}_{0}+{m}_{1})$, ${\beta }_{2}=({m}_{0}+{m}_{1}){m}_{2}/({m}_{0}+{m}_{1}+{m}_{2})$, ${\mu }_{1}=G({m}_{0}+{m}_{1})$, and ${\mu }_{2}=G({m}_{0}+{m}_{1}+{m}_{2})$; and, finally, the Lenz vector:

$$\begin{eqnarray}&&{{\boldsymbol{e}}}_{1}=\displaystyle \frac{({\dot{{\boldsymbol{r}}}}_{1}\times {{\boldsymbol{G}}}_{1})}{{\beta }_{1}{\mu }_{1}}-\displaystyle \frac{{{\boldsymbol{r}}}_{1}}{{r}_{1}}.\end{eqnarray} \tag{ 9 }$$

We define direction angles as

$$\begin{eqnarray}&&\mathrm{cos}{\theta }_{i}={\hat{{\boldsymbol{s}}}}_{i}\cdot {\hat{{\boldsymbol{k}}}}_{1},\qquad \mathrm{cos}{\epsilon }_{i}={\hat{{\boldsymbol{s}}}}_{i}\cdot {\hat{{\boldsymbol{k}}}}_{2},\qquad \mathrm{cos}{i}_{12}={\hat{{\boldsymbol{k}}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2},\end{eqnarray} \tag{ 10 }$$

where θ_i is the angle between the spin of the ith body (in the main text we use i_s1 to denote θ₀), ${\hat{{\boldsymbol{s}}}}_{i}$, and inner orbit, ${\hat{{\boldsymbol{k}}}}_{1}$, _i is the angle between the spin of the ith body, ${\hat{{\boldsymbol{s}}}}_{i}$, and outer orbit, ${\hat{{\boldsymbol{k}}}}_{2}$, and i₁₂ is the inclination between two orbits.

Averaging the equations of motion over the mean anomalies of the inner and outer bodies, we obtain the following equations for the conservative motion:

$$\begin{eqnarray}{\dot{{\boldsymbol{G}}}}_{1} & = & -\gamma (1-{e}_{1}^{2})\mathrm{cos}{i}_{12}{\hat{{\boldsymbol{k}}}}_{2}\times {\hat{{\boldsymbol{k}}}}_{1}+5\gamma ({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2}){\hat{{\boldsymbol{k}}}}_{2}\times {{\boldsymbol{e}}}_{1}\\ & & -\displaystyle \frac{25}{16}{\epsilon }_{{\rm{oct}}}\gamma \phantom{\rule{0ex}{3.2ex}}\{[2(1-{e}_{1}^{2})[({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{e}}}}_{2})\mathrm{cos}{i}_{12}\\ & & +({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})(\hat{{{\boldsymbol{k}}}_{1}}\cdot {\hat{{\boldsymbol{e}}}}_{2})]{\hat{{\boldsymbol{k}}}}_{1}\\ & & +2[(1-{e}_{1}^{2})(\hat{{{\boldsymbol{k}}}_{1}}\cdot {\hat{{\boldsymbol{e}}}}_{2})\mathrm{cos}{i}_{12}\\ & & -7({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})]{{\boldsymbol{e}}}_{1}]\times \ {\hat{{\boldsymbol{k}}}}_{2}\\ & & +\left[2(1-{e}_{1}^{2})({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})\mathrm{cos}{i}_{12}{\hat{{\boldsymbol{k}}}}_{1}+\left[\displaystyle \frac{8}{5}{{\boldsymbol{e}}}_{1}^{2}-\displaystyle \frac{1}{5}\right.\right.\\ & & -7{({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})}^{2}+(1-{e}_{1}^{2}){\mathrm{cos}}^{2}{i}_{12}\phantom{\rule{0ex}{3ex}}]{{\boldsymbol{e}}}_{1}]\times \hat{{{\boldsymbol{e}}}_{2}}\}\\ & & -\displaystyle \sum _{i}{\alpha }_{1i}\mathrm{cos}{\theta }_{i}{\hat{{\boldsymbol{s}}}}_{i}\times {\hat{{\boldsymbol{k}}}}_{1},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}{\dot{{\boldsymbol{G}}}}_{2} & = & -\gamma (1-{e}_{1}^{2})\mathrm{cos}{i}_{12}{\hat{{\boldsymbol{k}}}}_{1}\times {\hat{{\boldsymbol{k}}}}_{2}+5\gamma ({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2}){{\boldsymbol{e}}}_{1}\times {\hat{{\boldsymbol{k}}}}_{2}\\ & & +\displaystyle \frac{25}{16}{\epsilon }_{{\rm{oct}}}\gamma \phantom{\rule{0ex}{3ex}}\{[2(1-{e}_{1}^{2})[({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{e}}}}_{2})\mathrm{cos}{i}_{12}\\ & & +({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})(\hat{{{\boldsymbol{k}}}_{1}}\cdot {\hat{{\boldsymbol{e}}}}_{2})]{\hat{{\boldsymbol{k}}}}_{1}\\ & & +2[(1-{e}_{1}^{2})(\hat{{{\boldsymbol{k}}}_{1}}\cdot {\hat{{\boldsymbol{e}}}}_{2})\mathrm{cos}{i}_{12}\\ & & -7({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})]{{\boldsymbol{e}}}_{1}]\times {\hat{{\boldsymbol{k}}}}_{2}\\ & & +\left[2(1-{e}_{1}^{2})({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})\mathrm{cos}{i}_{12}{\hat{{\boldsymbol{k}}}}_{1}+\left[\displaystyle \frac{8}{5}{{\boldsymbol{e}}}_{1}^{2}-\displaystyle \frac{1}{5}\right.\right.\\ & & -7{({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})}^{2}+(1-{e}_{1}^{2}){\mathrm{cos}}^{2}{i}_{12}\phantom{\rule{0ex}{3ex}}]{{\boldsymbol{e}}}_{1}]\times \hat{{{\boldsymbol{e}}}_{2}}\}\\ & & -\displaystyle \sum _{i}{\alpha }_{2i}\mathrm{cos}{\epsilon }_{i}{\hat{{\boldsymbol{s}}}}_{i}\times {\hat{{\boldsymbol{k}}}}_{2},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}{\dot{{\boldsymbol{e}}}}_{1} & = & -\displaystyle \frac{\gamma (1-{e}_{1}^{2})}{| | {{\boldsymbol{G}}}_{1}| | }[\mathrm{cos}{i}_{12}{\hat{{\boldsymbol{k}}}}_{2}\times {{\boldsymbol{e}}}_{1}-2{\hat{k}}_{1}\\ & & \times {{\boldsymbol{e}}}_{1}-5({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2}){\hat{{\boldsymbol{k}}}}_{2}\times {\hat{{\boldsymbol{k}}}}_{1}]\\ & & -\displaystyle \frac{25}{16}{\epsilon }_{{\rm{oct}}}\gamma \phantom{\rule{0ex}{3ex}}\{2\sqrt{(1-{e}_{1}^{2})}[({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})\mathrm{cos}{i}_{12}{{\boldsymbol{e}}}_{1}\\ & & +\left.\left[\displaystyle \frac{8}{5}{{\boldsymbol{e}}}_{1}^{2}-\displaystyle \frac{1}{5}-7{({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})}^{2}+(1-{e}_{1}^{2}){\mathrm{cos}}^{2}{i}_{12}\right]{\hat{{\boldsymbol{k}}}}_{1}\right]\times \hat{{{\boldsymbol{e}}}_{2}}\\ & & +2\sqrt{(1-{e}_{1}^{2})}[[({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{e}}}}_{2})\mathrm{cos}{i}_{12}+({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})(\hat{{{\boldsymbol{k}}}_{1}}\cdot {\hat{{\boldsymbol{e}}}}_{2})]{{\boldsymbol{e}}}_{1}\\ & & +[(\hat{{{\boldsymbol{k}}}_{1}}\cdot {\hat{{\boldsymbol{e}}}}_{2})\mathrm{cos}{i}_{12}-7({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{e}}}}_{2})]{\hat{{\boldsymbol{k}}}}_{1}]\\ & & \times \left.{\hat{{\boldsymbol{k}}}}_{2}+\displaystyle \frac{16}{5}({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{e}}}}_{2})\sqrt{(1-{e}_{1}^{2})}{{\boldsymbol{k}}}_{1}\times {{\boldsymbol{e}}}_{1}\right\}\\ & & -\displaystyle \sum _{i}\displaystyle \frac{{\alpha }_{1i}}{| | {{\boldsymbol{G}}}_{1}| | }[\mathrm{cos}{\theta }_{i}{\hat{{\boldsymbol{s}}}}_{i}\times {{\boldsymbol{e}}}_{1}\\ & & +\left.\displaystyle \frac{1}{2}(1-5{\mathrm{cos}}^{2}{\theta }_{i}){\hat{{\boldsymbol{k}}}}_{1}\times {{\boldsymbol{e}}}_{1}\right],\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}{\dot{{\boldsymbol{e}}}}_{2} & = & \displaystyle \frac{\gamma }{\sqrt{1-{e}_{2}^{2}}}\displaystyle \frac{{\beta }_{1}\sqrt{{\mu }_{1}{a}_{1}}}{{\beta }_{2}\sqrt{{\mu }_{2}{a}_{2}}}\phantom{\rule{0ex}{3.3ex}}[(1-{e}_{1}^{2})\mathrm{cos}{i}_{12}{{\boldsymbol{e}}}_{2}\\ & & \times {\hat{{\boldsymbol{k}}}}_{1}-5({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2}){{\boldsymbol{e}}}_{2}\times {{\boldsymbol{e}}}_{1}\\ & & -\left[\displaystyle \frac{1}{2}-3{e}_{1}^{2}+\displaystyle \frac{25}{2}{({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})}^{2}\right.\\ & & -\left.\left.\displaystyle \frac{5}{2}(1-{e}_{1}^{2}){\mathrm{cos}}^{2}{i}_{12}\right]{\hat{{\boldsymbol{k}}}}_{2}\times {{\boldsymbol{e}}}_{2}\right]\\ & & -\displaystyle \frac{25}{16}{\epsilon }_{{\rm{oct}}}\gamma \sqrt{1-{e}_{2}^{2}}\displaystyle \frac{{\beta }_{1}\sqrt{{\mu }_{1}{a}_{1}}}{{\beta }_{2}\sqrt{{\mu }_{2}{a}_{2}}}\\ & & \times \phantom{\rule{0ex}{3ex}}\{2(1-{e}_{1}^{2})\phantom{\rule{0ex}{3ex}}[({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})({\hat{{\boldsymbol{k}}}}_{1}\cdot {{\boldsymbol{e}}}_{2}){\hat{{\boldsymbol{e}}}}_{2}\\ & & +\displaystyle \frac{1-{e}_{2}^{2}}{{e}_{2}}\left[\displaystyle \frac{8}{5}{{\boldsymbol{e}}}_{1}^{2}-\displaystyle \frac{1}{5}-7{({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})}^{2}\right.\\ & & +(1-{e}_{1}^{2}){\mathrm{cos}}^{2}{i}_{12}\phantom{\rule{0ex}{3ex}}]{\hat{{\boldsymbol{k}}}}_{2}]\times {{\boldsymbol{e}}}_{1}\\ & & -\left[2\left(\displaystyle \frac{1}{5}-\displaystyle \frac{8}{5}{{\boldsymbol{e}}}_{1}^{2}\right)({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{e}}}}_{2}){{\boldsymbol{e}}}_{2}\right.\\ & & +14(1-{e}_{1}^{2})({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})({\hat{{\boldsymbol{k}}}}_{1}\cdot {\hat{{\boldsymbol{e}}}}_{2})({\hat{{\boldsymbol{k}}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2}){{\boldsymbol{e}}}_{2}\\ & & +7({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{e}}}}_{2})\left[\displaystyle \frac{8}{5}{{\boldsymbol{e}}}_{1}^{2}-\displaystyle \frac{1}{5}-7{({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{k}}}}_{2})}^{2}\right.\\ & & +(1-{e}_{1}^{2}){\mathrm{cos}}^{2}{i}_{12}\phantom{\rule{0ex}{3ex}}]{{\boldsymbol{e}}}_{2}]\times {\hat{{\boldsymbol{k}}}}_{2})\}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\dot{{\boldsymbol{L}}}}_{i}={\alpha }_{1i}\mathrm{cos}{\theta }_{i}\ {\hat{{\boldsymbol{s}}}}_{i}\times {\hat{{\boldsymbol{k}}}}_{1}+{\alpha }_{2i}\mathrm{cos}{\epsilon }_{i}{\hat{{\boldsymbol{s}}}}_{i}\times {\hat{{\boldsymbol{k}}}}_{2},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{\epsilon }_{{\rm{oct}}}=\displaystyle \frac{{m}_{0}-{m}_{1}}{{m}_{0}+{m}_{1}}\displaystyle \frac{{a}_{1}}{{a}_{2}}\displaystyle \frac{{e}_{2}}{1-{e}_{2}^{2}},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\alpha }_{1i}=\displaystyle \frac{3{{Gm}}_{0}{m}_{1}{J}_{2i}{R}_{i}^{2}}{2{a}_{1}^{3}{(1-{e}_{1}^{2})}^{3/2}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\alpha }_{2i}=\displaystyle \frac{3{{Gm}}_{2}{m}_{i}{J}_{2i}{R}_{i}^{2}}{2{a}_{2}^{3}{(1-{e}_{2}^{2})}^{3/2}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{3{{Gm}}_{2}{\beta }_{1}{a}_{1}^{2}}{4{a}_{2}^{3}{(1-{e}_{2}^{2})}^{3/2}}.\end{eqnarray} \tag{ 19 }$$

In the above expressions, the parameter _oct quantifies the importance of the octupole term relative to the quadrupole term.

The magnetic braking as a spin-down process of the central star is modeled as

$$\begin{eqnarray}&&{\dot{{\boldsymbol{L}}}}_{0,\mathrm{mb}}=-{\alpha }_{{\rm{mb}}}{C}_{0}{\omega }_{0}^{3}{\hat{{\boldsymbol{s}}}}_{0},\end{eqnarray} \tag{ 20 }$$

where the spin-down rate α_mb is set to be $1.66\times {10}^{-13}$ years according to Barker & Ogilvie (2009). Incidentally, the same magnetic braking effect was incorporated in Xue et al. (2014), although it was not noted explicitly.

The correction due to GR induces the precession of the pericenter:

$$\begin{eqnarray}&&{\dot{{\boldsymbol{e}}}}_{1,\mathrm{GR}}=\displaystyle \frac{3{\mu }_{1}{n}_{1}}{{c}^{2}{a}_{1}(1-{e}_{1}^{2})}{\hat{{\boldsymbol{k}}}}_{1}\times {{\boldsymbol{e}}}_{1},\end{eqnarray} \tag{ 21 }$$

where c is the speed of light and n₁ is mean motion of the inner orbit.

For the tidal effect, we adopt the equilibrium tidal model with constant delay time ${\rm{\Delta }}{t}_{i}$ (Mignard 1979). Similarly, the averaged equations are

$$\begin{eqnarray}&&{\dot{{\boldsymbol{G}}}}_{2,\mathrm{tide}}=0,\qquad {\dot{{\boldsymbol{G}}}}_{1,\mathrm{tide}}=-{\dot{{\boldsymbol{L}}}}_{0}-{\dot{{\boldsymbol{L}}}}_{1},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}{\dot{{\boldsymbol{e}}}}_{1,\mathrm{tide}} & = & \displaystyle \sum _{i}\displaystyle \frac{15}{2}{k}_{{2}_{i}}{n}_{1}\left(\displaystyle \frac{{m}_{(1-i)}}{{m}_{i}}\right){\left(\displaystyle \frac{{R}_{i}}{{a}_{1}}\right)}^{5}{f}_{4}({e}_{1}){\hat{{\boldsymbol{k}}}}_{1}\times {{\boldsymbol{e}}}_{1}\\ & & -\displaystyle \sum _{i}\displaystyle \frac{{K}_{i}}{{\beta }_{1}{a}_{1}^{2}}\left[{f}_{4}({e}_{1})\displaystyle \frac{{\omega }_{i}}{2{n}_{1}}({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{s}}}}_{i}){\hat{{\boldsymbol{k}}}}_{1}\right.\\ & & -\left.\left(\displaystyle \frac{11}{2}{f}_{4}({e}_{1})\mathrm{cos}{\theta }_{i}\displaystyle \frac{{\omega }_{i}}{{n}_{1}}-9{f}_{5}({e}_{1})\right){{\boldsymbol{e}}}_{1}\right],\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}{\dot{{\boldsymbol{L}}}}_{i,{tide}} & = & {K}_{i}{n}_{1}\left[{f}_{4}({e}_{1})\sqrt{1-{e}_{1}^{2}}\displaystyle \frac{{\omega }_{i}}{2{n}_{1}}({\hat{{\boldsymbol{s}}}}_{i}-\mathrm{cos}{\theta }_{i}{\hat{{\boldsymbol{k}}}}_{1})\right.\\ & & -\left.{f}_{1}({e}_{1})\displaystyle \frac{{\omega }_{i}}{{n}_{1}}{\hat{{\boldsymbol{s}}}}_{i}+{f}_{2}({e}_{1}){\hat{{\boldsymbol{k}}}}_{1}+\displaystyle \frac{({{\boldsymbol{e}}}_{1}\cdot {\hat{{\boldsymbol{s}}}}_{i})(6+{e}_{1}^{2})}{4{(1-{e}_{1}^{2})}^{9/2}}\displaystyle \frac{{\omega }_{i}}{{n}_{1}}{{\boldsymbol{e}}}_{1}\right],\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&{K}_{i}={\rm{\Delta }}{t}_{i}\displaystyle \frac{3{k}_{{2}_{i}}{{Gm}}_{(1-i)}^{2}{R}_{i}^{5}}{{a}_{1}^{6}},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{f}_{1}(e)=\displaystyle \frac{1+3{e}^{2}+3{e}^{4}/8}{{(1-{e}^{2})}^{9/2}},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{f}_{2}(e)=\displaystyle \frac{1+15{e}^{2}/2+45{e}^{4}/8+5{e}^{6}/16}{{(1-{e}^{2})}^{6}},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{f}_{3}(e)=\displaystyle \frac{1+31{e}^{2}/2+255{e}^{4}/8+185{e}^{6}/16+25{e}^{8}/64}{{(1-{e}^{2})}^{15/2}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{f}_{4}(e)=\displaystyle \frac{1+3{e}^{2}/2+{e}^{4}/8}{{(1-{e}^{2})}^{5}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{f}_{5}(e)=\displaystyle \frac{1+15{e}^{2}/4+15{e}^{4}/8+5{e}^{6}/64}{{(1-{e}^{2})}^{13/2}}.\end{eqnarray} \tag{ 30 }$$

The three main short-range forces (GR, PT, and rotational distortion) modify ${\hat{{\boldsymbol{e}}}}_{1}$ and induce additional precession of ${\hat{{\boldsymbol{e}}}}_{1}$ around ${\hat{{\boldsymbol{k}}}}_{1}$:

$$\begin{eqnarray}&&{\dot{{\boldsymbol{e}}}}_{1}={\omega }_{{\rm{pre}}}{\hat{{\boldsymbol{k}}}}_{1}\times {{\boldsymbol{e}}}_{1}.\end{eqnarray} \tag{ 31 }$$

The precession rate, ω_pre, for the three main short-range forces can be read from the evolution equations in Appendix A in a straightfoward manner as

$$\begin{eqnarray}&&{\omega }_{{\rm{GR}}}=\displaystyle \frac{3{\mu }_{1}{n}_{1}}{{c}^{2}{a}_{1}(1-{e}_{1}^{2})}\propto \displaystyle \frac{1}{{a}_{1}^{5/2}(1-{e}_{1}^{2})},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\omega }_{{\rm{PT}}}=\displaystyle \frac{15}{2}{k}_{{2}_{1}}{n}_{1}\left(\displaystyle \frac{{m}_{0}}{{m}_{1}}\right){\left(\displaystyle \frac{{R}_{1}}{{a}_{1}}\right)}^{5}{f}_{4}({e}_{1})\propto \displaystyle \frac{1}{{a}_{1}^{13/2}{(1-{e}_{1}^{2})}^{5}},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\omega }_{{\rm{PRD}}}=\displaystyle \frac{{\alpha }_{11}}{| | {{\boldsymbol{G}}}_{1}| | }\displaystyle \frac{1}{2}(1-5{\mathrm{cos}}^{2}{\theta }_{1})\propto \displaystyle \frac{1}{{a}_{1}^{7/2}{(1-{e}_{1}^{2})}^{2}},\end{eqnarray} \tag{ 34 }$$

where ω_GR, ω_PT, and ω_PRD are the precession rates induced by GR, PT, and PRD, respectively.

Note that the above expressions are consistent with those of Liu et al. (2015) if the tidal Love number is set to be twice the deformation love number ${k}_{{2}_{i}}$, and if the spin and orbit of the inner planet are aligned (θ₁ = 0°).

Throughout the present analysis, we have adopted 10 days as the spin rotation period of the inner planet. If one considers Jupiter as a typical planet, then 10 hr, instead of 10 days, may be more relevant. Therefore, we repeat our fiducial run using the 10 hr period while keeping all of the other parameters unchanged. Figure 15 shows the result, which is basically identical to Figure 2. For a more quantitative comparison, we show the branching ratios of the final outcomes; PHJ 8.5%, RHJ 0.4%, NM 2.1%, and TD 89.0%. Thus, we conclude that the final result is very insensitive to the choice of the planetary spin period in this range.

Zoom In Zoom Out Reset image size