ORBITAL CIRCULARIZATION OF A PLANET ACCRETING DISK GAS: THE FORMATION OF DISTANT JUPITERS IN CIRCULAR ORBITS BASED ON A CORE ACCRETION MODEL

Distant extrasolar gaseous giant planets in nearly circular orbits have been detected by direct imaging observations in several systems (e.g., Kalas et al. 2008; Marois et al. 2008; Kuzuhara et al. 2013). In the conventional core accretion model, it is difficult to form cores that are massive enough to undergo runaway gas accretion at ≳ 30 AU within disk lifetime (∼ a few million years),³ because the core growth timescale is roughly proportional to a cube of the distance from the central star (e.g., Ida & Lin 2004a). Although the cores cannot be formed in such distant regions, gas giant planets formed interior to 30 AU can be scattered by other giant planets to attain semimajor axes a ≳ 100 AU (e.g., Marzari & Weidenschilling 2002; Nagasawa et al. 2008). Because the dynamical energy is lower in outer regions, a of scattered planets are broadly distributed up to ∼1000 AU. However, due to the conservation of the total angular momentum, the eccentricities e of the scattered orbits must be excited to be close to unity. While disk–planet interactions tend to damp such high values of e, they may not be efficient enough to account for the observed low eccentricities of distant gas giant planets, because the local protoplanetary disk mass may not be massive enough in the distant regions (Ida et al. 2013).

This difficulty has raised the possibility of the formation of gas giants by disk gravitational instability (e.g., Boss 2001; Helled et al. 2014 and references therein). Kratter et al. (2010) showed that if disk instability forms planetary mass clumps, it would form more abundant brown dwarfs and M-star companions. A population synthesis simulation based on the disk instability model (Forgan & Rice 2013) showed that most of such brown dwarfs and M-star companions may be retained in outer regions. However, it is not consistent with direct imaging surveys done so far. Furthermore, the observationally clear correlation between the fraction of stars with gas giants and stellar metallicity (Fischer & Valenti 2005) is not easily explained by the disk instability model, while the correlation is consistent with the core accretion model (e.g., Ida & Lin 2004b).

Based on the core accretion model, Crida et al. (2009) proposed outward type II migration for the origin of the distant gas giants in nearly circular orbits. However, the outward type II migration requires a pair of giant planets in a common gap, with the inner planet more massive than the outer one, and appropriate disk conditions.

Ida et al. (2013) found another path to form the distant gas giants in nearly circular orbits, based on the core accretion model: outward scattering of cores by a nearby gas giant followed by accretion of gas in outer regions. Because the orbital circularization through gas accretion was also shown by a hybrid N-body and two-dimensional (2D) hydrodynamical simulation (E. Thommes 2010, private communication), this path is a promising mechanism. Details of the path they found are as follows. Oligarchic growth produces similar-sized multiple cores (Kokubo & Ida 1998). Once some core starts runaway gas accretion, the planet's mass rapidly increases (e.g., Bodenheimer & Pollack 1986; Ikoma et al. 2000). When the mass increase is fast enough, the planet undergoes close encounters with nearby cores to strongly scatter them, rather than shepherd them (Zhou et al. 2007; Shiraishi & Ida 2008). Some cores are scattered to large distances, where the surface densities of both residual planetesimals and gas are relatively low. If the scattered planet is a gas giant, its mass could be comparable to or larger than the local disk mass, because the scattering may occur in a late disk evolution stage after the formation of the gas giant. Then, eccentricity damping due to dynamical friction from a local disk gas is inefficient.

On the other hand, in the core scattering model, core masses are usually well below the local disk mass. However, since the scattered cores have highly eccentric orbits and the relative velocity between the cores and the disk gas should be highly supersonic, dynamical friction from the local disk gas is less efficient (Ostriker 1999; Papaloizou & Larwood 2000; Muto et al. 2011). Muto et al. (2011) showed that the dynamical friction timescale in supersonic regime for a planet with mass M_p, orbital eccentricity e, and a semimajor axis a in a gas disk with surface density Σ is

$$\begin{eqnarray} \tau _{\rm DF} &\sim& \frac{1}{8\pi }\left(\frac{M_\ast }{\Sigma\, a^2}\right) \left(\frac{M_\ast }{M_p}\right) \left(\frac{c_s}{v_{\rm K}}\right) \left(\frac{v}{v_{\rm K}}\right)^{3}T_{\rm K}\nonumber\\ &\sim& 10^4 \left(\frac{\Sigma\, a^2}{0.01\,M_\ast }\right)^{-1} \left(\frac{M_p}{10\, M_{\oplus }}\right)^{-1} \left(\frac{c_s/v_{\rm K}}{0.1}\right) e^3 T_{\rm K},\quad \end{eqnarray} \tag{ 1 }$$

where M_* is the host star's mass, c_s is local sound velocity, and v_K and T_K are the Kepler velocity and its orbital period at a, respectively. For a planet with mass 10 M_⊕ in a highly eccentric orbit (e ∼ 1) with a ∼ 100 AU, τ_DF is as long as ∼10⁷ yr, which is longer than an observationally inferred disk lifetime ∼a few × 10⁶ yr. Thereby, the dynamical friction from local disk gas may be less effective than the orbital circularization via planetary gas accretion, which we discuss in this paper (see Section 5). Furthermore, since the dynamical friction also damps the semimajor axis efficiently in the course of eccentricity damping from values close to unity, we may not be able to retain the cores in outer regions.

Reduction in the planetesimal accretion rate decreases the critical core mass for the onset of gas accretion (e.g., Ikoma et al. 2000). The scattered cores in highly eccentric orbits spend most of the time at large distances, where the planetesimal accretion rate is significantly low, so it is possible for the scattering to trigger gas accretion onto the cores. In the case of highly eccentric orbits, since the cores spend most of the time near the apocenters at large distances, their orbits are circularized through the accretion of local gas with higher specific angular momentum. As a result, it is expected that the cores' orbits are circularized, keeping their apocenters almost fixed, in the course of gas accretion. We will show that the apocenter shrinks because the semimajor axis shrinks slightly due to energy dissipation by collision between the disk gas and the planet, and that the energy dissipation accelerates the orbital circularization (in the case of moderate eccentricity, the energy dissipation is more important for the orbital circularization than accretion of high angular momentum gas).

