Hot Jupiters driven by high-eccentricity migration in globular clusters

Hot Jupiters driven by high-eccentricity migration in globular clusters

Adrian S. Hamers and Scott Tremaine Institute for Advanced Study, School of Natural Sciences, Einstein Drive, Princeton, NJ 08540, USA hamers@ias.edu
Abstract

Hot Jupiters (HJs) are short-period giant planets that are observed around of solar-type field stars. One possible formation scenario for HJs is high-eccentricity (high-) migration, in which the planet forms at much larger radii, is excited to high eccentricity by some mechanism, and migrates to its current orbit due to tidal dissipation occurring near periapsis. We consider high- migration in dense stellar systems such as the cores of globular clusters (GCs), in which encounters with passing stars can excite planets to the high eccentricities needed to initiate migration. We study this process via Monte Carlo simulations of encounters with a star+planet system including the effects of tidal dissipation, using an efficient regularized restricted three-body code. HJs are produced in our simulations over a significant range of the stellar number density . Assuming the planet is initially on a low-eccentricity orbit with semimajor axis , for the encounter rate is too low to induce orbital migration, whereas for HJ formation is suppressed because the planet is more likely ejected from its host star, tidally disrupted, or transferred to a perturbing star. The fraction of planets that are converted to HJs peaks at for intermediate number densities of . Warm Jupiters, giant planets with periods between 10 and 100 days, are produced in our simulations with an efficiency of up to . Our results suggest that HJs can form through high- migration induced by stellar encounters in the centers of of dense GCs, but not in their outskirts where the densities are lower.

Subject headings:
planets and satellites: dynamical evolution and stability – globular clusters: general – gravitation – scattering

1. Introduction

Globular clusters (GCs) are among the oldest and densest stellar systems known. Recent findings have challenged our understanding of these systems. Whereas GCs were traditionally thought to have formed in a single starburst, in the last decade multiple stellar populations have been observed in a large fraction of GCs, and the origin of these populations is still actively debated (see, e.g., Gratton et al. 2012 for a review). Another puzzling find is that searches for planets around stars in GCs have been unsuccessful. In particular, no planets were found in extensive Hubble Space Telescope observations of 47 Tuc. This failure was originally interpreted to imply that the occurrence rate of short-period planets around stars in 47 Tuc is at least an order of magnitude lower than for field stars (Gilliland et al., 2000). However, a recent study by Masuda & Winn (2017) has shown that this estimate should be revised, given the now better-known radius and period relations of giant planets around field stars; they find that the expected number of planets is only , so the null result is marginally consistent with the abundance of planets around field stars (Masuda & Winn, 2017). Similar arguments apply to other surveys of planets in GCs, in particular the ground-based surveys of 47 Tuc by Weldrake et al. (2005) and of Cen by Weldrake et al. (2008). Therefore, the existence of short-period giant planets in GCs is still an open question.

A deficit of short-period planets in GCs would suggest that planet formation in GCs is inefficient, and/or that dynamical interactions in these dense environments destroy such planets after they form. Inhibited planet formation might, for example, be due to radiation from nearby massive stars (Armitage, 2000; Adams et al., 2004; Thompson, 2013). Alternatively, GCs have low metallicities and the giant-planet occurrence rate is known to correlate with metallicity in field stars (Santos et al., 2001; Fischer & Valenti, 2005). However, it is unclear whether this relation also applies to GCs, whose formation history is not well-understood.

Dynamical interactions can disrupt planetary systems (e.g., Sigurdsson 1992), but they can also enhance the numbers of short-period planets by exciting high planetary eccentricities: in particular, if the periapsis of a giant-planet orbit becomes as small as a few stellar radii, tidal dissipation in the planet excited by interactions with the host star may become strong enough to drive migration of the planet to a tight circular orbit with a period of a few days, creating a hot Jupiter (HJ) — the easiest class of planet to detect in transit surveys. Once the planet has migrated to such a tightly bound orbit, it becomes immune to further perturbations from passing stars111In some cases, further orbital decay driven by tidal dissipation in the star could shrink the orbit until the planet is disrupted by the star. However, the relatively high occurrence rate of HJs around field stars indicates that stellar tidal dissipation is typically inefficient in solar-type stars, and we shall ignore this process in the present paper..

High-eccentricity (high-) migration around field stars has been widely studied. The possible eccentricity excitation mechanisms include planet-planet scattering (Rasio & Ford, 1996; Chatterjee et al., 2008; Ford & Rasio, 2008; Jurić & Tremaine, 2008; Nagasawa et al., 2008; Beaugé & Nesvorný, 2012); Lidov-Kozai (LK) oscillations (Lidov, 1962; Kozai, 1962) in binary-star systems (Wu & Murray, 2003; Fabrycky & Tremaine, 2007; Naoz et al., 2012; Petrovich, 2015a; Anderson et al., 2016; Petrovich & Tremaine, 2016; Hamers, 2017b), triple-star systems (Hamers, 2017a; Grishin et al., 2017), and multiplanet systems (Petrovich, 2015b; Xue & Suto, 2016); and secular chaos in multiplanet systems (Wu & Lithwick, 2011; Lithwick & Wu, 2011, 2014; Hamers et al., 2017). Although many authors have studied the dynamics of planets in open clusters and GCs (e.g., Sigurdsson 1992; de La Fuente Marcos & de La Fuente Marcos 1997; Bonnell et al. 2001; Davies & Sigurdsson 2001; Fregeau et al. 2006; Spurzem et al. 2009; Malmberg et al. 2011; Boley et al. 2012; Chatterjee et al. 2012; Parker & Quanz 2012; Hao et al. 2013; Liu et al. 2013; Li & Adams 2015; Wang et al. 2015; Zheng et al. 2015; Cai et al. 2017), the interplay between perturbations from passing stars and dissipative planetary tides has scarcely been investigated. Shara et al. (2016) considered the formation of HJs in two-planet systems in open clusters and found that HJ formation occurs in of the planetary systems. To our knowledge, no study has focused on similar processes in the much denser GCs.

In this paper, we study the formation of HJs in GCs through high- migration induced by passing stars. A computational challenge in this problem is the wide range of timescales: from a few days for the orbital period of an HJ, to 30 yr for the encounter time of a star with impact parameter and relative velocity , to for the lifetime of a GC. We approach this problem by using an efficient regularized restricted three-body code that includes tides and general relativistic corrections. This method allows us to simulate the cumulative effect of encounters over the lifetime of a GC, which is a prohibitive endeavor using general-purpose direct -body integrators.

