CROWDING-OUT OF GIANTS BY DWARFS: AN ORIGIN FOR THE LACK OF COMPANION PLANETS IN HOT JUPITER SYSTEMS

Recent observations have found multiple-planet systems whose architectures differ strikingly from that of our Solar System. One of the most distinctive features is a lack of companion planets in hot Jupiter (HJ) systems (e.g., Latham et al. 2011; Steffen et al. 2012). Although some exceptions are being discovered by the Kepler survey, this trend is clearly seen in recent observational results. Several ideas for the inhibition of planet formation in the vicinity of HJs have been presented (Fogg & Nelson 2007; Raymond et al. 2011). We propose a new way of explaining this observational property.

The combination of radial velocity and transit detection observations provides information on the average densities of planets (e.g., Borucki et al. 2010), which shows that a significant fraction of HJs are gaseous giant planets. This clearly means that HJs are surrounded by gaseous components in protoplanetary disks, at least during their formation stages. A natural question in this respect is what happened in the protoplanetary disks in the presence of HJs. It might be interesting to theoretically probe the condition of planet formation in the presence gaseous giants.

We investigate the formation of terrestrial planets in the presence of an HJ. For the origin of HJs, although type II migration (e.g., Lin et al. 1996) and tidal circularization of high-eccentricity planets (e.g., Ford & Rasio 2006; Nagasawa et al. 2008) were proposed, we here introduce a hybrid scenario in which HJs can be formed through gravitational instability prior to the growth of terrestrial planets (Inutsuka 2009). Recent progress in our understanding of protostar formation provides an explanation regarding the formation of protoplanetary disks and their accretion onto the central protostar in their early evolutionary phases (see review by Inutsuka 2012). One of the findings from high-resolution numerical simulations that describe the gravitational collapse of molecular cloud cores indicates that new-born protoplanetary disks are gravitationally unstable and produce multiple planetary mass objects that tend to migrate chaotically in the massive disks (Vorobyov & Basu 2010; Machida et al. 2011). Although the final fate of those planetary mass objects is not entirely clear, the objects that were formed relatively late may survive and remain as gaseous giants in various locations in the disks. Therefore it is natural to consider the planet formation process under the influence of gaseous giants in protoplanetary disks.

This Letter proceeds as follows. In Section 2, we describe our numerical model; in Section 3, we present the results of our simulations. In Section 4, we derive the condition for migration of HJ. In Section 5 we present a discussion.

A HJ with Jovian mass is initially placed at a_HJ = 0.05 AU on a circular orbit, where a_HJ is the semimajor axis of the HJ. Under the influence of the HJ, formation of solid planets from planetary embryos are investigated using numerical simulations that combine the N-body code and the semianalytical code.

The simulation domain of N-body calculations is set between 0.02 and 0.5 AU and the gravitational forces between all of the bodies are calculated (in the same manner as in previous works by Ogihara et al. 2007; Ogihara & Ida 2009, 2012). For the initial solid disk, the solid surface density is assumed as

$$\begin{eqnarray} &&\Sigma _{\rm d} =10 f_{\rm d}\left(\frac{r}{1 {\rm \,AU}}\right)^{-3/2}{\rm \,g \,cm^{-2}},\end{eqnarray} \tag{ 1 }$$

where f_d is a scaling factor and r is the radial distance from the central star. According to this density distribution, planetary embryos are placed with the isolation mass of M_iso = 2πaΔaΣ_d, where Δa is the width of the feeding zone and assumed to be ≃ 10 Hill radii. Note that protoplanets reach the isolation masses for the parameters we are concerned with, unless significant loss of surface density caused by collisional fragmentation occurs (Kobayashi et al. 2010).

The effect of type I migration is considered; thus, protoplanets migrate into the simulation domain from beyond 0.5 AU. We also simulate the growth and migration of protoplanets beyond 0.5 AU not by the N-body –body code but by a semianalytical code in the same manner as in the population synthesis model (e.g., Ida & Lin 2008), where protoplanets migrate after growth. When protoplanets reach the boundary (a = 0.5 AU), the bodies are added to the N-body code.

We incorporate the effects of eccentricity and inclination damping and orbital migration due to gravitational interaction with the gas disk (e.g., Ogihara & Ida 2009; Ogihara et al. 2010). The decay rate of the semimajor axis is (Tanaka et al. 2002)

