Hot Jupiter and Ultra-cold Saturn Formation in Dense Star Clusters

2.1.1. Planetary System Setups

2.1.2. Critical Velocity for Multi-hierarchy Scattering

The critical velocity v_c is an important variable in multi-body scattering. It is defined such that, for a relative velocity at infinity between two objects v_∞ > v_c (objects can be isolated bodies or they can have complex internal orbital architectures as is the case for binaries, triples, etc.), the total energy of the system in the center of mass reference frame is positive (e.g., Hut & Bahcall 1983; Fregeau et al. 2004; Leigh & Geller 2012; Leigh et al. 2016a, 2016b). Thus, in the regime where the total interaction energy is positive, full ionization is a possible outcome scenario for objects containing internal orbital architectures. In this regime, the incident object usually does a flyby relative to the scattered object, and interacts with it only once. For v_∞ < v_c, the total energy of the scattering system is negative. The corresponding scatterings are likely to be resonant in this regime, such that a bound chaotic interaction ensues, with the incident object undergoing repeated flybys relative to the scattered object. These interactions typically progress until disruption, forming at least one tightly bound system (other orbits need not survive to satisfy the criteria for energy and angular momentum conservation).

We label the masses of the k components of the scattered system to be m₀₀, m₀₁, ⋯ m_0k and the masses of the l components of the incident system to be m₁₀, m₁₁, ⋯ m_1l . If the relative velocity between the centers of mass of the two objects is v_∞, the total energy of the scattering system in the center of mass reference frame is

$$\begin{eqnarray}&&\begin{array}{l}{E}_{\mathrm{tot}}=\displaystyle \frac{1}{2}\displaystyle \frac{{M}_{0}{M}_{1}}{{M}_{0}+{M}_{1}}{v}_{\infty }^{2}+\displaystyle \sum _{i}^{n}\displaystyle \frac{{m}_{0i}{\left({{\boldsymbol{v}}}_{0i}-{{\boldsymbol{V}}}_{0}\right)}^{2}}{2}\\ -\displaystyle \sum _{i\lt j\leqslant k}\displaystyle \frac{{{Gm}}_{0i}{m}_{0j}}{| {{\boldsymbol{r}}}_{0i}-{{\boldsymbol{r}}}_{0j}| }+\displaystyle \sum _{i}^{m}\displaystyle \frac{{m}_{1i}{({{\boldsymbol{v}}}_{1i}-{{\boldsymbol{V}}}_{1})}^{2}}{2}\\ -\displaystyle \sum _{i\lt j\leqslant l}\displaystyle \frac{{{Gm}}_{1i}{m}_{1j}}{| {{\boldsymbol{r}}}_{1i}-{{\boldsymbol{r}}}_{1j}| },\end{array}\end{eqnarray} \tag{ 1 }$$

where M₀ = ∑_i m_0i is the total mass of the scattered object and M₁ = ∑_i m_1i is the total mass of the incident object. V₀ and V₁ are the initial center of mass velocities of the scattered and incident objects, respectively. The critical velocity between two objects corresponding to a total energy of zero for the system is calculated using the scattering facilities of SpaceHub, specifically via the following equation:

$$\begin{eqnarray}&&\begin{array}{l}{v}_{{\rm{c}}}^{2}=\displaystyle \frac{{M}_{0}+{M}_{1}}{{M}_{0}{M}_{1}}\left[\displaystyle \sum _{i\lt j\leqslant k}\displaystyle \frac{2{{Gm}}_{0i}{m}_{0j}}{| {{\boldsymbol{r}}}_{0i}-{{\boldsymbol{r}}}_{0j}| }+\displaystyle \sum _{i\lt j\leqslant l}\displaystyle \frac{2{{Gm}}_{1i}{m}_{1j}}{| {{\boldsymbol{r}}}_{1i}-{{\boldsymbol{r}}}_{1j}| }\right.\\ \left.-\displaystyle \sum _{i}^{k}{m}_{0i}{\left({{\boldsymbol{v}}}_{0i}-{{\boldsymbol{V}}}_{0}\right)}^{2}-\displaystyle \sum _{i}^{l}{m}_{1i}{\left({{\boldsymbol{v}}}_{1i}-{{\boldsymbol{V}}}_{1}\right)}^{2}\right].\end{array}\end{eqnarray} \tag{ 2 }$$

2.1.3. Impact Parameter for High Mass Ratio Scattering

Due to gravitational focusing, incident objects at spatial infinity with impact parameter b will pass the scattered object at closest approach with a corresponding critical pericenter distance for the hyperbolic/parabolic orbit:

$$\begin{eqnarray}&&\begin{array}{l}{Q}_{{\rm{c}}}={a}_{\mathrm{0,1}}(1-{e}_{\mathrm{0,1}})\\ =\ {a}_{\mathrm{0,1}}\left(1-\sqrt{1+\displaystyle \frac{{b}^{2}}{{a}_{0,1}^{2}}}\right)\\ =\ -\displaystyle \frac{G({M}_{0}+{M}_{1})}{{v}_{\infty }^{2}}\left(1-\sqrt{1+\displaystyle \frac{{b}^{2}}{{\left(-\tfrac{G({M}_{0}+{M}_{1})}{{v}_{\infty }^{2}}\right)}^{2}}}\right).\end{array}\end{eqnarray} \tag{ 3 }$$

We only consider interactions satisfying Q < Q_c, where p_c is the critical pericenter distance. Therefore, to cover all parameter space in which the desired outcome is causally and physically permitted, the scattering experiments must be sampled with an impact parameter distribution of b², uniformly distributed in the range $[0,{b}_{\max }^{2}]$, where b_max is the corresponding impact parameter for p_c. We can rewrite Equation (3) in the more convenient form:

$$\begin{eqnarray}&&\begin{array}{l}{b}_{\max }={Q}_{{\rm{c}}}\sqrt{1-2\displaystyle \frac{{a}_{\mathrm{0,1}}}{{Q}_{{\rm{c}}}}}\\ =\ {Q}_{{\rm{c}}}\sqrt{1+2\displaystyle \frac{G({M}_{0}+{M}_{1})}{{v}_{\infty }^{2}{Q}_{{\rm{c}}}}}.\end{array}\end{eqnarray} \tag{ 4 }$$

For higher values of v_∞, ${b}_{\max }\to {Q}_{{\rm{c}}}$, while the lower limit gives ${b}_{\max }\sim \tfrac{1}{{v}_{\infty }}$. In Hut & Bahcall (1983), the two limits are combined to construct the following equation:

$$\begin{eqnarray}&&{b}_{\max }({v}_{\infty })=\left(\displaystyle \frac{C}{{v}_{\infty }/{v}_{{\rm{c}}}}+D\right)R,\end{eqnarray} \tag{ 5 }$$

where R is the contact distance between the two scattering objects (e.g., single-binary scattering is the SMA of the binary; binary–binary scattering is the sum of the semimajor axes of the binaries) and C and D are dimensionless constants that depend on the desired outcome (DR ∼ Q_c in the high v_∞ limit and CR ∼ Q_c in the low v_∞ limit). This equation is widely adopted in nearly equal mass single-binary and binary–binary scatterings. However, in this work, the planet mass is much smaller than the star mass, and hence special consideration should be given to this interesting limit. The velocity term v_∞ in the original equation is made dimensionless by the critical velocity v_c. If we take the limit of high v_∞ in both Equations (4) and (5), we have

