On the formation of hot and warm Jupiters via secular higheccentricity migration in stellar triples
Abstract
Hot Jupiters (HJs) are Jupiterlike planets orbiting their host star in tight orbits of a few days. They are commonly believed not to have formed in situ, requiring inwards migration towards the host star. One of the proposed migration scenarios is secular higheccentricity or high migration, in which the orbit of the planet is perturbed to high eccentricity by secular processes, triggering strong tidal evolution and orbital migration. Previous theoretical studies have considered secular excitation in stellar binaries. Recently, a number of HJs have been observed in stellar triple systems. In the latter, the secular dynamics are much richer compared to stellar binaries, and HJs could potentially be formed more efficiently. Here, we investigate this possibility by modeling the secular dynamical and tidal evolution of planets in two hierarchical configurations in stellar triple systems. We find that the HJ formation efficiency is higher compared to stellar binaries, but only by at most a few tens of per cent. The orbital properties of the HJs formed in the simulations are very similar to HJs formed in stellar binaries, and similarly to studies of the latter we find no significant number of warm Jupiters. HJs are only formed in our simulations for triples with specific orbital configurations, and our constraints are approximately consistent with current observations. In future, this allows to rule out high migration in stellar triples if a HJ is detected in a triple grossly violating these constraints.
keywords:
planets and satellites: dynamical evolution and stability – planetstar interactions – gravitation1 Introduction
Hot Jupiters (HJs) are Jupiterlike planets orbiting their host star on tight orbits, downward of 10 d. Despite the first detection of a HJ over two decades ago (Mayor & Queloz, 1995), it is still unclear how these planets could have formed. It is commonly believed that planet formation in the protoplanetary disk phase is not efficient so close to the star (e.g. Lin et al. 1996), although some recent studies have considered the possibility of in situ formation through core accretion (Lee et al., 2014; Boley et al., 2016; Batygin et al., 2016). Discounting the latter possibility, then the protoHJ must have formed in regions further away from the star, i.e. beyond the ice line of one to several AU, and subsequently migrated inwards.
Two main migration scenarios have been considered: (1) disk migration e.g. Goldreich & Tremaine 1980; Lin & Papaloizou 1986; Bodenheimer et al. 2000; Tanaka et al. 2002), and (2) migration induced by tidal dissipation in the HJ, requiring high orbital eccentricity, also known as ‘high’ migration. Several mechanisms have been proposed to drive the high eccentricity required for high migration, including
The orbital periods of HJs peak around . In the case of disk migration, migration in principle proceeds until the planet is engulfed by the star, and truncation of the disk needs to be invoked in order to explain the observed orbital period distribution. In contrast, high migration induced by secular processes naturally predicts a ‘stalling’ orbital period of . Also, for disk migration the obliquity (the angle between the stellar spin and orbit of the planet) is expected to be (close to) zero, unless the stellar spin axis was initially misaligned with respect to the protoplanetary disk (e.g. Bate et al. 2010; Foucart & Lai 2011; Batygin 2012; Lai 2014). In contrast, large obliquities, even retrograde, have been observed (e.g. Winn et al. 2011; Mazeh et al. 2015), and these are more naturally explained in the case of secular high migration, in which the orbit of the planet is continuously changing its orientation with respect to the star before the evolution is dominated by tidal dissipation.
However, secular high migration generally faces the problem that the theoretical occurrence rates are about an order of magnitude too low compared to observations. Also, a high efficiency of tidal dissipation in the HJ is required, and the predicted orbital periods are on the short end of the observations unless tidal dissipation is very efficient (e.g. Petrovich 2015a; Anderson et al. 2016).
Recently, HJs have been found orbiting stars in stellar triple systems. The three HJs are WASP12b (Hebb et al., 2009; Bergfors et al., 2013; Bechter et al., 2014), HATP8b (Latham et al., 2009; Bergfors et al., 2013; Bechter et al., 2014), and KELT4Ab (Eastman et al., 2016). In these systems, the HJs orbit the tertiary star in a stellar triple system, i.e. the HJ host star is orbited by a more distant pair of stars. From a dynamical point of view, such a system is a hierarchical fourbody system in the ‘2+2’ or ‘binarybinary’ configuration, where one of the binaries is composed of a single star + planet, and the other binary is composed of the two other stars. The binarity of the distant pair can affect the orbital evolution of the planet around its host star. In particular, Pejcha et al. (2013) showed that in these systems the parameter space for exciting high eccentricities triggered by flips of the orbital planes is larger compared to the equivalent triple systems, potentially enhancing the formation rate of HJs compared to the case of a single stellar companion.
The observations of HJs in stellar triples mentioned above constitute only one example of the configurations in which planets can orbit stars in stellar triple systems. For hierarchical fourbody systems, there exist two classes of longterm stable configurations: the ‘2+2’ and the ‘3+1’ configurations^{1}^{1}1Here, we do not consider the possibility of nonhierarchical orbits such as trojan orbits.. In the ‘2+2’ or ‘binarybinary’ configuration (cf. panel 1 of Fig. 1), two binary pairs are bound in a wider orbit. In the ‘3+1’ or ‘triplesingle’ configuration (cf. panels 24 of Fig. 1), a hierarchical triple, composed of an inner and outer orbit, is orbited by a more distant fourth body.
In the case of ‘2+2’ systems, the planet orbits one of the stars, and only one physically unique configuration is possible (cf. panel 1 of Fig. 1). For ‘3+1’ systems, three physically distinct configurations exist. The planet can orbit

one of the stars in the innermost binary (i.e. an ‘Stype’ or ‘circumprimary’ orbit; cf. panel 2 of Fig. 1);

the (center of mass of the) innermost binary (i.e. a ‘Ptype’ or ‘circumbinary’ orbit; cf. panel 3 of Fig. 1);

the (center of mass of the) triple system (i.e. a ‘circumtriple’ orbit; cf. panel 4 of Fig. 1).
Note that the nomenclature ‘Stype’ and ‘Ptype’, introduced by Dvorak (1982), strictly only applies to binary star systems.
In case (i), the planet is orbiting a single star; it can approach this star at a short distance (i.e. a fraction of an AU) while maintaining dynamical stability, i.e. not destroying the hierarchy of the system. However, in cases (ii) and (iii), the orbit of the planet is very likely unstable if the planet were to approach one of the stars at such a short distance because of the perturbations by the other stars. Therefore, for ‘3+1’ systems, likely only case (i) is relevant for HJ formation through secular high migration because the latter requires repeated close passages to a star. Nevertheless, it might be possible to form HJs through tidal capture of dynamically unstable planets by stars, triggered e.g. by secular evolution. We do not consider the latter possibility here. In case (i), similarly to ‘2+2’ fourbody systems, coupled LK oscillations between the orbits can result in extremely high eccentricities and secular chaotic motion (Hamers et al., 2015).
The above considerations suggest a higher formation efficiency of HJs in stellar triples compared to stellar binaries^{2}^{2}2Not taking into account the relative frequencies of binaries and triples., and this motivates study of the efficiency of HJ formation through secular high migration in stellar triples. In this paper, we investigate this through numerical simulations, focussing on the ‘2+2’ and ‘3+1’ (i) configurations (i.e. panels 1 and 2 of Fig. 1). We pay particular emphasis on comparing the efficiency of HJ production to the case of a stellar binary (i.e. panels 1 or 2 of Fig. 1 with removed).
Furthermore, we consider the conditions necessary for HJ formation through high migration in triples, and show that these conditions exclude the formation of HJs through high migration in triples with certain configurations. Although not presently the case, a future detection of a HJ violating these conditions would be strong indication for another formation mechanism at work such as disk migration.
The structure of this paper is as follows. In Section 2, we describe the secular method and other assumptions. We give an example of the rich secular dynamics of planets in stellar triples in Section 3. In Section 4, we present the results from the population synthesis study. We discuss our results in Section 5, and conclude in Section 6.
2 Methods and assumptions
Symbol  Description  (Range of) value(s) 

