EXTRASOLAR BINARY PLANETS. I. FORMATION BY TIDAL CAPTURE DURING PLANET–PLANET SCATTERING

We consider planetary systems that are composed of a central star and three gas giant planets. We take the origin of the coordinate at the central star with mass M_*. The equation of motion of planet i is

$$\begin{eqnarray} \frac{d^2\mathbf {r}_i}{dt^2} &=& -G\frac{M_\ast +M_i}{r_{i}^3}\mathbf {r}_{i}-GM_j\left(\frac{\mathbf {r}_j}{r_{j}^3}+\frac{\mathbf {r}_{ij}}{r_{ij}^3}\right) \nonumber\\ && -GM_k\left(\frac{\mathbf {r}_k}{r_{k}^3}+\frac{\mathbf {r}_{ik}}{r_{ik}^3}\right), \end{eqnarray} \tag{ 1 }$$

where M_i and r_i are masses and position vectors of planet i (= 1, 2, and 3), respectively, and r_ij ≡ r_i − r_j.

When a planet passes the pericenter to the star or another planet, we impulsively dissipate the orbital energy of the passing bodies according to tidal interactions, following Nagasawa & Ida (2011). We impose the tidal interactions with the host star only for encounters with the stellarcentric pericenter distance q < 0.04 AU to save computational time.

When a planet repeatedly undergoes close encounters with the central star and suffers the tidal dissipation, its semimajor axis and orbital eccentricity shrink, keeping q almost constant (e.g., Nagasawa et al. 2008; Nagasawa & Ida 2011). As a result, the planet becomes a hot Jupiter.

If a planet encounters another planet closely enough, the planets can be trapped through planet–planet tidal dissipation to form a gravitationally bound pair. At the initial phase after the trapping, eccentricity of the binary orbit is generally close to unity. Through repeated encounters, the binary orbit is circularized in a similar way to the formation of hot Jupiters. Since the relative motion is described by hyperbolic orbits at the trapping, and by highly eccentric orbits in most of the time during the circularization, the tidal interaction is dominated by dynamical tide, which is described in the following.

For tidal trapping and circularization of binary orbits, we use the formula for energy loss due to dissipation as a result of dynamical tide between two objects derived by Portegies Zwart & Meinen (1993), following Podsiadlowski et al. (2010). When planets i and j undergo tidal interactions, the tidal energy loss caused by a single close encounter between the planets is

$$\begin{eqnarray} E_{{\rm tide}} &=& \frac{GM_j^2}{R_i}\left[\left(\frac{R_i}{q_{ij}}\right)^6 T_2(\eta _i)+\left(\frac{R_i}{q_{ij}}\right)^8 T_3(\eta _i)\right] \nonumber \\ && + \frac{GM_i^2}{R_j}\left[\left(\frac{R_j}{q_{ij}}\right)^6 T_2(\eta _j)+\left(\frac{R_j}{q_{ij}}\right)^8 T_3(\eta _j)\right], \end{eqnarray} \tag{ 2 }$$

where q_ij is the pericenter distance between planets i and j, R_i is planetary physical radius, η_i ≡ {M_i/(M_i + M_j)}^1/2(q_{i, j}/R_i)^3/2, and T_{2, 3}(η_i) are given by Portegies Zwart & Meinen (1993) as a fifth-degree polynomial function.

For encounters between the planets with the pericenter distance less than 10(R_i + R_j), we subtract tidal dissipation energy from their orbital energy at the pericenter passage using impulse approximation. Because tidal interactions rapidly weaken as the distance between the bodies increases, even if we switch on the tidal interaction at larger distance, the results hardly change (see Section 3.2).

Strictly speaking, it may not be correct to use the impulse approximation and Equation (2) after the eccentricity of the binary orbit is significantly damped, because the oscillation of the planetary bodies raised by the tides is not dissipated before the next encounter. In such situations, how the semimajor axis and eccentricity of the binary orbit are changed by dynamical tides depends on the phase of oscillation, and a and e change chaotically (e.g., Press & Teukolsky 1977; Mardling 1995). However, we are most concerned with early phase of the tidal trapping to evaluate the formation probability of binary planets. So, we neglect the chaotic evolution due to incomplete oscillation damping.