$$\begin{eqnarray}&&\displaystyle \frac{{CR}}{{v}_{\infty }/{v}_{{\rm{c}}}}={b}_{\max }={Q}_{{\rm{c}}}\sqrt{\displaystyle \frac{2G({M}_{0}+{M}_{1})}{{Q}_{{\rm{c}}}}}/{v}_{\infty },\end{eqnarray} \tag{ 6 }$$

since CR ∼ Q_c. The equation above holds only if

$$\begin{eqnarray}&&{v}_{{\rm{c}}}\sim \sqrt{\displaystyle \frac{2G({M}_{0}+{M}_{1})}{{Q}_{{\rm{c}}}}}=\sqrt{\displaystyle \frac{{M}_{0}+{M}_{1}}{{M}_{0}{M}_{1}}\displaystyle \frac{2{{GM}}_{0}{M}_{1}}{{CR}}}.\end{eqnarray} \tag{ 7 }$$

Equation (2) then requires

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{2{{GM}}_{0}{M}_{1}}{{CR}}\sim \displaystyle \sum _{i\lt j\leqslant k}\displaystyle \frac{2{{Gm}}_{0i}{m}_{0j}}{| {{\boldsymbol{r}}}_{0i}-{{\boldsymbol{r}}}_{0j}| }+\displaystyle \sum _{i\lt j\leqslant l}\displaystyle \frac{2{{Gm}}_{1i}{m}_{1j}}{| {{\boldsymbol{r}}}_{1i}-{{\boldsymbol{r}}}_{1j}| }\\ -\displaystyle \sum _{i}^{k}{m}_{0i}{\left({{\boldsymbol{v}}}_{0i}-{{\boldsymbol{V}}}_{0}\right)}^{2}-\displaystyle \sum _{i}^{l}{m}_{1i}{\left({{\boldsymbol{v}}}_{1i}-{{\boldsymbol{V}}}_{1}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 8 }$$

For single-binary scatterings with component masses (m₀₀, m₀₁)–m₁, this becomes

$$\begin{eqnarray}&&\displaystyle \frac{G({m}_{00}+{m}_{01}){m}_{1}}{{Ca}}\sim \displaystyle \frac{{{Gm}}_{00}{m}_{01}}{2a}.\end{eqnarray} \tag{ 9 }$$

If m₀₀ ∼ m₀₁ ∼ m₁, this yields C ∼ 4, which is the frequently adopted value. However, if m₀₁ is a Jupiter mass planet, or m₀₀ ∼ m₁ ∼ 1000m₀₁, we obtain C ∼ 2000. Indeed,

$$\begin{eqnarray}&&C\propto \displaystyle \frac{{M}_{0}{M}_{1}}{{\displaystyle \sum }_{i\lt j\leqslant k}{m}_{0i}{m}_{0j}+{\displaystyle \sum }_{i\lt j\leqslant l}{m}_{1i}{m}_{1j}}.\end{eqnarray} \tag{ 10 }$$

Therefore, for high mass ratio systems, we cannot use Equation (5) by normalizing v_∞ by v_c and letting C ∼ [4–10] (Hut & Bahcall 1983; Fregeau et al. 2004). The scattering facilities of SpaceHub provide both choices given in Equations (4) and (5).

2.1.4. Tidal Factor Between Two Objects

The tidal factor is widely used in scattering experiments to estimate the tidal perturbation between two objects. It is defined as

$$\begin{eqnarray}&&{\delta }_{0\to 1}={F}_{\mathrm{tid},0\to 1}/{F}_{\mathrm{rel},1},\end{eqnarray} \tag{ 11 }$$

where F_tid,0→1 is the tidal force that object 0 exerts on object 1 and F_rel,1 is the internal gravitational force of object 1. If δ → 0, then object 0 and object 1 are isolated. This tidal tolerance factor is often used to save on computational runtime, by switching to integrating Keplerian orbits in the limit of very low tidal tolerance factors (i.e., this option is turned on whenever the tidal tolerance parameter drops below a critical user-defined value). In this paper, we only use this parameter to generate the initial position of the interloper but we do not invoke this approximation afterwards, and explicitly integrate the N-body problem without any such simplifying approximations.

In binary–binary and binary-single scatterings, as in Fregeau et al. (2004), we have

$$\begin{eqnarray}&&{F}_{\mathrm{tid},0\to 1}=\displaystyle \frac{2{{GM}}_{0}{M}_{1}}{| {{\boldsymbol{r}}}_{0}-{{\boldsymbol{r}}}_{1}{| }^{3}}{a}_{1}(1+{e}_{1})\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{rel},1}=\displaystyle \frac{{{Gm}}_{10}{m}_{11}}{| {a}_{1}(1+{e}_{1}){| }^{2}}.\end{eqnarray} \tag{ 13 }$$

Indeed, a more general equation for scatterings between any two generic objects is provided by the tidal radius

$$\begin{eqnarray}&&{r}_{\mathrm{tid},0\to 1}({\delta }_{0\to 1})=| {{\boldsymbol{r}}}_{0}-{{\boldsymbol{r}}}_{1}| ={\left(\displaystyle \frac{{M}_{0}}{{\delta }_{0\to 1}{M}_{1}}\right)}^{1/3}{D}_{1},\end{eqnarray} \tag{ 14 }$$

where D₁ is the geometric size of object 1 (i.e., the stellar radius for a single star, the orbital separation for a binary star, etc.). For binary-single scatterings we have D₁ = a₁(1 + e₁). The typical value δ ∼ 0.02 yields the classical tidal disruption equation

$$\begin{eqnarray}&&| {{\boldsymbol{r}}}_{0}-{{\boldsymbol{r}}}_{1}| =3.7{\left(\displaystyle \frac{{M}_{0}}{{M}_{1}}\right)}^{1/3}{D}_{1}.\end{eqnarray} \tag{ 15 }$$

The scattering facilities of SpaceHub take the maximum value of either r_tid,0→1(δ_0→1) or r_tid,1→0(δ_1→0) as the final tidal radius r_tid(δ).

For scattering experiments, strong interactions occur if p < p_c ∼ r_tid(δ = 0.02), while perturbative interactions occur if r_tid(δ = 0.02) < p < r_tid(δ = 10⁻⁴). The tidal radius is also used to calculate the initial relative distance between the scattered and incident objects. In this project, we take r_tid(δ = 10⁻⁵) as the initial distance separating the centers of mass for the two scattering objects. The trajectories from infinity to this initial distance are calculated analytically by regarding the scattering objects as isolated up until this critical point.

2.1.5. Initial Conditions for Orbital Angles

2.1.6. Termination Time of Flyby

The interloper is dropped at the point where the distance between the center of mass of the interloper and the center of mass of the planetary system is r_tid(δ = 10⁻⁵). In the host cluster environment, we adopt the 1D velocity dispersions σ ∼ 0.2 and 1 km s⁻¹, following Shara et al. (2016), which are typical of open star clusters (e.g., Geller & Leigh 2015; Leigh et al. 2016b). The corresponding relative velocities at infinity are, respectively, v_∞ ∼ 0.73 and ∼3.6 km s⁻¹. The velocity is not in the resonant interaction regime, i.e., v_∞ > v_c. Therefore, we terminate the flyby process at t = 2t_drop, where the drop-in time t_drop is the time that the interloper takes to travel from the starting point to the point of closest approach. In the hyperbolic asymptotic limit,