Assuming that the orbits of the scattered cores are quickly circularized to a degree that depends on the ratio between the planet mass and the local disk mass, keeping their semimajor axes fixed, Ida et al. (2013) performed a population synthesis simulation to predict statistical distributions of distant gas giants formed by this mechanism. They showed that the fraction of systems with the distant gas giants is a ∼0.1%–1% and most have low eccentricities (e ≲ 0.1). Although the fraction is lower, systems with multiple distant gas giants are also formed, because a single gas giant can scatter multiple cores in the inner regions. The HR8799 system has four distant gas giants and the outer three planets could be in 4:2:1 resonance. In the scattering core model, the formation of four distant gas giants is extremely rare and the probability for capture into the resonances during the inward migration associated with eccentricity damping is not clear. However, the inner two planets have semimajor axes ∼15 and 27 AU, which could be formed in situ without the scattering process. The formation of the HR8799 system by a core accretion scenario is a very interesting problem that should be addressed in the future. In this paper, we focus on a fundamental process of orbital circularization of an isolated planet through gas accretion.

In the core scattering model, a positive correlation between the semimajor axis and the mass of the distant gas giants is predicted. The critical planet mass for the gap opening is higher in the larger orbital radius (e.g., Ida & Lin 2004a). If the gap opening halts the growth of gas giants, a correlation is established, as shown in the population synthesis calculation in Ida et al. (2013).

Note that Ida et al. (2013) assumed that the eccentricities of the scattered cores are efficiently damped without any decrease in semimajor axes, implicitly assuming very efficient damping due to gas accretion, although they did not incorporate detailed orbital evolution by the gas accretion. As we will show, the eccentricity damping due to gas accretion is indeed efficient, while the degree of damping depends on how much the planets grow by accreting gas, and semimajor axes are also damped by a factor of ∼2. The population synthesis simulation must be improved by incorporating the damping formulas due to gas accretion derived in this paper, in order to discuss the distribution of distant gas giants in comparison with observations when the number of discovered planets becomes large enough for statistical argument.

Here, we investigate the orbital circularization of the scattered planets during gas accretion through detailed analytical calculations. Section 2 describes the assumptions of gas accretion onto the cores. In Section 3, we analytically derive the formulas for the orbital evolution in the course of gas accretion. In Section 4, we describe the orbital evolution by numerically solving the formulas. In Section 5, we show some results of the population synthesis calculations, by incorporating the prescriptions of orbital evolution through the mass growth due to the accretion of disk gas. Section 6 is devoted to our summary.

We start our calculation from the stage at which a core has already been scattered outward by a gas giant to attain eccentricity close to unity (e_ini) and a semimajor axis (a_ini) that is much larger than the original (a_ori). Note that the pericenter distance of the scattered planet's orbit must be close to a_ori: a_ori ≃ q_ini = a_ini(1 − e_ini). Since a_ori should be close to the gas giant's orbital radius, we can regard that q_ini ∼ 1–10 AU, based on a core accretion model (e.g., Ida & Lin 2004a). As we will show, the pericenter distance of the scattered planet quickly increases. Accordingly, the planet immediately becomes isolated from the perturbing gas giant, so that we neglect its further perturbations. We do not calculate the initial scattering process by the gas giant, but instead study the evolution of e and a of the scattered planet due to accretion of gas, for given e_ini, a_ini, and disk radius r_d, to derive general formulas for the orbital evolution. In the following, we explain our prescriptions for accretion of gas onto the scattered planet.

After a core mass exceeds a critical core mass, the pressure gradient no longer supports the gas envelope of the planet against the planetary gravitational force and a quasi-static contraction of the gas envelope starts (e.g., Mizuno 1980; Bodenheimer & Pollack 1986). The critical core mass is given by (Ikoma et al. 2000)

$$\begin{equation} M_{\rm c,crit} \simeq 10 \left(\frac{\dot{M}_{\rm c}}{M_\oplus /10^6\, {\rm yr}} \right)^{(0.2-0.3)} \left(\frac{\kappa }{\kappa _{0}} \right)^{(0.2-0.3)} M_{\oplus }, \end{equation} \tag{ 2 }$$

where κ is the opacity of the gas envelope and κ₀ is that of the minimum-mass solar nebula model (Hayashi 1981). Since the planetesimal accretion rate, $\dot{M}_{\rm c}$, determines the heat energy source to support the envelope, a lower value of $\dot{M}_{\rm c}$ leads to a smaller value of M_{c, crit}. In general, a planetesimal accretion rate rapidly decreases with orbital radius (e.g., Ida & Lin 2004a). After a core is scattered outward to acquire high orbital eccentricity, the core spends most of its time at much larger orbital radii than near the original location. Thereby, an orbit-averaged value of $\dot{M}_{\rm c}$ is significantly lowered and it is likely that a quasi-static contraction of the gas envelope is initiated by the outward scattering.

According to the quasi-static contraction, disk gas can be supplied to a Hill radius or a Bondi radius of the planet. This means that the gas accretion rate onto the planet is regulated by heat transfer through the envelope, rather than environmental disk gas conditions, except in the final stage in which the contraction is very fast. The (Kelvin–Helmholtz) timescale of the envelope contraction is given by (Ikoma et al. 2000; Ikoma & Genda 2006)

$$\begin{equation} \tau _{\rm KH} \simeq 10^{10} \left(\frac{M}{M_\oplus } \right)^{-(3-4)} \left(\frac{\kappa }{\kappa _{\rm ini}}\right)\,{\rm yr}. \end{equation} \tag{ 3 }$$

From these arguments, we assume that the gas accretion rate does not depend on the position of the eccentric orbit during an orbital period, although environmental disk conditions change considerably during one orbital cycle for highly eccentric orbits.

In general, the dependence of τ_KH on orbital radius r is weak for the radiation-dominated envelope (e.g., Ikoma et al. 2000). If a convective envelope develops, the envelope contraction rate can be affected by disk temperature and density (Ikoma et al. 2001; Piso & Youdin 2014). However, even if τ_KH has the r-dependence, the dependence may be smoothed out when τ_KH is longer than the orbital period, that is, when M ≲ 100 M_⊕, because the response time of the envelope structure is given by τ_KH. As we show in the following, the gas accretion rate onto the planet may be regulated by disk gas supply, rather than by envelope contraction for M ≳ 100 M_⊕.

Note that in the highly eccentric orbit, the Bondi and Hill radii significantly change during one orbital circulation. The change might also induce oscillation of the gas envelope, which could affect gas accretion and heat generation/cooling. An investigation of this effect is left for a future work. We will only assume the constant accretion rate during one orbit inside the disk, but will not adopt any particular form of τ_KH.