The plan of the paper is as follows. In Section 2, we present our regularized restricted three-body code and test it using three-body integrations with a highly accurate -body code. In Section 3, we apply our method to a population-synthesis study of planets in the centers of dense GCs, and we describe the properties of the migrating planets. In addition, we present an analytic model that approximately describes the period distribution of the HJs that are formed in the simulations. We discuss our results in Section 4, and we conclude in Section 5.

2. Methodology

2.1. Numerical integration method

Consider a star of mass orbited by a planet with semimajor axis and mass . The orbit of the planet, also referred to as the “binary”, is perturbed by a passing star with mass . We assume that , in which case the planet can be interpreted as a test particle that does not affect the stellar motion. The effect of the encounter on the planetary orbit can be calculated numerically using direct -body integration. However, rather than using standard integration methods, we take advantage of the fact that we are dealing with a restricted three-body problem. In particular, we assume that the perturber moves on a hyperbolic orbit having separation with respect to the planet-hosting star, and we regularize the motion of the planet with respect to its host star. We use Kustaanheimo-Stiefel (KS) regularization (Kustaanheimo & Stiefel, 1965; Stiefel & Scheifele, 1971) with the time transformation , where and are the physical and fictitious times, respectively, and is the planet-host star separation. This approach allows us to compute the effect of the perturber on the planetary orbit with a factor of 20-100 performance increase compared to a direct -body code, without much loss of accuracy (see Section 2.2 below).

In addition to the gravitational force from the point-mass perturber, we also include general relativistic corrections to the first post-Newtonian (PN) order and tidal effects induced in the planet by its host star. The latter are required to be able to model the formation of HJs. The perturbing acceleration to the regularized motion of the planet is given by

(1)

Here, is given by

(2)

where and are, respectively, the relative separation vectors of the planet and the perturber with respect to the host star; the first and second terms in equation (2) are the ‘direct’ and ‘indirect’ terms, respectively.

We take into account the planetary tidal evolution with the equilibrium tide model. Here, we assume that the planetary rotational spin vector is always aligned with the normal to the planetary orbit. Generally, the expression for the tidal perturbing acceleration is then (cf. Eq. 8 of Hut 1981)

(3)

where , , and are the planetary radius, apsidal motion constant, and tidal time lag; the overhead dot denotes time derivatives, , is the true anomaly, is the azimuthal unit vector in the orbital plane, pointing in the direction of increasing , and is the planetary spin frequency. The planet’s rotational angular momentum is much smaller than its orbital angular momentum; therefore, the tidal torque is small as well. This justifies the assumption that the orbit-averaged tidal torque is zero, i.e., , which is equivalent to or , where is a function of eccentricity given by

(4)

with the mean motion. Equation (4) is equivalent to pseudosynchronization, i.e., (cf. equation 42 of Hut 1981). In the numerical integrations, we set , where is the instantaneous eccentricity. We note that this treatment of tides is slightly different from Wu & Lithwick (2011; cf. Eq. A7) and Antonini et al. (2016; cf. Eq. 6), who set the term in proportional to to zero.

Lastly, the relativistic perturbing acceleration is (Einstein et al., 1938)

(5)

2.2. Validation

Figure 1.— Comparison of the regularized restricted 3-body code (RR3; red dashed lines) with 3-body integrations carried out with ARCHAIN (green solid lines). The left and middle panels show the planetary semimajor axis and eccentricity, respectively, as a function of time. The right panels show the absolute value of the relative error in the semimajor axis (black solid lines) and the eccentricity (black dashed lines) in RR3 relative to ARCHAIN. Each row corresponds to a different value of the perturber’s closest approach , indicated above each panel.

We test our regularized restricted three-body code (RR3) with integrations carried out with ARCHAIN (Mikkola & Merritt, 2006, 2008), a high-accuracy -body code that uses chain regularization. We use ARCHAIN implemented within AMUSE (Portegies Zwart et al., 2013; Pelupessy et al., 2013). In all validation integrations shown in this Section, we only include Newtonian point-mass dynamics, i.e., we set in equation (1), and the speed of light in ARCHAIN is set to (i.e., effectively infinity) to eliminate relativistic effects. In contrast, in the simulations of Section 3 below, both tidal and general relativistic effects are taken into account.

The following initial conditions are assumed, representing some of the strong and weak encounters in the population-synthesis calculations of Section 3. The masses are set to , (the Jupiter mass ), and . The initial binary semimajor axis and eccentricity are and . The perturber’s velocity at infinity is , and its periapsis distance is taken to be either 1, 3, or 10 . The binary’s inclination with respect to the perturber is zero (i.e., the binary orbit is prograde relative to the perturber), and the longitudes of periapsis of the binary and perturber initially differ by . The initial true anomaly of the binary is zero. The duration of the integrations is with the periapsis passage occurring at ; for the adopted values of , , 1.4 and 4.7 yr, respectively. Note that these tests are particularly challenging because (i) the orbits are coplanar and (ii) when the unperturbed orbits of the planet and the perturber can collide (in the validation tests, we assume that the radii of all three bodies are zero and hence do not check for the occurrence of collisions).

In the left and middle panels of Fig. 1, we show the semimajor axis and eccentricity as a function of time. The right panels show the relative errors in the semimajor axis and the eccentricity in RR3 as a function of time, assuming the ARCHAIN code gives the exact result. Each row corresponds to a different value of , indicated above each panel.

The encounter with results in a destructive perturbation to the planetary orbit, i.e., the planet becomes unbound from its host star ( and ). Despite the large perturbation, RR3 computes the final semimajor axis and eccentricity with error compared to ARCHAIN. The performance increase with respect to ARCHAIN is a factor of .

The encounter with is less destructive, leaving the planet bound to its host star. The eccentricity increases from 0.2 to and the semimajor axis decreases by . The error with respect to ARCHAIN is again , and the performance increase is a factor of .

Lastly, the case corresponds to a ‘secular’ encounter in which the angular speed of the perturber at periapsis is much lower than the binary mean motion, i.e., the associated ratio of these quantities,

(6)

satisfies (in this formula, we assumed ). For , , indicating that the encounter is of the secular type (in contrast, and for and 3 , respectively). Secular encounters produce a permanent change in the eccentricity but no permanent change in semimajor axis. This is indeed the case in the bottom row of Fig. 1: the semimajor axis returns to its original value after being perturbed by , whereas the eccentricity is changed by after the encounter. The error made with respect to ARCHAIN is less than , whereas the performance gain is a factor of . Generally, the performance gain with RR3 tends to increase for more secular encounters.