$$\begin{eqnarray}&&{t}_{\mathrm{drop}}=M\sqrt{\displaystyle \frac{{a}^{3}}{{{GM}}_{\mathrm{tot}}}},\end{eqnarray} \tag{ 16 }$$

where M is the relative mean anomaly between the interloper dropping point and the distance of closest approach, a is the corresponding SMA of the incident hyperbolic orbit and M_tot is the total mass of the system.

After the flyby, the planetary system is perturbed, and high eccentricity excitations can occur due to several channels, including planet–planet scatterings (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Ford & Rasio 2006; Chatterjee et al. 2008), LK oscillations (Wu & Murray 2003; Fabrycky & Tremaine 2007; Katz et al. 2011; Naoz et al. 2012, 2011; Teyssandier et al. 2013; Petrovich 2015), and coplanar secular interactions (Li et al. 2014). The (excited) high eccentricity effectively increases the rate of tidal dissipation due to the tidal force at pericenter, which makes planet migration possible. We integrate the after-flyby planetary systems up to 1 Gyr, and discuss the results of the flybys at 0.01, 0.1, and 1 Gyr.

2.2.1. Tidal Dissipation

We adopt the static tidal mode by implementing the tidal perturbation acceleration between stellar particles and planet particles in the simulation. The tidal acceleration that the star exerts on the planet is

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{\mathrm{tid}}=-3{k}_{\mathrm{AM},{\rm{p}}}\displaystyle \frac{{{GM}}_{* }^{2}}{{M}_{{\rm{p}}}{r}^{2}}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{r}\right)}^{5}\left[\left(1+{\tau }_{{\rm{p}}}\displaystyle \frac{\dot{r}}{r}\right)\hat{{\boldsymbol{r}}}-({{\rm{\Omega }}}_{{\rm{p}}}-\dot{\theta }){\tau }_{{\rm{p}}}\hat{{\boldsymbol{\theta }}}\right],\end{eqnarray} \tag{ 17 }$$

where r, R_p, k_AM,p and τ_p are the particle distance, the planetary radius, the apsidal motion constant and the tidal time lag, respectively. Ω_p is the planetary spin frequency, where we use the orbit-averaged assumption (Hut 1981; Hamers et al. 2017),

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{p}}}=n\displaystyle \frac{1+\tfrac{15}{2}{e}^{2}+\tfrac{45}{8}{e}^{4}+\tfrac{5}{16}{e}^{6}}{(1+3{e}^{2}+\tfrac{3}{8}{e}^{4}){\left(1-{e}^{2}\right)}^{3/2}},\end{eqnarray} \tag{ 18 }$$

where $n=\sqrt{{{GM}}_{* }/{a}^{3}}$ is the mean motion of the planetary orbit and a and e are the SMA and the eccentricity, respectively. In this work, we use k_AM,p = 0.25 and τ_p = 0.66 s (Hamers et al. 2017). Typically, the timescale for stellar tides to operate is much longer than for planetary tides. Therefore, we do not take it into consideration.

In the N-body code, the variables a and e are calculated from the instant relative positions and velocities between particles, that is,

$$\begin{eqnarray}&&a=\displaystyle \frac{\mu }{\mu /| d{\boldsymbol{r}}| -| d{\boldsymbol{v}}{| }^{2}}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\boldsymbol{e}}=\left(\displaystyle \frac{| d{\boldsymbol{v}}{| }^{2}}{\mu }-\displaystyle \frac{1}{| d{\boldsymbol{r}}| }\right)d{\boldsymbol{r}}-\displaystyle \frac{d{\boldsymbol{r}}\cdot d{\boldsymbol{v}}}{\mu }d{\boldsymbol{dv}},\end{eqnarray} \tag{ 20 }$$

where μ = GM_tot is the gravitational parameter of the planet–star pair and d r and d v are the relative position and velocity between the planet–star pair, respectively. Since the accelerations are calculated between particle pairs in the code, the variables d r and d v are always calculated between particles. However, for systems with more than two particles, once a hierarchical triple system is formed, the outer object will orbit around the center of mass of the inner binary. Therefore, to calculate a and e of the outer orbit, the d r and d v must be the relative position and velocity between the outer object and the center of mass of the inner binary. In our case, the inner binary is a star–planet system, and hence the center of mass of the inner binary is close to the star. Therefore, it is a reasonably good approximation to simply use position and velocity of the star to calculate a and e of the outer orbit.

2.2.2. General Relativistic (GR) Effects

In secular eccentricity excitation processes (including the standard LK and coplanar secular processes), the periapsis of the inner orbit precesses around its center of mass. However, the general relativity effect precesses the orbit in the opposite direction of the LK oscillations, and can thus significantly suppress the eccentricity excitation (Liu et al. 2015). In this work, we implement GR precession via the first order post-Newtonian (PN) method. The corresponding pair-wise acceleration is

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{1\mathrm{PN}}=\displaystyle \frac{{{GM}}_{* }}{{c}^{2}{r}^{3}}\left[4({\boldsymbol{r}}\cdot \dot{{\boldsymbol{r}}})\dot{{\boldsymbol{r}}}+4{{GM}}_{* }\hat{{\boldsymbol{r}}}-(\dot{{\boldsymbol{r}}}\cdot \dot{{\boldsymbol{r}}}){\boldsymbol{r}}\right].\end{eqnarray} \tag{ 21 }$$

The second order PN term contributes very little to the GR precession, which can be regarded as a correction to the first order PN term with the efficiency being lower by an order of magnitude in v/c. Thus, it can be safely ignored in our calculations. The 2.5 PN gravitational wave radiation term in the planetary system can also be ignored due to its very long timescale. The timescale for merger due to gravitational radiation in a circular isolated binary with initial SMA a₀ is (Peters 1964)

$$\begin{eqnarray}&&\begin{array}{l}{\tau }_{\mathrm{GW}}(e=0)=\displaystyle \frac{5{c}^{5}{a}_{0}^{2}}{256{G}^{3}{m}_{12}{m}_{1}{m}_{2}}\\ =\ 3.4\times {10}^{19}\left(\displaystyle \frac{{M}_{\odot }}{{m}_{12}}\right)\left(\displaystyle \frac{{M}_{J}}{{m}_{1}{m}_{2}/{m}_{12}}\right){\left(\displaystyle \frac{{a}_{0}}{\mathrm{au}}\right)}^{4}\,\mathrm{yr},\end{array}\end{eqnarray} \tag{ 22 }$$

where m₁₂ is the total mass of the isolated binary and M_J is Jupiter mass. For an eccentric isolated binary, the corresponding timescale is

$$\begin{eqnarray}&&{\tau }_{\mathrm{GW}}(e\sim 1)=\displaystyle \frac{768}{425}{\left(1-{e}^{2}\right)}^{7/2}{\tau }_{\mathrm{GW}}(e=0).\end{eqnarray} \tag{ 23 }$$

In order for the timescale τ_GW(e ∼ 1) to be smaller than the age of the universe, e needs to be larger than ∼0.999, which happens in only a tiny fraction of our simulations even when eccentricity excitation mechanisms are important.