When the envelope contraction is faster than the supply of gas and the supply has the r-dependence, the assumption of the constant accretion rate is violated. The supply can be limited by global disk accretion and Bondi accretion. The limit by global disk accretion becomes important for M ≳ 100 M_⊕, because the quasi-static contraction rate is given by $\dot{M}_{\rm KH} \sim M/\tau _{\rm KH} \sim 10^{-10}(M/M_{\oplus })^{(4-5)}(\kappa /\kappa _{\rm ini})^{-1}\,M_{\oplus }\,{\rm yr^{-1}}$ (Equation (3)) and the observationally inferred typical value of disk accretion rates onto T Tauri stars is $\dot{M}_{\rm disk} \sim 10^{-8}M_{\odot }\,{\rm yr^{-1}} \sim 3 \times 10^{-3}\,M_{\oplus }\,{\rm yr^{-1}}$. However, the r-dependence of the disk accretion rate is very weak in the regions of r ≪ r_d, where a steady-accretion-disk approximation is valid.

The Bondi accretion rate is given by $\dot{M}_{\rm Bondi} \sim \pi \rho _g (GM/v^2)^2 v$, where ρ_g is disk gas density and v is the relative velocity between the planet and disk gas. Both ρ_g and v sensitively depend on r. In general, $\dot{M}_{\rm Bondi} < \dot{M}_{\rm KH}$ for high e and large M_p. However, our calculations start from small M_p, and e is already damped when M_p becomes large. We found that in most orbital evolutions we consider, $\dot{M}_{\rm Bondi} \gtrsim \dot{M}_{\rm KH}$ and the supply limit by Bondi accretion does not occur.

Thus, our assumption of constant gas accretion may be justified. As a result of the time-independent gas accretion rate, the planet accretes gas preferentially in outer regions where the planet spends most of the time. The specific angular momentum of the planetary orbit is given by $\ell _{\rm p} = \sqrt{GM_*a(1-e^2)} \sim \sqrt{2GM_*q}$, where q = a(1 − e) is the pericenter distance and e ∼ 1 is assumed. That of the local gas near the apocenter (Q = a(1 + e) ≃ 2a) is given by $\ell _{\rm g} \sim \sqrt{2GM_*a}$ (we assume circular Keplerian motion for the disk gas). Bcause q = a(1 − e) ≪ a for e ∼ 1, the planet's specific angular momentum is increased by the accretion of local gas. Thus, the planetary orbit tends to be circularized with the apocenter distance (Q) fixed.

If the orbit deviates from the disk with a finite size r_d, beyond which gas density is significantly declined, we halt gas accretion at r > r_d. In the following derivations, we consider two cases: (1) the apocenter is inside the disk (Q < r_d) and (2) it is outside the disk (Q > r_d). We will refer to cases (1) and (2) as an embedded case and deviated case, respectively. In the deviated case, the planet mostly accretes gas at r ∼ r_d and planetary orbits tend to be fitted to circular orbits at r_d rather than those at Q.

The specific orbital energies of the planet and local gas near the apocenter are _p = −GM_*/2a and _g = −GM_*/2a(1 + e) ∼ −GM_*/4a, respectively. Near the apocenter, the planet's specific orbital energy is increased by the accretion of local gas near the apocenter. However, the accretion of lower specific energy near the pericenter is significant due to a deep potential near the pericenter, in spite of the fast passage of the pericenter. As we show in Section 4.1, in the embedded case, the orbit-averaged specific orbital energy of accreting gas is exactly the same as that of the planet, irrespective of orbital eccentricity. The planet's specific orbital energy actually decreases if collisional dissipation between the planet and disk gas is taken into account. It also contributes to eccentricity damping, with slight decay of the semimajor axis.

In the embedded case, the gas accretion rate onto the planet is independent of the phase of the orbits. In the deviated case, we assume a constant gas accretion rate at r < r_d and a zero accretion rate at r > r_d. We do not assume the value of the constant accretion rate, because we will derive orbital evolution as a function of planetary mass M_p, but not as a function of time.

The relative velocity between the planet and the disk gas would be supersonic almost everywhere for highly eccentric orbits with e ≳ h/r ∼ 0.1, where h is the disk scale height. For incident supersonic gas flow to stay in the Hill radius or Bondi radius, we need some energy dissipation. Bow shock in front of the planet may provide the energy dissipation. We will leave the full hydrodynamic simulations on a bow shock for future work and assume that the planet accretes disk gas in unperturbed flow and the accretion rate is independent of orbital phase in most of calculations. Note, however, that the relative velocity between the gas flow and the planet is smaller in the post-shock flow than in the unperturbed flow, which may make the eccentricity damping less efficient. In Section 4.3, we perform calculations, taking into account the effect of the shock with a simple one-dimensional (1D) model, and show that the eccentricity damping is indeed slowed down but does not significantly change our conclusion.

When cores are scattered by a gas giant in the early stage, eccentricities are preferentially pumped up compared with inclinations. However, if the core undergoes repeated close encounters before its orbit is circularized, orbital inclinations are also excited. We also calculated e and a evolution with non-zero inclinations, and found that the final values of e and a change by less than 5% if the inclination is smaller than 30 degrees. So, we here show the results with zero inclinations.

In summary, the assumptions we use in most of runs are as follows.

• 1.  
  The gas disk is in Keplerian rotation. Because the relative velocity between the gas and the planet is generally supersonic, the gas is hardly perturbed by the planetary gravitational perturbations.
• 2.  
  The motions of the planet and the gas disk are coplanar.
• 3.  
  The gas accretion rate onto the planet is consistent with time during one orbit, which means that orbit-averaging can be done. The planet captures the local unperturbed disk gas, conserving mass and angular momentum (energy is not conserved).
• 4.  
  If the gas disk has the finite size, we truncate gas accretion during the period in which the planet goes out of the disk.

In Section 3 through Section 4.2, we adopt these assumptions and derive analytical formulas to describe the orbital circularization process through the gas accretion. Even if we adopt assumption 4, analytical formulas are derived, because the constant gas accretion rate is still applied at r < r_d and analytical orbit averaging can be done. If we include the effect of shock dissipation, analytical integration is not possible, so we show the orbital evolution obtained by numerical integration (Section 4.3).

With the assumptions 1 to 4 described in the previous section, we analytically derive formulas to calculate the orbital evolution in the form of differential equations. Numerically integrating the differential equations, we will show the evolution paths of e and a that are uniquely determined by their initial values and r_d.

According to discussions in Section 2, we first calculate changes in the angular momentum and energy of the planet, ΔL and ΔE, during one orbital period, assuming that the mass accretion rate is constant with time during one orbital period. We also assume that the changes in orbital elements are small enough over one orbital period; in other words, the mass of the captured gas during one orbit (ΔM) is much smaller than the instantaneous planetary mass (M).

The changes ΔL and ΔE are then given by

