Formation of Wide-orbit Gas Giants Near the Stability Limit in Multi-stellar Systems

Most of the stars in our Galaxy are members of binary or triple star systems. Many exoplanets have been discovered in circumstellar orbits in multiple star systems by transit and radial velocity surveys. Although the circumstellar orbits of planets can be destabilized by secular perturbations from a binary companion, they are stable if their orbital radii are less than a critical value. Holman & Wiegert (1999) derived a fitting formula for the critical separation for stability as

$$\begin{eqnarray}&&{a}_{{\rm{c}}}\simeq (0.46-0.38\mu -0.63{e}_{{\rm{b}}}+0.59\mu {e}_{{\rm{b}}}){a}_{{\rm{b}}},\end{eqnarray} \tag{ 1 }$$

where μ is the binary companion mass scaled by the sum of the host and companion stars, e_b and a_b are the eccentricity and semimajor axis, respectively, of the binary companion orbit, and circular orbits are assumed for the planets. For an equal-mass binary pair (μ = 0.5) with a circular binary orbit, a_c ≃ 0.27a_b.

Wide-orbit extrasolar gaseous giant planets in nearly circular orbits have been detected by direct imaging observations for several systems (e.g., Kalas et al. 2008; Marois et al. 2008; Kuzuhara et al. 2013). One of the latest announcements is a (4 ± 1) M_Jup planet in the HD131399 triple star system (Wagner et al. 2016). The primary star, HD131399A, is an A-type star with mass ≃1.82 M_⊙, which the discovered planet (HD131399Ab) orbits. The close binary of a G-type star with mass ≃0.96 M_⊙ (HD131399B) and a K-type star with mass ≃0.6 M_⊙ (HD131399C) is orbiting HD131399A with semimajor axis of ∼350 au and eccentricity of ∼0.13 (Wagner et al. 2016). With these values, Equation (1) shows a_c ≃ 0.24 a_b ≃ 84 au. The estimated semimajor axis and eccentricity of the planet HD131399Ab are ${82}_{-27}^{+23}$ au and 0.35 ± 0.25 (Wagner et al. 2016). Although Equation (1) is for planets in circular orbits, the estimated planetary semimajor axis is close to a_c and the planet's orbit would be marginally stable or possibly unstable on long timescales (Wagner et al. 2016; Veras et al. 2017). Although the planet could be a background star (Nielsen et al. 2017), if it is a planet, how to form such a marginally stable world in a multiple planetary system is a challenge for planet formation theory.

Because in situ formation of gas giant planets is difficult at ∼80 au with the conventional core accretion scenario (e.g., Ida & Lin 2004), Wagner et al. (2016) discussed the following three possibilities: (i) planet–planet scattering in the star A system, (ii) formation of a circumbinary planet around the B/C pair followed by transfer to the A system, and (iii) formation of the planet around either of the stars and evolution in the course of stellar orbital evolution of the triple stellar system before settling down to the current stellar configuration. For model (i), an additional unseen massive planet is required in the inner region. Because the eccentricity after the scattering must be close to unity, an eccentricity damping mechanism for the scattered planet is also required. Although models (ii) and (iii) are not ruled out, it is not clear what evolutionary path leads to the current planetary and stellar orbital configurations.

Here, we are focused on the scattering model and discuss how a gas planet with baffling features (a wide-separation moderate-eccentric orbit close to the stability limit by perturbations from a companion(s)) such as HD131399Ab can be reproduced. Planet scattering by a more massive planet can significantly increase the apoastron distance, but not the periastron distance. To explain a wide-orbit gas giant with moderate orbital eccentricity, a mechanism to lift the periastron distance, equivalently, to damp the eccentricity without significant semimajor axis damping is required.