Note that the formulas for dynamical tide could include large uncertainty, because the effect of dynamical tide depends on what wave modes are excited and how they dissipate. However, the tidal dissipation energy is inversely proportional to several powers of mutual distance between the planets. In the case of Equation (2), $E_{\rm tide} \propto q_{ij}^{-6}$. Even if E_tide changes by a factor of 10, the capture separation changes only by 50%.

After the circularization, quasi-static tides become predominant instead of dynamical tides. Although quasi-static tides are weaker by orders of magnitude than dynamical tides for high orbital eccentricity, quasi-static tides work also for circular binary orbits and last until a spin-orbit synchronous state is established. Thus, the cumulative effect in orbital evolution due to quasi-static tide is important when we consider observational detectability of binary planets. In Section 3.4, we calculate such longer-term tidal evolution due to quasi-static tides between planets, taking into account planetary spins and also tides between the planets and the host star. The formulation is based on tidal dissipation functions Q, instead of Equation (2), as explained in the Appendix. We discuss orbital stability against tidal evolution due to quasi-static tides during the main-sequence lifetime of solar-type host stars (∼10 Gyr). We will show that a binary system is stable if stellarcentric distance is ≳0.3–0.4 AU. Also note that as long as the tidal capture occurs beyond ∼0.3 AU from the central star, envelope removal due to the tidal dissipation does not take place (see Section 3.2).

To evaluate the formation rate of the binary planets, we set three giant planets. In early studies, planet–planet scattering simulations started from two-planet systems, but recent simulations are all done using three planets or more (e.g., Beaugé & Nesvorný 2012).

However, since the pioneering work, Podsiadlowski et al. (2010), carried out two-planet simulations as one set of runs, we also performed the same two-planet simulations to confirm that our simulation code reproduces Podsiadlowski et al. (2010)'s results. For each semimajor axis, we carried out 20–30 runs. We found that in the two-planet cases, the formation rates of binary planets are 1 out of 20 at 0.2 AU, 5 out of 20 at 1.0 AU, and 11 out of 30 at 5.0 AU, respectively. The results are consistent with Podsiadlowski et al.'s (2010), although the numbers of our runs are smaller than theirs.

Our main results in this paper are obtained by three-planet simulations. We use initial conditions as follows: the semimajor axes of three planets are 1.0a₁, 1.45a₁, and 1.9a₁, and we test four different initial semimajor axes of the innermost planet, a₁ = 1, 3, 5, and 10 AU. Initial orbital eccentricities are zero for all the planets, but we set small orbital inclinations as I₁ = 05, I₂ = 10, and I₃ = 15 to ensure three-dimensional motions, following Marzari & Weidenschiling (2002). We set the other angle variables (the longitudes of ascending node, those of pericenter, and the mean longitude) randomly. For different runs with the same a₁ and planetary radii (R_j), we use different seeds for the random number generation. The mass and radius of the central star are 1 M_☉ and 1 R_☉, and those of three planets are 1M_J and 2R_J, respectively. Since the orbital crossing and formation of binary planets may occur at the timings before the envelopes of gas giants fully contract, we use relatively large R_i. We also performed an additional set of calculations with a₁ = 0.5 AU and R_j = 1R_J. Unless we particularly note that R_j = 1R_J, we use R_j = 2R_J as a nominal parameter (also see Table 1).

For each a₁, we carried out 100 runs to follow the orbital evolution over 10 Myr with a fourth-order Hermite scheme. We stopped calculations when a pair of planets collide (r_ij ⩽ R_i + R_j) or when they form a binary system and the binary eccentricity e_bi < 0.01). We regard that a planet is ejected from the system when the instantaneous distance of the planet from the star becomes more than 10,000 AU and the planet's eccentricity exceeds unity. We also check the collision between the central star and a planet, but when a planet approaches the central star, it tends to become a hot Jupiter due to the tidal interaction with the central star rather than a collision. We neglect spins of the planets and the central star in the N-body simulations for simplicity. (The spin evolution is calculated in the long-term evolution due to quasi-static tide shown in Section 3.4.)