$$\begin{eqnarray} \Delta L & \simeq \int l_{\rm gas} \; {dM} = \Delta M \; \frac{1}{t_{\rm d}} \int _{-t_{\rm d}/2}^{t_{d}/2} l_{\rm gas}\, {dt}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \Delta E &\simeq& \int (\epsilon _{\rm gas} - \epsilon _{\rm coll}){dM}\nonumber\\ &=& \Delta M \; \frac{1}{t_{d}} \int _{-t_{\rm d}/2}^{t_{\rm d}/2} (\epsilon _{\rm gas} - \epsilon _{\rm coll} ){dt}, \end{eqnarray} \tag{ 5 }$$

where the integral is during one orbit, t_d is a duration at r < r_d (t_d ⩽ T_K), t = 0 is a pericenter passage, l_gas and _gas are the specific angular momentum and the orbital energy of accreting gas, and _coll is energy dissipation by collision between the planet and accreting gas.

Through ΔL and ΔE during the mass growth of ΔM, the specific angular momentum and the orbital energy of the planet, ℓ_p and _p, are changed. Since ΔL = (Δℓ_p + ℓ_p)(M + ΔM) − Mℓ_p ≃ Δℓ_p · M + ℓ_pΔM and ΔE ≃ Δ_p · M + _pΔM, the change rate of ℓ_p and _p of the planet in one orbital period are expressed as

$$\begin{eqnarray} \frac{\Delta \ell _{\rm p}}{\ell _{\rm p}} & \simeq \frac{\Delta L}{M\ell _{\rm p}} - \frac{\Delta M}{M} = \frac{\Delta M}{M} \; f_\ell, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \frac{\Delta \epsilon _{\rm p}}{\epsilon _{\rm p}} & \simeq \frac{\Delta E}{M\epsilon _{\rm p}} - \frac{\Delta M}{M} = \frac{\Delta M}{M} \; f_\epsilon, \end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray} f_\ell &=\frac{1}{t_{\rm d}} \int _{-t_{\rm d}/2}^{t_{\rm d}/2} \left(\frac{l_{\rm gas}}{l_{\rm p}} -1 \right) {dt}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} f_\epsilon &=\frac{1}{t_{\rm d}} \int _{-t_{\rm d}/2}^{t_{\rm d}/2} \left(\frac{\epsilon _{\rm gas}}{\epsilon _{\rm p}} - \frac{\epsilon _{\rm coll}}{\epsilon _{\rm p}} -1 \right){dt}. \end{eqnarray} \tag{ 9 }$$

Since $\ell _{\rm p} = \sqrt{GM_\ast a(1-e^2)}$ and _p = −GM_*/2a, where M_* is the host star's mass and G is the gravitational constant, the changes Δℓ_p and Δ_p are related with the changes of the eccentricity and semimajor axis (Δe and Δa) as

$$\begin{eqnarray} \frac{\Delta \ell _{\rm p}}{\ell _{\rm p}} &= \frac{\Delta (\sqrt{a})}{\sqrt{a}} + \frac{\Delta (\sqrt{1-e^2})}{\sqrt{1-e^2}} \simeq \frac{\Delta a}{2a} - \frac{e}{1-e^2}\Delta e, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \frac{\Delta \epsilon _{\rm p}}{\epsilon _{\rm p}} &= \frac{\Delta (a^{-1})}{a^{-1}} \simeq - \frac{\Delta a}{a}. \end{eqnarray} \tag{ 11 }$$

So, Δe and Δa are given by

$$\begin{eqnarray} \Delta e &\simeq - \frac{1-e^2}{e} \left(\frac{\Delta \ell _{\rm p}}{\ell _{\rm p}} + \frac{1}{2}\frac{\Delta \epsilon _{\rm p}}{\epsilon _{\rm p}} \right)= - \frac{\Delta M}{M} f_e, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \frac{\Delta a}{a} &\simeq - \frac{\Delta \epsilon _{\rm p}}{\epsilon _{\rm p}} = - \frac{\Delta M}{M} f_a, \end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray} f_e &= \frac{1-e^2}{e} \left(f_\ell + \frac{1}{2} f_\epsilon \right), \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} f_a &= f_\epsilon. \end{eqnarray} \tag{ 15 }$$

We derive the differential equations for the orbit-averaged evolution of e and a in terms of planetary mass M:

$$\begin{eqnarray} \frac{d e}{d \log M} &\simeq \frac{\Delta e}{\Delta M/M} = - f_e, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \frac{d \log a}{d \log M} &\simeq \frac{\Delta a/a}{\Delta M/M}= - f_a. \end{eqnarray} \tag{ 17 }$$

From these equations, we also obtain

$$\begin{equation} \frac{d \log a}{d e}\simeq \frac{f_a}{f_e}. \end{equation} \tag{ 18 }$$

So far, we have not assumed any specific forms for ℓ_gas and _gas. Here we assume that the planet captures gas in circular Keplerian motion to analytically derive formulas f_e and f_a. (In Section 4.3, we calculate f_e and f_a for post-shocked gas flow using a simple 1D model.) Note that the analytical formulas f_e and f_a are derived even for deviated cases where Q > r_d.

For unperturbed gas flow (circular Keplerian flow),

$$\begin{eqnarray} \ell _{\rm gas} & = \sqrt{GM_\ast r}, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} \epsilon _{\rm gas} & = - \frac{GM_\ast }{2r}, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} \epsilon _{\rm coll} & = \frac{v_{\rm {rel}}(r)^2}{2}, \end{eqnarray} \tag{ 21 }$$

where v_rel(r) is the relative velocity between the planet and the local gas and r is the instantaneous distance of the planet from the central star. The radial and tangential components of instantaneous velocity of an eccentric Keplerian orbit of the planet at r are given by

$$\begin{eqnarray} v_r & = v_{\rm K} \sqrt{2-\frac{r}{a}-\frac{a}{r}(1-e^2)}, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} v_\phi & = v_{\rm K} \sqrt{\frac{a}{r}(1-e^2)}, \end{eqnarray} \tag{ 23 }$$

where a and e are the planet's semimajor axis and eccentricity. Because the local Keplerian velocity is given by $v_{\rm K}=\sqrt{GM_\ast /r}$, the square of relative velocity is

$$\begin{eqnarray} v_{\rm {rel}}(r)^2 &= v_r^2 + (v_\phi - v_{\rm K})^2\qquad\qquad\qquad \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} &\quad\quad\;\,\,= \frac{GM_\ast }{r} \left(3 - \frac{r}{a} - 2\sqrt{\frac{a}{r} (1-e^2) } \right). \end{eqnarray} \tag{ 25 }$$

For integrating Equations (8) and (9), we convert time to eccentric anomaly using the Kepler equation. The time average of the powers of r^α(α = 1/2, −1, −3/2) are analytically integrated using the conversion:

$$\begin{eqnarray} &&{\frac{1}{t_{\rm d}} \int _{-t_{\rm d}/2}^{t_{\rm d}/2} \left(\frac{r}{a} \right)^{1/2} dt =\frac{2\sqrt{1+e}}{3} f_{1/2}(e,u_\mathrm{d});} \nonumber \\ &&{f_{1/2}(e,u_\mathrm{d})=\frac{4(E(k)-E(y,k)) - (1-e)(K(k)-F(y,k)) - e\sin 2y \sqrt{1-k^2\sin ^2y} }{u_\mathrm{d}-e\sin u_\mathrm{d}},\quad }\nonumber\\ && \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} &&\frac{1}{t_{\rm d}} \int _{-t_{\rm d}/2}^{t_{\rm d}/2} \left(\frac{r}{a} \right)^{-1} dt =f_{-1}(e,u_\mathrm{d}); \nonumber\\ && f_{-1}(e,u_\mathrm{d})=\frac{u_\mathrm{d}}{u_\mathrm{d}-e\sin u_\mathrm{d}}, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} &&\frac{1}{t_{\rm d}} \int _{-t_{\rm d}/2}^{t_{\rm d}/2} \left(\frac{r}{a} \right)^{-3/2} dt =\frac{2}{\sqrt{1+e}}f_{-3/2}(e,u_\mathrm{d});\nonumber\\ && f_{-3/2}(e,u_\mathrm{d})=\frac{K(k)-F(y,k)}{u_\mathrm{d}-e\sin u_\mathrm{d}}, \end{eqnarray} \tag{ 28 }$$

where $k \equiv \sqrt{2e/(1+e)}$, y ≡ (π − u_d)/2, K(k) is the complete elliptic integral of the first kind, E(k) is the complete elliptic integral of the second kind, F(y,k) is the elliptic integral of the first kind, E(y, k) is the elliptic integral of the second kind, and u_d is the maximum eccentric anomaly (0 < u_d < π) within the disk (r < r_d), which is given by

$$\begin{eqnarray} u_\mathrm{d}\equiv \left\lbrace \begin{array}{@{}ll@{}} \displaystyle \cos ^{-1} \left[ \frac{1}{e} \left(1- \frac{r_\mathrm{d}}{a} \right)\right]& [{\rm for}\ Q>r_\mathrm{d}], \\ \pi & [{\rm for} \ Q <r_\mathrm{d}]. \\ \end{array} \right. \end{eqnarray} \tag{ 29 }$$

In the embedded case, u_d = π and t_d = T_K, while u_d < π and t_d < T_K, depending on r_d, in the deviated cse.

With Equations (26)–(28), Equations (8) and (9) are written as

$$\begin{eqnarray} f_\ell (e,u_\mathrm{d}) &=\frac{2}{3\sqrt{1-e}}f_{1/2}(e,u_\mathrm{d}) -1, \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} f_\epsilon (e,u_\mathrm{d}) &=4 f_{-1}(e,u_\mathrm{d}) - 4\sqrt{1-e}f_{-3/2}(e,u_\mathrm{d}) -2.\quad \end{eqnarray} \tag{ 31 }$$

Note that f_ℓ and f are functions of only e, independent of a, in the embedded case. Even in the deviated case, the a-dependence enters f_ℓ and f only through the scaled quantity a/r_d in u_d.

From Equations (14) and (15),

$$\begin{eqnarray} \frac{d e}{d \log M} & =& -\!f_e(e,u_\mathrm{d})\nonumber\\ &=& -\!\frac{1-e^2}{e} \left(f_\ell (e,u_\mathrm{d}) + \frac{1}{2} f_\epsilon (e,u_\mathrm{d}) \right), \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray} \frac{d \log a}{d \log M} & = -\!f_a(e,u_\mathrm{d}) = -\!f_\epsilon (e,u_\mathrm{d}). \end{eqnarray} \tag{ 33 }$$

These equations show that while the semimajor axis is damped only by Δ_p, the orbital eccentricity is damped by both Δℓ_p and Δ_p.

By numerically integrating Equations (32) and (33), we obtain the evolution of the orbital elements according to the planetary mass growth. In Figure 1, we plot f_ℓ(e, u_d), f_a(e, u_d)[ = f(e, u_d)], and f_e(e, u_d) as a function of e for u_d = π, π/2, and π/4. Because f_e(e, u_d), f_a(e, u_d) > 0 for any values of e and u_d, e and a monotonically decrease as the planet grows through the accretion of disk gas. For u_d = π (embedded case), f_ℓ dominates the eccentricity damping. That is, the accretion of disk gas with high specific angular momentum near the apocenter is responsible for the eccentricity damping. On the other hand, for u_d = π/2 and π/4 (deviated cases), f_ℓ is small or negative except for high e. In these cases, the energy dissipation by collision between incident gas and the planet is responsible for the eccentricity damping (see Section 4.2).

Zoom In Zoom Out Reset image size

When e becomes small enough, dlog a/dlog M quickly approaches zero (f_a → 0). Thus, the asymptotic values of a are uniquely determined by the initial values of e, a, and r_d. In the next section, we show the numerically obtained evolution paths.

4.1. The Embedded Case

First, we consider the embedded case; that is, entire parts of a planetary orbit are embedded in the disk. Because Equation (16) is independent of multiplication of M by a constant factor, we can adopt a scaled quantity M/M_ini as a variable, where M_ini is the initial value of M. Then, Equation (16) includes only e and M/M_ini, so the evolution of e is uniquely given as a function of M/M_ini for any initial values of eccentricity (e_ini). The evolutional paths for representative values of e_ini that were obtained by a numerical integration of Equation (16) are shown in Figure 2(a). This figure shows that e decreases to values below 0.2 e_ini when M attains 10M_ini. The typical core mass required to start runaway gas accretion is ∼10 M_⊕, which means that e is reduced to values smaller than 0.2 when the planet acquires Saturnian mass (∼100 M_⊕), even if its initial orbit was close to a parabolic orbit (e_ini ∼ 1). If the mass is scaled by that at e = 0.1, denoted by M₀₁, e is uniquely determined by M/M₀₁. Figure 2(b) shows the self-similar solution of e as a function of M/M₀₁ (the solid curve). The fitting formula given by Equation (41), which is presented in Section 5, is also plotted with the dashed curve.

Zoom In Zoom Out Reset image size