Stellar mass  
Planetary mass  1  
Stellar mass  
Stellar mass  
Stellar radius  
Planetary radius  
Stellar radius  
Stellar radius  
Tidal disruption factor  2.7  
Stellar viscous timescale  
Planetary viscous timescale  
Stellar apsidal motion constant  0.014  
Planetary apsidal motion constant  0.25  
Stellar gyration radius  0.08  
Planetary gyration radius  0.25  
Stellar spin period  
Planetary spin period  
Stellar obliquity (stellar spinplanetary orbit angle)  
Planetary orbital semimajor axis  14 AU  
Period of orbit 2  1 d  
Period of orbit 3  1 d  
Planetary orbital eccentricity  0.01  
Orbit eccentricity ()  0.01  0.9  
Planetary orbital inclination  0  
Orbit inclination ()  0180  
Orbit Argument of pericenter  0360  
Orbit Longitude of ascending node  0360 

Mass ratio sampled from a flat distribution between 0.1 and 1.0.

Kippenhahn & Weigert (1994).

Lognormal distribution (Raghavan et al., 2010).

Rayleigh distribution with an rms width of 0.33 (Raghavan et al., 2010).
2.1 Configurations and notation
In panels 1 and 2 of Fig. 1, the hierarchical configurations of planets in triples considered in this work are depicted schematically. In addition to these two configurations, we included the case of a stellar binary, i.e. panels 1 or 2 of Fig. 1 with removed. Invariably, we use to denote the mass of the star that the planet is orbiting, and to denote the planetary mass. The other two stellar masses are denoted with and . Depending on the configuration, there are two to three orbits in the system, labeled 1 to 3. For example, for ‘2+2’ systems, the semimajor axis of the orbit of the planet around its star is denoted with , the semimajor axis of the orbit of stars and is denoted with , and the semimajor axis of the ‘superorbit’ of the two bound pairs is denoted with .
2.2 Secular dynamics
In Hamers et al. (2015), the equations describing the secular evolution of hierarchical fourbody systems were given for both the ‘2+2’ and ‘3+1’ configurations. Here, we used the updated code SecularMultiple (Hamers & Portegies Zwart, 2016) interfaced in the AMUSE framework (Portegies Zwart et al., 2013; Pelupessy et al., 2013). The SecularMultiple code applies to arbitrary hierarchies and numbers of bodies; evidently, its underlying algorithm reduces to the same equations for four bodies. The method is based on an expansion of the Hamiltonian in terms of ratios of the orbital separations (ratios of all combinations of small to large separations), which are assumed to be small. Subsequently, the Hamiltonian is averaged assuming Keplerian motion, and the resulting equations of motion for the orbital vectors and for all orbits are integrated numerically. Note that unlike in Hamers et al. (2015), SecularMultiple includes the ‘cross’ terms at the octupole order that depend on all three orbits simultaneously, and which were not included in the integrations of Hamers et al. (2015). Here, in addition to including all terms at the octupole order, we also included terms corresponding to pairwise binary interactions up and including the fifth order in the separation ratios. Relativistic corrections were included to the first postNewtonian (PN) order, neglecting terms in the PN potential associated with interactions between binaries.
2.3 Tidal evolution
We included tidal evolution for the planet and stars directly orbiting other bodies. For the ‘2+2’ configuration, tides were taken into account in bodies 0 and 1 in orbit 1, and bodies 2 and 3 in orbit 2 (cf. Fig. 1). For the ‘3+1’ configuration, tides were taken into account only for bodies 0 and 1 in orbit 1. For tidal evolution, we only considered twobody interactions, i.e. we neglected perturbations of other bodies on the direct tidal evolution. Evidently, indirect coupling is taken into account by the changing of the eccentricities due to secular evolution.
We adopted the equilibrium tide model of Eggleton et al. (1998), also taking into account the spin evolution (magnitude and direction) of all bodies involved in tidal evolution. For the planethosting star, we assumed zero initial obliquities. We included the effect of precession of the orbits due to tidal bulges and rotation.
The equilibrium tide model is described in terms of the viscous timescale (or equivalent related timescales), the apsidal motion constant , the gyration radius and the initial spin period . Our assumed values for the planet and stars, if applicable, are given in Table 1. Most of the values are adopted from Fabrycky & Tremaine (2007). For all tidal interactions, we assumed a constant tidal viscous timescale during the simulations. For high migration, the viscous timescale of the planet, , is the most important quantity. Apart from its simplicity, a temporally constant during high migration follows from the equations of motion with a number of physicallymotivated assumptions (Socrates & Katz, 2012).
3 Example of differences in secular evolution in stellar triples versus stellar binaries
Before proceeding with the population synthesis study in Section 4, we here briefly demonstrate some of the differences in the secular evolution of planets in stellar triple systems compared to stellar binaries. In particular, we consider a planet in a triple in the ‘2+2’ configuration in which the orbit of the planet is excited to high eccentricity, and subsequently a HJ is formed. In contrast, in the equivalent case of a stellar binary, i.e. with orbit 2 replaced by a point mass, the eccentricity of the orbit of the planet is less excited, and no HJ is formed (within 10 Gyr).
For the example system, we selected a system forming a HJ in the ‘2+2’ triple configuration in the population synthesis of Section 4. The parameters are , , , ; , , ; , , ; , , ; , , ; , , and . The other parameters are given in Table 1. Note that for these initial parameters, the initial inclination between orbits 1 and 3 is , and for orbits 2 and 3.
In the equivalent binary configuration, all orbital parameters associated with orbit 2 are removed, and they are replaced by the parameters of orbit 3 in the triple configuration. The binary companion mass, (here the tilde indicates the binary case), is now given by the sum of the masses within orbit 2 in the triple configuration, i.e. . Effectively, orbit 2 in the triple configuration is then reduced to a point mass.
In Fig. 2, we show the time evolution of the eccentricities, relative inclinations and semimajor axes for the triple and binary cases. In the binary case, the initial relative inclination of the orbit of the planet with respect to its ‘outer’ orbit 2, is . The maximum eccentricity induced in the former orbit, , is consistent with the simplest, lowestorder expectation, . This maximum eccentricity is not high enough to induce strong tidal dissipation in orbit 1 (the semimajor axis remains constant), nor trigger a tidal disruption. As expected, the relative inclination oscillates between the initial value and .
If the stellar companion to is replaced by a binary, then the evolution of orbit 1 is very different. Instead of regular LK oscillations, oscillates much more irregularly. Whereas the maximum inclination between orbit 1 and its outer orbit was limited to the initial value in the binary case, here this maximum inclination is much higher, reaching values of up to . The maximum eccentricities reached in orbit 1 are much higher, exceeding 0.9. At an age of , becomes high enough for tidal dissipation to become important, and a HJ is rapidly formed.
These differences between the binary and triple cases can be ascribed to the nodal precession of orbit 3 induced by LK oscillations acting between orbits 2 and 3. This nodal precession can enhance the LK oscillations between orbits 1 and 3 if the ratio of the LK timescales of associated with the orbital pairs (1,3) and (2,3) are commensurate. The latter condition can be quantified by considering the ratio of LK timescales as in Hamers & Portegies Zwart (2016), i.e.
(1) 
where we neglected factors of order unity in estimating the LK timescales. If , orbit 2 can be considered effectively a point mass, i.e. the binary nature of this orbit can be neglected (see also Section 2.4 of Hamers & Portegies Zwart 2016). In the intermediate regime of , the evolution is generally complex and chaotic as demonstrated in Fig. 2, and potentially very high eccentricities can be attained in orbit 1. For the system in Fig. 2, is indeed close to unity, i.e. .
The rich dynamics of the regime will be explored in more detail in a later paper.
4 Population synthesis
Using the methods described in Section 2, we carried out a population synthesis of planets in hierarchical stellar triple systems, considering two distinct hierarchical configurations (cf. panels 1 and 2 in Fig. 1). For reference, and in order to investigate the differences, we also carried out simulations of a planet in a stellar binary mimicking previous studies of high migration in stellar binaries (Wu & Murray, 2003; Fabrycky & Tremaine, 2007; Naoz et al., 2012; Petrovich, 2015a; Anderson et al., 2016; Petrovich & Tremaine, 2016). With some approximations, high migration fractions in stellar binaries can also be calculated analytically (Muñoz et al., 2016).
Consistent with formation beyond the ice line, the planet was initially placed on a nearly circular orbit () with a semimajor axis ranging between 1 and 4 AU around either the tertiary star (in the case of the ‘2+2’ configuration) or the primary star (in the case of the ‘3+1’ configuration). In the case of a stellar binary, the planet was placed on an orbit around the primary star (Stype orbit). Subsequently, we followed the secular dynamical evolution of the system, also taking into account tidal evolution.
4.1 Initial conditions
For each of the three configurations, we selected three different values of the viscous timescale of the planet and for each combination, we generated systems through Monte Carlo sampling. The viscous timescales considered were and 1.4 yr, respectively. For gas giant planets and high migration, Socrates et al. (2012) provided the constraint , by requiring that a HJ at 5 d is circularized in less than 10 Gyr. Our largest value of corresponds to this value; the other values correspond to 10 and 100 more efficient tidal dissipation. Note that corresponds to a tidal quality factor of (cf. equation 37 from Socrates et al. 2012), and .
We made the following assumptions in the Monte Carlo sampling. The mass and radius of the planet were set to and , respectively. The mass and radius of the star hosting the planet were set to and , respectively. The masses of the companion stars, and , were computed from and , where and were sampled from flat distributions between 0.1 and 1.0. The radii of the companion stars, and (of importance for stellar tides and collisions), were computed using the approximate massradius relation for (Kippenhahn & Weigert, 1994).
The initial orbit of the planet around its host star was assumed to have an eccentricity and a semimajor axis between 1 and 4 AU, sampled from a flat distribution. The orbital periods and were sampled from a normal distribution in with a mean of 5.03 and width of 2.28 (Raghavan et al., 2010), and lower and upper limits of 1 and days, respectively. The semimajor axes and were computed from these orbital periods and the sampled masses using Kepler’s law. The eccentricities and were computed from a Rayleigh distribution with an rms width of 0.33 between 0.01 and 0.9, approximating the distributions found by Raghavan et al. (2010).
In the case of ‘2+2’ systems, sampled orbital combinations were rejected if was smaller than the critical value for dynamical stability according to the criterion of Mardling & Aarseth (2001). Here, we treated orbits 2 and 3 as an isolated triple system with the mass of the ‘tertiary’ given by , which is likely a good approximation given the small relative mass of the planet. Also, we required that be larger than the critical value for dynamical stability according to the criterion of Holman & Wiegert (1999), computed treating the stars orbiting outside the planet as a single body with mass .
Similarly, in the case of ‘3+1’ systems, stability of orbits 2 and 3 was assessed using the criterion of Mardling & Aarseth (2001), treating orbit 1 as a point mass with mass . For the planetary orbit, we again applied the criterion of Holman & Wiegert (1999), this time treating orbits 1 and 2 as an isolated system (i.e. neglecting the mass ).
For the reference case of a stellar binary, only the Holman & Wiegert (1999) criterion was applied to orbits 1 and 2.
Without loss of generality, the inclination of orbit 1 was set to ; the other inclinations and were sampled from flat distributions in . The arguments of pericenter and the longitudes of the ascending nodes were sampled from flat distributions between and . These choices correspond to completely random orbital orientations.
For the other fixed initial parameters, we refer to Table 1. Regarding the stars, we adopted a constant viscous timescale of , following Fabrycky & Tremaine (2007). Assuming a tidal frequency of 4 d, this corresponds to a tidal quality factor of or , which is typical for Solartype stars (Ogilvie & Lin, 2007).
In Fig. 3, we show the initial distributions of the semimajor axes for the three configurations. Thick blue (thin green) solid lines indicate the sampled distributions of (); each row corresponds to a different configuration (note that for a stellar binary, or 3 bodies, only applies). Dashed lines show distributions from observed F and Gtype stellar triples from the sample of Tokovinin (2014), rejecting systems in which either the masses or orbital periods are unknown.
For ‘2+2’ systems, the distributions generated for the simulations approximate the observed distributions from Tokovinin (2014), except for small values of , i.e. . This is because there is an enhancement of short orbital periods in the observations which is believed to be due to LK oscillations induced by the tertiary star combined with tidal dissipation (Fabrycky & Tremaine, 2007). Note that for the ‘2+2’ configuration, tidal dissipation is included the stellar binary (orbit 2), hence the latter effect is taken into account in the simulations.
For the ‘3+1’ systems, the sampled distributions are distinctly different from the observations, especially regarding . This is due to the imposed requirement of stability of the system with the added planet orbiting one of the stars in the inner binary of the stellar triple, pushing to larger values, and, subsequently, as well. Presumably, this is also reflected in the true distributions of the semimajor axes of stellar triples with embedded planets, which are so far very poorly constrained given the small number of observed systems.
Lastly, in the case of stellar binary systems (cf. the bottom row of Fig. 3), there is a similar cutoff in imposed by dynamical stability.
4.2 Stopping conditions
The integrations were stopped if one of the following conditions was met.