We have performed more than 1 million scattering experiments with fixed incident angles, closest approach distances, SMA ratios, and cluster velocity dispersions as described above. The orbits of all planets that remain bound to their host stars, described by their eccentricity—SMA distributions, are shown in Figures 2–5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Several features of Figures 2–5 encapsulate some of the key findings of this paper. They illustrate the time evolution of the orbits of Jupiters and Saturns that survive an encounter, for the cases a_s/a_j = 2, 4, and 8. The numbers of Jupiters reaching Saturn's orbit, and vice versa, is largest for a_s/a_j = 2. Perhaps most striking is that hot Jupiters are never immediately formed as an aftermath of a close flyby. In contrast, both warm Jupiters (with a_j ∼ 1 au) and "ultra-cold" Saturns (with a_s ∼ 100–10,000 au are immediately formed, and they persist for longer than 100 Myr. By 10⁶ yr hot Jupiters have been formed, and increase in number at 10⁷ and 10⁸ yr due to tidal circularization. Ultra-cold Saturns decrease in number as these weakly bound objects occasionally acquire enough energy to escape their host stars.

The different orbital outcomes for a range of incident angles and closest approaches by the intruder star, immediately after the flyby, and at 10⁶, 10⁷, and 10⁸ yr, respectively, are shown in Figures 6–9. Both Saturn and Jupiter's initial orbits are negligibly perturbed when the intruder's closest approach is further than twice Saturn's SMA, and these are omitted from the figures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As Figure 6 shows, a flyby star approaching as close as Saturn strongly perturbs that planet on the flyby timescale. It sends some Saturns inwards as close to their host stars as Jupiter, and flings others outwards as far as 10,000 au, beginning their ultra-cold Saturn existence. Strongly preferred and avoided regions in a–e phase space are obvious in this figure. On timescales of 10⁶–10⁸ yr both planets migrate enough, under their mutual gravitational influence and tides, to smear out the sharp a–e features of Figure 6, as shown in Figures 7–9.

Figure 7 shows that hot Jupiters are produced within 10⁶ yr of the encounter. Figure 7 suggests, and Figures 8 and 9 confirm, that high inclination flybys with closest approaches inside Jupiter's orbit are by far the most effective way to make hot Jupiters.

In each subplot of Figures 6–9, with fixed i and Q, the initial phases of the two planets are uniformly generated in time with eccentricity equal to zero. We see that for perturbations with Q > 2a_p, where a_p is the SMA of the outer planet (Saturn), the flyby star hardly changes the SMA of the Jupiter, regardless of the incident inclinations. With Q < 2a_p and even Q < a_p, where the interloper penetrates the planetary system, there are no directly formed hot Jupiters with ${a}_{{\rm{j}}}^{{\prime} }\lt 0.09\,\mathrm{au}$. This is because the flyby star must extract most of the energy from the Jupiter orbit to make it shrink below ∼0.09 au, which is almost impossible within a single flyby. From this plot we also see that, for a given i and Q, the cluster velocity dispersion has little effect on the distribution of the SMA of the Jupiter right after the flyby. This indicates that, in the strong interaction regime Q < a_p, the energy change of the inner planet is insensitive to the velocity dispersion of the environment. This is because the velocity of the interloper at the closest approach is

$$\begin{eqnarray}&&{v}_{Q}=\sqrt{{v}_{\infty }^{2}+\displaystyle \frac{2{{GM}}_{\mathrm{tot}}}{Q}},\end{eqnarray} \tag{ 24 }$$

where M_tot is the total mass of the scattering system. For a velocity dispersion σ = 0.2 or 1 km s⁻¹, v_∞ is negligible compared to 2GM_tot/Q with small Q. Therefore, the term v_Q is not sensitive to the environmental velocity dispersion in the strong interaction regime.

In each subplot, we also see clear differences between the prograde and retrograde orbits. For scatterings with fixed Q but different i, we see that in the prograde scatterings, the Jupiter shrinks its SMA more than in the case of retrograde scatterings. Since in the prograde scattering the interloper has more time to interact with the Jupiter at the point of closest approach due to the low relative velocity between the interloper and the planets, the interloper is able to extract more energy from Jupiter's orbit. This trend begins to disappear as the inclination i increases to 90°. In the 90° case, with the interloper traveling perpendicular to the planet orbital plane, the interaction time between the interloper and the planet is the same in the prograde and retrograde scatterings.

As shown in each row of Figure 10, as i increases, the SMA changes become smaller, due to the shorter interaction time between the interloper and planet at the point of closest approach. As i increases, the maximum energy that can be extracted from the planetary system becomes smaller.

Zoom In Zoom Out Reset image size

We also see in many of the subplots in Figure 10 that there are two peaks in both prograde and retrograde scatterings. In total, we have four peaks from different relative positions at closest approach: (i) prograde scattering with the interloper and planet at the same side of the star (the interloper and planet move in the same direction with small relative distance, thus contributing to the dominant peak), (ii) retrograde scattering with the interloper and planet at the opposite side of the star (the interloper and planet move in the same direction with large relative distance), which contributes to the secondary dominant peak, (iii) prograde scattering with the interloper and planet at opposite sides of the star (the interloper and planet move in opposite directions with large relative distance), which contribute to the third prominent peak, and (iv) retrograde scattering with the interloper and planet at the same side of the star (the interloper and planet move in the opposite direction with small relative distance), which contributes to the weakest peak that is almost invisible.

The distance between the peaks becomes smaller as Q increases. The reason is that, for large Q, the distance between the interloper and planet, Q ± a_j, is dominated by Q instead of by ±a_j. Thus, no matter which side of the star the planet is located on, i.e., +a_j or −a_j, the distance between the interloper and the planet in either case becomes similar.

Figure 11 shows the distribution of the SMA of the Saturn right after the flyby. Because the outer planet orbit is more weakly bound, it is easier to modify the orbit of the Saturn relative to that of the Jupiter. In general, the orbit of the Saturn right after the flyby shows properties similar to those of the Jupiter. However, we see more differences arising when adopting different velocity dispersions, especially in prograde scatterings. This is because, with fixed Q, the relative velocity between the interloper and the planets determines the interaction time between the interloper and the planet, which affects the properties of the planet right after the flyby. The relative velocity between the interloper and planet at the distance of closest approach is determined by the planet's orbital velocity and the velocity of the interloper at closest approach v_Q. In the case of Jupiter, due to the tight orbit, the orbital velocity is dominant over the velocity of the interloper, thus the velocity dispersion does not affect much the relative velocity between the interloper and the planet at the point of closest approach. However, for the Saturn with smaller orbital velocity, the orbital velocity becomes less dominant over the velocity of the interloper. Thus, the velocity dispersion has a larger impact on the relative velocity between the interloper and the Saturn at the point of closest approach.

Zoom In Zoom Out Reset image size

The SMA change of the planets can be estimated by

$$\begin{eqnarray}&&{dE}={E}_{{\rm{p}}}-{E}_{{\rm{p}},0}=-\displaystyle \frac{{{GM}}_{* }{M}_{{\rm{p}}}}{2{a}_{{\rm{p}}}^{{\prime} }}+\displaystyle \frac{{{GM}}_{* }{M}_{{\rm{p}}}}{2{a}_{{\rm{p}}}},\end{eqnarray} \tag{ 25 }$$