Figure 1 shows an example of orbital evolution to form binary planets. The left panel shows the time evolution of the semimajor axes of the three planets. Two thin lines represent the pericenter distance a_i(1 − e_i) and apocenter distance a_i(1 + e_i). Orbital crossing begins at t ∼ 23, 000 yr. Immediately after that, planets 2 (dashed light green line) and 3 (dash-dotted blue line) form a binary.

Zoom In Zoom Out Reset image size

The right panel represents the path of approaching planets in the local Hill coordinate. The inner planet (planet 2) is set at the origin. In this panel, the X-axis lies along the line from the central star to the inner planet, and the Y-axis is perpendicular to X. The axes are normalized by Hill radius, defined by

$$\begin{eqnarray} &&r_{\rm Hill} = \left(\frac{M_i+M_j}{3M_\ast }\right)^{1/3}\frac{M_ia_i+M_ja_j}{M_i+M_j}. \end{eqnarray} \tag{ 3 }$$

The outer planet (denoted by the solid red line) enters the Hill sphere (denoted by the dash-dotted green line) at t ∼ 23, 000 yr. After that, two planets undergo repeated close encounters to suffer tidal interactions (dotted black line indicates the region within which tidal interaction is taken into account) for ∼300 yr, and they form a gravitationally bound pair (binary). We found that when the planet enters the Hill sphere from the upper left (lower left) direction in this plot, the planets tend to form a prograde (retrograde) binary.

We present the distribution of the stellarcentric semimajor axes and the eccentricities of the binary's barycenters (a_G and e_G) in Figure 2. Four different symbols represent the different initial semimajor axes a₁ = 1 AU (plus), 3 AU (cross), 5 AU (asterisk), and 10 AU (square). This figure shows that the binary planets are formed near their initial orbits. This is because the binary planets are formed in the early stage of orbital instability before they are significantly diffused by many scatterings. The (stellarcentric) eccentricities of the barycenters of the binary pairs are distributed in the range of 0.01–0.6 with the mean value ∼0.15. The mean value corresponds to the eccentricity required for close encounters from initial orbital separation of planets ∼4r_Hill. Since the resultant eccentricities e_G depend on the phase angle of the encounters, they are distributed in a broad range.

Zoom In Zoom Out Reset image size

Figure 3 shows the stellarcentric orbital inclination of formed binary planets, $I_{\rm bi}=\arccos (h_{ij,Z}/h_{ij})$, where h_{ij, Z} is the component of h_ij = r_ij × v_ij that is perpendicular to the orbital plane, and v_ij is the relative velocity between planets i and j. Because we do not see any significant semimajor-axis dependence of I_bi, we superpose the results from different initial a₁. If tidal trapping occurs isotropically, I_bi becomes a sine distribution. The K-S test suggests that this distribution is not the isotropic distribution. The distribution is slightly skewed to small I_bi, which might reflect the fact that the trapping occurs in an early phase before significant orbital excitations. However, retrograde binary planets are formed in a non-negligible fraction of runs (18/40). Note that since we do not take into account the spins of the planets, the retrograde orbits do not necessarily imply the tidally unstable configuration.

Zoom In Zoom Out Reset image size

Figure 4 shows the final semimajor axis of the binary orbits, a_bi. Because we plot the values after the binary eccentricities have been significantly damped (e_bi < 0.01), a_bi is equivalent to the binary separation. Here we also superpose the results from different initial stellarcentric a₁. The distribution of a_bi is peaked at 2R_tot − 4R_tot, where R_tot = R_i + R_j. This is about twice as large as the pericenter distances just after binary capture (see Section 3.2). The factor of two is attributed to angular momentum conservation during the tidal circularization as below.

Zoom In Zoom Out Reset image size