Since Equation (18) has a similar structure to Equation (16), the evolution of e is uniquely given as a function of a/a_ini. Figure 3(a) shows the evolutionary paths on the a-e plane for representative values of e_ini. Because both e and a keep decreasing, the evolution starts at the right end and moves leftward. Damping of e is dominated over that of a except for e ≃ 1. Even if e_ini = 0.99, the asymptotic value of a for e → 0, which we denote as a_final, is as much as ∼0.48 a_ini. As in the case of the e-M relation, if the semimajor axis is scaled by that at a specific value of e, the evolution of e is expressed by a single curve, regardless of e_ini and a_ini. Figure 3(b) shows the self-similar relation, a/a₀, as a function of e (the solid curve), where a₀ is a at e = 0. The fitting formula given by Equation (43), which is presented in Section 5, is also plotted with the dashed curve. This plot clearly shows that damping of a is much smaller than that of e: for the damping of e from ∼1 to 0, a is decreased by ∼50%, and for that from 0.8 to 0, the decrease in a is only ∼30%.

Zoom In Zoom Out Reset image size

Figure 4(a) shows specific orbital angular momentum and energy of local gas (ℓ_gas and _gas) scaled by those of the planet (ℓ_p and _p) for e = 0.9 as functions of t/T_K. The pericenter passage is at t/T_K = 0 and t/T_K = 1, and the apocenter passage is at t/T_K = 0.5, respectively. Because the orbit is highly eccentric, for most of an orbital period, ℓ_gas > ℓ_p and _gas > _p (because the energy is negative, |_gas| < |_p|) except in the regions close to pericenter (t/T_K = 0 and t/T_K = 1). As a result, an orbit-averaged value of ℓ_gas/ℓ_p is considerably larger than unity (in this case, it is 〈ℓ_gas/ℓ_p〉 = f_ℓ(0.9, π) + 1 = 2.66) and the specific angular momentum of the planet increases through these accretion of disk gas.

Zoom In Zoom Out Reset image size

On the other hand, it is analytically shown that an orbit-averaged value of _gas/_p is unity, because

$$\begin{eqnarray} \left \langle \frac{\epsilon _{\rm gas}}{\epsilon _{\rm p}} \right \rangle &=&\frac{1}{T_{\rm K}} \int _{-T_{\rm K}/2}^{T_{\rm K}/2} \left(\frac{\epsilon _{\rm gas}}{\epsilon _{\rm p}}\right)d t=\frac{1}{T_{\rm K}} \int _{-T_{\rm K}/2}^{T_{\rm K}/2} \left(\frac{r}{a} \right)^{-1} d t \nonumber\\ &=& f_{-1}(e,\pi) = \frac{\pi }{\pi - e \sin \pi } = 1. \end{eqnarray} \tag{ 34 }$$

Although |_gas| < |_p| most of time, |_gas| is much larger than |_p| near the pericenter passage (Figure 4(a)). The contribution of large |_gas| compensates for the excess energy accretion in outer regions. However, the orbital energy of the planet decreases because it is also contributed by the collisional energy dissipation $\epsilon _{\rm coll}({=}v_{\rm {rel}}^2/2)$. Equation (12) shows that the collisional energy dissipation also damps e. If the collisional energy dissipation is neglected, a is conserved and e is damped by the accretion of higher specific angular momentum gas. With the effect of the collisional energy dissipation, damping of e is faster and a is also damped while the a-damping is slower than the e-damping.

From these relations, together with _p = −GM_*/2a and $\ell _{\rm p} = \sqrt{GM_\ast a(1-e^2)}$, it is readily found that both a and e always decrease in the case of constant gas accretion rate. Figure 4(b) shows ℓ_gas/ℓ_p and _gas/_p at e = 0.2. In this case, the integrals are more symmetric about ℓ_gas/ℓ_p = 1 and _gas/_p = 1. Thereby, when e is reduced to ≲ 0.2, the orbit-averaged values of ℓ_gas/ℓ_p and _gas/_p are nearly unity and the decrease in e and a due to planetary mass growth slows down.

The initial pericenter distance q_ini of the core's orbit before the e-damping process would correspond to the original semimajor axis of the core (a_ori) before the core was scattered by a gas giant, which may be ∼1–10 AU. Figure 5 shows asymptotic semimajor axis a_final scaled by q_ini. Because a_final is the final semimajor axis of a gas giant formed from a scattered core after e is damped, a_final/q_ini indicates an efficiency to send a planet to outer regions. In this figure, we find that a core originally at the inner region (a_ori ∼ q_ini ∼ 10 AU) becomes a gas giant with a large radius (a_final ≳ 30 AU when e_ini ≳ 0.73, and a_final ≳ 100 AU when e_ini ≳ 0.94).

Zoom In Zoom Out Reset image size

Figure 6 shows the evolution of the pericenter distance scaled by the initial one, q/q_ini, due to planetary mass growth, for representative values of e_ini. It is shown that q/q_ini quickly increases, which justifies our assumption that the scattered planet becomes quickly isolated and the further perturbations from the gas giant in the inner region are neglected.

Zoom In Zoom Out Reset image size

4.2. The Deviated Case

Next, we consider the deviated case in which Q > r_d. In this case, gas accretion is halted at r > r_d, so the planet cannot accrete gas with higher specific angular momentum and energy, which results in smaller f_ℓ and larger f(= f_a) (see Figure 1). The increase of f(= f_a) is more effective than the decrease of f_ℓ, so f_e is larger. Thus, both the e and a dampings in the deviated case are more efficient than in the embedded case.

Since there is a characteristic length r_d, a self-similar solution like Figure 3(b) does not exist. However, it is clear that the evolutions of Q/r_d and q/r_d should be the same for the same initial values. Figure 7 shows the evolutions of Q/r_d and q/r_d. The evolutions are in the direction of the bottom right. The evolutions in the deviated case correspond to those in the region of Q > r_d. We also added the following embedded evolutions in the region of Q < r_d. Because the evolution paths do not cross each other, we can parameterize the evolution paths with one parameter. In Figure 7, we used the value of q/r_d at the time when Q is reduced to be r_d, as the parameter.

Zoom In Zoom Out Reset image size

The evolutions of e and a corresponding to the solutions in Figure 7 are plotted in Figure 8. The orbital evolution is toward the direction of the bottom left. This figure shows that in the early phase of Q > r_d, the semimajor axis is predominantly damped. Note that even in this phase where a is rapidly reduced, Figure 7 shows that q increases so quickly that the planet becomes isolated from the perturbing gas giant. For e_ini ∼ 1, a is damped by orders of magnitude until Q ∼ a(1 + e) ∼ 2a is reduced to ∼r_d. In the following embedded phase of Q < r_d, however, a is reduced at most by a factor of two, as shown. Thus, in this case, a_final ∼ r_d/4, independent of the values of a_ini, as long as e_ini ∼ 1. In other words, we can infer the values of r_d from a_final.

Zoom In Zoom Out Reset image size