where a_p is the initial SMA and ${a}_{{\rm{p}}}^{{\prime} }$ is the SMA immediately after the flyby. It is straightforward to note that, if (E_p − E_p,0)/∣E_p,0∣ < 0, the corresponding orbit shrinks, while if 0 < (E_p − E_p,0)/∣E_p,0∣ < 1, the orbit expands. For (E_p − E_p,0)/∣E_p,0∣ > 1, the planet is ejected from its host star. As the incident angle increases, the maximum relative binding energy shift becomes smaller as the interaction time at the point of closest approach decreases. Therefore, the binding energy shift becomes weaker as Q increases.

The prograde scatterings are responsible for the majority of the planet ejections, where (E_p − E_p,0)/∣E_p,0∣ > 1. The prograde scatterings with the interloper and planets at the same side of the host star have the longest interaction times (and hence impart the largest impulse) needed to remove the binding energy from the planet. Since the retrograde encounter is not as efficient as the prograde case in removing the binding energy of the planet due to its shorter interaction time, the retrograde scattering can hardly change the SMA of the Jupiter at Q ∼ 4a_j and the SMA of the Saturn at Q ∼ 4a_s, whereas prograde scatterings do change the SMA of the planet significantly (by up to about one order of magnitude).

Figures 12 and 13 show the eccentricity of the Jupiter and Saturn right after the flyby. Since the initial value of the eccentricity of the planets is zero, this also indicates the eccentricity change right after the first flyby. For Q < 4a_j and Q < 4a_s (i.e., the strong interaction regime), a flyby can significantly kick the eccentricity of a planet to a large value. For Q < a_j or Q < a_s, where the interloper penetrates within the orbit of the planetary system, the majority of the planets are perturbed into the "near ejection" state where e ∼ 1. As for Q > 4a_j and Q > 4a_s, where scattering happens in the perturbative regime, the eccentricity shift due to the flyby rapidly drops below 0.1. These results indicate that just one flyby does not efficiently increase the eccentricity of a planet, which is essential to enable the subsequent high eccentricity tidal migration process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figures 12 and 13, we can also see that as the inclination i increases, retrograde scatterings become similar to the prograde scatterings in removing the angular momentum of the planets.

Figure 14 shows the inclination between the Jupiter and Saturn orbits right after the flyby if both the Jupiter and Saturn still remain bound to their original host star. When the closest approach Q is smaller than a_p, we see that it is possible for the interloper to flip the orbital orientation of the Saturn to the opposite direction by transferring sufficient angular momentum to its orbit. Also, as the incident inclination i increases, we see that the probability of creating an inclined Saturn orbit with respect to the Jupiter orbit increases as well, especially in the prograde flyby case. It is difficult to get a highly inclined orbit when Q > a_p, such that LK oscillations are unlikely to operate on the system after the flyby to help make hot Jupiters.

Zoom In Zoom Out Reset image size

In Figures 10–14, we explore the properties of the planetary system right after the flyby with two different values of the environmental velocity dispersion, σ = 0.2 or 1 km s⁻¹. We find that close flybys create an ideal situation for tidal migration to produce hot Jupiters. Close flybys can excite the eccentricity of the planets by extracting significant angular momentum from planetary systems and kick the original coplanar planetary system into an inclined planetary system. These two processes increase the efficiency of the tidal migration process.

We find that the orbital properties are insensitive to the low velocity dispersion. The velocity of the interloper at infinity contributes very little to the velocity of the interloper at the distance of closest approach for small Q due to the strong gravitational focusing, as indicated in Equation (24).

Besides exploring how the cluster velocity dispersion changes the orbital properties of the planetary system after the flyby, we also explore the flyby scattering with different initial SMA ratios. We fix the velocity dispersion to be 1 km s⁻¹ and perform scatterings with initial SMA ratios equal to 2, 4, and 8.

Before proceeding, we briefly comment on the possible impact of orbital resonances on our results. The naive expectation here is that, for smaller initial SMA ratios, such resonances should be more important as the initial SMA ratio is closer to a strong resonance (e.g., the 2:1 resonance for an initial SMA ratio of 2). As the initial ratio of SMAs increases, we expect resonances to become less important, as these resonances would tend to operate at very high ratios of the orbital periods, where we expect resonances to be very weak and have little to no dynamical effects. In our simulations, we see only mild evidence for systems becoming locked into orbital resonances for initial SMA ratios of 2, with at most of order a few percent of our total simulations ending up stabilized in a resonant state. For higher initial SMA ratios, these effects disappear. This is expected since the perturbed planetary systems will be scattered closest to mainly higher-order resonances, which tend to get weaker as the ratio of periods increases. In the end, we find that orbital resonances have a very negligible effect on our overall results, statistics, and conclusions, at least for the initial conditions considered here. Thus, these effects can safely be ignored in this section when trying to understand how the product distributions from our simulations depend on the initial SMA ratio.

As in Figure 10, Figure 15 shows the SMA of the Jupiter right after the flyby for different initial planet SMA ratios (see also Figures 16–19). It is indicated in the figure that, as we double the SMA ratio, the distribution of a_j is almost identical to the distribution of a_j with the original ratio but with Q doubled, e.g., the a_s/a_j = 8 with Q = 0.125a_p is almost identical to a_s/a_j = 4 with Q = 0.25a_p. It is not surprising that the absolute closest approach determines the post-flyby SMA and eccentricity (as shown in Figure 17) of the Jupiter. The reason is that the inter-planet interaction in the close flyby is negligible due to the short interaction flyby time. Thus, we can treat the perturbation exerted on each planet almost independently.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As already noted and discussed, orbital inclinations are created by highly inclined flybys. The changing z-direction to the interloper changes the direction and magnitude of the angular momentum of the planets. For fixed Q, as the SMA ratio decreases, the SMA of the Jupiter decreases while keeping the SMA of the Saturn constant. Since the angular momentum of the planet is

$$\begin{eqnarray}&&{L}_{{\rm{p}}}\propto \sqrt{a(1-{e}^{2})},\end{eqnarray} \tag{ 26 }$$

then, for a smaller SMA, it is easier to change the direction of the inner planet by importing/extracting the same angular momentum to/from the planet's orbit. Hence, we can see in Figure 19 that, with a smaller SMA ratio, the planetary system tends to be more inclined after the flyby. This argument is also valid in coplanar flybys: smaller SMA ratios also create more retrograde orbits in all those interactions corresponding to the closest approaches.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can further analyze the SMA ratio after the flyby with the same method by treating the perturbation exerted on each planet independently. The binding energy of the planet is

$$\begin{eqnarray}&&{E}_{{\rm{p}}}\propto 1/a.\end{eqnarray} \tag{ 27 }$$

Therefore, it is harder to change the SMA of the planet with smaller initial SMA by importing/extracting the same energy to/from the planet's orbit. With fixed Q, a smaller initial SMA ratio gives a smaller absolute SMA for the Jupiter. It becomes more difficult to change the SMA of the Jupiter. Thus, the final SMA ratio right after the flyby depends more on the change in SMA of the Saturn.