Because tidal force is a strong function of a separation distance between two planets, the two planets just after the trapping usually have a highly eccentric binary orbit with q_{bi, 0} ∼ R_tot − 2R_tot. The initial angular momentum of the binary orbit is $h_{\rm bi,0} \simeq [a_{\rm bi,0}(1-e_{\rm bi,0}^2)]^{1/2}\simeq (2q_{\rm bi,0})^{1/2}$, while the final angular momentum is $h_{\rm bi} \simeq a_{\rm bi}^{1/2}$. The angular momentum conservation indicates that a_bi ∼ 2q_{bi, 0}. This explains the peak at 2 ≲ a_bi/R_tot ≲ 4 and deficit of binary planets at 1 ≲ a_bi/R_tot ≲ 2 in Figure 4. We found that the peaked value of a_bi is much smaller than the Hill radius ∼40R_tot at 1 AU, and the distribution of a_bi does not depend on a₁. This indicates that the binary is hardly affected by interactions with the star or a third planet outside the Hill sphere, unless a_G is very small.

In this subsection we show the process of tidal trapping in more detail. We assumed that the relative velocity of the planets is impulsively decreased at their closest approach. We expect that the binary planets are formed when the relative velocity (v_ij) of two planets immediately after the impulsive tidal dissipation is smaller than their escape velocity,

$$\begin{eqnarray} &&v_{ij} < v_{\rm {esc}}=\sqrt{\frac{2G(M_i+M_j)}{r_{ij}}}. \end{eqnarray} \tag{ 4 }$$

In Figure 5, we present the relative distances (q_ij) and the relative velocities (v_peri) after the tidal dissipation of two planets at each pericenter passage, in the case of a₁ = 0.5 AU with R_{i, j} = 1R_J. The velocities are normalized by the relative velocity of the contact binary, v_{ij, contact} = [G(M_i + M_j)/R_tot]^1/2, and the pericenter distances are normalized by R_tot. When q_ij/R_tot becomes less than unity, two planets collide. The solid red line represents v_peri = v_esc. The passages below this line correspond to those of bound orbits. The dash-dotted (light blue) lines represent the binary's circular orbital velocity v_{K, bi} = [G(M_i + M_j)/r_ij]^1/2. When v_peri reaches this line, their tidal circularization has been completed. Open blue circles, magenta squares, and green crosses represent a series of the pericenter passages of binary, escaping, and colliding pairs, respectively. At points along the line of v_peri = v_esc, a binary (a gravitationally bound pair) is first formed. As the binary bodies repeat close approaches, v_peri is decreased by tidal dissipation and q_ij is increased by the angular momentum conservation, which are represented by a chain of blue open circles from the upper-left to lower-right direction in the figure. We found that the first tidal captures of binary planets occur when q_ij ∼ R_tot − 2R_tot, that is, q_ij is small enough for the tidal force to be strong enough but larger than R_tot to avoid a collision. In all of our simulations, no tidal capture was found from non-bound orbits with q_ij > 3R_tot (the blue open circles outside of 4R_tot in the figure are wandering passages during the circularization that start from q_ij ∼ R_tot − 2R_tot).

Zoom In Zoom Out Reset image size

As already mentioned, we included the tidal dissipation when the closest approach occurs within 10R_tot. As a stellarcentric distance of binary planets increases, the Hill radius becomes larger, but we did not change the threshold distance of 10R_tot. We carried out extra calculations at a₁ = 5 AU with the threshold distance of 50R_tot, to check its effect (set 5a in Table 1). The formation rate of binary planets is 8% in these calculations, and binary formation occurred at q_ij > 10R_tot in only 1 of 100 runs. So, the results hardly change even if the tidal force is incorporated from more distant encounters than q_ij = 10R_tot.

Note that tidal destruction of planets hardly occurs in this capture process. The relative velocity at grazing approach (q_ij ∼ R_tot − 2R_tot) that causes tidal capture is $v \simeq [ v_{\rm gr}^2 + (e v_{\rm Kep})^2]^{1/2}$, where ev_Kep is relative velocity when the planets are sufficiently separated and v_gr is a contribution by the gravitational acceleration between interacting planets, which is given by ∼(0.5–1)v_esc for q_ij ∼ 2R_tot − 1R_tot. For nominal parameters, 1M_J and 2R_J, the surface escape velocity is v_esc ≃ 44 km s⁻¹. The Keplerian velocity is v_Kep ≃ 30(a/1AU)^−1/2 km s⁻¹. Since typical stellarcentric eccentricity is e ≲ 0.3 at the capture (Figure 4), v_gr is dominated and v ∼ (0.5–1)v_esc for a ≳ 0.3 AU, where the orbits of binary planets are stable on timescales of ∼10 Gyr (see Section 3.4). For collision cases, smoothed particle hydrodynamics (SPH) simulation (e.g., Ikoma et al. 2006) showed that significant envelope loss occurs only for v ≳ 2v_esc. Therefore, the tidal destruction is unlikely.

