Introduction
More than 100 of hot Jupiters (HJs) are presently known around main sequence stars. Although there is no precise definition, for our purposes we will include in this group those planets with observed masses m>0.8m Jup and orbital periods P⩽10 days. The lower limit for the mass may appear arbitrary, and is dictated more by dynamical considerations than by physical properties of planetary bodies. The upper limit on orbital periods, however, is more easily justified, but also dynamical in nature. For instance, around solar-type stars, this period corresponds to a limit which giant planets in circular orbits suffer significant tidal effects.
The origin of these planets is still a matter of debate. It appears very unlikely that they formed in situ (e.g., Lin et al. Reference Lin, Bodenheimer and Richardson1996), so their present location must have been achieved after a significant orbital decay from outside the ice line. Although several evolutionary mechanisms were proposed, including disc–planet interactions (e.g., Lin et al. Reference Lin, Bodenheimer and Richardson1996; Benítez-Llambay et al. Reference Benítez-Llambay, Masset and Beaugé2011) and planet–planet scattering (e.g. Rasio & Ford Reference Rasio and Ford1996; Juric & Tremaine Reference Juric and Tremaine2008), a smooth planetary migration due to disc–planet interactions appeared as the best candidate.
Until fairly recently, all detected HJs were consistent with (the assumption of) circular orbits and, more importantly, with values for the misalignment angle consistent with aligned systems (see Winn et al. Reference Winn, Fabrycky, Albrecht and Johnson2010). These orbital characteristics are expected from disc-induced migration, which led further credibility to this scenario. In the past few years, however, the picture changed. A larger population of HJs was analysed for the so-called Rossiter–McLaughlin effect, leading to a large portion of bodies displaying significant values of the misalignment angle (currently ~40%), including about ~15% of planets in retrograde orbits with respect to the stellar spin. So, instead of having a rather simple and ‘cold’ population of HJs, we are now faced with a more complex dynamics.
Since significant misalignment angles are not consistent with smooth planetary migration, their origin must lie elsewhere. At this point we are faced with two questions: (i) what other driving mechanism could explain highly inclined and even retrograde planets and (ii) are all HJs consistent with this new scenario, or must we assume two separate populations of HJs?
Tidal effects tend to align orbits, so any observed misalignment of the orbits must have been caused by the migration mechanism itself. This speaks of a high excitation mechanism which must have affected both the inclination and eccentricity, although the latter may have been later damped by tides. Two scenarios have been proposed for such a mechanism: Lidov–Kozai trapping with a binary companion (e.g., Naoz et al. Reference Naoz, Farr, Lithwick, Rasio and Teyssandier2011, Reference Naoz, Farr, Lithwick and Rasio2012) and planet–planet scattering within an initially cold but dynamically unstable planetary system (Nagasawa et al. Reference Nagasawa, Ida and Bessho2008; Nagasawa & Ida Reference Nagasawa and Ida2011; Beaugé & Nesvorný Reference Beaugé and Nesvorný2012).
Although both scenarios are completely different, the end result is the same. An initially circular orbit of a giant planet beyond the HJ region is excited to high eccentricities (usually close to parabolic orbits) in such a way that the pericentric distance is so close to the star that tidal effects are not only significant but dominant over the gravitational perturbation that generated the excitation. If this process also affected the inclination, then the subsequent orbital evolution of the planet would damp the eccentricity and semimajor axis, leaving as a final product a hot planet with highly inclined but near-circular orbit.
In Beaugé & Nesvorný (Reference Beaugé and Nesvorný2012) we showed that planet–planet scattering may explain many of the observed orbital characteristics of HJs, including the eccentricity-semimajor axis distribution, the so-called 3-day pile-up, and the distribution of misalignment angles. Obviously the result depends on the tidal model and the adopted values for the tidal parameters, but this may actually serve as observational constraints on these little-known parameters.
Although Lidov–Kozai trapping also proved an efficient mechanism, and also explains the existence of highly misalignment HJs, there has been a significantly less comparison with the observed HJ population. For example, it is not clear that this scenario explains the 3-day pile-up or the eccentricity distribution. In short, which of the observed orbital characteristics are due to the excitation mechanism and which to the subsequent tidal evolution?
In this work, we wish to address precisely these issues. The main idea is to explore the Lidov–Kozai model using similar tools and dynamical models as developed in Beaugé & Nesvorný (Reference Beaugé and Nesvorný2012), changing planet–planet scattering by Kozai resonance. Our aim is then to present a consistent comparison between the predictions of both scenarios, and try to deduce which observed characteristics of the planets are robust and which are model-dependent.
Numerical simulations
Let m p be a planet orbiting a star m A and perturbed by its secondary companion m B. We will denote by a the semimajor axis, e the eccentricity, I the inclination, M the mean anomaly, ω the argument of pericenter and Ω the longitude of the ascending node. As with the masses, the orbital elements of the planet will be identified by a subscript p, while those of the binary companion by B. We will assume Jacobi coordinates with m A as the central mass, and the reference plane perpendicular to the total orbital angular momentum.
The equator of m A will coincide with the orbital reference plane, so the inclination I p will be equal to the misalignment angle. Finally, we will denote by I mut the mutual inclination between the orbits of m p and m B, and fundamental to the orbital evolution of the planet due to the gravitational perturbations of the secondary star.
Initial conditions
The simulations presented by Naoz et al. (Reference Naoz, Farr, Lithwick, Rasio and Teyssandier2011, Reference Naoz, Farr, Lithwick and Rasio2012) started with a planet in circular orbit and aligned with the central star, while the mass and orbit of the secondary was chosen randomly from a predefined distribution. Here we adopted a different route.
In Martí & Beaugé (Reference Martí and Beaugé2012) we presented a formation scenario for planets in tight binary systems based on the concept of stellar scattering. Initially we assumed a planet with m p=0.002 M⊙ in a circular orbit with a p=2 AU around a central star of mass m A=1.59 M⊙, plus a secondary star of mass m B=0.4 M⊙ with a B=18.5 AU and e B=0.36. The system was then perturbed by a hyperbolic fly-by with a third star m C with different initial conditions. This generated a sequence of new binary systems in which both the planet and the secondary star acquired new orbital elements.
Although in that work we focused our interest in reproducing the observed γ-Cephei system, the result of our simulations generated a series of synthetic binary systems with different orbital characteristics. We will use those end-products of our stellar scattering experiments as initial conditions for our tidal trapping experiments.
We chose a total of 21 226 scattering experiments using an impactor mass of m C=0.4 m ⊙, a value similar to the mean of the stellar initial mass function for stars (see Kroupa 2001). However, we also checked other values and found no appreciable change in the results. The minimum distance between the impactor m Cand the binary was taken randomly between 0.1 and 4.0 times the apocentric distance between both primary stars. Figure 1 shows the post-scattering distribution of both the planet and the binary companion which will be taken as initial conditions for our tidal simulations. While most of the planets and binaries retained semimajor axes near their original values, there is a considerable spread in both orbits. The same is also noted in the eccentricities, where e B has a very narrow distribution around ~0.36 and practically all values of e p are restricted to <0.2.
Post-scattering distribution of the orbits of both the planet m p and binary companion m B. The vertical dashed lines in the bottom frame correspond roughly to the limits of the Lidov–Kozai resonance domain.
Finally, while most values of I p remained close to zero, the mutual inclinations I mut show a wider range. Approximately 53% have values I mut<40° and thus below the approximate limit necessary for the appearance of the Lidov–Kozai resonance (e.g., Lidov & Kozai Reference Lidov1961; Kozai Reference Kozai1962; Harrington Reference Harrington1968; Innanen et al. Reference Innanen, Zheng, Mikkola and Valtonen1997; Libert & Henrard Reference Libert and Henrard2007; Haghighipour Reference Haghighipour2010). However, the exact limiting inclination for the resonance is a function of the semimajor axis ratio (Libert & Henrard Reference Libert and Henrard2007), so the value adopted here is illustrative but not rigorous. Close to 40% of the initial conditions have 40°<I mut<140° and may be affected by the secular resonance. Lastly, near 7% have highly retrograde orbits, with I mut>140°. Notice that there seems to be two peaks in the distribution of the mutual inclinations (identified by vertical dashed lines), each associated with the limiting values for the appearance of the Lidov–Kozai libration domain.
Equations of motion
This set of 21 226 initial conditions were numerically integrated for the total time span of 109 years. The gravitational effects were modelled with an octupole secular Hamiltonian following the description by Laskar & Boué (Reference Laskar and Boué2010). We also included tidal effects according to the modified model presented in Beaugé & Nesvorný (Reference Beaugé and Nesvorný2012), and adopted the same values for the tidal parameters: Q′p=5×106 and Q′*=107. The tidal equations contained both the dissipative terms and the precessional effects due to the deformation of the extended body. Finally, for completeness, we also included general relativity effects and stellar oblateness. The expression for both may be found in Beaugé & Nesvorný (Reference Beaugé and Nesvorný2012).
Together with the orbital elements, we also integrated the spin equations, using the same descriptions as in Correia et al. (Reference Correia, Laskar, Farago and Boué2011). The initial rotational periods where chosen equal to 28 days (for the star) and 10 days (for the planet). This last value is not really sensitive, since a pseudo-synchronization of the planet was achieved early in the simulation. The initial obliquity of the star was taken equal to zero, and thus the equator coincided with the pre-scattering orbits of both the planet and binary companion. The obliquity of the planet, however, was chosen equal to its inclination with respect to the central star.
All simulations were stopped if the planet reached the Roche radius of the central star, which corresponds approximately to r Roche≃0.009 AU for the cases studied here.
Results
Figure 2 shows the final distribution of planetary orbits at the end of the simulations. The top frame shows the distribution of eccentricities as a function of the semimajor axis, where the colour code indicates the initial value of I mut. Most of the runs lead to very close encounters with the central star and to pericentric distances closer than the Roche radius. Only a few of these still survived in the system after 1 Gyr, and are shown in the plot. The bottom frame shows the final inclination I with respect to the stellar equator, where we only show those runs with final a p>r Roche. Additional information is shown in Table 1, which indicates the total number of initial conditions in each interval of I mut, as well as their outcomes.
Final distribution of 21 226 initial conditions evolved for 1 Gyr including gravitational perturbations, tidal effects, GR and stellar oblateness. Colour code indicates the original mutual inclination between the planet and the secondary star. The Roche radius for the central star is located at ~0.009 AU. Most of the tidally trapped bodies do not survive and are engulfed by m A.
Outcomes of the simulations, as function of the initial mutual inclination
Practically all runs with initial inclinations I mut<40° (black circles) were not tidally trapped and remained with semimajor axes close to their original values. About ~3%, however, were lost due to engulfment by the central star and only ~0.05% gave origin to HJs. For 40°<I mut<90°, most of the initial conditions were susceptible to excitation by the Lidov–Kozai resonance. The end results show precisely this effect, where only ≃14% of the bodies remained close to their original semimajor axes, while almost 85% were tidally trapped. However, practically all had pericentric distances smaller than the Roche radius and were disrupted, so the overall formation efficiency of HJs was actually very low, of the order of ≃0.7%.
A similar result was found for orbits initially retrograde with respect to the perturber. In this case, the tidal trapping appears to have almost the same efficiency; only ≃7% remained near their original semimajor axis) while almost the same percentage as in the prograde case of the initial conditions lead to the formation of stable HJs. Finally, initial conditions retrograde with respect to the perturber but outside the Lidov–Kozai libration domain suffered a similar fate as those in direct orbits, with most of them avoiding tidal trapping. However, the effectivity of the formation of stable HJs is damped to 0.2% as expected, because of the planet being outside the Lidov–Kozai resonance. So, it appears that initial conditions with I mut>90° seemed to be slightly more effective in forming stable HJs than those in direct orbits.
Although the Lidov–Kozai mechanism proved to be inefficient in generating stable HJs, the final distributions of these bodies in the (e p, a p) plane are very similar to that presented in Beaugé & Nesvorný (Reference Beaugé and Nesvorný2012). Practically all form a kind of ‘cascade’ or correlated distribution with e p ~0.8 at a p=0.2 AU and ending with circular orbits at a p ≃0.03 AU. The location of this cascade depends on the tidal parameter Q′p as well as the integration time, shifting towards larger semimajor axis for longer timescales or smaller values of Q′p.
The final distribution in the (I p, a p) plane, however, presents significant differences with that obtained from planet–planet scattering. There is a noticeable lack of planets in low inclination orbits (both direct and retrograde), and the distribution appears almost symmetric with respect to polar orbits. Although Naoz et al. (Reference Naoz, Farr, Lithwick and Rasio2012) found a more uniform distribution, both distributions are in contrast with the observed distribution obtained from data of real HJs. Although there are a few cases in which the initial conditions flip from direct to retrograde orbits (or vice versa), most of the final distribution retains the original direction of the orbital angular momentum.
As shown in Fig. 1, most of our stellar companions have semimajor axis a B ~10−30 AU, generating a very compact binary system (although observations of planets in tight binary systems seem to show a minimum separation of 20 AU for planetary formation, we carried out our simulations for separations smaller than 20 AU as a matter of completeness). To analyse how the final distributions change with the stellar separation, we repeated the simulations applying a multiplicative factor N to the initial value of a B in our initial conditions. Results are shown in Fig. 3 for N=2, 5 and 10, while the statistics of the outcomes are summarized in Table 2. In these runs we have only considered initial conditions within the Lidov–Kozai resonance domain (I mut∈[40°, 140°]).
Final distribution of initial conditions with I mut∈[40°, 140°] where the semimajor axis of the binary was increased by a factor N (i.e., a B→N aB). Top frame shows results for N=2, middle for N=5 and bottom plot for N=10. Total integration time was 1 Gyr. Colour code is the same as previous figure.
Outcomes of the simulations, as function of the separation factor N of the binary companion
On one hand, we note a significant growth in the proportion of stable HJs as function of N, which increase from 1.5% for N=2 to almost 30% for wide binaries (N=10). Much of this new population is fed from planets that were tidally disrupted for lower values of N. Thus, wider binaries lower the number of planets engulfed by the star, many of whom are able to become tidally trapped at larger orbital distances and thus remain dynamically stable. This result is in agreement with those of Naoz et al. (Reference Naoz, Farr, Lithwick and Rasio2012), although our simulations cover a much greater range of semimajor axes.
The distribution in the (e p, a p) and (I p, a p) planes, however, shows little change as function of N. HJs still show a cascade-type distribution in eccentricity as function of the semimajor axis, although there appears to be a stronger correlation between both elements for wider binaries. The distribution of final inclinations still seems to favour values near polar orbits, and most of the initial orbital orientations are preserved after the planets become HJs.
In conclusion, although larger values of a B increase the efficiency of the formation process of hot planets, the final distribution of their orbits seems fairly invariant. A more detailed comparison is shown in the left-hand frame of Fig. 4, where we plot the distribution of HJs in semimajor axis for different values of N. All show a peak near a p≃0.03 AU, which corresponds roughly to the 3-day pile-up in the case of solar-type central stars. Thus, the accumulation of observed giant planets around this value appears to be caused by the tidal trapping itself, and not by the driving mechanism (planet–planet scattering or Lidov–Kozai excitation). This fact can be observed in the right-hand frame of the same figure. Here we show the value of −da p/dt, as a function of a p, as calculated with our tidal model, initial values for a p and e p were chosen for illustration purposes only. The continuous line corresponds to initial conditions equal to a p(t=0)=1 AU and e p(t=0)=0.99, while the dashed line corresponds to e p(t=0)=0.1. As the semimajor axis decays with time, the evolution of the orbit occurs from right to left. The location of the cascade in the (a p, e p) plane is actually determined by the tidal parameter. In this example, we choose its value such that the limit for e p→0 occurs close to the 3-day pile-up.
Left: Final distribution of simulated HJs as function of semimajor axis, for different values of N (black: N=1, blue: N=2, red: N=5, green: N=10). The peak near 0.03 AU corresponds to the 3-day pile-up in solar-type stars. Right: Value of the time derivative of the semimajor axis (sign changed), as function of a p, for initial conditions at 1 AU and e p=0.99 (continuous line) and e p=0.1 (dashed line).
For initial quasi-parabolic orbits of the planet (corresponding to e p close to unity), during the first stage the planet's eccentricity remains finite and the planetary tide are much larger than the stellar counterpart. The magnitude of |da p/dt| decreases for smaller values of a p, until it reaches a minimum for e p~0 (see Rodríguez & Ferraz-Mello Reference Rodríguez and Ferraz-Mello2010). From this point onwards, orbital decay continues, but now fuelled by stellar tides, whose effect increases for the smaller semimajor axis. The value of the semimajor axis for which e p≃0 is a (very pronounced) minimum of |da p/dt|, and implies that bodies will remain near this spot for very long timescales.
However, this minimum is only noted when tidal trapping occurs at high initial eccentricities. As noted in the dashed line, for low initial eccentricities (e p(t=0)=0.1), e p is very quickly damped and the planet reaches the hot-planet region already in almost circular orbit and suffers a monotonic increase in the magnitude of orbital decay for all values of a p.
Thus, the so-called 3-day pile-up appears to be both a consequence of tidal evolution and initial conditions with very high eccentricities, and constitutes observable evidence of tidal trapping. Moreover, the value of a p associated to this accumulation may be used to set the tidal parameter Q′p, as shown by Beaugé & Nesvorný (Reference Beaugé and Nesvorný2012).
It should be pointed out that through-out the whole set of simulations the values of the masses of the bodies were held fixed at the masses of the bodies which form the actual γ-Cephei system. It is important to state that our results are valid only for the mass ratios used in this work. Although, there is no reason to believe that this particular selection of mass ratios are preferred in order to reproduce the same orbital characteristics that were obtained from planet–planet scattering and also some of the properties of the observed population of HJs. We suppose that these values of the masses of the bodies are sufficiently representative for the study of the Lidov–Kozai trapping as a mechanism for the HJ formation.
The cascade
As with the case of planet–planet scattering, we found that the final distribution of HJs in the (e p, a p) plane shows a strong correlation between those parameters and roughly resembles a cascade. Although its shape is similar to a curve of constant orbital angular momentum, both are not equivalent. This raises the questions as to what mechanism causes the appearance of the cascade and on which parameters it depends.
We will concentrate on the dynamical evolution of the planet after it has been tidally trapped, and can therefore disregard the gravitational perturbations from the binary. This simplifies the discussion, since we can study only the tidal effects in the two-body problem. However, it is important to stress that this approximation is only employed here; all our previous simulations included the gravitational effects of the binary. Let us define a set of initial conditions in the (e p, a p) plane defined by the intervals a p∈[a 0−Δa p, a 0+Δa p] and e p∈[e 0−Δe p, e 0+Δe p].
If we choose Δa p=Δe p=0, all initial conditions will define a single point and evolve simultaneously along a curve of constant angular momentum 
$G \propto \sqrt {a_{\rm p} (1 - e_{\rm p}^2 )} $. If we assume that the initial time is different for each initial condition, then the initial configuration will disperse, but still along a curve corresponding to a single value of G. This, however, does not correspond to the ‘cascade’ type shape discussed above, which encompasses a finite interval of values of G.
Let us now consider Δe p=0 but Δa p=0.2 AU, generating an ensemble of initial conditions with a finite spread in semimajor axis but with the same eccentricity. This will show us how planets with different semimajor axis will cover the cascade. The evolution of such a set is depicted in the left-hand frame of Fig. 5, where each colour corresponds to a different time (its value given in years). For comparison, curves of constant G are shown in grey. Even though we chose a fairly significant spread Δa p, the ensemble still evolves as a close pack with little relative dispersion. Although its shape resembles a segment of the cascade, it never completes the figure, and again does not seem to explain the observed distribution.
Simulations of tidal evolution of ensembles of initial conditions in the two-body problem. Colours show these sets for different times (values in years). The grey curves correspond to constant orbital angular momentum. Left: Initial semimajor axes a p∈[1.8, 2.4] AU, and all eccentricities e p=0.99. Right: All initial semimajor axis chosen equal to a p=2 AU and eccentricities in the interval e p∈[0.952, 0.992].
The right-hand frame shows the results of a new set of initial conditions, this time with Δa p=0 and e p∈[0.952, 0.992], which will now illustrate how planets with different eccentricities cover the cascade. Again the colour code shows the shape of the ensemble for different times, and the individual values are the same as those in the other graph. Even though we adopted a small value for Δe, since we are in the realm of almost parabolic orbits it generates a large dispersion in pericentric distance. This has two consequences: first, a large spread in G and second, very different tidal decay timescales. So, while those initial conditions with the lower eccentricities suffer almost no change in the semimajor axis throughout the simulation (even for times ~1010 years), those with higher values of e p suffer a much quicker orbital decay. The result is a cascade covering a significant interval in values of G and e p, and very similar to that obtained from both our planet–planet scattering and Lidov–Kozai trapping experiments.
Summarizing, the cascade-like structure in the (e p, a p) plane obtained for HJs in the tidal trapping simulations seems primarily caused by small differences in the planetary eccentricity at the moment of tidal capture. Differences in semimajor axis cause little effect, while differences in time of capture will also be slight. For example, there is not much change in the location of the cascade from t=109 to 1010 years, and any difference in the initial time will be constrained by these two curves.
Conclusions
In this paper, we have analysed the formation scenario of HJs based on tidal trapping from quasi-parabolic orbits excited by Lidov–Kozai resonances with a binary stellar component. We have used the same tidal and gravitational model as in Beaugé & Nesvorný (Reference Beaugé and Nesvorný2012), where we discussed the same problem but assuming planet–planet scattering as the catalyst. Our aim has been twofold. First, compare the final distribution of HJs predicted by both mechanisms using similar dynamical tools, models and parameters. Second, try to understand the origin of some of the observed characteristics of the simulations and the observed planets.
We have found that several of the final orbital properties of the simulated HJs (as well as the real bodies) are mainly caused by the process of tidal trapping, and independent of the excitation mechanism. This includes both the distribution in the (e p, a p) and the 3-day pile-up. Both may be considered as observational evidence that a significant portion (or even most) of the observed HJs originated from tidal trapping and not from smooth disc-induced migration.
However, we have also found that the final distribution of the inclinations is different with respect to those obtained from planet–planet scattering experiments. Figure 6 shows a comparison of both, together with the observed distribution of the sky-projected misalignment angles currently available for HJs. While the Lidov–Kozai mechanism shows a notorious absence of low-inclination orbits, the distribution obtained from scattering and the observed values show a greater proportion of orbits nearly coplanar with the stellar equator.
Distribution of the final inclination between the planet and the stellar equator, obtained from our Lidov–Kozai simulations (left plot), as compared with the results from planet–planet scattering experiments (Beaugé & Nesvorný (Reference Beaugé and Nesvorný2012)) (centre plot) and with the distribution of the sky-projected misalignment angle λ for observed HJs (right plot).
The reason behind this difference may lie in the excitation mechanism itself. In planet–planet scattering there is no direct correlation between the excitation in eccentricity and that in inclination, and it is possible (albeit likely) that high-eccentricity orbits remain with low inclinations (e.g., below ~45°). Thus, in this scenario, most of the HJs will preferably remain aligned with the equator of the central star, as observed in real planets. In Lidov–Kozai trapping, formation of HJs requires a high inclination with respect with the perturbing mass. Unless the binary is conveniently located in a polar orbit with respect to the equator of m A, the final inclinations of the hot planets will show a more random distribution, with no preference for almost aligned orbits.
Naoz et al. (Reference Naoz, Farr, Lithwick and Rasio2012) proposed to solve this issue arguing in favour of the existence of a second population of HJs generated by smooth disc-induced planetary migration and, consequently, containing low inclinations. However, as noted from Fig. 4, this second population should not show a 3-day pile-up in orbital periods, and thus cannot explain the fact that most of the real HJs in this pile-up show small misalignment angles. Planet–planet scattering suffers from none of these limitations, and thus appears to be more consistent with the distribution of real planets. Even so, it is indeed possible that all three proposed mechanisms (Lidov–Kozai, scattering and smooth planetary migration) could have contributed to the complete sample of known HJs. Only future work, including combined scenarios and additional data, will allow us to speculate as to the effective role of each of them.
Acknowledgements
This work has been partially funded by CONICET and SECYT/UNC. The authors would like to express their gratitude to N. Haghighipour for the invitation us to participate in this special issue, and to both referees for valuable suggestions.