$$\begin{eqnarray} t_a &=& \frac{1}{C_{\rm I}}\frac{1}{2.7+1.1q(r)} \left(\frac{M}{M_*}\right)^{-1}\left(\frac{\Sigma _{\rm g} r^2}{M_*}\right)^{-1}\left(\frac{c_{\rm s}}{v_{\rm K}}\right)^{2} \Omega ^{-1}, \nonumber\\ \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \, \, &=& 1.6 \times 10^3 C_{\rm I}^{-1} f_{\rm g}^{-1}\left(\frac{2.7+1.1q(r)}{4.35}\right)^{-1}\left(\frac{r}{0.1 {\rm \,AU}}\right)^{3/2}\nonumber \\ &&\times \left(\frac{M}{M_\oplus }\right)^{-1}\left(\frac{M_*}{M_\odot }\right)^{-1/2}\left(\frac{L_*}{L_\odot }\right)^{1/4}{\rm yr}, \end{eqnarray} \tag{ 3 }$$

where C_I and −q(r) denote the type I migration efficiency factor (Ida & Lin 2008) and the gas surface density gradient (q(r) = −dln Σ_g/dln r), respectively.¹ The migration can be halted when the planets are captured in a mean motion resonance with the HJ (e.g., Ogihara & Kobayashi 2013) or they reach a region with a positive density gradient. The latter corresponds to locations on the disk inner edge (e.g., the magnetospheric cavity) or the outer edge of a gap that is opened by the HJ; terrestrial planets in these locations can gain a positive torque from the disk from the steep surface density gradient (Masset et al. 2006). Planets stop their migration at density gaps or resonance locations, whichever are encountered first by the planets (i.e., whichever is at the more distant location from the star).

The unperturbed gas surface density is considered to be of the form

$$\begin{equation} \Sigma _{\rm g} = 2400 f_{\rm \,g} \left(\frac{r}{1 {\rm \,AU}}\right)^{-3/2}\exp \left(\frac{-t}{t_{\rm dep}}\right){\rm \,g \,cm^{-2}}, \end{equation} \tag{ 4 }$$

where f_g is a scaling factor (Hayashi 1981). We usually assume f_g = 1 in our simulations. For the case of f_d/f_g = 1, the solar metallicity is considered. Note that there is a vast range in the metallicity of stars; we vary the values of f_d and/or f_g in order to change the metallicity. The gas disk exponentially dissipates on the depletion timescale of t_dep = 10⁶ yr ≃ 3 × 10⁷T_K, where T_K is the orbital period at 0.1 AU. The unperturbed gas surface density profile is shown by the dashed line in Figure 1(a). We consider a disk inner edge, r_in, at 0.05 AU from the central star and inside of which there is a cavity. The HJ opens up a gap near its orbit; for gap profile computations we use an analytical description developed by Crida et al. (2006), which depends on the disk viscosity. In this study, a relatively high value for the turbulent viscosity is assumed (α = 10⁻²), where α denotes the disk viscosity. We will discuss the dependence of the disk viscosity in M. Ogihara et al. (in preparation). The gas surface density profile, which includes the density gap, is shown by the solid line in Figure 1(a). At the gap edge, q(r) becomes smaller than −2.7/1.1, which causes a positive torque on the planet. The location of q(r) = −2.7/1.1 for the gap edge r_g is shown by the dotted line in Figure 1. The location of the 2:1 resonance with the HJ is shown by the dot-dashed line, which is further from the central star than the gap assuming a_HJ = 0.05 AU. When the HJ comes very close to the central star (≲ 0.03 AU), the locations of the 2:1 resonance and the gap edge move inside r_in (Figure 1(b)).

Zoom In Zoom Out Reset image size

If HJs move to very close-in orbits, the tidal interaction with the star would be important. The migration timescale induced by the tidal torque is (Goldreich & Soter 1966)

$$\begin{eqnarray} t_{a,{\rm tide}} &=& \left|\frac{a}{\dot{a}}\right| =\frac{Q_*}{3 k_{2,*}}\left(\frac{M}{M_*}\right)^{-1}\left(\frac{a}{R_*}\right)^{5} \Omega ^{-1}, \end{eqnarray} \tag{ 5 }$$

where Q_* and k_{2, *} are the tidal dissipation function and the Love number of the central star, respectively. Assuming M_* = M_☉ and R_* = R_☉, the timescale is reduced to