3. Population-synthesis study

3.1. Gravitational dynamics and tidal evolution

We apply the regularized restricted three-body code RR3 described in Section 2 to model numerically the effects of encounters on planetary orbits in dense stellar systems such as GCs. The planet is assumed to have an initial semimajor axis , 2, or 4 . These values approximately span the range in which most known giant exoplanets are found, although many as-yet undetected giant planets are likely to be present at larger semimajor axes. The initial planetary eccentricity is assumed to follow from a Rayleigh distribution with an rms value of 0.33 (Jurić & Tremaine, 2008), cut off at a maximum eccentricity of 0.6. The resulting distribution has an rms value of , which is close to the rms value of for planets with periods above 10 days around stars in the solar neighborhood222Rms value obtained from www.exoplanets.org on October 30 2017.. We checked that there is little to no dependence of our results on the initial eccentricity.

In addition to the Newtonian accelerations due to the host and perturbing stars, we include the lowest-order precession of the planetary apsides due to general relativity and tidal evolution of the planetary orbit. These effects are implemented with the assumptions described in Section 2.1, and the associated parameters , , and are given in the top part of Table 1. The tidal time lag is assumed to be constant (Socrates & Katz, 2012). We set , for which an HJ with a 5-day orbital period is circularized in less than 10 Gyr (Socrates et al., 2012). This time lag corresponds to a tidal quality factor (cf. equation 37 from Socrates et al. 2012), and . We assume that stellar tides are negligible.

3.2. Generating encounters

3.2.1 Flux of perturbing stars

In our simulations, encounters are generated continuously until the current age of the GC is reached (assumed to be 10 Gyr). To sample the encounters, we assume a locally homogeneous stellar background with stellar number density and one-dimensional velocity dispersion independent of stellar mass. We assume a Maxwellian stellar velocity distribution at large distances from the host star, such that the distribution function is

(7)

where is the fraction of stars with masses in the interval and is the relative velocity dispersion (Binney & Tremaine, 2008). When the gravitational attraction of the host star is included, the distribution function must still be a function of the energy , where and is the distance from the host star. Therefore, the distribution function is modified to

(8)

for , and zero otherwise.

We now introduce an imaginary ‘encounter sphere’ centered at the host star with a radius . We also choose , where is the length scale beyond which and are no longer (approximately) constant. Stars impinging on the encounter sphere are considered to be perturbers. Using equation (8), we find that the number density of perturbers at the encounter sphere within a mass range and with velocities between and is

(9)

where is the Heaviside step function. Integration of equation (3.2.1) over all perturber masses and velocities gives

(10)

where

(11)

and is the complementary error function. Therefore, the fraction of perturbers at the encounter sphere with mass is proportional to ; we will use this result in our sampling procedure described below (Section 3.2.2). Equation (10) shows that the stellar number density at the encounter sphere, , is larger than due to gravitational focusing (note that for ). Henceforth, when we use the term “number density” we will always be referring to .

Consider a point on the encounter sphere with position vector relative to the host star. Next, define a local coordinate system centered on this point in which the axis is directed toward the host star, i.e., , and the and axes lie on the tangent plane of on the encounter sphere. The differential flux of stars into the encounter sphere is given by (Hénon, 1972), i.e.,

(12)

Integrating the differential flux over all perturber masses, velocities, and the entire encounter sphere, we obtain a total encounter rate of

(13)

In the limit of large (, i.e., weak or negligible gravitational focusing), equation (3.2.1) reduces to

(14)

independent of the perturber mass function.

3.2.2 Numerical sampling procedure

We generate encounters using the following procedure.

  1. The perturber mass is assumed to follow a Salpeter distribution (Salpeter, 1955), modified to account for the finite lifetime of stars and the gravitational focusing implied by equation (10).

    Specifically, an initial mass is sampled from a Salpeter distribution, , with lower and upper limits 0.1 and 100 . Using the SSE stellar evolution code (Hurley et al., 2000) as implemented in AMUSE (Portegies Zwart et al., 2013; Pelupessy et al., 2013) and assuming a metallicity , this initial mass is replaced by the mass after 5 Gyr of stellar evolution, . Subsequently, we compute the associated value of and (eq. 11), and reject the sampled mass if , where is a random number between 0 and 1, and is the maximum value of over the allowed range of .

    We do not account for the possibility that stars may be ejected from the cluster due to asymmetric mass loss or other effects, nor do we account for binaries among the perturbers.

  2. A random location on the encounter sphere is chosen. The velocities , , and are then sampled from the distribution implied in equation (3.2.1). From these velocities, the periapsis distance and the speed at infinity, , are computed. From the velocities and the orientation of the hyperbolic orbit is determined.

  3. The next encounter is generated assuming that the probability for that the time delay between encounters exceeding is , with given by equation (3.2.1).

3.2.3 Planetary perturbations and the encounter sphere radius

The gravitational effect of each encounter on the planetary orbit is followed from the time the perturber enters the encounter sphere until the time that it again impinges on the encounter sphere on the opposite side of the hyperbolic orbit. The effects of the perturber on the planet when the perturber is outside the encounter sphere are neglected. Note that we do take into account the attraction of the host star on the perturber at all distances by including the effects of gravitational focusing as described in Section 3.2.1.

In principle, the radius of the encounter sphere could be taken to be close to such that the gravitational effect of each encounter on the planetary orbit is fully accounted for. However, this approach is computationally impractical. Fortunately, gravitational forces from perturbers at large distances contribute negligibly to orbital changes because tidal forces fall off as (see equation 2). Therefore, it is sufficient to restrict to relatively small values. In the simulations below, we vary between 25 and 100 and show that these encounter radii are large enough for our purposes (e.g., the outcome fractions are largely independent of ).

We assume that there is at most one perturber within the encounter sphere at a given time. This is justified because in our simulations, the typical timescale for a perturber to pass through the encounter sphere is short compared to the timescale for the next perturber to enter the encounter sphere. More quantitatively, let the encounter passage time be estimated by333This analysis assumes that the trajectory of the perturber relative to the binary is not strongly perturbed by the encounter, that is where is the impact parameter. This assumption is not correct for the closest encounters but should be adequate for this argument. , and the time to the next encounter by , with estimated from equation (3.2.1). Then,

(15)

where we substituted numerical values corresponding to the largest ratio in the simulations. This shows that our assumption of at most one perturber in the encounter sphere is justified for our simulations. We also assume that the orbital phase of the planet is randomized each time a new perturber is introduced in the encounter sphere, which is justified by a similar argument:

(16)

where is the orbital period of the planet.

3.2.4 Isolated tidal evolution

We assume that the orbit of the planet evolves due to tides only when there are no perturbing stars within the encounter sphere. In this process, the semimajor axis and eccentricity of the planet are evolved according to the orbit-averaged version of equation (3),

(17a)
(17b)

where

(18)

Note that equations (17a) and (17b) conserve the semilatus rectum, , and consequently, the orbital angular momentum. Also, note that equations (17a) and (17b) can be obtained from equations (9) and (10), respectively of Hut (1981), by replacing in the latter equations by its pseudosynchronous value, equation (4).

Symbol Description (Range of) Value(s)
Planetary System
Mass of planet-hosting star
Planetary mass
Planetary radius
Planetary tidal time lag
Planetary apsidal motion constant 0.25
Initial stellar obliquity (stellar spin-planetary orbit angle)
Initial planetary orbital semimajor axis 1, 2, 4
Initial planetary orbital eccentricity 0.01–0.6
Cluster and Encounter Properties
Perturber mass 0.1-100
Stellar number density
Stellar (not relative) velocity dispersion
Cluster age
Cluster metallicity 0.001
Encounter sphere radius 25, 50, 75, 100
Derived Encounter Quantities
Stellar mass density
Stellar luminosity density
Estimate of the smallest periapsis distance (0.0068 – 5.1)
Estimate of the total number of encounters 10 –
  • Rayleigh distribution with rms 0.33; subsequently cut off at 0.6.

  • Salpeter initial mass function (Salpeter, 1955), corrected for a stellar age of 5 Gyr (metallicity ), and for gravitational focusing (see Section 3.2).

  • 20 values with logarithmic spacing.

  • Computed from the final mass function assumed for the perturbers.

  • Computed from the final mass and luminosity functions assumed for the perturbers.

  • Minimum and maximum values for the assumed parameters; computed using equation (3.4).

  • Minimum and maximum values for the assumed parameters; computed according to with given by equation (3.2.1).

Table 1Description of symbols used and assumed values in the simulations. Top part: properties of the planet and planet-hosting star. Middle and bottom parts: cluster and encounter properties, and some derived quantities.

3.3. Cluster parameters

For various combinations of and (considered as fixed grid parameters), we simulate systems for a time , with different initial planetary eccentricities and random seeds for the encounters. The adopted values of are 25, 50, 75, and 100 . An overview of the initial conditions is given in Table 1.

Figure 2.— Distribution of central luminosity densities for the MW GCs in the Harris catalog (Harris 1996, 2010 Edition). The assumed range of the luminosity densities in our simulations (Table 1) is indicated by the two vertical dashed lines.

In particular, we choose the one-dimensional velocity dispersion to be , comparable to the typical velocity dispersion of Milky Way (MW) GCs (the mean velocity dispersion in the Harris catalog is ; Harris 1996, 2010 Edition). The number densities in our simulations range between and . Using the SSE stellar evolution code (Hurley et al., 2000) to compute the final masses and luminosities, assuming an age of 5 Gyr and a metallicity of , we find that these number densities correspond to mass densities between and , and luminosity densities between and .

In Fig. 2, we show the distribution of the central luminosity densities, , in the MW GCs (Harris 1996; 2010 Edition). The assumed range in our simulations is indicated with the two vertical dashed lines. The MW GC central luminosity densities peak around , which is near the logarithmic midpoint of our range of simulated densities. We exclude number densities lower than even though these densities occur in some MW GCs, because number densities lower than yield few to no migrating planets in our simulations (see Section 3.6.2 below).

3.4. Encounter properties

Figure 3.— Encounter properties. Solid colored lines: the cumulative distributions of the periapsis distances (top panel), the perturber masses (middle panel), and the speeds at infinity (bottom panel), according to the method described in Section 3.2. Each color corresponds to a different encounter radius , indicated in the legends. In the top panel, the vertical solid black line shows the minimum estimated (eq. 3.4). In the top and bottom panels, the colored dashed lines show the expected analytic distributions, equations (20) and (3.4), respectively. Also in the top and bottom panels, the upper and lower sets of colored horizontal dotted lines show (the Poisson noise in the simulations) and for each value of , where is the number of sampled encounters.

Before presenting our main results in Section 3.6, we briefly discuss properties of the encounters generated via the procedure outlined in Section 3.2 and compare to analytic results. In Fig. 3, we show the cumulative distributions of the periapsis distances (top panel), the perturber masses (middle panel), and the speeds at infinity (bottom panel), obtained numerically through the method described in Section 3.2. In each panel, distributions are shown for the four different encounter sphere radii adopted in the simulations (see Table 1) and in all panels .

The top panel shows the distributions of . The vertical black solid line shows the expected minimum to occur, which can be estimated by setting the collision time for encounters with equal to the age of the cluster, . Following a similar derivation as in Binney & Tremaine (2008, S7.5.8) we find

(19)

where is the mean perturber mass (in the simulations, ). Note that in the absence of gravitational focusing, this reduces to , which is also easily obtained from equation (3.2.1) by setting with . The simulated encounters satisfy . Note that is determined by a single encounter, and therefore it fluctuates among the realizations with different . The upper colored horizontal dotted lines show an estimate of the noise level in the CDF, where is the number of sampled encounters; the lower colored horizontal dotted lines show .

The colored dashed lines in the top panel of Fig. 3 show the expected cumulative distribution,

(20)

(This can be derived from eq. 3.4 by solving for and noting that the CDF out to some periapsis distance must be proportional to for .) Equation (20) is normalized by the condition that the CDF is unity at , since this is the largest possible sampled . The simulated encounters (solid colored lines) are consistent with equation (20), although there is noticeable noise for because of the small number of close encounters.

The middle panel of Fig. 3 shows the distribution of the perturber masses. The median perturber mass is , and the mean is . Despite gravitational focusing, there is little observable dependence of the mass distributions on in these plots, largely because the effects of focusing are strongest for the small fraction of stars with the largest masses.

Finally, the bottom panel of Fig. 3 shows the distribution of the perturber speeds at infinity. The colored dashed lines show the expected cumulative distribution,

(21)

The sampled ’s are consistent with equation (3.4), although there is noticeable noise for . The median sampled depends weakly on , varying from for to for .

3.5. Stopping conditions