Since the absolute SMA of the planets determines if the interaction between the planet and the interloper is in the fast/slow regime, it is difficult to make a simple argument on how the final SMA ratio depends only on the initial SMA ratio. For instance, if the initial SMA ratio is 4, where both Saturn and Jupiter reside in the slow regime where the flyby tends to increase the SMA, then the final SMA tends to become larger. The reason is that, while the SMA of both planets tends to become larger, the change in Saturn's SMA is more significant. If the initial SMA ratio is 2, Jupiter's orbit may enter into the faster regime, where the flyby tends to tighten Jupiter's orbit while keeping the Saturn in the slow regime. Thus, although smaller initial SMA ratios make the final SMA ratio more weakly dependent on Jupiter's orbit, the final SMA ratio tends to become larger than in the case where the initial SMA ratio is equal to 4 because the two planets migrate in opposite directions.

To make the planets interact more efficiently, the SMA ratio of the planets right after the flyby should not be too large nor too small. Two close orbits make planet–planet scattering possible during the subsequent long-term evolution after the flyby and shorten the timescale for secular processes to operate (under the condition that the system is stable).

Because of the short flyby time, the interaction between the planets during the flyby at the point of closest approach can be safely ignored. Hence, we can treat the perturbation exerted by the interloper on the two planets independently. Therefore, the scattering with different initial SMA ratios affects the properties of the planets independently, as shown in Figures 15–18.

However, the initial SMA ratio affects the inclination and final SMA of the planetary system right after the flyby, thus influencing the probability and timescale for subsequent tidal migration. We will discuss the corresponding "sweet spot" for the initial SMA ratio in the next section.

As listed in Section 3, the outcomes after the flyby are classified into "planet–planet collision," "star–planet collision," "star–star collision," "Jupiter ejection," "Saturn ejection," "both ejection," and "both remain." In this subsection we will discuss the fraction of each outcome for different environmental velocity dispersions and initial SMA ratios. We pay special attention to the both-remain systems since these are obviously required for post-interaction tidal migration.

Figure 20 shows the standard case where a_s/a_j = 2 and σ = 1 km s⁻¹. We can see that for Q > 4a_p, the flyby will never change the architecture of the planetary system, and therefore all the outcomes are both remain. As the closest approach comes close to the planetary system with Q ∼ 2a_p, the flyby begins to eject the outer planet Saturn and then eject the inner planet Jupiter. As shown in the row with Q = 2a_p, high inclination incident angles where the interaction time is shorter decrease the ejection fraction for Saturn alone.

Zoom In Zoom Out Reset image size

The fraction of the both-remain systems in the close flyby regime where Q < a_p is roughly 10%. We also find that most of the collisions are planet–star collisions that only appear in coplanar scatterings. The reason for this is related to the volume-filling factor of the objects, which is much smaller for coplanar configurations relative to the isotropic case (Leigh et al. 2017, 2018). In other words, there is more free space in the isotropic case, increasing the timescale for two objects to coincide in both time and space.

As an aside we note that planet–star collisions can give rise to a rich set of possible outcomes (De Marco & Soker 2011; Schuler et al. 2011; Ginsburg et al. 2012; Kratter & Perets 2012; Veras et al. 2013; Deal et al. 2015). In particular, planets that are engulfed may shape the ejecta of their host stars, while swallowed and disrupted planets may enrich the outer layers of stars in detectable lithium, calcium, magnesium, and iron.

Figure 21 shows the results for the same initial conditions as in Figure 20, but with σ equal to 0.2 km s⁻¹. As discussed before, v_∞ contributes only minimally to the velocity of the interloper at the distance of closet approach v_Q. Therefore, the fractions are very similar to the standard case with σ equal to 1 km s⁻¹. Hence, we may conclude that in stellar clusters with low velocity dispersions where gravitational focusing is significant, the fraction of planetary systems that may form hot Jupiters after a close flyby is not particularly sensitive to the velocity dispersion.

Zoom In Zoom Out Reset image size

Figures 22 and 23 show the cases where a_s/a_j = 4 and 8, respectively. With a larger SMA ratio, we see a larger fraction in the ratio of Saturn ejections, and a much lower both-ejection fraction. Since planet ejection requires a close encounter with the interloper, it is more difficult for the interloper to have a close encounter with both planets if the two planets are far apart.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In a hierarchical triple system, the inner orbit (here, the star-Jupiter binary) can undergo LK oscillations. In LK oscillations, the energy (∝ − M_tot,in/a) of each orbit is conserved but the inner and outer orbits continuously exchange angular momentum ($\propto \sqrt{a(1-{e}^{2})}$). Therefore, the eccentricity of the inner orbit can be excited to very high values as its angular momentum reaches a minimum, which is balanced to conserve angular momentum by reducing the orbital inclination between the inner and outer orbits. The timescale for the quadrupole LK cycle to operate is

$$\begin{eqnarray}&&{\tau }_{\mathrm{LK}}=\displaystyle \frac{1}{{n}_{\mathrm{in}}}\left(\displaystyle \frac{{M}_{\mathrm{tot},\mathrm{in}}}{{M}_{\mathrm{out}}}\right){\left(\displaystyle \frac{{a}_{\mathrm{out}}}{{a}_{\mathrm{in}}}\right)}^{3}{\left(1-{e}_{\mathrm{out}}^{2}\right)}^{3/2},\end{eqnarray} \tag{ 31 }$$

where n_in is the mean motion of the Jupiter orbit, M_tot, in is the total mass of the inner orbit, and M_out is the mass of the perturber. LK oscillations require the triple system to be hierarchical and stable. These oscillations operate in the range where $-\sqrt{3/5}\lt \cos (I)\lt \sqrt{3/5}$, with I being the inclination between the two planets' orbits.

In this work, the flyby perturbation generates the necessary inclination, and the two giant planet orbits ensure that the timescale τ_LK is shorter than the age of the universe. This combination of factors allows the eccentricity of the Jupiter to be increased to a high enough value that tidal dissipation becomes efficient enough to make a hot Jupiter.

An example of hot Jupiter formation from the LK effect is shown in Figure 24. It can be seen that the effect starts to operate on a timescale of about 1–1.5 Myr.

Zoom In Zoom Out Reset image size

Close planet–planet scatterings can convert Keplerian shear into an angular momentum deficit, triggering high eccentricity migration (e.g., Rasio & Ford 1996; Weidenschilling & Marzari 1996; Ford & Rasio 2006; Chatterjee et al. 2008). In planet–planet scattering, the planets typically undergo several close encounters to grow their eccentricities to high values. However, if the orbital velocity of a planet is larger than the escape velocity at the planet's surface, the cross section for collision would be larger than the cross section for scattering. This may lead to a planet collision rather than just eccentricity growth. Goldreich et al. (2004), Ida et al. (2013) and Petrovich et al. (2014) give the upper limit of the eccentricity from planet–planet scattering as

$$\begin{eqnarray}&&{e}_{{\rm{p}}-{\rm{p}}}\sim \displaystyle \frac{\sqrt{{{GM}}_{{\rm{p}}}/{R}_{{\rm{p}}}}}{2\pi {a}_{{\rm{p}}}/P},\end{eqnarray} \tag{ 32 }$$

where M_p is the mass of planet, R_p is the radius of the planet, a_p is the SMA of the planet, and P is the orbital period of the planet. A Jupiter mass planet with a_p ∼ 1 au has roughly e_p−p ∼ 1, which makes high eccentricity tidal migration possible.