$$\begin{eqnarray} t_{a,{\rm tide}} &=& 6\times 10^{9}\left(\frac{Q_*}{10^6}\right)\left(\frac{k_{2,*}}{0.3}\right)^{-1} \nonumber\\ &&\times \left(\frac{M_{\rm HJ}}{M_{\rm J}}\right)^{-1}\left(\frac{a}{0.03 {\rm \,AU}}\right)^{6.5} \; {\rm \,yr}, \end{eqnarray} \tag{ 6 }$$

which is a strong function of a. We also incorporate this damping in our simulations, adopting a slightly smaller valuer of Q_* = 10⁵ to reduce the computational cost, which does not change our conclusions for Q_* ≲ 10⁶. Note that if HJs reach 0.02 AU, they can fall onto the star within several 10⁸ yr even assuming conservative values of Q_* = 10⁶ and R_* = R_☉. Thus, we set the inner boundary of the computational region at 0.02 AU.

In Figure 2(a), we show the time evolution of the semimajor axis for our fiducial model parameters, where f_d = f_g = C_I = 1 is assumed. The time unit is expressed by the orbital period at 0.1 AU (lower axis) and by yr (upper axis). All bodies are plotted as filled circles connected with solid lines, whereas the HJ is represented by open circles. The dotted line represents the location of q(r) = −2.7/1.1; the inner disk profile is determined by the gap opened by the HJ (t ≲ 9 × 10⁷T_K), whereas it is determined by the inner cavity in the late stage of the simulation. Planetary embryos that are initially placed inside 0.5 AU undergo inward migration. The innermost planet stops its migration when it is captured into the 2:1 mean motion resonance with the HJ at t ≃ 2 × 10⁵T_K. Other migrating protoplanets approach the inner planet, which is captured in the 2:1 resonance, and interact with each other, leading to mergers or captures into mean motion resonances. As a result, several planets relax to a quasi-steady state captured in mutual mean motion resonances, which is called a chain of resonant planets or a "resonant chain."

Zoom In Zoom Out Reset image size

After that, several planetary embryos, whose growth is calculated using the semianalytical code, migrate inward from outside 0.5 AU.² These protoplanets also interact with the planets in the resonant chain, and eventually participate in the resonant chain. Because the innermost planet is located outside the gap edge, the planets in the resonant chain reside in the region of negative density gradient in the gas disk and keep losing the orbital angular momentum by the negative type I migration torque. Through gravitational interaction between the HJ and the planets in the resonant chain, the HJ migrate inward by being pushed in by the resonant chain as a whole, although the HJ itself is in the cavity of the gas disk. We call this phenomenon a "crowding-out" effect. Note that the crowding-out may also occur due to the interaction between an HJ and a massive swarm of planetesimals feeling strong aerodynamical gas drag (Raymond et al. 2006; Mandell et al. 2007).

When the HJ comes to a very close-in orbit (a ≃ 0.03 AU), the location of the 2:1 resonance with the HJ is comparable to the location of the disk inner edge at 0.05 AU (Figure 1(b)). The innermost planet in the resonant chain reaches the disk inner edge and gains positive torque from the disk, which compensates for the net negative torque exerted on the planets in the resonant chain; the inward migration of the HJ due to the crowding-out effect is halted at t ≃ 9 × 10⁷T_K.³

The disk gas is depleted on a timescale of 10⁶ yr ≃ 3 × 10⁷T_K. Then, only the HJ moves inward because of the tidal torque from the host star, which results in widening separation between the HJ and the innermost planet. This leads to a deviation from the 2:1 resonance, and subsequently some planets exhibit orbit crossing. In this run, two giant impact events between planets are observed at t ≃ 2.1 × 10⁸T_K, which result in losses of commensurate relationships. The HJ is lost to the central star,⁴ and two terrestrial planets remain at the end of the simulation. The largest planet has M = 2.3 M_⊕.

Figures 2(b) and (c) show the results with decreased (f_d = 0.1) and increased (f_d = 10) solid surface density, respectively. Although the HJ interacts gravitationally with the resonant chain, the inward migration of the HJ caused by the crowding-out is not seen in the result of decreased solid density (f_d = 0.1) because small planets formed in the less massive disk, feeling weak type I migration torques, cannot effectively reduce the angular momentum of the HJ through resonant interactions. The largest planet has a mass of 0.15 M_⊕ at the end of the simulation. The migration of the HJ due to the tidal torque from the star is also not observed because the HJ has a larger semimajor axis (a ≃ 0.05 AU). Several planets have commensurate relations in the final state. On the other hand, the crowding-out of the HJ is clearly observed in the result of increased density (f_d = 10). In fact, the timescale of inward migration of the HJ is shortened. Relatively large planets with a maximum mass of 24 M_⊕ remain in the final state.