In the simulations, we distinguish among the following six outcomes (see Section 3.6.1 below for examples):

  1. HJ formation: the planetary orbit shrinks to an orbital period that lies between 1 and 10 d, due to encounters, tidal dissipation, or a combination of the two. If the orbit is circularized (), we immediately stop the integration — there could be further perturbations due to tidal forces from encounters, but these are extremely weak because the semimajor axis of the orbit is so small. As stated earlier, we neglect stellar tides, so once the orbit is circularized, tidal evolution ceases. In some cases, after the orbital period lies between 1 and 10 d but the orbit is not yet fully circularized. In this case, we also consider the planet to have become an HJ, unless (see below).

  2. Warm Jupiter (WJ) formation: the planetary orbit has shrunk after to an orbital period between 10 and 100 d. The planetary orbit may still be (moderately) eccentric.

  3. Tidal disruption: the planet-host star separation has reached an instantaneous value , where the tidal disruption radius is assumed to be

    (22)

    Here, is a dimensionless parameter, which we assume to be (Guillochon et al., 2011). For a solar-mass host star and a planet of Jupiter’s mass and radius, which are assumed here, . We do not consider collisions with the perturber as a separate outcome, although a fraction of the “tidal disruption” outcomes will in fact be collisions with the host star.

  4. Unbinding: after an encounter, the semimajor axis of the planet with respect to its host star is negative or larger than . In the latter case, the planet is technically still bound to its host star. However, the wide orbit makes the planet extremely susceptible to future encounters even for the lowest densities that we consider, very likely leading to ejection.

  5. Transfer: the planet is captured by a perturber, which we consider to be the case if the semimajor axis of the planet with respect to the perturber after the encounter is positive and less than . In principle, a new simulation could be started with the captured planet around the last perturber, subject to further perturbations from encounters that could, e.g., lead to HJ formation. However, this further investigation is beyond the ambitions of this paper.

  6. No migration: none of the above occurred. Nevertheless, the semimajor axis and/or eccentricity may still have been affected significantly.

3.6. Simulation results

3.6.1 Examples

Figure 4.— Six examples to illustrate the typical eccentricity and semimajor axis evolution for the outcomes defined in Section 3.5. In these examples, the initial planetary semimajor axis is , the encounter radius is , and the stellar density is . The vertical axes show the semimajor axis (black dotted lines), the periapsis distance (black solid lines), and the encounter periapsis distances (red solid lines). Top row from left to right: HJ formation, WJ formation, tidal disruption; bottom row: unbound planet, transferred planet, and no migration. In each panel, the yellow horizontal dashed line shows the tidal disruption radius (eq. 22). Note that each change in the height of the red solid lines corresponds to a new encounter with a certain ; the horizontal sections indicate the time between encounters, not the passage time of the encounters themselves. The ratio of the passage time to the time between encounters for the chosen parameters is (see equation 3.2.3).

In Fig. 4, we show a number of examples to illustrate the typical eccentricity and semimajor axis evolution for the outcomes defined in Section 3.5. In these examples, the initial planetary semimajor axis is , the encounter radius is , and the stellar density is .

In the top-left panel, an HJ is formed within of evolution. After , a number of strong perturbers with a of a few gradually decrease the periapsis distance, while leaving the semimajor axis relatively unchanged. After a strong encounter at , the periapsis distance is reduced to a small enough value that strong tidal evolution is triggered. The orbit circularizes at a semimajor axis of . In the top-middle panel, a WJ is formed through a series of encounters in which the semimajor axis decreases to and the periapsis distance decreases to . The planet survives for 10 Gyr, and has final orbital properties characteristic of WJs (; ). Tidal dissipation played no role in this process. The top-right panel shows an example of tidal disruption; the semimajor axis is not much affected by encounters, whereas the periapsis distance is gradually decreased. A single encounter at increases the semimajor axis to , making the planet more susceptible to further perturbations, and it is disrupted shortly thereafter.

In the bottom-left panel, a strong encounter with unbinds the planet at . Another strong encounter occurs at in the bottom-middle panel, and the planet is transferred to the perturber. Lastly, in the bottom-right panel, no destructive or migration-inducing encounters occur, and the planet survives without significant tidal dissipation, although the eccentricity has been slightly excited and the semimajor axis slightly increased.

3.6.2 Simulation outcome fractions