The Lidov–Kozai mechanism (Kozai 1962; Lidov 1962) is one of the possible eccentricity variation mechanisms. The secular perturbation from the binary companions orbiting around HD131399A can induce large amplitude oscillations of the eccentricity e and inclination i of HD131399Ab without changing the semimajor axis due to conservation of the vertical component of the orbital angular momentum ($\propto \sqrt{1-{e}^{2}}\cos i$) around A, if the initial argument of periastron is proper. In this case, however, the proper range of the initial argument of periastron is narrow. Furthermore, even if the narrow range of initial conditions is satisfied, the periastron distance repeatedly becomes smaller according to the oscillation, and HD131399Ab's orbit can also be repeatedly perturbed by the massive planet in the inner region that scattered HD131399Ab. It is not clear if HD131399Ab's orbit can be stable during the age of the host star.

Another possibility is dynamical friction from the protoplanetary disk gas if exists after the scattering. Bromley & Kenyon (2014) investigated the orbital evolution of eccentric planets by dynamical friction. They found that the periastron distance is increased by a factor of at most four times in the case of uniform disk gas depletion. In the case of inside-out depletion, dynamical friction from the outer part of the disk contributes to significant lift-up of the periastron distance, but it is unlikely that the inside-out disk gas depletion continues up to the ∼50 au that is required to reproduce the orbit of HD131399Ab.

The model proposed by Kikuchi et al. (2014) is also based on planet–planet scattering, but considers scattering of a solid core and the eccentricity damping during gas accretion onto the core to form a gas giant. The gas accretion is most efficient when the core is near the apoastron in the highly eccentric orbit after the scattering. Even if the disk gas surface density of $\propto {r}^{-3/2}-{r}^{-1}$, the periastron distance is significantly increased. This model does not require a gas giant in an inner region more massive than the wide-orbit gas giant in the final state. This model is a promising option to reproduce the current orbital configuration of HD131399Ab.

The model also naturally explains why the orbit is close to the stability limit. Kikuchi et al. (2014) considered the formation of a wide-orbit gas giant around a single star and found that the final orbit is regulated by the radius of the disk's outer edge, because gas accretion onto the planet does not occur beyond the disk's outer edge. When a stellar companion exists, the circumstellar disk is truncated at ∼1/3 of the binary separation by the perturbations from the companion with e_b ∼ 0.1, and the truncation radius becomes larger for a smaller-mass companion and for smaller e_b (e.g., Artymowicz & Lubow 1994). If the model by Kikuchi et al. (2014) is adopted, the final semimajor axis of the planet settles down to the location just inside the disk truncation radius, which is close to the stability limit. This result is not specific to the HD131399 system, but can be applied generally to circumstellar gas giants in multiple stellar systems. We use the HD131399 system as a reference system. Note that circumstellar gas giant planets, γ Cephei Ab, HD196885Ab, and HD41004Ab, in closer binary systems (a_b ∼ 20 au) also have orbits close to the stability limit (e.g., Thebault & Haghighipour 2015). Because these gas giants are at ∼2 au, they can be in principle formed by the standard core accretion model. However, it is not clear if planet formation proceeds in such highly perturbed environments (Thebault & Haghighipour 2015). These planets could also be formed by the model here, if another planet with ≳20–30 Earth masses existed at ≲1 au, as discussed later.

We briefly summarize the model of Kikuchi et al. (2014) in Section 2.1 and describe the assumptions and initial conditions for the discussions on HD131399 system in Section 2.2. We present the results of numerical simulations in Section 3. Section 4 is devoted to the summary and discussion.

2.1. The Model by Kikuchi et al. (2014)

The scenario proposed by Kikuchi et al. (2014) is based on the conventional core accretion model as follows: (i) an icy core (or a core with some amount of gaseous envelope) accretes from planetesimals in inner regions at semimajor axis ≲30 au, (ii) it is scattered outward by a nearby gas giant or more massive core to acquire a highly eccentric orbit with a periastron distance close to the original semimajor axis, (iii) its orbit is circularized with significant lift-up of the periastron distance through accretion of local protoplanetary disk gas along the eccentric orbit, and (iv) through the local gas accretion, the planet becomes a gas giant. For a given initial high eccentricity and semimajor axis, Kikuchi et al. (2014)investigated the details of the process in step (iii), assuming that the gas disk motion is Keplerian, the motions of the planet and the gas disk are coplanar (they found that if the inclination is smaller than 30°, the final values of eccentricity and semimajor axis change by less than 5%), the gas accretion rate onto the planet from local disk regions is regulated by its envelope contraction rather than supply from the disk and accordingly, it is independent of a radial distribution of disk gas surface density except that it vanishes beyond the diskʼs outer edge.