We now derive the condition for inward migration of the HJ due to the crowding-out within the disk lifetime. When the HJ moves close to the central star, it can subsequently be lost to the star because of tidal evolution. Therefore, the condition for crowding-out can simply correspond to the condition for the loss of HJ.

The timescale of the HJ's induced migration due to the crowding-out effect, t_{a, HJ}, which depends on the total mass in planets in the resonant chain (M_chain) and on the migration timescale of planets in the chain (t_{a, chain}), is given by

$$\begin{eqnarray} &&t_{a, {\rm HJ}} = -\frac{J_{\rm HJ}}{2 T_{\rm chain}}\simeq t_{a,{\rm chain}} \left(\frac{a_{\rm HJ}}{a_{\rm chain}}\right)^{1/2}\frac{M_{\rm HJ}}{M_{\rm chain}}, \end{eqnarray} \tag{ 7 }$$

where J_HJ and T_chain are, respectively, the angular momentum of the HJ and the total migration torque exerted on the planets in the resonant chain. a_chain and t_{a, chain} are the typical values for the semimajor axis and the migration timescale of the planets in the resonant chain.

Assuming that all the bodies that migrate to near 0.1 AU are included in the resonant chain, we derive the total mass of the migrating planets M_chain. Because the type I migration speed increases with increasing mass, bodies grow to a certain mass and then start migration. The mass of migrating planets is expressed by the critical mass for migration, which is determined by a balance between the accretion timescale and the migration timescale, or the isolation mass, whichever is smaller. The isolation mass is smaller in the inner region; thus, we can define a boundary a_trans, inside which all the migrating planets are with the isolation mass.

The total mass of the migrating planets is given by

$$\begin{eqnarray} &&M_{\rm chain}\simeq \int ^{a_{\rm max}}_{a_{\rm in}}\Sigma _{\rm d} 2\pi r dr, \end{eqnarray} \tag{ 8 }$$

where a_in is the radius of the inner edge of the initial solid disk, and a_max denotes the maximum semimajor axis, inside which solid bodies can start migration before the disk gas is depleted. We assume a_in = 0 for simplicity. When input parameters (f_d, C_I, and f_g) are specified, a unique value for a_max is determined.

We first consider the case where planets with the isolation mass can only migrate inward; the condition for a_trans > a_max is given by

$$\begin{eqnarray} &&2 \left(\frac{t_{\rm dep}}{10^6 {\rm \,yr}}\right)^{4/3} f_{\rm d}^{27/11} C_{\rm I}^{59/33} f_{\rm g}^{53/33}\lesssim 1, \end{eqnarray} \tag{ 9 }$$

where M_* = M_☉ and L_* = L_☉ are assumed hereafter. Then the maximum semimajor axis a_max is deduced by equating the migration timescale $t_a|_{M=M_{\rm iso}}$ with the disk depletion timescale t_dep. Substituting a_max into Equation (8), we obtain

$$\begin{eqnarray} M_{\rm chain, iso} &\simeq & 0.1\left(\frac{t_{\rm dep}}{10^6 {\rm \,yr}}\right)^{2/3}\left(\frac{f_{\rm d}}{0.1}\right)^{3} C_{\rm I}^{2/3} f_{\rm \,g}^{2/3} \,M_\oplus. \end{eqnarray} \tag{ 10 }$$

Note that this should be overestimated by about a factor of two to five because the integral in Equation (8) is performed from the central star for simplicity. Then the migration timescale of the HJ is approximated by

$$\begin{eqnarray} t_{a, {\rm HJ, iso}} &\simeq & 4\times 10^7\zeta ^{-1} \left(\frac{a_{\rm HJ}}{0.05 {\rm \,AU}}\right)^{1/2}\left(\frac{M_{\rm HJ}}{M_{\rm J}}\right)\left(\frac{t_{\rm dep}}{10^6 {\rm \,yr}}\right)^{-4/3}\nonumber \\ &&\times \left(\frac{f_{\rm d}}{0.1}\right)^{-4} C_{\rm I}^{-7/3} f_{\rm \,g}^{-7/3}\left(\frac{a_{\rm chain}}{0.1 {\rm \,AU}}\right){\rm \,yr}, \end{eqnarray} \tag{ 11 }$$

where a correction factor ζ is introduced. In addition to the overestimation of M_chain, the gas density profile also accounts for this correction. Although we adopt Equation (4) and derive the migration timescale, both Σ_g and q(r) become smaller around r = 0.1 AU. We put all of these corrections into the factor ζ.