Figure 5.— Fractions of the six outcomes (Section 3.5) as a function of the stellar number density . The initial planetary semimajor axis is () in the left (right) set of panels. Red solid lines correspond to HJ formation, blue solid lines correspond to WJ formation, yellow dashed lines correspond to tidal disruption events, dark blue dashed lines correspond to unbound planets, green dotted-dashed lines correspond to transferred planets, and solid black lines (‘no migration’) correspond to planets that remain bound to their host star with periods above the assumed WJ threshold, 100 d (see Section 3.5). Results are shown for four values of the encounter sphere radius ; line thickness increases with increasing . The bottom panels are zoomed-in versions of the top panels, showing in more detail the less common outcomes. The data (for ) are also given in Table 2.
3.0 25 0.00 0.00 0.00 0.05 0.00 0.94
3.0 50 0.00 0.00 0.00 0.06 0.00 0.94
3.0 75 0.00 0.00 0.00 0.06 0.00 0.93
3.0 100 0.00 0.00 0.00 0.07 0.01 0.92
3.16 25 0.00 0.00 0.00 0.07 0.01 0.92
3.16 50 0.00 0.00 0.00 0.07 0.01 0.92
3.16 75 0.00 0.00 0.00 0.09 0.00 0.90
3.16 100 0.00 0.00 0.00 0.10 0.01 0.89
3.32 25 0.00 0.00 0.00 0.09 0.01 0.89
3.32 50 0.00 0.00 0.00 0.11 0.01 0.87
3.32 75 0.00 0.00 0.00 0.13 0.01 0.86
3.32 100 0.00 0.00 0.00 0.15 0.01 0.84
3.47 25 0.00 0.00 0.00 0.13 0.01 0.85
3.47 50 0.00 0.00 0.01 0.15 0.02 0.82
3.47 75 0.00 0.00 0.00 0.17 0.01 0.80
3.47 100 0.00 0.00 0.01 0.20 0.01 0.78
3.63 25 0.00 0.00 0.01 0.19 0.02 0.77
3.63 50 0.01 0.00 0.01 0.21 0.02 0.74
3.63 75 0.00 0.00 0.01 0.24 0.02 0.73
3.63 100 0.01 0.00 0.01 0.26 0.01 0.71
3.79 25 0.01 0.00 0.01 0.27 0.03 0.68
3.79 50 0.01 0.01 0.01 0.30 0.03 0.65
3.79 75 0.01 0.00 0.01 0.33 0.02 0.63
3.79 100 0.01 0.00 0.01 0.37 0.02 0.59
3.95 25 0.01 0.01 0.01 0.38 0.04 0.55
3.95 50 0.01 0.00 0.02 0.41 0.04 0.52
3.95 75 0.01 0.00 0.01 0.44 0.03 0.50
3.95 100 0.01 0.00 0.01 0.49 0.03 0.47
4.11 25 0.02 0.01 0.03 0.50 0.05 0.40
4.11 50 0.02 0.01 0.02 0.57 0.05 0.34
4.11 75 0.02 0.00 0.02 0.58 0.03 0.34
4.11 100 0.01 0.00 0.02 0.63 0.03 0.31
4.26 25 0.03 0.01 0.03 0.63 0.06 0.24
4.26 50 0.02 0.01 0.03 0.68 0.06 0.20
4.26 75 0.03 0.01 0.03 0.72 0.04 0.18
4.26 100 0.02 0.01 0.03 0.76 0.03 0.16
4.42 25 0.04 0.01 0.03 0.72 0.08 0.12
4.42 50 0.03 0.01 0.04 0.78 0.06 0.08
4.42 75 0.02 0.00 0.03 0.82 0.04 0.08
4.42 100 0.02 0.00 0.03 0.85 0.03 0.07
4.58 25 0.03 0.01 0.04 0.80 0.08 0.04
4.58 50 0.03 0.01 0.05 0.82 0.07 0.03
4.58 75 0.03 0.00 0.04 0.86 0.05 0.02
4.58 100 0.02 0.00 0.03 0.88 0.04 0.03
4.74 25 0.03 0.00 0.04 0.83 0.08 0.01
4.74 50 0.03 0.00 0.05 0.85 0.06 0.01
4.74 75 0.03 0.00 0.04 0.88 0.05 0.01
4.74 100 0.02 0.00 0.03 0.91 0.04 0.00
4.89 25 0.03 0.00 0.04 0.84 0.08 0.00
4.89 50 0.03 0.00 0.05 0.85 0.06 0.00
4.89 75 0.03 0.00 0.05 0.87 0.05 0.00
4.89 100 0.02 0.00 0.04 0.90 0.04 0.00
5.05 25 0.03 0.00 0.04 0.85 0.08 0.00
5.05 50 0.02 0.00 0.05 0.87 0.06 0.00
5.05 75 0.02 0.00 0.04 0.89 0.05 0.00
5.05 100 0.02 0.00 0.04 0.91 0.04 0.00
5.21 25 0.02 0.00 0.04 0.85 0.09 0.00
5.21 50 0.02 0.00 0.05 0.86 0.07 0.00
5.21 75 0.01 0.00 0.05 0.89 0.05 0.00
5.21 100 0.01 0.00 0.05 0.91 0.04 0.00
5.37 25 0.02 0.00 0.06 0.83 0.09 0.00
5.37 50 0.02 0.00 0.06 0.86 0.07 0.00
5.37 75 0.01 0.00 0.04 0.90 0.04 0.00
5.37 100 0.01 0.00 0.04 0.91 0.04 0.00
5.53 25 0.01 0.00 0.05 0.85 0.08 0.00
5.53 50 0.02 0.00 0.05 0.87 0.07 0.00
5.53 75 0.01 0.00 0.05 0.90 0.04 0.00
5.53 100 0.01 0.00 0.04 0.90 0.04 0.00
5.68 25 0.02 0.00 0.06 0.85 0.07 0.00
5.68 50 0.01 0.00 0.05 0.87 0.07 0.00
5.68 75 0.01 0.00 0.05 0.89 0.05 0.00
5.68 100 0.01 0.00 0.04 0.91 0.04 0.00
5.84 25 0.01 0.00 0.06 0.85 0.09 0.00
5.84 50 0.01 0.00 0.06 0.86 0.07 0.00
5.84 75 0.01 0.00 0.05 0.90 0.05 0.00
5.84 100 0.00 0.00 0.05 0.90 0.04 0.00
6.0 25 0.01 0.00 0.06 0.85 0.08 0.00
6.0 50 0.01 0.00 0.06 0.87 0.07 0.00
6.0 75 0.01 0.00 0.05 0.91 0.04 0.00
6.0 100 0.01 0.00 0.05 0.92 0.03 0.00
Table 2 Simulated outcome fractions for the six outcomes listed in Section 3.5 (: unbound; : transferred; : no migration). The same data are also shown in Fig. 5. We define . The statistical uncertainty in the outcome fractions is about 0.02.

Fig. 5 shows the fractions of the six outcomes (see Section 3.5) as a function of the stellar number density . The initial planetary semimajor axis is () in the left (right) set of panels. In Fig. 5 and in subsequent figures, red solid lines correspond to HJ formation, blue solid lines correspond to WJ formation, yellow dashed lines correspond to tidal disruption events, dark blue dashed lines correspond to unbound planets, green dotted-dashed lines correspond to transferred planets, and solid black lines correspond to nonmigrating planets. The bottom panels are zoomed-in versions of the top panels, showing in more detail the less common outcomes. The expected Poisson fluctuation in the outcome fractions from our realizations is . The data (for ) are also given in Table 2.

We first discuss the case . For the lowest densities, the encounter rate is too low to strongly affect the planetary orbit, and the fraction of nonmigrating planets is near 100%. As is increased, the nonmigrating fraction decreases smoothly to zero near , while the fraction of unbound planets increases to at the same density and stays near that value for even higher densities. The migration fraction, i.e., the sum of , and , is nearly zero for , increases to near , and remains approximately constant at higher densities. Note that the sum of the migration fraction and the nonmigration fraction is not unity, because both categories exclude the unbinding and transfer outcomes.

The fraction of HJs peaks at at (the fraction approaches for the smallest encounter radii but these values are less reliable). For higher densities, the HJ fraction decreases, whereas the fraction of tidal disruptions and transferred planets increases. These trends can be understood qualitatively by noting that the density needs to be sufficiently high for strong encounters to decrease the periapsis distance of the planet and trigger tidal migration. However, as the density increases, strong encounters become increasingly likely, and therefore more catastrophic outcomes, like tidal disruption or a transfer to the perturber, start to dominate.

The fraction of unbound planets is weakly dependent on . Nevertheless, the fractions shown in Fig. 5 indicate convergence with respect to ; in particular, there are only small differences between and 100 .