Figure 9 shows the evolution of e as a function of the planetary mass M/M_Q, where M_Q is the planetary mass at Q = r_d. In the deviated phase, the e-damping is more efficient than in the embedded phase, although it is slightly slower because of high eccentricity.

Zoom In Zoom Out Reset image size

4.3. The Effect of Shock

We have considered the energy dissipation by collision between the incident gas flow and the planet. The dissipation is needed to bind the gas around the planet and it accelerates eccentricity damping as shown. The dissipation should occur through bow shock in front of the planet. The shock not only causes the energy dissipation, but also makes the relative velocity lower. So far, we have neglected the relative velocity damping by shock, which should weaken the eccentricity and semimajor axis damping. Here we evaluate the effect of the shock using a simple 1D model. Full 2D or three-dimensional hydrodynamical simulations will be done in a separate paper.

The simple 1D model we use is as follows. The ratio of post-shock velocity (sv_rel) to pre-shock velocity (v_rel) is

$$\begin{equation} s =\frac{(\gamma -1){\cal M}^2+2}{(\gamma +1){\cal M}^2} \simeq \frac{1}{4} \left(1+ \frac{3}{{\cal M}^2} \right), \end{equation} \tag{ 35 }$$

where γ is the specific heat ratio (γ = 5/3 in the monatomic molecule) and ${\cal M}$ is the Mach number for pre-shock gas flow, which is given by

$$\begin{equation} {\cal M}=\frac{v_{\rm {rel}}}{c_s} = 30 \left(\frac{r}{1\,{\rm AU}} \right)^{-1/4} \left(3 - \frac{r}{a} - 2\sqrt{\frac{a}{r} (1-e^2) } \right)^{1/2}, \end{equation} \tag{ 36 }$$

where we used an optically thin disk temperature, T = 280(r/1 AU)^−1/2 K, for an evaluation of sound velocity c_s. For the subsonic case (${\cal M} <1$), s = 1. The radial and tangential components of the velocity of the post-shock gas (u_r, u_ϕ) and those of the pre-shock gas (0, v_K) are related to those of a planet (v_r, v_ϕ), which are given by Equations (22) and (23), as

$$\begin{equation} u_r - v_r = s (0 - v_r), \end{equation} \tag{ 37 }$$

$$\begin{equation} u_\phi - v_\phi = s (v_{\rm K} - v_\phi). \end{equation} \tag{ 38 }$$

Then, the integrants of Equations (8) and (9) for f_ℓ and f are

$$\begin{eqnarray} \frac{l_{\rm gas}}{l_{\rm p}} -1 &= \frac{r u_\phi }{r v_\phi } -1 = s \left(\frac{r v_{\rm K}}{r v_\phi } -1 \right) = s \left(\frac{l_{\rm gas,0}}{l_{\rm p}} -1 \right),\qquad \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} \frac{\epsilon _{\rm gas}}{\epsilon _{\rm p}} - \frac{\epsilon _{\rm coll}}{\epsilon _{\rm p}} -1 &= \frac{\left(u_r^2+u_\phi ^2\right)/2 - GM_\ast /r - s^2 v_{\rm {rel}}^2/2}{\left(v_r^2+v_\phi ^2\right)/2 - GM_\ast /r} -1 \nonumber \\ &= s \left(\frac{v_{\rm K}^2/2 - GM_\ast /r - v_{\rm {rel}}^2/2}{\left(v_r^2+v_\phi ^2\right)/2 - GM_\ast /r} -1 \right)\nonumber\\ &= s \left(\frac{\epsilon _{\rm gas,0}}{\epsilon _{\rm p}} - \frac{\epsilon _{\rm coll,0}}{\epsilon _{\rm p}} -1 \right), \end{eqnarray} \tag{ 40 }$$

where quantities with the subscript "0" are those for unperturbed gas flow neglecting shock ($l_{\rm gas,0}=\sqrt{GM_\ast r}$, _{gas, 0} = −GM_*/2r, and $\epsilon _{\rm coll,0}=v_{\rm {rel}}^2/2$).

As we showed in the previous subsections, the integrations for f_ℓ and f (Equations (8) and (9)) can be analytically done in the case neglecting the damping of the relative velocity (equivalently, s = 1). However, since in the present case s varies along the orbit (Figure 10), we integrate Equations (8) and (9) numerically. Here we consider the embedded case. f_ℓ and f depend on a through s in the present case, so we assume a = 100 AU (a is variable when we consider the evolution paths). Note that evolution paths in the a–e plane are the same as those in the case without the relative velocity damping. In Equation (18), dlog a/de, is given approximately by orbit-averaged f_a and f_e. But, more exactly, dlog a/de must be integrated every time. From Equations (39) and (40), it is apparent that s completely cancels and dlog a/de is the same.

Zoom In Zoom Out Reset image size

The functions f_ℓ, f_a(= f), and f_e in the case with shock are compared with those without shock in Figure 11. The effect of shock lowers all the functions. Accordingly, both e and a dampings are slowed down, whereas the evolution paths on the a–e plane do not change. The evolution paths on the e–M, a–M, and q–M planes are shown in Figure 12. The initial conditions are q_ini = 10 AU, e_ini = 0.9, and M_ini = 10 M_⊕, respectively. While e declines to values <0.2 at M_p ∼ 50 M_⊕ in a non-shock case, it does not become <0.2 until M_p ∼ 3000 M_⊕ in a shock case. However, because the masses of direct-imaged planets are relatively large (∼10 M_J), the orbital circularization is still effective. The a damping is also slowed down, but the semimajor axis at M_p ∼ 10 M_J is not significantly larger than that in a non-shock case. The pericenter distance q is still quickly increased, so the assumption that the scattered planet becomes quickly isolated and the further perturbations from the gas giant in inner region are neglected is justified. Thus, although the eccentricity damping is less efficient, the formation of distant gas giants in nearly circular orbits is not significantly inhibited by the effect of shock.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The self-similar solution in Figure 2(b) can be approximately fitted by

$$\begin{equation} e \simeq 0.1\left[1- \log _{10}\left(\frac{M}{M_{01}}\right)\right]^3. \end{equation} \tag{ 41 }$$

where M₀₁ is M at e = 0.1. Then, for any given e_ini and M_ini, e for M(> M_ini) are evaluated by

$$\begin{eqnarray} e & \simeq 0.1\left[1 - \log _{10}\left(\frac{M_{\rm ini}}{M_{01}}\right) - \log _{10}\left(\frac{M}{M_{\rm ini}}\right)\right]^3 \nonumber \\ & = 0.1\left[\left(\frac{e_{\rm ini}}{0.1}\right)^{1/3} - \log _{10}\left(\frac{M}{M_{\rm ini}}\right) \right]^3. \end{eqnarray} \tag{ 42 }$$