When the planets migrate inward from the outer region where the mass is determined by the critical mass for migration (a_trans < a_max), the total mass in the resonant chain is estimated by

$$\begin{eqnarray} M_{\rm chain, crit} &\simeq & 5.6\left(\frac{t_{\rm dep}}{10^6 {\rm \,yr}}\right)^{5/24} f_{\rm d}^{37/32} C_{\rm I}^{5/96} f_{\rm \,g}^{11/96} \,M_\oplus. \quad \end{eqnarray} \tag{ 12 }$$

Note again that this is somewhat overestimated. The migration timescale of the HJ is given by

$$\begin{eqnarray} t_{a, {\rm HJ, crit}} &\simeq & 1\times 10^4\zeta ^{-1}\left(\frac{a_{\rm HJ}}{0.05 {\rm \,AU}}\right)^{1/2}\left(\frac{M_{\rm HJ}}{M_{\rm J}}\right)\left(\frac{t_{\rm dep}}{10^6 {\rm \,yr}}\right)^{-5/12}\nonumber\\ &&\times f_{\rm d}^{-37/16} C_{\rm I}^{-53/48} f_{\rm \,g}^{-59/48}\left(\frac{a_{\rm chain}}{0.1 {\rm \,AU}}\right){\rm \,yr}, \end{eqnarray} \tag{ 13 }$$

where ζ is the correction factor as mentioned above. By comparing this expression with the actual migration timescale observed in simulations, we find ζ ≃ 0.01.

In Figure 3, we plot the condition for migration (t_{a, HJ} < t_dep) due to the crowding out on the f_d–f_g plane for ζ = 1 (gray-hatched region) and 0.01 (dark-hatched region). Figure 3(a) is the case where the full type I migration rate is considered (C_I = 1), whereas Figure 3(b) shows the condition for C_I = 0.1. We find that the parametric region in which the crowding-out is effective is not so limited; in fact, the effect becomes rather common when the gas or solid surface densities are larger than the minimum mass solar nebular model (f_gf_d > 1). The condition for crowding-out is not expected to depend on the metallicity (f_d/f_g). There is a correlation between the stellar metallicity and the probability of a HJ (e.g., Fischer & Valenti 2005), which can be explained if HJs are more likely to form around high-metallicity stars.

Zoom In Zoom Out Reset image size

The origin of the lack of companion planets close to HJs can be naturally explained using our model. When the crowding-out of HJ is effective (e.g., the disk mass (f_df_g) is large), the HJ can be pushed inward by a chain of resonant planets, leading to a collision with the central star. Therefore, the HJ is never observed, whereas other terrestrial planets remain around 0.1 AU. On the other hand, when the crowding-out is not effective (e.g., in a low-mass disk), the HJ and companion planets may coexist. However, their masses and/or sizes are below the detection limit of current surveys. In both cases, systems with HJ and close-in terrestrial planets have little chance of being detected.

Our model provides a theoretical prediction for future observations. That is, if additional planets are found just outside the HJ, the masses of planets are likely to be small. In those systems, the maximum mass for solid planets can be estimated by equating t_{a, HJ} and t_dep as

$$\begin{eqnarray} M_{\rm chain, max} &\simeq & 0.6 \, \zeta ^{-1/2}\left(\frac{a_{\rm HJ}}{0.05 {\rm \,AU}}\right)^{1/4}\left(\frac{M_{\rm HJ}}{M_{\rm J}}\right)^{1/2}\left(\frac{t_{\rm dep}}{10^6 {\rm \,yr}}\right)^{-1/2}\nonumber \\ &&\times C_{\rm I}^{-1/2} f_{\rm \,g}^{-1/2}\left(\frac{a_{\rm chain}}{0.1 {\rm \,AU}}\right)^{1/2} \,M_\oplus. \end{eqnarray} \tag{ 14 }$$

In addition, if many exoplanets that contradict our prediction are detected, it means that some of the assumptions made in our model are invalid. One possibility is that the inner edge of the gaseous disk is much larger than our choice at the formation stage of terrestrial planets, which makes the crowding-out inefficient. Whatever observational results we obtain, the implications of our results may impose some constraints on planet formation theory.

We thank the anonymous referee for useful comments. We also thank Shoichi Oshino, Eiichiro Kokubo, and Yasunori Hori for fruitful discussions. Numerical computations were in part conducted on the general-purpose PC farm at CfCA of NAOJ.