The tidal dissipation could inflate the planetary envelope, which accelerates tidal circularization. However, while this may change total circularization timescale, it would not significantly change binary orbital separations after the circularization because the separations are regulated by q_ij at the trapping when the inflation has not been caused.

We summarize the results of four sets of 100 runs with initial stellarcentric semimajor axes a₁ = 1, 3, 5, and 10 AU (set 1, 3, 5, and 10) for R_i = 2R_J in Figure 6 and Table 1. The ejection rate increases and the collision one decreases as stellarcentric semimajor axis increases, because R_i/r_Hill decreases with the semimajor axis. The important result is that the formation rate of the binary planets is ∼10%, almost independent of the semimajor axis (in other words, almost independent of the value of R_i/r_Hill) as long as a₁ = 1–10 AU.

Zoom In Zoom Out Reset image size

Compared with the results of two-planet systems, the binary formation rate is lower, because the initial conditions of two-planet simulations cause tidal capture before stellarcentric eccentricity is excited. In the three-planet cases, the orbital behavior leading to tidal capture is much more complicated than in the two-planet systems, and eccentricity is excited enough before the capture that the frequency of close encounters with relatively low relative velocity is diminished. Nevertheless, the formation probability of the binary planets in three-planet systems is still as large as ∼10%.

In order to check how the formation rate of the binary planets depends on other parameters, we carried out four additional calculations (sets 0.5a, 3a, 3b, and 5a in Table 1). In these additional calculations, we changed one of the parameters (planetary radius, initial planet–planet distances, or the tidal limit) and kept the other parameters, including planetary masses, the same. In set 0.5a, we used half-sized planets (R_i = 1R_J), keeping the planetary masses the same. Since a₁ = 0.5 AU, this set has the same value of R_i/r_Hill as set 1. As a result, the ejection/collision ratio is similar, as shown in Table 1. Although the effect of tidal dissipation is weaker than in set 1, we found a similar (9/100) binary formation rate. The binary planets are formed through grazing encounters. When we use the smaller planetary radius, while the tidal dissipation becomes weaker, some fraction of the close encounters that lead to collisions for R_i = 2R_J result in binary formation. These two effects tend to cancel, so that we did not find a large difference between set 1 and set 05a in our small number of statistics.

In set 3a (a₁ = 3 AU), while we keep R₁ = 2R_J, R₂ and R₃ were reduced to 1R_J. The distribution of semimajor axes of half-sized binary planets is peaked at a_bi ∼ 2.5R_tot, which is similar to Figure 4. However, the formation rate of the binary planets is increased to be ∼17%. The increase comes from formation of 2R_J − 1R_J binaries (11% of the 17%). This may be because the 2R_J can suffer stronger tidal forces without collisions.

In set 3b (a₁ = 3 AU), the planet–planet initial orbital separation is set to be ∼6r_Hill at a₁ = 3 AU, which is 1.5 times larger than the standard cases. Table 1 shows that the formation probability of binary planets is slightly lower (∼6%). This is because the planets have to build up larger eccentricities for close encounters to occur, and accordingly the relative velocities are higher.

In set 5a (a₁ = 5 AU), we increased the threshold distance inside which the tidal force is incorporated, by five times at a₁ = 5 AU. The formation probability of the binary planets does not change, as we mentioned in Section 3.2.