The system reached an age of 10 Gyr.

A HJ was formed and no further evolution is expected within a Hubble time, i.e. and .

The innermost planet was tidally disrupted by its host star, i.e. , where is given by
(2) Here, is a dimensionless parameter; throughout, we assumed (Guillochon et al., 2011).

The orbit of the planet around its host star became dynamically unstable because of perturbations by the other stars, evaluated using the criterion of Holman & Wiegert (1999).

The stellar orbits became dynamically unstable according to the criterion of Mardling & Aarseth (2001). Note that this only applies to stellar triples.

In the ‘2+2’ case, stars collided in orbit 2. In principle, the collision remnant could subsequently secularly excite the orbit of the planet producing a HJ. However, the probability of stellar collision is small (see below).

The run time of the simulation exceeded 12 CPU hrs (imposed for practical reasons).
0.782  0.771  0.787  0.040  0.046  0.038  0.155  0.158  0.153  0.000  0.000  0.000  0.015  0.015  0.015  0.002  0.004  0.003  
0.324  0.300  0.338  0.039  0.045  0.036  0.399  0.403  0.393  0.226  0.228  0.222  0.000  0.000  0.000  0.010  0.023  0.009  
0.817  0.805  0.820  0.033  0.037  0.031  0.145  0.146  0.143  0.001  0.001  0.001  0.000  0.000  0.000  0.005  0.013  0.006  
0.788  0.772  0.793  0.015  0.021  0.016  0.174  0.182  0.169  0.000  0.000  0.000  0.013  0.013  0.013  0.002  0.005  0.003  
0.326  0.303  0.340  0.019  0.023  0.016  0.423  0.429  0.415  0.217  0.218  0.213  0.000  0.000  0.000  0.014  0.026  0.014  
0.839  0.827  0.840  0.013  0.015  0.012  0.143  0.144  0.141  0.000  0.000  0.000  0.000  0.000  0.000  0.005  0.014  0.006  
0.796  0.781  0.802  0.003  0.005  0.003  0.179  0.186  0.172  0.000  0.000  0.000  0.015  0.015  0.015  0.002  0.007  0.003  
0.332  0.305  0.346  0.002  0.004  0.002  0.444  0.450  0.433  0.208  0.210  0.202  0.000  0.000  0.000  0.012  0.028  0.014  
0.828  0.820  0.831  0.002  0.003  0.002  0.164  0.167  0.162  0.000  0.000  0.000  0.000  0.000  0.000  0.005  0.011  0.005 
4.3 Results
4.3.1 Outcome fractions
In Table 2, fractions are given of several pathways in the simulations. We distinguish between no migration, i.e. the planet’s semimajor axis not changing significantly, the formation of a HJ, the planet being tidally disrupted (‘TD’), the planet becoming dynamically unstable (‘DI, planet’), stars colliding in orbit 2 (this applies only to the ‘2+2’ cases), and the exceeding of the CPU run time. Each row corresponds to a specific hierarchical configuration (‘HC’) and viscous timescale of the planet, . The fractions are given after 5 and 10 Gyr of evolution, and by picking a random time between 0.1 and 10 Gyr for each of the simulations.
Table 2 shows that the HJ fractions between the various configurations are very similar. Nevertheless, there is a systematic trend, for all viscous timescales, that the HJ fractions for the stellar triple ‘2+2’ and ‘3+1’ configurations are higher compared to the stellar binary configuration. This confirms our expectation of Section 1, although the enhancement is only small, not exceeding a few per cent relative to the sampled systems.
A striking difference between the ‘3+1’ and the other configurations is the high fraction of tidally disrupted systems. This can be attributed to smaller values of in simulations with the ‘3+1’ configuration compared to the ‘outer’ orbits of orbit 1 in other configurations (cf. Fig. 3), which results in a stronger excitation of the eccentricity of the orbit of the planet, orbit 1, by orbit 2 (see also Section 4.3.3). Another difference of the ‘3+1’ configuration compared to the others, is a large fraction of systems in which the planet becomes dynamically unstable, of up to . This can be attributed to enhanced eccentricity of orbit 2 excited by LK oscillations of orbit 3. There are a number of outcomes of such a dynamical instability; this was investigated for a similar situation (in this case with a Ptype circumbinary planet, rather than an Stype planet) using body integrations by Hamers et al. 2016b. In the latter work, it was found that ejection of the planet, which is much less massive than the stars, is most likely. In principle, planets could also be tidally captured by stars and forming HJs in this manner (tidal effects were not taken into account in Hamers et al. 2016b), although we expect that this is not very likely. Even if a fraction of 0.1 of unstable systems would lead to tidal capture, the enhancement of the HJ fraction would at most be a few per cent. Nevertheless, calculating the contribution of HJs from dynamical instability triggered by secular evolution in these systems is an interesting avenue for future work.
Collisions between stars occur only in the ‘2+2’ configuration, in which the stars labeled and can approach each other at short distances as orbit 2 is driven to high eccentricity by orbit 3. As mentioned above, the collision remnant could subsequently secularly excite the orbit of the planet, producing a HJ. However, the probability of stellar collision is small and therefore also the contribution to forming HJs.
4.3.2 HJ orbital period distribution
In Fig. 4, we show the orbital period distributions of the shortperiod planets formed in our simulations for each configuration. In the top three panels, results are combined from the different values of . Red solid, blue dashed and green dotted lines show the distributions from the simulations after 5 Gyr, 10 Gyr and after a random time (‘’). Black crosses with error bars indicate the observations from Santerne et al. (2016), and the orbital periods of the three HJs in stellar triples known to date, WASP12, HAT8P and KELT4Ab, are indicated with the vertical black dotted lines. For reference, we also include with black open circles the distribution from Fig. 23 of Anderson et al. (2016) for and , where and is the tidal time lag of the planet (cf. Table 1 of the latter paper). With (cf. Table 1), and , or corresponds to a viscous timescale of . Distributions are normalized to unit total area. In the bottom panel, the cumulative distributions are compared for the different configurations at 10 Gyr, and for different viscous timesscales shown with different line thicknesses.
The distributions of the orbital periods are very similar for the different configurations. Compared to the case of a stellar binary and combining the different , the KS and values are and for the ‘2+2’ configuration and and for the ‘3+1’ configuration. These distributions are also similar to those of Anderson et al. (2016), although the latter are truncated at a smaller orbital period.
The dependence on the viscous timescale , shown in the bottom panel of Fig. 4, is easily understood by noting that stronger tides result in larger stalling separations of the HJ (e.g. Petrovich 2015a). Another HJ property that depends significantly on is the formation time (cf. Section 4.3.4). Other properties such as the initial semimajor axes and eccentricities are weakly dependent on , and below the corresponding results are shown for all three viscous timescales combined.
Companions of HJs in the observed sample of Santerne et al. (2016) are poorly constrained, and so it is unclear to which extent the HJ orbital period distribution is different for single stars compared to higherorder systems. The orbital periods of the three HJs observed in triple systems so far, indicated with the black vertical dotted lines in Fig. 4, are consistent with the simulations. The period of WASP12 () lies at the low end of the simulations, whereas the periods of HAT8P and KELT4Ab ( and d, respectively) lie at the peak of the simulations around 3 d.
The simulations fail to produce WJs, i.e. planets with periods between 10 and 100 d, in any meaningful numbers. The number of WJs (out of 4000 systems) at 5 Gyr, 10 Gyr and a random time are 0, 0 and 0 for the ‘2+2’ configuration, 1, 1 and 2 for the ‘3+1’ configuration, and 1, 1, and 2 for the binary configuration, respectively. These few systems appear as low bumps in the orbital period distribution between 10 and 100 days. This is in contrast to observations, which show a much larger number of planets in this period regime compared to HJs. This result is similar to previous studies of high migration in stellar binaries (e.g. Petrovich & Tremaine 2016; Antonini et al. 2016), and in multiplanet systems (Hamers et al., 2016a).
4.3.3 Orbital properties of the stellar triple arranged by outcome
In Fig 5, we show the distributions of the initial semimajor axes, arranged by various outcomes in the simulations. As mentioned above, these distributions for the HJs are weakly dependent on and the results in Fig 5 are combined for the three viscous timescales. Thick (thin) lines apply to orbit 2 (3). We distinguish between no migration (black dotted), HJ formation (red solid), planetary tidal disruption (yellow dashed) and dynamical instability of the planet (blue dotdashed).
Several interesting constraints on the semimajor axes are revealed in Fig. 5. For ‘2+2’ systems, HJs are only formed if , and needs to be and . These limits can be understood as follows. The lower limit on coincides with the lowest sampled value of (cf. the thin red solid and black dotted lines in the top panel of Fig. 5), and which do not lead to tidal disruption. The lowest sampled value of arises from the requirement of dynamical stability of the system (cf. Section 4.1 and Fig. 3). The upper limit on for HJ systems arises from the wellknown phenomenon that LK cycles are suppressed if the ‘outer’ orbit (here, orbit 3) is wide compared to the ‘inner’ orbit (here, orbit 1), in which case precession due to GR, tidal bulges and/or rotation quenches LK oscillations (e.g. Naoz et al. 2013).
This upper limit on , , translates to an upper limit on because of the requirement of dynamical stability. A typical ratio of required for dynamical stability in our simulations is (assuming the criterion of Mardling & Aarseth 2001). This implies an upper limit on of , and which is indeed roughly consistent with the upper limit found in the simulations.
Similar limits apply to the tidal disruption systems in the ‘2+2’ configuration. Most of these systems have an initial smaller compared to the HJ systems: for very small , close to dynamical instability with respect the planet (orbit 1), strong secular dynamical interactions excite very high eccentricities in orbit 1, leading to tidal disruption of the planet. The upper limit is lower compared to HJs systems, implying a lower value of the upper limit on , in this case . Similarly to the HJ systems, the upper limits on and are roughly consistent with requiring dynamical stability of orbits 2 and 3, respectively.
Considering ‘3+1’ systems, for HJ systems is larger compared to the ‘2+2’ configuration, which can be ascribed to the different hierarchy, i.e. is always required for dynamical stability for ‘3+1’ systems, whereas this is not the case for ‘2+2’ systems (for the latter, systems with similar and are typically dynamically stable unless is small). The typical values of for HJs in ‘3+1’ systems are very similar to those of for HJs in ‘2+2’ systems. These similar limits can again be explained by the requirement of dynamical stability and the suppression of LK oscillations for large ‘outer’ orbit semimajor axes (i.e. in the ‘3+1’ case, the outer orbit is orbit 2, and in the ‘2+2’ case, the outer orbit is orbit 3).
For ‘3+1’ planetary tidal disruption systems, is smallest, with . For systems in which the planet becomes dynamically unstable, lies roughly in between the typical values for tidal disruption and HJ formation. These instabilities are driven by secular eccentricity excitation of orbit 2 by the torque of orbit 3, i.e. instability would not have occurred in the absence of the fourth body.
The implications of these constraints on and in triple systems are discussed in Section 5.3.
In addition to the semimajor axes, a quantity of interest is the ratio of LK timescales applied to different orbits, with its definition depending on the hierarchy. For ‘2+2’ systems, the relevant ratio, , was introduced in Section 3 (cf. equation 1). As mentioned there, particularly high eccentricities, and therefore HJs and tidal disruptions, are expected in the regime . In the top panel of Fig. 6, the distribution of is shown for ‘2+2’ systems, arranged by the outcomes in the simulations (combining results for the three values of ). HJs systems indeed show a preference for compared to nonmigrating systems. For tidal disruption systems, the distribution is even more concentrated towards , consistent with the notion that more violent eccentricity excitation is more likely to result in tidal disruption.
For ‘3+1’ systems, the relevant LK timescale ratio is associated with the orbital pairs (1,2) and (2,3) (Hamers et al., 2015), and is given by
(3) 
As shown by Hamers et al. (2015), complex eccentricity oscillations and potentially high eccentricities are expected if . In the bottom panel of Fig. 6, we show the distribution of for the ‘3+1’ configuration, arranged by the outcomes. The HJ systems indeed display a preference for larger values of close to unity (median value of ) compared to the nonmigrating systems (median value of ).
The tidally disrupted systems, on the other hand, show a distribution of which is biased toward small values of , seemingly contradicting the above expectation. However, as already mentioned in Section 4.3.1, the fraction of tidally disrupted planets in the ‘3+1’ configuration is large owing to the small value of compared to in the ‘2+2’ configuration, and in the binary configuration (cf. Fig. 3). As shown in Fig. 5, for the tidal disruption systems is indeed strongly biased to small values, , in which case a very high eccentricity can be induced in orbit 1 regardless of orbit 3, leading to tidal disruption. Such small values of correspond to small values of (cf. equation 3).
In Fig 7, the initial distributions of the eccentricities are shown arranged by the various outcomes, similarly to Fig. 5. As might be expected, there is some preference for HJ and tidal disruption systems to have larger initial (for ‘2+2’ systems) and (for ‘3+1’ systems). The planetary dynamical instability outcome shows a strong preference for high values of . If is high initially, only a small change of due to the secular torque of orbit 3 is sufficient to drive instability of the planet.
4.3.4 HJ formation times
The distributions of the ‘end times’ associated with various outcomes of the simulations are shown in Fig. 8, each row corresponding to a different configuration. Here, ‘end time’ corresponds to the time of HJ formation or time of disruption for the tidally disrupted planets. Results are shown separately for the three assumed values of with different line thicknesses (decreasing with increasing line thickness).
As expected, the shorter values of (stronger tides) lead to shorter HJ formation times. Similarly to the distribution of the HJ orbital periods, the HJ formation times are not significantly different between the three configurations. The red open circles show data for simulations in stellar binaries from the second panel of Fig. 22 () by Anderson et al. (2016), and which are roughly consistent with our simulations. Tidal disruptions, however, occur earlier in our simulations of stellar binaries compared to Anderson et al. (2016) (cf. the yellow open circles in Fig. 8), and this may be because of different assumptions on the orbital parameters. There are no large differences in the tidal disruption time distributions between the different configurations.
4.3.5 Stellar obliquities
In Fig. 9, we show the distributions of the stellar obliquity, i.e. the angle between the spin of the planethosting star and the orbit of the planet, for the nonmigrating (black dotted lines) and HJ systems (red solid lines). Results are combined for the three values of . Data from Fig. 24 of Anderson et al. (2016) for and are shown with black open circles. Additionally, we show with black stars observational data adopted from Lithwick & Wu (2014). The initial obliquities in the simulations were assumed to be zero.
In the case of a stellar binary, there are two distinct, but broad, peaks in the obliquity distribution near and . Such peaks in the obliquity distribution are well known to arise in stellar binaries (Fabrycky & Tremaine 2007; Naoz & Fabrycky 2014; Anderson et al. 2016; see Storch et al. 2016 for a detailed study on the origin of the bimodel distribution). There appears to be a discrepancy with Anderson et al. (2016), who find peaks near values of and . This may be due various reasons, including the initially shorter spin period of rather than 10 d for the host star, magnetic braking, which was included in Anderson et al. (2016) but not in this work (magnetic braking results in the spin down of the star, which, due to spinorbit coupling, can significantly change the final obliquity, Storch et al. 2016), or different .
For stellar triples, the obliquity distributions seem marginally less peaked compared to the case of a stellar binary. In the bottom panel of Fig. 9, the cumulative distributions between the three configurations are compared directly. Compared to the case of a stellar binary, the KS and values are and for the ‘2+2’ configuration and and for the ‘3+1’ configuration. The differences in the obliquity distributions in stellar triples compared to binaries might be ascribed to the more chaotic nature of the eccentricity and inclination oscillations in the regime compared to binaries, in which the regime for chaotic evolution is much smaller (i.e. only for marginally hierarchical systems, e.g. Li et al. 2014). This implies larger possible variations in inclinations at times of maximum eccentricity and therefore less peaked distributions of the obliquity after tidal dissipation has shrunk the planetary orbit.
4.3.6 Initial inclinations
Lastly, we consider the initial inclinations. In Fig. 10, the distributions of the initial mutual inclination between the planetary orbit and its parent orbit are shown, arranged by the simulation outcomes (combining results for the three ). This mutual inclination is for ‘3+1’ and ‘binary’ systems; for the ‘2+2’ configuration, it is (cf. Fig. 1). We recall that the initial orbital orientations were assumed to be random (cf. Section 4.1). As expected, high mutual inclinations, near , are required to excite high eccentricities leading to tidal disruptions or HJs. There are some differences between the triple and binary star cases: for the former, the range of inclinations for tidal disruption and HJ systems around is larger, even extending to inclinations lower than or higher than . The enlargement of the inclination window for high eccentricities is most pronounced for the ‘3+1’ configuration, for which HJs can be formed even for nearly coplanar systems (prograde or retrograde). For ‘3+1’ configurations, this can be understood by noting that orbits 1 and 2 can be made highly inclined by changes of the absolute inclination of orbit 2 driven by the torque of orbit 3 (Hamers et al., 2015).
In Fig. 11, we show a similar figure with the distributions of the initial mutual inclination between orbits 2 and 3 (stellar orbits only). This inclination only applies to the triple star cases. For the ‘2+2’ configuration, the tidal disruption and HJ systems show no strong dependence on , although there is a preference for tidal disruption systems and a lack of HJ systems near ; the origin for this is unclear. Regarding the ‘3+1’ configuration, systems in which the planet becomes dynamically unstable show a strong preference for high inclinations of near . For such high inclinations, the maximum eccentricity excited in orbit 2 is high, making a dynamical instability with respect to orbit 1 likely. This effect prevents a large fraction of highly inclined systems to become tidal disruption or HJ systems. The latter, instead, show peaks in the initial distribution of at lower inclinations near and higher inclinations near .
5 Discussion
5.1 Caveats of the simulations
The systems in our simulations were sampled assuming a planet initially orbiting its host star between and and assuming lognormal distributions for the stellar orbits, rejecting systems that are dynamically unstable according to (approximate) analytic stability criteria (cf. Section 4.1). Although this method excludes systems that are almost certainly dynamically unstable, it may not reflect the intrinsic population of stellar triples with planets, for which the orbital distributions are currently essentially unconstrained. In other words, aside from dynamical stability constraints, it is unclear if the planet formation process(es) that occur in stellar triples with one (or more) planet(s) have other effects on the properties of triples. We believe that the method adopted in this work is warranted, given the lack of observational constraints and better alternatives. Nevertheless, it should be taken into account that the our results depend on the assumed orbital distributions, most importantly the distributions of the semimajor axes and eccentricities, and the assumed viscous timescales. The latter dependence is shown explicitly in Table 2, and in Figs. 4 and 8.
It was mentioned in Section 4.3.1 that the imposed run time of 12 CPU hrs was exceeded in a few per cent of systems, but that these systems are unlikely to produce HJs if the integrations had not been stopped before reaching an age of 10 Gyr. In these systems, the LK timescales are extremely short compared to the integration time of 10 Gyr, such that not all cycles could be completed within the CPU time limit. These systems have distinct orbital properties compared to the HJ systems, and are therefore unlikely to produce HJs if the integrations had not been stopped before completion. This is illustrated in Figures 6 and 8, where the distributions of and and of the end times, respectively, are shown for the run time exceeded systems with the thick black dashed lines. The ratios and are typically very small for the run time exceeded systems, implying very short LK timescales for the innermost pairs (i.e. orbits 13 for the ‘2+2’, and orbits 12 for the ‘3+1’ configurations). Also, the ages at which the simulations were terminated (cf. Fig. 8) are typically longer than the HJ formation times.
5.2 Comparisons of the HJ fractions to other high migration scenarios and observations
The highest HJ fraction in stellar triples found in our simulations was , compared to in our fiducial binary star simulations (cf. Table 2). The enhancement due to the richer (secular) dynamics therefore turn out to be only small. We note, however, that the fraction of tidal disruptions in ‘3+1’ configurations is markedly higher compared to binaries (up to ), mainly due to the planet+star subsystem being orbited by a relatively tight outer orbit. This finding of apparently low HJ fractions is similar to the result of Hamers et al. (2016a) for secular eccentricity excitation in multiplanet systems (between 3 and 5 planets orbiting a single star). In the latter work, it was found that, although potentially very high eccentricities can be attained in orbit of the innermost planet, these high eccentricities typically lead to a ‘fast’ tidal disruption, rather than a slower process of orbital tidal decay.
To compare our simulated HJ fractions to observations, we assume a stellar triple fraction of (e.g. Tokovinin 2014) and a giant planet fraction in triples of . Adopting our highest fraction of , the expected HJ fraction around Solartype MS stars due to our scenario of secular high ‘triple’ migration is . In other words, given a Solartype MS star (either single or multiple), the fraction of these stars around which a HJ is expected to form through our scenario is . The observed HJ fraction around Solartype MS stars is (e.g. Wright et al. 2012), i.e. 20 times larger. Therefore, our ‘triple’ scenario cannot resolve the existing discrepancy that the predicted HJ formation rates of high migration in stellar binaries are too low compared to observations. This is because the formation efficiency is only marginally higher, whereas the fraction of stellar triple systems is lower compared to binaries.
5.3 Constraints on stellar triples hosting HJs and comparisons to observations
As discussed in Section 4.3.3, we find that HJs can only be formed through high migration in stellar triples with specific orbital configurations (cf. Fig. 5). The three HJs observed in triples so far, WASP12b (Hebb et al., 2009; Bergfors et al., 2013; Bechter et al., 2014), HATP8b (Latham et al., 2009; Bergfors et al., 2013; Bechter et al., 2014), and KELT4Ab (Eastman et al., 2016), orbit the tertiary star of the triple, i.e. the ‘2+2’ configuration applies. For this configuration, we found the constraints and .
The projected separations of the stellar binaries in WASP12b and HATP8 were found by Bechter et al. (2014) to be and , respectively. For simplicity, we adopt these separations as semimajor axes, i.e. and for these systems, respectively. For KELT4Ab, Eastman et al. (2016) found a projected separation of for the stellar binary, which we adopt as . In all three cases, the semimajor axis satisfies .
Assuming an angular separation of between WASP12b and its stellar binary companion (Bergfors et al., 2013) and a distance of 427 pc (Chan et al., 2011), the projected separation is . For HATP8, an angular separation of between HATP8 and its stellar binary companion (Bergfors et al., 2013) and a distance of 230 pc (Latham et al., 2009) imply a projected separation of . Eastman et al. (2016) found a projected separation for KELT4Ab of . Again adopting the projected separations as semimajor axes, the semimajor axes of these systems satisfy .
This shows that, at least approximately given that only the projected separations are known, the observed HJs are consistent with the constraints given by the simulations. Moreover, if a HJ were to be found in a stellar triple with distinctly different orbital parameters (e.g. for the ‘2+2’ configuration), this would indicate HJ formation through another mechanism such as disk migration.
6 Conclusions
We have considered the formation of HJs through high migration in stellar triple systems. In this scenario, the orbit of a Jupiterlike planet is excited to high eccentricity by secular perturbations from the stars. The subsequent close pericenter passages lead to tidal decay of the orbit, potentially producing a HJ. We have considered two hierarchical configurations of the planet in the stellar triple system: the planet orbiting the tertiary star (‘2+2’ configuration; cf. panel 1 of Fig. 1), and the planet orbiting one of the stars in the inner binary (‘3+1’ configuration; cf. panel 2 of Fig. 1). For reference, we also carried out simulations in stellar binaries, in which case the planet is orbiting one of the stars in an Stype orbit. Our main conclusions are as follows.
1. Although the HJ fraction in stellar triple systems is larger compared to stellar binaries, the enhancement is small, at most a few tens of per cent. In our simulations, the HJ fraction (at an age of 10 Gyr) is largest for the ‘2+2’ systems, i.e. assuming a viscous timescale of for the planet. The corresponding fraction for ‘3+1’ systems is , compared to for stellar binaries. For longer viscous timescales, the fractions are lower but the trend of slightly larger HJ fractions for triples persists.
2. The orbital period distributions of the HJs in our simulations are very similar for the different configurations, peaking around 3 d. The orbital periods of the three HJs in stellar triples observed so far are consistent with the orbital period distributions from the simulations. No significant number of WJs is produced in the simulations for any of the configurations.
3. In our simulations, HJs are formed only for specific ranges of the semimajor axes of the orbits of the triple system (cf. Section 4.3.3). Assuming the HJ initially orbited its host star at a separation between 1 and 4 AU, for ‘2+2’ systems the semimajor axis of the stellar binary orbit should satisfy . Also, the semimajor axis of the ‘outer’ orbit, , should satisfy .
These constraints are approximately consistent with the three HJs observed so far in hierarchical triple systems, which are all of the ‘2+2’ type. Moreover, a future detection of a HJ violating these conditions would be strong indication for another formation mechanism of HJs in these triples, such as disk migration.
For HJhosting ‘3+1’ configurations, the semimajor axis of the first orbit around the star+planet subsystem should satisfy . For the the second orbit, the semimajor axis should be .
4. For similar ranges of semimajor axes (for details, see Fig. 5), we also expect planets to be tidally disrupted by the star. This suggests possible enrichment of stars by planets in hierarchical triple systems satisfying specific orbital configurations. Alternatively, if the planet is not completely destroyed during the disruption event it may be partially or completely stripped of its envelope (Liu et al., 2013), producing hot Neptunes or (shortperiod) superEarths in these triple systems.
Acknowledgements
We thank Dong Lai for stimulating discussions and comments on the manuscript, and the anonymous referee for a helpful report. ASH gratefully acknowledges support from the Institute for Advanced Study.
References
 Anderson et al. (2016) Anderson K. R., Storch N. I., Lai D., 2016, MNRAS, 456, 3671
 Antonini et al. (2016) Antonini F., Hamers A. S., Lithwick Y., 2016, preprint, (arXiv:1604.01781)
 Bate et al. (2010) Bate M. R., Lodato G., Pringle J. E., 2010, MNRAS, 401, 1505
 Batygin (2012) Batygin K., 2012, Nature, 491, 418
 Batygin et al. (2016) Batygin K., Bodenheimer P. H., Laughlin G. P., 2016, ApJ, 829, 114
 Beaugé & Nesvorný (2012) Beaugé C., Nesvorný D., 2012, ApJ, 751, 119
 Bechter et al. (2014) Bechter E. B., et al., 2014, ApJ, 788, 2
 Bergfors et al. (2013) Bergfors C., et al., 2013, MNRAS, 428, 182
 Bodenheimer et al. (2000) Bodenheimer P., Hubickyj O., Lissauer J. J., 2000, Icarus, 143, 2
 Boley et al. (2016) Boley A. C., Granados Contreras A. P., Gladman B., 2016, ApJ, 817, L17
 Chan et al. (2011) Chan T., Ingemyr M., Winn J. N., Holman M. J., SanchisOjeda R., Esquerdo G., Everett M., 2011, AJ, 141, 179
 Chatterjee et al. (2008) Chatterjee S., Ford E. B., Matsumura S., Rasio F. A., 2008, ApJ, 686, 580
 Dvorak (1982) Dvorak R., 1982, Oesterreichische Akademie Wissenschaften Mathematisch naturwissenschaftliche Klasse Sitzungsberichte Abteilung, 191, 423
 Eastman et al. (2016) Eastman J. D., et al., 2016, AJ, 151, 45
 Eggleton et al. (1998) Eggleton P. P., Kiseleva L. G., Hut P., 1998, ApJ, 499, 853
 Evans (1968) Evans D. S., 1968, QJRAS, 9, 388
 Fabrycky & Tremaine (2007) Fabrycky D., Tremaine S., 2007, ApJ, 669, 1298
 Ford & Rasio (2008) Ford E. B., Rasio F. A., 2008, ApJ, 686, 621
 Foucart & Lai (2011) Foucart F., Lai D., 2011, MNRAS, 412, 2799
 Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
 Guillochon et al. (2011) Guillochon J., RamirezRuiz E., Lin D., 2011, ApJ, 732, 74
 Hamers & Portegies Zwart (2016) Hamers A. S., Portegies Zwart S. F., 2016, MNRAS, 459, 2827
 Hamers et al. (2015) Hamers A. S., Perets H. B., Antonini F., Portegies Zwart S. F., 2015, MNRAS, 449, 4221
 Hamers et al. (2016a) Hamers A. S., Antonini F., Lithwick Y., Perets H. B., Portegies Zwart S. F., 2016a, MNRAS,
 Hamers et al. (2016b) Hamers A. S., Perets H. B., Portegies Zwart S. F., 2016b, MNRAS, 455, 3180
 Hebb et al. (2009) Hebb L., et al., 2009, ApJ, 693, 1920
 Holman & Wiegert (1999) Holman M. J., Wiegert P. A., 1999, AJ, 117, 621
 Jurić & Tremaine (2008) Jurić M., Tremaine S., 2008, ApJ, 686, 603
 Kippenhahn & Weigert (1994) Kippenhahn R., Weigert A., 1994, Stellar Structure and Evolution
 Kozai (1962) Kozai Y., 1962, AJ, 67, 591
 Lai (2014) Lai D., 2014, MNRAS, 440, 3532
 Latham et al. (2009) Latham D. W., et al., 2009, ApJ, 704, 1107
 Lee et al. (2014) Lee E. J., Chiang E., Ormel C. W., 2014, ApJ, 797, 95
 Li et al. (2014) Li G., Naoz S., Holman M., Loeb A., 2014, ApJ, 791, 86
 Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
 Lin & Papaloizou (1986) Lin D. N. C., Papaloizou J., 1986, ApJ, 309, 846
 Lin et al. (1996) Lin D. N. C., Bodenheimer P., Richardson D. C., 1996, Nature, 380, 606
 Lithwick & Wu (2011) Lithwick Y., Wu Y., 2011, ApJ, 739, 31
 Lithwick & Wu (2014) Lithwick Y., Wu Y., 2014, Proceedings of the National Academy of Science, 111, 12610
 Liu et al. (2013) Liu S.F., Guillochon J., Lin D. N. C., RamirezRuiz E., 2013, ApJ, 762, 37
 Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
 Mayor & Queloz (1995) Mayor M., Queloz D., 1995, Nature, 378, 355
 Mazeh et al. (2015) Mazeh T., Perets H. B., McQuillan A., Goldstein E. S., 2015, ApJ, 801, 3
 Muñoz et al. (2016) Muñoz D. J., Lai D., Liu B., 2016, MNRAS, 460, 1086
 Nagasawa et al. (2008) Nagasawa M., Ida S., Bessho T., 2008, ApJ, 678, 498
 Naoz & Fabrycky (2014) Naoz S., Fabrycky D. C., 2014, ApJ, 793, 137
 Naoz et al. (2012) Naoz S., Farr W. M., Rasio F. A., 2012, ApJ, 754, L36
 Naoz et al. (2013) Naoz S., Kocsis B., Loeb A., Yunes N., 2013, ApJ, 773, 187
 Ogilvie & Lin (2007) Ogilvie G. I., Lin D. N. C., 2007, ApJ, 661, 1180
 Pejcha et al. (2013) Pejcha O., Antognini J. M., Shappee B. J., Thompson T. A., 2013, MNRAS, 435, 943
 Pelupessy et al. (2013) Pelupessy F. I., van Elteren A., de Vries N., McMillan S. L. W., Drost N., Portegies Zwart S. F., 2013, A&A, 557, A84
 Petrovich (2015a) Petrovich C., 2015a, ApJ, 799, 27
 Petrovich (2015b) Petrovich C., 2015b, ApJ, 805, 75
 Petrovich & Tremaine (2016) Petrovich C., Tremaine S., 2016, ApJ, 829, 132
 Portegies Zwart et al. (2013) Portegies Zwart S., McMillan S. L. W., van Elteren E., Pelupessy I., de Vries N., 2013, Computer Physics Communications, 183, 456
 Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
 Rasio & Ford (1996) Rasio F. A., Ford E. B., 1996, Science, 274, 954
 Santerne et al. (2016) Santerne A., et al., 2016, A&A, 587, A64
 Socrates & Katz (2012) Socrates A., Katz B., 2012, preprint, (arXiv:1209.5723)
 Socrates et al. (2012) Socrates A., Katz B., Dong S., 2012, preprint, (arXiv:1209.5724)
 Storch et al. (2016) Storch N. I., Lai D., Anderson K. R., 2016, preprint, (arXiv:1607.03937)
 Tanaka et al. (2002) Tanaka H., Takeuchi T., Ward W. R., 2002, ApJ, 565, 1257
 Tokovinin (2014) Tokovinin A., 2014, AJ, 147, 87
 Winn et al. (2011) Winn J. N., et al., 2011, AJ, 141, 63
 Wright et al. (2012) Wright J. T., Marcy G. W., Howard A. W., Johnson J. A., Morton T. D., Fischer D. A., 2012, ApJ, 753, 160
 Wu & Lithwick (2011) Wu Y., Lithwick Y., 2011, ApJ, 735, 109
 Wu & Murray (2003) Wu Y., Murray N., 2003, ApJ, 589, 605
 Xue & Suto (2016) Xue Y., Suto Y., 2016, ApJ, 820, 55