An example of hot Jupiter formation via planet–planet scattering is shown in Figure 25. This occurs rather continuously up to a time of about 13 Myr, when the orbit of Jupiter achieves circularization. After that time, Jupiter and Saturn no longer interact strongly.

Zoom In Zoom Out Reset image size

LK oscillations operate only when the orbital inclination is within a certain range. However, Li et al. (2014) found that a perturber with a highly elliptical orbit, where the octupole LK effect is strong, can also excite the eccentricity of an initially coplanar system. Thus, tidal migration may also occur in this regime. Indeed, we find that usually the flyby star extracts the angular momentum of the outer planet rather than that of the inner planet. This preferentially forms an elliptical outer orbit for the Saturn, which drives the further evolution of the system. The coplanar octupole effect comes into play when

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

where e_in is the initial eccentricity of the inner binary, ω is the argument of the periapsis of the outer orbit, and is the LK octupole parameter

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}}\displaystyle \frac{{e}_{\mathrm{out}}}{1-{e}_{\mathrm{out}}^{2}}.\end{eqnarray} \tag{ 34 }$$

No hot Jupiters were found to be formed via this channel in our simulations.

There are other mechanisms that can excite the eccentricity of a planet, including LK oscillations from a stellar companion (Naoz et al. 2012; Anderson et al. 2016) and secular chaos due to multiple LK timescales. Stellar LK oscillations demand a companion star, either from binary star formation as an outcome of the flyby (which is very unlikely) or a flyby binary interloper, which is beyond the scope of this work. Secular chaos generally demands more than two planets, also beyond the scope of this paper. Thus, we do not consider eccentricity excitation from these mechanisms.

After the flyby, the modified planetary system may fall into the eccentricity excitation regime, where hot Jupiter formation is possible. For all the both-remain systems, as discussed in more detail in the preceding section, we divide the outcomes into several different regimes based on the criteria below:

LK: The inclination between the two planets' orbits satisfies $-\sqrt{3/5}\lt \cos (I)\lt \sqrt{3/5}$ and the quadrupole timescale is τ_LK < 1 Gyr.

Planet–planet scattering: The pericenter of the outer planet/the apocenter of the inner planet satisfies $\tfrac{{a}_{{\rm{s}}}(1-{e}_{{\rm{s}}})}{{a}_{{\rm{j}}}(1+{e}_{{\rm{j}}})}\lt 1$, where the planets could undergo close encounters.

Coplanar octupole: The inclination between the two planet orbits is within 5° and Equation (33) is satisfied.

After the flyby, we perform long-time integrations of the perturbed planetary systems up to 1 Gyr, including tidal dissipation and GR effects as described in Section 2.2. Since for Q > 2a_p the perturbation becomes too weak to transition the planetary system into the high eccentricity tidal migration regime, we only perform long integrations for systems with Q < 2a_p.

During the long-time evolution, tidal dissipation and GR effects may break the planetary system by ejecting either the Jupiter or Saturn through a close planet–planet encounter or collision. If the eccentricity of the remaining planet is not high enough to make tidal dissipation efficient, then we will not observe it to form a hot Jupiter within the age of the universe.

Figure 26 shows the distribution of the disruption time (i.e., the timescale over which the planetary system is disrupted) in log scale. For the coplanar case, this timescale is much shorter than for the other cases. The planetary system tends to disrupt in ∼10⁴ yr, while the other cases are more likely disrupt at around 10⁶ yr. This is because for the coplanar case, there is a much higher probability for the two planets to undergo close encounters. The distribution is truncated by the flyby time on the left side and by the integration upper limit on the right side.

Zoom In Zoom Out Reset image size

Figure 27 shows the SMA of the Jupiter right after the flyby (gray), 10⁶ yr after the flyby (blue), 10⁷ yr after the flyby (purple), and 10⁸ yr after the flyby (orange). We see that, for Q < 0.5a_p, a non-negligible fraction of the Jupiters evolve into hot Jupiters with a_j between 10⁻² and 10⁻¹ au due to tidal dissipation and the aforementioned high eccentricity excitation mechanisms. For Q > 1 a_p, it is extremely difficult to form a hot Jupiter progenitor. We do not see any hot Jupiters even 10⁸ yr after the flyby.

Zoom In Zoom Out Reset image size

The evolution of the orbital SMA of Jupiter with time, at several representative times up to 1 Gyr, is shown in Figure 28. For most configurations, hot Jupiter formation happens promptly, already present at 1 Myr.

Zoom In Zoom Out Reset image size

Figure 29 shows the same distribution as Figure 27 but with an initial SMA ratio equal to 4. A large fraction of the Jupiters evolve into hot Jupiters with Q = 0.125a_p and i = 90°. This is the ideal parameter space where the flyby has a large probability to create a highly inclined system with an appropriate SMA ratio that boosts the LK effect. For an initial SMA ratio equal to 2, as indicated in Figure 27, the two giant planets are too close to maintain the architecture of the planetary system. Thus, there is a higher chance for the planetary system to undergo close planet–planet interactions rather than secular LK excitations. For an initial SMA ratio equal to 8, as shown in Figure 30, we see that the fraction of the hot Jupiters drops again in highly inclined flybys. The timescale of the quadrupole LK effect is given by Equation (31). Unlike the LK mechanism in binary star systems, where the perturber is a star, the perturber in the planetary system is only a giant planet. Thus, the system requires a much smaller SMA ratio to make the LK effect operate within the age of the cluster. Therefore, with the larger SMA ratio equal to 8, the LK effect becomes less efficient than for an SMA ratio equal to 4. For a Jupiter mass planet-solar mass host star binary with a_j = 5 au and with a Saturn mass perturber at 10 au, the typical timescale given by Equation (31) that is sensitive to the SMA ratio is ∼2.5 × 10⁶ yr. With an initial SMA ratio equal to 8, the timescale is about one order of magnitude longer. This explains why the hot Jupiter formation fraction drops quickly for a_j/a_s = 8. However, the SMA ratio cannot be too small in the LK regime, otherwise it breaks the stability of the planetary system. Therefore, for a_j/a_s = 2, the main source of hot Jupiters is from planet–planet scatterings.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 31–33 show the distribution of the eccentricity of the Jupiter right after the flyby (gray), 10⁶ yr after the flyby (blue), 10⁷ yr after the flyby (purple), and 10⁸ yr after the flyby (orange). In the close flyby regime where Q < 2a_j, the general trend of the Jupiter eccentricity after long interaction times is to decrease, while in the perturbative regime where Q > 2a_j the general trend of the Jupiter eccentricity is to increase. Indeed, in the close flyby regime, the interloper increases the eccentricity of the Jupiter by extracting angular momentum from the planet. There is a large probability of creating a Jupiter orbit with a small pericenter. However, in the perturbative regime, the flyby interloper barely creates a small pericenter Jupiter orbit. Therefore, in the close flyby regime, the tidal dissipation will dominate in the long-time integration, while in the perturbative regime, planet–planet interactions become the dominant effect. The tidal dissipation process will circularize the Jupiter orbit but planet–planet interactions will increase the eccentricity of the Jupiter orbit. Therefore, we see opposite general trends for the eccentricity evolution in the close flyby and perturbative regimes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 34 shows the hot Jupiter formation fraction of all both-remain systems after the flyby. Different colors represent the fractions from different channels and the sum of the fractions indicates the total formation fraction. We find that there is no hot Jupiter formed from the coplanar secular effect due to its very narrow parameter space, such that flybys cannot create the ideal initial conditions in this regime. For a_s/a_j = 2, most of the hot Jupiters formed from planet–planet interactions. For a_s/a_j = 4, the LK channel contributes most of the hot Jupiters, especially the perpendicular flybys with a formation fraction one order of magnitude higher. For a_s/a_j = 8, we see that most of the hot Jupiters formed from the LK channel. Indeed, with a_s/a_j = 2, the flybys more easily create closer planet orbits to activate planet–planet scatterings. However, with a_s/a_j = 8, where the two orbits are distant from each other, it is harder for the flybys to bring the two orbits close enough together. Therefore, we see that hot Jupiter formation from planet–planet scattering dominates in the small SMA ratio regime, while hot Jupiter formation from secular LK effects dominates in the large SMA ratio regime. Figure 35 shows the formation fraction with lower velocity dispersion, where σ = 0.2  km s⁻¹. The formation fraction is slightly higher than for σ = 1 km s⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To estimate the hot Jupiter formation rate from different channels in different environments, we need the hot Jupiter formation fractions of each channel from isotropic scatterings. Then, the hot Jupiter formation rate in each channel can be estimated by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{HJ}}={{\rm{\Gamma }}}_{\mathrm{strong}}f={\sigma }_{{{\rm{a}}}_{p}}{{nv}}_{\mathrm{rms}}f,\end{eqnarray} \tag{ 35 }$$