In previous sections, we showed that the binary gas giant planets are formed through the scattering of three planets orbiting around a central star and planet–planet dynamical tide with non-negligible probability (∼10%). After the orbital circularization, the dynamical tide diminishes and the orbital evolution of the binary system is governed by long-term quasi-static tides. In this section we calculate long-term tidal evolution of the formed binary systems through the planet–planet and planet–star quasi-static tidal interactions, instead of the dynamical tide that has been considered in previous sections. The evolution on the main-sequence lifetime of solar-type stars (∼10 Gyr) is very important for detectability of binary planets in extrasolar planetary systems by, e.g., transit observations (see Paper II).

We calculate this tidal evolution process, basically following Sasaki et al. (2012). Sasaki et al. (2012) neglected the tidal interactions between the central star and the binary companion and their spins. Because we consider binary planets with relatively small stellarcentric radius where the detectability by transit observation is not too small, we include all the tidal interactions and the spin angular momenta in the planet-companion-star systems.

Although we consider a pair of comparable planets, we call one planet "primary" and the other "companion" for convenience. We use the subscripts "*" for the central star, "pp" for the primary planet, and "cp" for the companion planet. In this subsection, we use the following assumptions.

• 1.  
  The total angular momentum, that is, the sum of the angular momenta of the stellar and the planets' spins and those of the binary and stellarcentric orbits, is conserved.
• 2.  
  All orbits in the system are circular and coplanar.
• 3.  
  All the spin angular momentum vectors are parallel to their orbital angular momentum vectors.
• 4.  
  The separation of two objects can be changed only by the tidal interactions between them and those with the host star. It is not affected by the other objects.
• 5.  
  The masses of the planets are negligible compared with that of the central star, i.e., M_* ≫ M_pp, M_cp.

The parameters we adopt in the calculations are shown in Table 2. We set t = 0 at the completion of tidal circularization of the binary. We assume that the initial planetary and stellar spin periods are 10 hr and 30 days, respectively. The binary orbital period is calculated by a_bi; it is about three days for a_bi = 2.5R_tot. The stellarcentric orbital period is calculated by its semimajor axis of the binary barycenter; it is one yr for a_G = 1 AU. Thus, it is reasonable to assume that Ω_pp = Ω_cp > n_bi > n_G at t = 0, where Ω_pp and Ω_cp are the planetary spin frequencies and n_bi and n_G are mean motions of binary and stellarcentric orbits, which are given by

$$\begin{eqnarray} n_{\rm bi} &=& \sqrt{\frac{G(M_{\rm pp}+M_{\rm cp})}{a_{\rm bi}^3}},\qquad \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} n_{\rm G} &=& \sqrt{\frac{G(M_\ast + M_{\rm pp}+M_{\rm cp})}{a_{\rm G}^3}}. \end{eqnarray} \tag{ 6 }$$

In the Appendix, the quasi-static tidal torque equations are integrated to give n_G(t), n_bi(t), Ω_*(t), Ω_pp(t) (= Ω_cp(t)) as explicit functions of time t:

$$\begin{eqnarray} n_{\rm G}(t) &=& \left[\frac{39}{2}\frac{1}{G(M_{\rm pp}+M_{\rm cp})(GM_\ast)^{2/3}}\left(\frac{k_{\rm 2pp}}{Q_{\rm pp}}R_{\rm pp}^5+\frac{k_{\rm 2cp}}{Q_{\rm cp}}R_{\rm cp}^5\right.\right. \nonumber\\ &&\left.\left. +\frac{M_{\rm pp}^2+M_{\rm cp}^2}{M_\ast ^2}\frac{k_{2\ast }}{Q_\ast }R_\ast ^5\right)t+n_{\rm G,0}^{-13/3}\right]^{-3/13}, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} n_{\rm bi}(t) &=& \left[\frac{39}{2}\frac{1}{\lbrace G(M_{\rm pp}+M_{\rm cp})\rbrace ^{5/3}}\frac{1}{M_{\rm pp}M_{\rm cp}} \left(\frac{k_{\rm 2pp}}{Q_{\rm pp}}R_{\rm pp}^5M_{\rm cp}^2\right.\right. \nonumber\\ &&\left.\left. +\frac{k_{\rm 2cp}}{Q_{\rm cp}}R_{\rm cp}^5M_{\rm pp}^2\right)t+n_{\rm bi,0}^{-13/3}\right]^{-3/13}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&{\Omega _\ast (t) = -\displaystyle\frac{k_{2\ast }R_\ast ^3}{\alpha _\ast Q_\ast } \displaystyle\frac{M_{\rm pp}+M_{\rm cp}}{M_\ast } \displaystyle\frac{(GM_{\rm pp})^2+(GM_{\rm cp})^2}{(GM_\ast)^{4/3}} }\nonumber \\ &&{\times \displaystyle\frac{ \left\lbrace n_{\rm G}^{-1/3}(t)-n_{\rm G,0}^{-1/3} \right\rbrace }{\left[k_{\rm 2pp}R_{\rm pp}^5/Q_{\rm pp}+k_{\rm 2cp}R_{\rm cp}^5/Q_{\rm cp}+\lbrace (M_{\rm pp}/M_\ast)^2+(M_{\rm cp}/M_\ast)^2\rbrace k_{2\ast }R_\ast ^5/Q_\ast \right]} }\nonumber\\ &&{ +\, \Omega _{\ast,0},} \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&{\Omega _{\rm pp}(t) =\Omega _{\rm cp}(t) = -\displaystyle\frac{k_{\rm 2pp}R_{\rm pp}^3}{\alpha _{\rm pp} Q_{\rm pp} M_{\rm pp}} }\nonumber \\ &&{\times\! \left[\!\displaystyle\frac{(GM_{\rm pp})M_{\rm cp}^3}{\lbrace G(M_{\rm pp}+M_{\rm cp})\rbrace ^{1/3}}\displaystyle\frac{n_{\rm bi}^{-1/3}(t)-n_{\rm bi,0}^{-1/3}}{\left(M_{\rm cp}^2k_{\rm 2pp}R_{\rm pp}^5/Q_{\rm pp}+M_{\rm pp}^2k_{\rm 2cp}R_{\rm cp}^5/Q_{\rm cp} \right)}\right. }\nonumber \\ &&{+\! \left.\displaystyle\frac{(M_{\rm pp}+M_{\rm cp})(GM_\ast)^{2/3} \left\lbrace n_{\rm G}^{-1/3}(t)-n_{\rm G,0}^{-1/3} \right\rbrace }{\left[ k_{\rm 2pp}R_{\rm pp}^5/Q_{\rm pp}\,{+}\,k_{\rm 2cp}R_{\rm cp}^5/Q_{\rm cp}\,{+}\,\right\lbrace (M_{\rm pp}/M_\ast)^2\,{+}\,(M_{\rm cp}/M_\ast)^2 \left\rbrace k_{2\ast }R_\ast ^5/Q_\ast \right]} \!\right]}\nonumber\\ &&{+\, \Omega _{\rm pp,0},} \end{eqnarray} \tag{ 10 }$$

where the subscripts, "0" represent the values at t = 0, k₂ are Love numbers, and Q are tidal dissipation functions. We adopt the estimate for the current solar and Jovian values of k₂ and Q as parameter values of the host star and the planets, which are shown in Table 2.

Using these equations, we show in Figure 7 the orbital evolutions of binaries at a_G = 0.2 AU (upper panel), 0.3 AU (middle panel), and 0.4 AU (lower panel). The dash-dotted green and dashed blue lines represent n_bi and Ω_pp. The vertical red lines represent the lifetime of the main-sequence phase of solar-type stars (t = 10¹⁰ yr), the upper horizontal black lines are the binary orbital angular velocities for a contact binary ([G(M_pp + M_cp)/(R_pp + R_cp)³]^1/2), and the lower one represents n_crit that is determined by the critical semimajor axis of the binary orbit, below which the binary separation is so large that the orbit is destabilized by stellar gravitational force (Equation (A17)). When the dash-dotted green lines stay in the region surrounded by two horizontal black lines (n_crit < n_bi < [G(M_pp + M_cp)/(R_pp + R_cp)³]^1/2), the binary planets are stable. The two planets collide for n_bi ⩾ [G(M_pp + M_cp)/(R_pp + R_cp)³]^1/2, and they escape from each other for n_bi ⩽ n_crit.

Zoom In Zoom Out Reset image size