Here, we summarize their results (for detailed derivations of formulas, refer to Kikuchi et al. 2014). The changes of the energy and angular momentum are derived analytically, although it must be solved numerically. The damping rates of eccentricity e and semimajor axis a as a function of the planetary mass M are given by

$$\begin{eqnarray}&&\displaystyle \frac{{de}}{d\mathrm{log}M}=-\displaystyle \frac{1-{e}^{2}}{e}\,\left({f}_{{\ell }}(e,{u}_{{\rm{d}}})+\displaystyle \frac{1}{2}{f}_{\epsilon }(e,{u}_{{\rm{d}}})\right),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{log}a}{d\mathrm{log}M}=-{f}_{\epsilon }(e,{u}_{{\rm{d}}}),\end{eqnarray} \tag{ 3 }$$

where u_d is the maximum eccentric anomaly ($0\lt {u}_{{\rm{d}}}\lt \pi$) within the gas disk of the outer edge at r_d, which is given by

$$\begin{eqnarray}{u}_{{\rm{d}}}\equiv \left\{\begin{array}{ll}{\cos }^{-1}\left[\tfrac{1}{e}\left(1-\tfrac{{r}_{{\rm{d}}}}{a}\right)\right] & \ \ {\rm{for}}\ \ [Q\gt {r}_{{\rm{d}}}]\\ \pi & \ \ {\rm{for}}\ \ [Q\lt {r}_{{\rm{d}}}]\end{array}\right.,\end{eqnarray} \tag{ 4 }$$

where Q is the apoastron distance of the planet. The functions f_ℓ and f are the change rates of the specific angular momentum and orbital energy over one orbital period expressed as

$$\begin{eqnarray}&&{f}_{{\ell }}={\int }_{-{t}_{{\rm{d}}}/2}^{{t}_{{\rm{d}}}/2}s\,\left(\displaystyle \frac{{l}_{\mathrm{gas}}}{{l}_{{\rm{p}}}}-1\right){dt},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{f}_{\epsilon }={\int }_{-{t}_{{\rm{d}}}/2}^{{t}_{{\rm{d}}}/2}s\left(\displaystyle \frac{{\epsilon }_{\mathrm{gas}}}{{l}_{{\rm{p}}}}-\displaystyle \frac{{\epsilon }_{\mathrm{coll}}}{{l}_{{\rm{p}}}}-1\right){dt},\end{eqnarray} \tag{ 6 }$$

where t_d is a duration at $r\lt {r}_{{\rm{d}}}$, r is the instantaneous distance of the planet from the central star, l_p and _p are the specific angular momentum and orbital energy of the planet, l_gas and _gas are those of accreting gas at instantaneous locations along the orbit, and _coll is energy dissipation by collision between the planet and accreting gas. The factor s is given as follows. The relative velocity between the planet and the unperturbed gas flow (circular Keplerian flow), is readily derived by

$$\begin{eqnarray}&&{v}_{\mathrm{rel}}{(r)}^{2}=\displaystyle \frac{{{GM}}_{* }}{r}\,\left(3-\displaystyle \frac{r}{a}-2\sqrt{\displaystyle \frac{a}{r}(1-{e}^{2})}\right),\end{eqnarray} \tag{ 7 }$$

where M_* is the host star mass and a and e are the semimajor axis and the eccentricity of the instantaneous planetary orbits. However, because we are concerned with damping of initially highly eccentric orbits, the reduction of the relative velocity due to bow shock that occurs in front of the planet must be taken into account. Kikuchi et al. (2014) derived the reduction factor s for post-shock relative velocity (${{sv}}_{\mathrm{rel}}$) with a 1D plane-parallel approximation as

$$\begin{eqnarray}&&s=\displaystyle \frac{(\gamma -1){{ \mathcal M }}^{2}+2}{(\gamma +1){{ \mathcal M }}^{2}}\simeq \displaystyle \frac{1}{4}\left(1+\displaystyle \frac{3}{{{ \mathcal M }}^{2}}\right),\end{eqnarray} \tag{ 8 }$$

where γ is the specific heat ratio (γ = 7/5 for H₂ molecules) and ${ \mathcal M }$ is Mach number of pre-shock,

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{{v}_{\mathrm{rel}}}{{c}_{s}}=30\,{\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{-1/4}{\left(3-\displaystyle \frac{r}{a}-2\sqrt{\displaystyle \frac{a}{r}(1-{e}^{2})}\right)}^{1/2},\end{eqnarray} \tag{ 9 }$$

where we used an optically thin disk temperature, $T=280{(r/1\mathrm{au})}^{-1/2}\,{\rm{K}}$, for evaluation of sound velocity c_s. Note that Equations (2) and (3) are independent of a disk model except the estimation of Mach number.

2.2. Initial Conditions or Settings

We apply Kikuchi et al.'s model to a reference planet, HD131399Ab. We first evaluate the initial orbital elements and mass of a precursor body for HD131399Ab and the gas disk that are substituted into Kikuchi et al.'s model, where "initial" means the moment just after the scattering by an unseen planet. We do not simulate the scattering process itself, but evaluate a reasonable range of the "initial" parameters after the scattering.

The initial periastron distance q_i must be around the unseen massive planet. In the conventional core accretion model, q_i may be ≲40 au. The upper limit of the initial apoastron distance Q_i may be limited by the existence of the binary companions. The best-fit of the orbit of the binary by Wagner et al. (2016) is a_b = 349 ± 28 au and e_b = 0.13 ± 0.05 where a_b and e_b are the semimajor axis and eccentricity of the binary with respect to HD131399A. We assume a_b = 350 au and e_b = 0.13 and consider the cases of Q_i ≲ 300 au. The initial semimajor axis a_i and eccentricity e_i of the planet are calculated by q_i and Q_i as

$$\begin{eqnarray}&&{a}_{{\rm{i}}}=\displaystyle \frac{{Q}_{{\rm{i}}}+{q}_{{\rm{i}}}}{2},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{e}_{{\rm{i}}}=\displaystyle \frac{{Q}_{{\rm{i}}}-{q}_{{\rm{i}}}}{{Q}_{{\rm{i}}}+{q}_{{\rm{i}}}}.\end{eqnarray} \tag{ 11 }$$

The inclination to the gas disk is assumed to be zero for simplicity. The initial mass of the planet (M_i) must be larger than the critical core mass for the onset of gas accretion. We use M_i = 20 M_E and 50 M_E, where M_E is an Earth mass. The final mass of the planet is given by the observational best-fit value, 4 M_J (Wagner et al. 2016). We stop the calculation when the mass of the planet reaches the final mass. In Kikuchi et al.'s model, orbital evolution is described by the planetary mass evolution but not by time evolution. We do not need to assume gas accretion rate onto the planet.

We use 100 au and 150 au for the gas disk radius—because the circumstellar disk is tidally truncated by the perturbations from the binary companions—at ∼1/3 of the distance between the primary star HD131399A and the binary (e.g., Artymowicz & Lubow 1994). The disk truncation radius is comparable to or slightly larger than the long-term orbital stability limit radius. We assume the gas disk exists until the end of the calculation.

Using the initial semimajor axis and eccentricity, we calculate the final semimajor axis and eccentricity for different values of disk size r_d and the initial core mass M_i, integrating Equations (2) and (3) up to M = 4 M_J. Figure 1 shows some examples of orbital evolution as a function of the planetary mass growth. According to the mass growth of the planet, its eccentricity and semimajor axis are damped and periastron q and apoastron Q converge. Periastrons are lifted up quickly and the scattered planet becomes isolated from the perturber, the more massive planet. Apoastrons Q settle down to ${{fr}}_{{\rm{d}}}$ with f ∼ 0.7–1.0. Kikuchi et al. (2014) showed that the final Q is smaller for smaller q_i/r_d and smaller final e. Because the final e is similar in each panel of Figure 1, the difference of the final values of Q reflects the different values of q_i/r_d.

Zoom In Zoom Out Reset image size

The damping of Q is associated with the planetary mass increase. This figure shows that the mass increase by 10 times from the initial value is enough to damp the apoastrons Q to ≲r_d. The mass increase factor to the "final" state is 4 M_J/M_i = 25 or 64. Thus, the planets' final orbits are inevitably close to the disk truncation radius. Orbital eccentricities are damped significantly from the initial values ∼1, but have not been completely damped in this case. The semimajor axes are slightly smaller than Q in the final state, corresponding to the retained e. Because the orbital stability limit is slightly closer to the primary star than the disk outer edge is, the final semimajor axes are comparable to the orbital stability limit. With $Q\sim {{fr}}_{{\rm{d}}}$, the final a would be

$$\begin{eqnarray}&&a=(1-e)Q\sim (1-e){{fr}}_{{\rm{d}}}\sim (1-e)(f/3)\left(\displaystyle \frac{{r}_{{\rm{d}}}}{{a}_{{\rm{b}}}/3}\right){a}_{{\rm{b}}}.\end{eqnarray} \tag{ 12 }$$

For e ∼ 0.35, f ∼ 0.7–1.0, and ${r}_{{\rm{d}}}\sim {a}_{{\rm{b}}}/3$, a ∼ (0.15–0.22) a.

Figure 2 shows the values of the semimajor axis a and eccentricity e of the planet when the planetary mass reaches 4 M_J, on the Q_i–q_i plane with available ranges of 50 au < Q_i < 300 au and q_i < 40 au for disk sizes, r_d = 100 au and 150 au, and the initial planetary masses, M_i = 20 M_E and 50 M_E. The parameter f is controlled by q_i/r_d, so that q_i/r_d and r_d/a_b (the binary separation is fixed to be a_b = 350 au) determine the final a, as shown in Equation (12). Because the mass increase factor determines the final e and the final mass is fixed to 4 M_J, M_i determines the final e in this case. Each square on the map corresponds to one calculation that we have performed. The observationally estimated semimajor axis and eccentricity of HD131399Ab are ${82}_{-27}^{+23}$ au and 0.35 ± 0.25, respectively (Wagner et al. 2016). Different colors for the squares represent how the calculated final values of a and e are consistent with the observational estimates, from yellow, light blue, orange, and black. The yellow squares show the best ones, where 77 au < a < 87 au and 0.3 < e < 0.4. Figure 2 shows that the current orbit of HD131399Ab is reproduced from the initial conditions with q_i ∼ 20–30 au, a relatively large initial mass (M_i = 50 M_E), and relatively large disk truncation radius (r_d = 150 au).

Zoom In Zoom Out Reset image size

If ${Q}_{{\rm{i}}}$ is close to the separation of the binary companion, there is a risk that the planet is again scattered by the binary companion pair until Q is sufficiently damped. For M_i = 50 M_E and r_d = 150 au, the yellow squares exist down to Q_i ∼ 150 au, where the risk may be low. The relatively massive initial mass, M_i = 50 M_E, corresponds to a planet on its way to becoming a gas giant. In this case, a gas giant more massive than Saturn must have existed at 20–30 au to scatter the planet.

We have investigated the formation of circumstellar wide-orbit gas giants in multiple star systems, motivated by the announcement of the discovery of planet HD131399Ab by direct imaging (Wagner et al. 2016). This planet is a circumstellar gas giant in a hierarchical triple star system and has perplexing orbital features. The semimajor axis is about 1/4 of the separation between A and the B/C binary pair, which means that the planet's orbit would be marginally stable or possibly unstable on long timescales. We have applied the formation model for wide-orbit gas giants by Kikuchi et al. (2014) and found that the orbit close to the stability limit is not a fortunate accident, but rather an inevitable result.

The suggested scenario is as follows: (i) a core with or without a gaseous envelope is scattered by a nearby more massive planet at 20–30 au so that it develops a highly eccentric orbit with the periastron at 20–30 au, (ii) the scattered planet spends most of its time near the apoastron on a highly eccentric orbit after the scattering, where the planetesimal accretion rate must be very small, (iii) rapid gas accretion starts because of the low planetesimal accretion rate, and the gas accretion is most efficient near the apoastron, (iv) through gas accretion onto the planet, the orbit is circularized such that the periastron distance is significantly increased and the apoastron distance is adjusted to the tidal truncation radius of the disk by the binary companion (r_d), which is comparable to the orbital stability limit. Therefore, this model is a promising option for reproducing the current orbital configuration of HD131399Ab. The best-fit results are obtained by the initial scattering location q_i ∼ 20–30 au, the initial planet mass M_i ∼ 50 M_E and the disk truncation radius r_d ∼ 150 au.

If the orbital data of HD131399Ab is improved, the predicted initial conditions are better constrained. Even if HD131399Ab turns out be a background star, our model can be generally applied to circumstellar gas giants in multiple stellar systems, and it predicts that the planets would often be located close to the stability limit against the perturbations of a companion(s). The results depend on the scaled parameters, q_i/r_d and r_d/a_b, where a_b is the binary separation and q_i corresponds to the original semimajor axis of the planet before scattering. The parameter r_d/a_b is usually ∼1/3, except for cases with high e_b or a very different mass companion from the host star, and the final Q does not change sensitively for q_i/r_d ≲ 0.1 (Kikuchi et al. 2014).

Other circumstellar gas giant planets close to the stability limit have already been discovered in more compact binary systems (a_b ∼ 20 au): γ Cephei Ab, HD196885Ab, and HD41004Ab, all of which are located at ∼2 au. The parameter r_d/a_b for these systems should be similar to HD131399Ab. Our model can reproduce their orbits close to the stability limit with reasonable ranges of values of q_i and M_i (see Equation (12)), while these planets could also be formed by the standard in situ core accretion model. Note that our model needs another planet with ≳20–30 M_E that scattered the core, which has not been detected in the γ Cephei A, HD196885A, and HD41004A systems. Furthermore, in compact binary systems, there is a greater risk that the scattered core will be perturbed by the binary companion before the scattered orbit shrinks. For very wide binary systems, our model may not work well either. In our model, Q_i must be between r_d ∼ a_b/3 and a_b, while it is likely that q_i ≲ 30 au. From Equation (11), $1\gt {e}_{{\rm{i}}}\gtrsim ({a}_{{\rm{b}}}-90)/({a}_{{\rm{b}}}+90)$. The range of e_i that satisfies this condition becomes narrower as a_b increases. Furthermore, the disk size cannot become infinitely large when a_b is very large. In this case, r_d is independent of the stability limit. Thus, our model is better applied to binary systems with separations of, say, 100–1000 au, while it could also be applied to more compact or more extended binary systems.

It is also pointed out that circumbinary gas giant planets tend to be located near the inner limit of the stable region around binary systems because the curcumbinary disks are truncated by an inner binary pair and planetary migrations stop at the inner edge (Kley & Haghighipour 2014, 2015). Here, we predict that circumplanetary wide-orbit gas giant planets tend to be located near the outer limit of the stable region around the host star. These predictions should be tested by more samples of direct imaging observations for wide-orbit gas giants in multiple stellar systems.

This paper was motivated by a group project at Lorentz center workshop "New Directions in Planet Formation," held in 2016 July, that S.I. participated in. S.I. thanks the team members of the group project, Nader Haghighipour, Judit Szulágyi, Steven Rieder, Katherine Kretke, and Jiwei Xie. We thank Akihiro Kikuchi for providing us the code developed for our previous paper, Kikuchi et al. (2014) and an anonymous referee for useful comments. Data analyses were in part carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was supported by JSPS KAKENHI grant 15H02065.