Comparing the two cases and 4 , we note that the curves showing the fraction of nonmigrating and unbound planets are shifted to the left (lower ) for , which can be understood qualitatively because planets with larger semimajor axes are more weakly bound. Similarly to the case , the migration fraction is approximately constant at high densities (), although the absolute value of the migration fraction, , is lower by a factor of . Furthermore, the HJ fraction peaks at a lower density, in this case around , with a value of which is a factor of lower than the case .

The fractions for the case (not presented here) show intermediate trends between those depicted in Fig. 5. We conclude that a larger tends to lower the HJ and WJ fractions and, more generally, the fractions of ‘interesting’ outcomes such as tidal disruption and planet transfer. In the remainder of Section 3.6, we focus mainly on the case . Most of the results shown in these sections do not depend strongly on .

3.6.3 Orbital changes

Figure 6.— Orbital elements in the plane at the end of the simulation, or after a stopping condition occurred, for the HJ (red stars), WJ (blue triangles), tidal disruption (yellow crosses), and no migration cases (black dots). Each panel corresponds to a different density, indicated therein (), and the initial semimajor axis is . The yellow dashed line shows the assumed threshold for tidal disruption, , with given by equation (22). Transfer outcomes are not shown; see Fig. 11 for a similar figure showing the orbital elements of the transferred planets.

To illustrate how encounters affect the planetary orbit, we show in Fig. 6 the orbital elements in the plane at the end of the simulation, or after a stopping condition occurred, for the outcomes HJ (red stars), WJ (blue triangles), tidal disruption (yellow crosses), and no migration (black dots). Each panel corresponds to a different density. In all panels . We recall that the initial eccentricities were sampled from a Rayleigh distribution with an rms of 0.33 and maximum value of 0.6 (see Section 3.1).

The yellow dashed line shows the assumed threshold for tidal disruption, , with given by equation (22). Obviously, the points corresponding to the tidal disruption outcome lie below this line. Most tidal disruptions occur with relatively large semimajor axes, , implying that the encounters leading to tidal disruptions mainly did so by exciting high eccentricities as opposed to strongly decreasing the semimajor axes.

The HJs populate the lower-right part of the plane. They tend to be either completely circular or have a small eccentricity. Their semimajor/orbital period distribution is discussed in more detail below. The WJs lie in between the HJs and nonmigrating planets, and are substantially more eccentric than the HJs (WJs around field stars in the solar neighborhood also tend to be eccentric, but probably for different reasons). The nonmigrating planets populate a large region in the diagram, even though they all began with a common semimajor axis of and small eccentricities. Encounters change the semimajor axes, typically by a factor of a few (note that we consider planets to be unbound if ), and can excite the eccentricities to values as high as (or even higher, at which point the planets become HJs or are tidally disrupted).

Figure 7.— Cumulative distributions of the semimajor axes separately for the three ‘nondisruptive’ outcomes, i.e., HJ and WJ formation and no migration, for four values of and various values of . The top (bottom) set of panels corresponds to (). Increasing line widths correspond to larger . Note that not all outcomes occur at all densities. In all panels, the black vertical dashed lines show predictions from the analytic model discussed in Section 3.7.
Figure 8.— Colored lines: the orbital period distributions for the nondisruptive outcomes in the simulations (HJ and WJ formation and no migration, combining data from all values of , and ). Black crosses and error bars: the observed distribution for giant planets in the solar neighborhood (Santerne et al., 2016). Black open circles: the distribution from simulations by Anderson et al. (2016) of high- HJ formation in stellar binaries in the field (data are shown from Fig. 23 of Anderson et al. 2016 for and , where and is the planetary tidal time lag; with , , and , or corresponds to a viscous timescale of ). The distributions are normalized to unit total area. The short black vertical dashed lines show the predictions from the analytic model (Section 3.7).

In Fig. 7, we show the cumulative distributions of the semimajor axes separately for the “nondisruptive” cases, i.e., HJ and WJ formation and no migration, for four values of and the various values of . The top (bottom) set of panels corresponds to (). We repeat that we did not include stellar tides in the simulations, which could affect the HJ semimajor axis distribution if the stellar tides are efficient (in particular, if the star has a convective layer). There is little dependence of the semimajor axis distributions on , which provides reassurance that our choices of are large enough. The HJ semimajor axes tend to be slightly smaller for higher densities and larger . These trends can be explained with an analytic model that is described below in Section 3.7. The predictions from that model are shown in Fig. 7 with the vertical black dashed lines.

In Fig. 8, we show the distributions of the final orbital periods for . The data have been combined for all three nondisruptive cases and all four values of . The different colored lines correspond to two values of . For comparison, we also show with black crosses and error bars the observed distribution for giant planets in the solar neighborhood by Santerne et al. (2016), and, with black open circles, the distribution from simulations by Anderson et al. (2016) of HJ formation by high-eccentricity migration in stellar binaries.

Our simulations predict a peak in the HJ period distribution in dense clusters at 2.5 day, which is similar to the predicted peak for stellar binaries in the field, but shorter than the observed peak for field HJs at 5 day. The shorter orbital periods in our simulations and those of Anderson et al. (2016) compared to the observations are typical for high- migration scenarios. Note that the final orbital period depends sensitively on the planetary radius (the final semimajor axis is approximately proportional to , see Section 3.7 below), and we assumed a single planetary radius (). Therefore, our simulations can be made to be more consistent with the observations by, e.g., assuming inflated HJs (e.g., Laughlin et al. 2011).

The occurrence of WJs in our simulations is a factor of a few lower than HJs (see Fig. 5). In contrast, the observations of Santerne et al. (2016) for the solar neighborhood indicate that a large fraction of planets have orbital periods characteristic of WJs. Simulations of high- migration in stellar binaries (e.g,. Petrovich & Tremaine 2016; Antonini et al. 2016), multiplanet systems (Hamers et al., 2017), and stellar triples (Hamers, 2017a) also produce few WJs, or a number of WJs inconsistent with observations. Nevertheless, the WJ fractions in our simulations are higher than, e.g., in Hamers et al. (2017) and (Hamers, 2017a). This is because the WJs in our simulations tend to be produced directly by encounters without substantial tidal evolution, whereas in purely secular high- migration, the semimajor axis can only change due to tidal interactions.

3.6.4 Obliquities

Figure 9.— Obliquity distributions for the HJ, WJ, tidal disruption, and no migration outcomes in the simulations and for four densities, combining data from the different . The ‘no migration’ and WJ curves are not plotted at the highest density, because these outcomes were too rare to produce good statistics. The top (bottom) panels show the probability (cumulative) distributions.