The spins of the individual planets are slowed down by the planet–planet tidal interaction, and accordingly the binary orbital angular momentum is increased. As a result, both Ω_pp(t) and n_bi(t) decrease. Since Ω_pp(t)'s deceleration is faster than that of n_bi(t), Ω_pp(t) catches up with n_bi(t) to establish a synchronous state. After that, the binary planets keep the synchronized state while n_G(t) is decelerated by the tidal torque from the star. In this synchronous state with Ω_pp = Ω_cp = n_bi, the total angular momentum is given by

$$\begin{eqnarray} L &=& L_{\rm bi}+L_{\rm G}+I_\ast \Omega _\ast +I_{\rm pp}\Omega _{\rm pp}+I_{\rm cp}\Omega _{\rm cp},\qquad\qquad\ \ \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} & = & \frac{M_{\rm cp}(GM_{\rm pp})}{\lbrace n_{\rm {bi}}G(M_{\rm pp}+M_{\rm cp})\rbrace ^{1/3}}+ \frac{(M_{\rm pp}+M_{\rm cp})(GM_\ast)^{2/3}}{n_{\rm {G}}^{1/3}} \nonumber \\ && + \alpha _\ast R_\ast ^2M_\ast \Omega _\ast + (\alpha _{\rm pp} R_{\rm pp}^2M_{\rm pp}+ \alpha _{\rm cp} R_{\rm cp}^2M_{\rm cp})n_{\rm bi}, \end{eqnarray} \tag{ 12 }$$

where α are moment of inertia ratios, the values of which we adopted are shown in Table 2. Since L_G (∝n_G(t)^−1/3) increases continuously, Ω_pp(t) and n_bi(t) increase to keep the synchronized state from the total angular momentum conservation (L_bi and I_*Ω_* decrease, but L_G, I_ppΩ_pp, and I_cpΩ_cp increase). Because n_bi(t) keeps increasing, a_bi keeps decreasing, and eventually the binary planets collide with each other.

The synchronous state among Ω_pp, Ω_cp, and n_bi is established in about 0.3 Myr in the cases of a_G = 0.3 AU and 0.4 AU. However, in the case of a_G = 0.2 AU, it becomes n_bi < n_crit, and two planets escape from each other in about 0.1 Myr before the synchronous state is established. The lifetime of the binary planets is longer than the main-sequence phase of solar-type stars (∼10 Gyr) for a_G = 0.4 AU. For a_G = 0.3 AU, the lifetime is about 7 Gyr. However, since the tidal torque is proportional to fifth order of the planetary radius (Equation (A1)), the lifetime for R_pp, R_cp = 1R_J with the same M_pp, M_cp = 1M_J is lengthened by about 10 times that in this plot. Because the gas envelope may fully contract in 0.1 Gyr, R_pp, R_cp = 1R_J may be more appropriate than 2R_J for the estimate of the binary lifetime. Therefore, the binary would survive also for a_G = 0.3 AU.

We adopted 10 hr as the initial planetary spin periods, 2π/Ω_{pp, 0} and 2π/Ω_{cp, 0}. For smaller Ω_{pp, 0} and Ω_{cp, 0}, the lifetime of the binary is shorter, but the lifetime is still ∼20 Gyr for a_G = 0.4 AU even if Ω_{pp, 0} = n_{bi, 0}. On the other hand, when the planetary spin period is ≲ 4 hr, the binary is separated more than a_crit before reaching the synchronous state.

The Q value may include large uncertainty. However, the binary stability condition for a_G comes from the tidal evolution timescale due to stellar tide. Since the timescale is proportional to $Q_{\rm pp} a_{\rm G}^{6.5}$ (Equation (A20)), the condition for a_G is not severely affected by the uncertainty in the Q value.

Note also that we neglect the spins in planet–planet scattering calculations. The planetary spin axis can be reversed in scattering. The binary planets approach and impact each other quickly by tidal evolution when binary orbit and spin are in retrograde. When the stellar centric orbit and binary orbit are in retrograde, the first, fourth, and fifth terms in Equation (11) change their signs and the binary is separated away. That means even if a_G ≳ 0.3 AU, a part of the binary planets cannot survive.