where ${\sigma }_{{{\rm{a}}}_{p}}$ is the corresponding cross section for Q = 1a_p for which a hot Jupiter can be formed (Q < a_p), n is the number density of the environment, v_rms is the root mean square velocity of the environment, and f is the hot Jupiter formation fraction of each channel after a close flyby. The cross section is given by

$$\begin{eqnarray}&&{\sigma }_{{{\rm{a}}}_{p}}=\pi {b}_{{{\rm{a}}}_{p}}^{2},\end{eqnarray} \tag{ 36 }$$

where

$$\begin{eqnarray}&&{b}_{{{\rm{a}}}_{p}}={a}_{{\rm{p}}}\sqrt{1+2\displaystyle \frac{G({M}_{0}+{M}_{1})}{{v}_{\infty }^{2}{a}_{{\rm{p}}}}}.\end{eqnarray} \tag{ 37 }$$

The number density of a virialized cluster is (e.g., Binney & Tremaine 1987; Spitzer 1987)

$$\begin{eqnarray}&&n\sim \displaystyle \frac{{M}_{{\rm{c}}}/\bar{m}}{4\pi {R}_{{\rm{c}}}^{3}/3}\sim \displaystyle \frac{6{v}_{\mathrm{rms}}^{6}}{\pi {G}^{3}{M}_{{\rm{c}}}^{2}\bar{m}},\end{eqnarray} \tag{ 38 }$$

where M_c and $\bar{m}$ are the total stellar mass and mean stellar mass of the cluster, respectively. The upper panel of Figure 36 shows the number density of the virialized open, young massive and globular clusters as a function of the velocity dispersion and the bottom panel of Figure 36 shows the number of the "close" flyby per star as a function of the velocity dispersion. Then, the hot Jupiter formation rate can be estimated from

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Gamma }}}_{\mathrm{HJ}}=\pi \left({a}_{{\rm{p}}}^{2}+2\displaystyle \frac{G({M}_{0}+{M}_{1})}{{v}_{\infty }^{2}}{a}_{{\rm{p}}}\right)\displaystyle \frac{6{v}_{\mathrm{rms}}^{7}}{\pi {G}^{3}{M}_{{\rm{c}}}^{2}\bar{m}}f\\ =\ \left(\displaystyle \frac{6{v}_{\mathrm{rms}}^{7}{a}_{{\rm{p}}}^{2}}{{G}^{3}{M}_{{\rm{c}}}^{2}\bar{m}}+\displaystyle \frac{12{v}_{\mathrm{rms}}^{5}{a}_{{\rm{p}}}}{{G}^{2}{M}_{{\rm{c}}}^{2}}\right)f\\ =\ \displaystyle \frac{6{v}_{\mathrm{rms}}^{5}{a}_{{\rm{p}}}}{{G}^{2}{M}_{{\rm{c}}}^{2}}\left(\displaystyle \frac{{v}_{\mathrm{rms}}^{2}}{{v}_{{\rm{p}}}^{2}}+2\right)f,\end{array}\end{eqnarray} \tag{ 39 }$$

with ${M}_{0}\sim {M}_{1}\sim \bar{m}$, ${v}_{\mathrm{rms}}^{2}=2{v}_{\infty }^{2}$, and $G\bar{m}/{a}_{{\rm{p}}}\sim {v}_{{\rm{p}}}^{2}$. For most clusters v_rms ≪ v_p, therefore,

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Gamma }}}_{\mathrm{HJ}}=\frac{12{v}_{\mathrm{rms}}^{5}{a}_{{\rm{p}}}}{{G}^{2}{M}_{{\rm{c}}}^{2}}f\\ =\ 1.6\times {10}^{-4}{\left(\frac{\sigma }{1\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{5}\left(\frac{{a}_{{\rm{p}}}}{20\ \mathrm{au}}\right){\left(\frac{{M}_{{\rm{c}}}}{1000{M}_{\odot }}\right)}^{-2}\,{\mathrm{Gyr}}^{-1}.\end{array}\end{eqnarray} \tag{ 40 }$$

Since the hot Jupiter formation timescale is typically much shorter than the age of its host star cluster, the hot Jupiter percentage per star in a cluster is

$$\begin{eqnarray}&&{F}_{\mathrm{HJ}}={{\rm{\Gamma }}}_{\mathrm{HJ}}{N}_{* }{\tau }_{{\rm{c}}}/{N}_{* }=\displaystyle \frac{12{v}_{\mathrm{rms}}^{5}{a}_{{\rm{p}}}}{{G}^{2}{M}_{{\rm{c}}}^{2}}f{\tau }_{{\rm{c}}},\end{eqnarray} \tag{ 41 }$$

where τ_c is the age of the cluster.

We performed isotropic scatterings with cluster velocity dispersions σ =  0.2 and 1 km s⁻¹. The impact parameter b of each scattering was randomly generated from the probability distribution function f(b) ∝ b within [0, b_max], where b_max is chosen to be large enough that the system is not significantly perturbed. Based on our scattering experiments described in the last section, the planetary systems after the flyby perturbation will be almost unchanged if Q > a_p. Therefore, we set the corresponding b with Q = a_p as the b_max in isotropic scatterings. All other setups are the same as described in Section 2.

Finally, it is worth emphasizing that the properties of star clusters change over time, primarily due to two-body relaxation and the evaporation of the cluster in a Galactic tidal field. We have not accounted for these effects in any systematic way in this paper, and instead defer properly accounting for these effects to future work. Briefly, the rate of stellar escape is inversely proportional to stellar mass, causing older clusters to become more depleted in low-mass stars and tending toward a top-heavy stellar mass function (relative to standard initial mass functions measured in young star-forming regions). In future work, a Monte Carlo approach will be needed for a given host galaxy to generate a realistic initial star cluster mass function, a realistic distribution of 3D star cluster galactocentric distances, and so on.