The stellar obliquity , the angle between the stellar spin and planetary orbit, is known to be broadly distributed for field HJs, with almost half of HJs having obliquities that exceed 10–20 and a handful having obliquities as high as (see Winn & Fabrycky 2015 for a review). In Fig. 9, we show the obliquity distributions for the HJ, WJ, tidal disruption, and no migration outcomes in the simulations with , combining data from the different . The top (bottom) panels show the probability density (cumulative) distributions. In the simulations, the spin of the star was not modeled; to compute the obliquity, we assume that the stellar spin and planetary orbit were initially aligned, and that the stellar spin direction is fixed.

For all of the four outcomes, the final obliquity distribution is quite broad and not far from isotropic (i.e., a flat distribution in ). For nonmigrating planets, the final distributions are not exactly flat, i.e., some memory of the initial zero obliquity is retained. Stronger interactions are involved in the other outcomes (in particular, HJ and WJ formation); consequently, the obliquity distributions for these cases tend to be more isotropic. The fraction of retrograde HJs ranges from to , higher than the few per cent observed in the field (Winn & Fabrycky, 2015). A high fraction of retrograde obliquities would be a telltale sign for an encounter-induced high- migration origin of HJs in dense clusters.

3.6.5 Stopping times

Figure 10.— Cumulative distributions of the stopping times, i.e., the times required to form an HJ, to tidally disrupt the planet, to unbind the planet, or to transfer the planet, combining different data for the different , and for .

The stopping times in the simulations are the times needed for encounters to form an HJ, to tidally disrupt the planet, to unbind the planet, or to transfer the planet. In Fig. 10, we show the cumulative distributions of these times for different densities, combining data for the different (with ). As expected, the stopping times are strongly dependent on the density: a higher density implies that the aforementioned stopping conditions occur earlier. Generally, the unbound and transfer outcomes occur earlier than the tidal disruption outcomes, which in turn occur earlier than the HJ outcomes. This can be understood by noting that unbinding or transfer generally occurs through a single strong encounter, whereas tidal disruption may involve a larger number of weaker encounters. Finally, the formation of HJs tends to occur later than tidal disruptions, given that some tidal evolution is typically involved to produce the HJ, adding to the formation time.

The stopping times are generally shorter for larger (results not shown here), which can be understood qualitatively from the stronger effect of perturbers if the planet is less bound to its host star.

3.6.6 Transferred planets

Figure 11.— Semimajor axis and eccentricity of the planet with respect to the perturber to which it was transferred, combining results from all densities and encounter sphere radii and assuming . To reduce clutter, only a fifth of the available data points are shown. The top and right panels show the marginalized eccentricity and semimajor axis distributions, respectively. The yellow line shows the assumed tidal disruption limit, and the red dashed line shows , approximately the periapsis distance for which tidal effects are important.

The planet is transferred to a perturbing star in up to of our simulations with (see Fig. 5). In Fig. 11, we show the semimajor axis and eccentricity of the planet with respect to the perturber to which it was transferred, combining results from all densities and encounter sphere radii, and assuming . The top and right panels show the marginalized eccentricity and semimajor axis distributions, respectively. The yellow line shows the assumed tidal disruption limit, and the red dashed line shows , approximately the periapsis distance for which tidal effects are important. The captured planets have semimajor axes ranging from roughly 0.2–5 (the upper limit is an artifact of the definition described in Section 3.5). The bulk of the transferred planets are not captured onto orbits for which tidal evolution is immediately important (of course, subsequent perturbations by encounters with other stars could drive tidal migration at later stages).

3.6.7 Perturber properties

Figure 12.— Cumulative distributions of for the last perturber associated with each outcome. Data are combined for the different , and . Note that there are no WJ or ‘no migration’ outcomes for the highest density shown.
Figure 13.— Similar to Fig. 12, now showing the cumulative distributions of the perturber masses, .
Figure 14.— Similar to Fig. 12, now showing the cumulative distributions of the perturber speeds at infinity, .

In Figs. 12, 13 and 14, we illustrate how the simulation outcomes are determined by the perturber properties. In particular, in these three figures we show the cumulative distributions of , and , respectively, for the last perturber associated with each outcome. For example, in the case of the unbound outcome, , , and correspond to the perturber that triggered the ejection of the planet. As may be expected intuitively, the ‘strong’ outcomes (tidal disruption, unbinding or transfer) tend to be associated with lower values of , higher values of , and lower values of . Note that for the HJs, the last encounter is usually not the encounter that drove tidal migration: typically, one or several strong encounters drive high eccentricity and tidal migration, and subsequently the planet is decoupled from further perturbations. The figures show the properties of the last perturber prior to the HJ stopping condition. The figures shown here apply only to the case , but the perturber properties are not much changed in the simulations with larger .

3.7. Analytic model for encounter-induced high- migration

Figure 15.— Colored solid lines: tracks in the diagram of several individual systems from the simulations that lead to HJs. For all tracks, , the density is , and is either 25, 50, 75 or 100 . The black dashed lines show analytic results from Section 3.7 (eq. 32; setting , approximately the mean perturber mass in the simulations), for three different values of the initial semimajor axis: 1, 2.5 and 5 (line thickness increases with increasing initial semimajor axis). In the lower-right part of the figure, the short black horizontal dotted lines show the simplest estimate of the stalling semimajor axis, equation (35); the short black horizontal dashed lines show the more accurate estimate, equation (3.7). The yellow dashed line is the tidal disruption line, with given by equation (22).

In this section, we discuss a simplified analytic model for the evolution of the planetary orbit including perturbations by encounters and tidal evolution (assuming equilibrium tides, as usual). The main aim is to obtain an estimate of the final ‘stalling’ semimajor axis of the HJs, which was observed in the simulations to be (Section 3.6.3 and Fig. 8).

First, we show in Fig. 15 with colored solid lines tracks in the plane of several individual systems from the simulations with . All systems start at some location on the right part of the horizontal line at (the initial eccentricity distribution was assumed to be a Rayleigh distribution; see Table 1). Due to encounters, the semimajor axis and eccentricity change randomly, and the planet random-walks in the plane. For the examples shown, eventually the periapsis distance becomes sufficiently small for tidal evolution to become important, after which the system evolves with approximately constant angular momentum or semilatus rectum until the orbit is circularized, at .

To model this type of evolution, we consider the coupled differential equations for and , where is the normalized angular momentum. We approximate the orbit-averaged tidal evolution equations (17) by replacing by (since we are mainly interested in the competition between tidal circularization and eccentricity excitation by encounters at high eccentricity), giving