Because this is an approximate formula, Equation (42) can be negative for large values of M/M_ini. In such cases, we set e = 0, because e has very small values in the exact solution. On the other hand, the self-similar solution in Figure 3(b) can be fitted by

$$\begin{equation} \frac{a}{a_0} \simeq 1 + 0.6\, e^3 + \frac{0.007}{1-e}, \end{equation} \tag{ 43 }$$

where a₀ is the asymptotic values of a at e → 0. Because a similar relation holds for a_ini and e_ini,

$$\begin{equation} \frac{a}{a_{\rm ini}} \simeq \frac{1 + 0.6\, e^3 + 0.007/(1-e)}{1 + 0.6\, e_{\rm ini}^3 + 0.007/(1-e_{\rm ini})}. \end{equation} \tag{ 44 }$$

In the population synthesis simulation, when a core with mass M_c closely encounters a gas giant, the eccentricity and semimajor axis that are excited by the scattering are evaluated with a Monte-Carlo procedure (see, e.g., Ida et al. 2013). We set the eccentricity, semimajor axis, and M_c as e_ini, a_ini, and M_ini in these equations, respectively. The mass growth of the planet due to gas accretion after the scattering is also calculated in the population synthesis simulation. The mass growth is truncated when the disk gas is severely depleted or a clear gap along the planetary orbit is opened (see, e.g., Ida & Lin 2004a; Ida et al. 2013). From Equation (42), the value of e when the planet mass increases to M can be derived from e_ini and M_ini. The semimajor axis a at M is derived from e_ini, e, and a_ini from Equation (44). Note that Ida et al. (2013) simply assumed e = 0 and a = a_ini.

Figures 13 show the e–a distributions of gas giant planets around solar-type stars, obtained by a population synthesis calculation with similar parameters to those of the results in Figure 7 of Ida et al. (2013). Note that rocky and icy planets with smaller masses are not plotted here. In panel (a), neither the dynamical friction to cores nor the damping via gas accretion is included. Most of giant planets at ≳ 30 AU have large eccentricities because they suffered strong gravitational scattering by other giants. In panel (b), only the dynamical friction is included. Eccentricities of a small number of planets are damped, but the effects are not significant. On the other hand, both effects are included in panel (c). Eccentricities are damped to values below 0.2 for ∼30% of giant planets. However, since the damping is not as efficient as the simple treatment in Ida et al. (2013), the fraction of systems that have gas giants with a > 30 AU and e < 0.2 is ∼0.1%, which is smaller than the probability (∼0.4%) in Figure 7 of Ida et al. (2013). It is a future problem to check if such a low fraction is consistent with direct imaging surveys. As already pointed out in Ida et al. (2013), the formation rate of high-eccentricity gas giants at ∼O(1) AU is lower in the theoretical prediction than that found by radial velocity surveys. If the theoretical prediction is improved so that a more frequent formation of high-eccentricity gas giants is reproduced, the theoretically predicted fraction of systems with distant gas giants in nearly circular orbits may also be increased.

Zoom In Zoom Out Reset image size

We have investigated orbital circularization due to planet growth through accreting disk gas. We have analytically derived the differential equations for evolutions of orbital eccentricity e and semimajor axis a, and numerically integrated them to discuss the solutions.

The motivation of these calculations is to examine our scenario for the formation of the distant gas giants in nearly circular orbits, which are being discovered by direct imaging surveys. Our scenario is based on the conventional core accretion model as follows: (1) Icy cores accrete from planetesimals in inner regions at a ≲ 30 AU, (2) they are scattered outward by a nearby gas giant to acquire highly eccentric orbits, (3) their orbits are circularized through the accretion of local protoplanetary disk gas, and (4) through the local gas accretion, the planets become gas giants. We started our calculations after step (2). For given initial e and a, we followed the process in step (3).

For highly eccentric orbits, a planet spends most of its time in the outer regions where disk gas has higher specific orbital angular momentum than the planet. Because the gas accretion rate from the disk is regulated by envelope contraction except for the final gas accretion phase, we assume that disk gas accretion rate is constant within one orbit. Even in the final phase when the accretion rate is limited by the supply of gas due to global disk accretion, the assumption is valid if steady disk accretion is established. Thus, the specific angular momentum of the planet increases with planet accretion, resulting in circularization of the planetary orbit. Energy dissipation by collision between the disk gas and the planet also induces the eccentricity damping.

Just after step (2), the core's pericenter distance must be close to it original location. We found that a pericenter distance is quickly raised by the orbital circularization, so that perturbations of the gas giant in the inner region can be neglected in the orbital circularization process. Thus, we investigated the orbital evolution of isolated planets accreting disk gas.

The orbital evolutions that we found include the following.

• 1.  
  The eccentricity is reduced to <0.2 before the planetary mass is increased by a factor of 10 (for example, if an icy core with ∼10 M_⊕ starts gas accretion, its orbit is circularized with e < 0.2 before it acquires a Saturnian-mass.)
• 2.  
  The eccentricity damping is dominated over the semimajor axis damping. While e is reduced from ∼1 to zero, a is decreased only by a factor of 2.
• 3.  
  These show that planetary growth and orbital circularization concurrently proceed and the orbital circularization is very efficient. If we take into account the effect of bow shock for supersonic incident gas flow, the orbital circularization becomes slower, but it is still efficient enough to account for the observed orbital properties of distant gas giants. The orbit is left in large orbital radii, which are about half of the semimajor axes that the scattered icy cores initially acquire.

We performed the population synthesis calculation by incorporating the fitting formulas for the eccentricity and semimajor axis damping by planet mass growth to show that the damping is efficient and giants with e ≲ 0.2 are left in distant regions at a ∼ 30–300 AU. However, with more detailed prescription using the formulas derived here, the fraction of systems that have such distant jupiters is as small as ∼0.1%, which is lower by a factor of 4 than that predicted in Ida et al. (2013) using simpler prescription.

We also consider the effect of the finite disk size. If the eccentric orbits of the scattered cores are deviated from the protoplanetary disk near their apocenters, their semimajor axes shrink to a quarter of the disk sizes. In other words, if observations show a concentration of distant gas giants at some orbital radius, it could reflect typical sizes of the protoplanetary disks, in a similar way that the pile-up location of hot jupiters could reflect the size of magnetospheric cavity (the size of the disk inner edge) where type II migration could be stalled.

We thank Professor Andrew Youdin for valuable and useful comments as a referee. We also thank Takayuki Muto and Taku Takeuchi for discussions. Our study was supported by JSPS KAKENHI grant No. 23103005.