A NEW HYBRID N-BODY-COAGULATION CODE FOR THE FORMATION OF GAS GIANT PLANETS

Disk evolution sets the context for our planet formation calculations. Viscous processes within the disk transport mass inward onto the central star and angular momentum outward into the surrounding molecular cloud. Heating from the central star modifies the temperature structure within the disk and removes material from the upper layers of the disk. All of these processes change disk structure on timescales comparable to the growth time for planetesimals. Thus, planets form in a rapidly evolving disk.

For a disk with surface density Σ and viscosity ν, conservation of angular momentum and energy leads to a nonlinear diffusion equation for the time evolution of Σ (e.g., Lynden-Bell & Pringle 1974; Pringle 1981),

$$\begin{equation} \frac{\partial \Sigma }{\partial t} = 3 R^{-1} \frac{\partial }{\partial R} \left(R^{1/2} \frac{\partial }{\partial R} \lbrace \nu \Sigma R^{1/2} \rbrace \right) - \left(\frac{\partial \Sigma }{\partial t} \right)_{{\rm ext}}, \end{equation} \tag{ 1 }$$

where R is the radial distance from the central star and t is the time. The first term is the change in Σ from viscous evolution; the second term is the change in Σ from other processes such as photoevaporation (e.g., Alexander et al. 2006; Owen et al. 2010) or planet formation (e.g., Alexander & Armitage 2009). The viscosity is ν = αc_sH, where c_s is the sound speed, H is the vertical scale height of the disk, and α is the viscosity parameter. The sound speed is c²_s = γkT_d/μm_H, where γ is the ratio of specific heats, k is Boltzman's constant, T_d is the midplane temperature of the disk, μ is the mean molecular weight, and m_H is the mass of a hydrogen atom. The scale height of the disk is H = c_sΩ⁻¹, where Ω = GM_⋆/R³ is the angular velocity.

There are two approaches to solving Equation (1). Analytic solutions adopt a constant mass flow rate $\dot{M} = 3 \pi \nu \Sigma$ through the disk, approximate a vertical structure, and solve directly for Σ(R, t), T_d(R, t), and other physical variables. Numerical solutions either adopt or solve for the vertical structure and then solve for the time variation of $\dot{M}$ and other physical variables (e.g., Hueso & Guillot 2005; Mitchell & Stewart 2010).

Here, we consider planet formation using both approaches. For the analytic solution, we follow Chambers (2009), who derives an elegant time-dependent model for a viscous disk irradiated by a central star. For the numerical solution, we assume that the midplane temperature is derived from the sum of the energy generated by viscous dissipation (subscript "V") and the energy from irradiation (subscript "I")

$$\begin{equation} T_d^4 = T_V^4 + T_I^4. \end{equation} \tag{ 2 }$$

The viscous temperature is

$$\begin{equation} T_V^4 = \frac{27 \kappa \nu \Sigma ^2 \Omega ^2}{64 \sigma }, \end{equation} \tag{ 3 }$$

where κ is the opacity and σ is the Stefan–Boltzmann constant (e.g., Ruden & Lin 1986; Ruden & Pollack 1991). With ν = αc²_sΩ⁻¹ and t₂ = (27α/64σ)(γk/μm_H)κΩΣ², the viscous temperature is

$$\begin{equation} T_V^4 = t_2 T_d. \end{equation} \tag{ 4 }$$

If the disk is vertically isothermal, the irradiation temperature is T⁴_I(R) = (θ/2)(R_⋆/R)²T_⋆, where R_⋆ and T_⋆ are the radius and effective temperature of the central star and

$$\begin{equation} \theta = \frac{4}{3 \pi } \left(\frac{R_{\star }}{R} \right)^3 + R \frac{\partial (H/R)}{\partial R} \end{equation} \tag{ 5 }$$

(Chiang & Goldreich 1997). Thus, the irradiation temperature is

$$\begin{equation} \left(\frac{T_I}{T_\star } \right)^4 = \frac{2}{3 \pi } \left(\frac{R_{\star }}{R} \right)^3 + \frac{H}{2R} \left(\frac{R_{\star }}{R} \right)^2 \left(\frac{\partial {\rm ln} H}{\partial {\rm ln} R} - 1 \right). \end{equation} \tag{ 6 }$$

Following Chiang & Goldreich (1997) and Hueso & Guillot (2005), we set ∂lnH/∂lnR = 9/7. With H = c_sΩ⁻¹, we set t₀ = (2T⁴_⋆/3π)(R_⋆/R)³ and t₁ = (R_⋆/R)²(7RΩ)⁻¹(γk/μm_H)^1/2T⁴_⋆. The irradiation temperature is then

$$\begin{equation} T_I^4 = t_0 + t_1 T_d^{1/2}. \end{equation} \tag{ 7 }$$

Although viscous disks are not vertically isothermal (Ruden & Pollack 1991; D'Alessio et al. 1998), this approach yields a reasonable approximation to the actual disk structure.

Because T_V and T_I are functions of the midplane temperature, we solve Equation (2) with a Newton–Raphson technique. Using Equations (4) and (7), we re-write Equation (2) as

$$\begin{equation} f(T_d) = T_d^4 - \big(t_0 + t_1 T_d^{1/2} + t_2 T_d\big) = 0. \end{equation} \tag{ 8 }$$

Adopting an initial T_d ≈ t^1/3₂ or T_d ≈ t^2/7₁, the derivative

$$\begin{equation} \frac{\partial f}{\partial T_d} = 4 T_d^3 - \frac{t_1}{2} T_d^{-1/2} - t_2 \end{equation} \tag{ 9 }$$

allows us to compute

$$\begin{equation} \delta T_d = f \left(\frac{\partial f}{\partial T_d} \right)^{-1} \end{equation} \tag{ 10 }$$

and yields a converged T_d to a part in 10⁸ in 2–3 iterations.

In the inner disk, the temperature is often hot enough to vaporize dust grains. To account for the change in opacity, we follow Chambers (2009) and assume an opacity

$$\begin{equation} \kappa = \kappa _0 \left(\frac{T_d}{T_e} \right)^n \end{equation} \tag{ 11 }$$

with n = −14 in regions with T_d > T_e= 1380 K (Ruden & Pollack 1991; Stepinski 1998). For simplicity, we assume κ = κ₀ when T_d < T_e.

To solve for the time evolution of Σ, we use an explicit technique with N annuli on a grid extending from x_in to x_out where x = 2 R^1/2 (Bath & Pringle 1981, 1982). We adopt an initial surface density

$$\begin{equation} \Sigma _0 = \frac{M_{d,0}}{2 \pi R R_0} e^{-R/R_0}, \end{equation} \tag{ 12 }$$

where M_d,0 is the initial disk mass and R₀ is the initial disk radius (Hartmann et al. 1998; Alexander & Armitage 2009). For the thermodynamic variables, we adopt γ = 7/5 and μ = 2.4.

Figure 1 compares analytic and numerical results for a disk with α = 10⁻², M_d,0 = 0.04 M_☉, and R₀ = 10 AU surrounding a star with M_⋆ = 1 M_☉. The numerical solution tracks the analytic model well. At early times, the surface density declines steeply in the inner disk (where dust grains evaporate) and more slowly in the outer disk (where viscous transport dominates). At late times, irradiation dominates the energy budget; the surface density then falls more steeply with radius.

Zoom In Zoom Out Reset image size

Figure 2 compares the evolution of the disk mass and accretion rate at the inner edge of the disk. In both solutions, the disk mass declines by a factor of roughly two in 0.1 Myr, a factor of roughly four in 1 Myr, and a factor of roughly 10 in 10 Myr. Over the same period, the mass accretion rate onto the central star declines by roughly four orders of magnitude.

Zoom In Zoom Out Reset image size

In addition to viscous evolution, we consider photoevaporation of disk material. Following Alexander & Armitage (2007), we assume high energy photons from the central star ionize the disk and drive a wind. As long as the inner disk is optically thick, recombination powers a "diffuse wind," where the change in surface density is

$$\begin{equation} \left(\frac{\partial \Sigma }{\partial t}\right)_{{\rm diffuse}} = -2 n_0(R) u_0(R) \mu m_H. \end{equation} \tag{ 13 }$$

In this expression,

$$\begin{equation} n_0(R) = 0.14 \left(\frac{3 \Phi _\star }{4 \pi \alpha _B R_g^3} \right)^{1/2} \left(\frac{2}{(R/R_g)^{15/2} + (R/R_g)^{25/2}} \right)^{1/5}, \end{equation} \tag{ 14 }$$

and

$$\begin{equation} u_0(R) = 0.3423 c_{s,w} \left(\frac{R}{R_g} - 0.1 \right)^{0.2457} e^{-0.3612 (R/R_g - 0.1)}. \end{equation} \tag{ 15 }$$

The recombination coefficient is α_B = 2.58 × 10⁻¹³ cm³ s⁻¹ (Cox 2000). The gravitational radius is R_g = GM_⋆/c²_s,w, where c_s,w = 10 $\rm km\, s^{-1}$ is the sound speed in the wind. For typical stars, the luminosity of H-ionizing photons, Φ_⋆, is ∼10⁴⁰–10⁴² s⁻¹. However, Owen et al. (2010) show that X-rays drive a more powerful wind than lower energy photons. Here, we assume that the mass-loss rate in the wind, $\dot{M}_{{\rm wind}} = \int 2 \pi R \partial \Sigma / \partial t dR$, is a free parameter and consider disk evolution for a range in $\dot{M}_{{\rm wind}}$. This approach is equivalent to assuming a range in Φ_⋆.

As the disk evolves, the diffuse wind removes more and more material throughout the disk. Eventually, the inner disk becomes optically thin to ionizing photons. Ionization then drives a "direct wind." The change in surface density with time is then much larger,

$$\begin{eqnarray} \left(\frac{\partial \Sigma }{\partial t}\right)_{{\rm direct}} &=& 0.47 \mu {\rm m}_H c_{s,w} \left(\frac{\Phi _\star }{4 \pi \alpha _B (H/R) R_{\rm in}^3(t)} \right)^{1/2}\nonumber\\ && \times \left(\frac{R}{R_{\rm in}(t)} \right)^{-2.42}, R > R_{\rm in}, \end{eqnarray} \tag{ 16 }$$

where R_in(t) is the inner edge of the optically thick portion of the disk. Inside R_in(t), the mass-loss rate from the disk is zero.

For mass-loss rates $\dot{M}_{{\rm wind}} \approx 10^{-10}\hbox{--}10^{-8}$ M_☉ yr^-1, the surface density evolution of a viscous disk with photoevaporation follows the evolution in Figure 1 for ∼10⁵–10⁶ yr. As the accretion rate through the disk drops, mass loss from the wind becomes more and more important. Once the inner disk becomes optically thin, the direct wind rapidly empties the inner disk of material. As the system evolves, the size of this inner "hole" grows from R_h ≈ R_g≈ 3 AU to R_h≈ 30–100 AU in 0.01–0.1 Myr (see also Alexander & Armitage 2009; Owen et al. 2010). A few thousand years later, the gaseous disk is gone.

Figure 3 shows the variation of disk mass and mass accretion rate for the disk in Figure 2 and several mass-loss rates. In general, the wind starts to empty the disk when the mass accretion rate through the disk falls below the mass-loss rate in the wind (see Alexander & Armitage 2009; Owen et al. 2010, and references therein). For very low mass-loss rates, $\dot{M}_{{\rm wind}} = 10^{-10}$ M_☉ yr^-1, the wind evaporates the disk on timescales of 10 Myr (see also Alexander & Armitage 2007, 2009). As the mass-loss rate from photoevaporation grows, the disk evaporates on shorter and shorter timescales. For $\dot{M}_{{\rm wind}} = 10^{-8}$ M_☉ yr^-1, the disk disappears on timescales shorter than 1 Myr.

Zoom In Zoom Out Reset image size

These evolutionary models capture the main observable features of disks around pre-main-sequence stars with ages of 1–10 Myr. For 1 Myr old young stars, disk masses of ∼0.001–0.1 M_☉ and mass accretion rates of 1–100 × 10⁻¹⁰ M_☉ yr^-1 are common (Gullbring et al. 1998; Hartmann et al. 1998; Andrews & Williams 2005, 2007). Typical disk lifetimes are 1–3 Myr (Hartmann et al. 1998; Haisch et al. 2001). Few pre-main-sequence stars have ages larger than ∼10 Myr (Currie et al. 2007; Mamajek 2009; Kennedy & Kenyon 2009). Thus, these calculations yield reasonable physical conditions for planet formation.

In our previous calculations, we set the surface density of the gaseous disk as $\Sigma = \Sigma _0 x_m a^{-n} e^{-t/t_g}$, with x_m as a scaling factor, n as a constant power law, and t_g as a constant gas depletion time. In our new approach, adopting an analytic (Chambers 2009) or a numerical (Equation (1)) solution to the radial diffusion equation has several important advantages for little computational effort.

• 1.  
  Following the evolution of the mass accretion rate onto the central star enables robust comparisons between observations of accretion in pre-main-sequence stars with the timing of planet formation throughout the disk.
• 2.  
  Treating the evolution of photoevaporation allows us to make links between the so-called transition disks and the formation of planets (e.g., Alexander & Armitage 2009).
• 3.  
  Tracking the time evolution of the radial position of the snow line (e.g., Kennedy & Kenyon 2008) in a physically self-consistent fashion gives us a way to include changes in the composition of planetesimals with time.
• 4.  
  Calculating the radial expansion of the disk with time provides a way to compare the observed properties of protostellar disks (e.g., Andrews & Williams 2005, 2007; Andrews et al. 2009) with the derived properties of planet-forming disks.

It is straightforward to derive the properties of power-law disks that yield results similar to those of our new calculations. In analytic and numerical solutions to the diffusion equation, the viscous time t_ν ≈ R²/3ν is roughly the gas depletion time t_g. For the Chambers (2009) analytic approach, t_ν≈ 4 Myr (α/10⁻³); for our numerical solution, t_ν≈ 1 Myr (α/10⁻³). In previous calculations, we adopt t_g = 10 Myr; thus, new calculations with α = 1–4 × 10⁻⁴ roughly correspond to the disks in Kenyon & Bromley (2008, 2009, 2010).

The new calculations generally yield shallower power-law slopes—n = 0.6 for Chambers (2009) and n = 1 for our numerical solution—than the n = 1.0–1.5 assumed in our previous studies. Because the timescale for planet formation is t ∝ a^n+1.5 (see Kenyon & Bromley 2008, and references therein), planets form relatively faster in the outer disk in these new calculations. However, the range of disk masses considered here, M_d,0 = 0.01–0.1 M_☉, overlaps the M_d,0 = 0.0003–0.25 M_☉ adopted in Kenyon & Bromley (2010). Thus, calculations with similar initial disk masses and similar initial size distributions of planetesimals will yield similar timescales for the formation of the first oligarchs.

As a concrete example, we compare several disk surface density distributions quantitatively. In Kenyon & Bromley (2008, 2009, 2010), we often adopted an initial surface density distribution of solid material, Σ = 30 g cm⁻² x_m (a/1 AU)^−3/2. At 5 AU, our numerical solutions for disks with M_d,0 = 0.1 M_☉ and R_d,0 = 30 AU have the same surface density as power-law disks with x_m = 3. For the same initial disk mass and radius, the Chambers (2009) analytic model has the same surface density as power-law disks with x_m≈ 2.5. For identical initial conditions and equivalent viscous timescales, these three surface density distributions yield similar evolutionary times for the growth of protoplanets.

Kenyon & Bromley (2008) describe the main details of the coagulation code we use to calculate the growth of small solid particles (planetesimals) into planets. Briefly, we perform calculations on a cylindrical grid with inner radius R_in and outer radius R_out. The model grid contains N concentric annuli with widths δR_i centered at radii R_i. Calculations begin with a mass distribution n(m_ik) of planetesimals with horizontal and vertical velocities h_ik(t) and v_ik(t) relative to a circular orbit. The horizontal velocity is related to the orbital eccentricity, e²_ik(t) = 1.6 (h_ik(t)/V_K,i)², where V_K,i is the circular orbital velocity in annulus i. The orbital inclination depends on the vertical velocity, i²_ik(t) = sin⁻¹(2(v_ik(t)/V_K,i)²).

The mass and velocity distributions of planetesimals evolve in time due to inelastic collisions, drag forces, and gravitational forces. The collision rate is nσvf_g, where n is the number density of objects, σ is the geometric cross-section, v is the relative velocity, and f_g is the gravitational focusing factor (Safronov 1969; Lissauer 1987; Spaute et al. 1991; Wetherill & Stewart 1993; Weidenschilling et al. 1997; Kenyon & Luu 1998; Krivov et al. 2006; Thébault & Augereau 2007; Löhne et al. 2008; Kenyon & Bromley 2008). The collision outcome depends on the ratio of the collision energy needed to eject half the mass of a pair of colliding planetesimals Q^⋆_d to the center-of-mass collision energy Q_c. If m₁ and m₂ are the masses of two colliding planetesimals, the mass of the merged planetesimal is

$$\begin{equation} m = m_1 + m_2 - m_{\rm {ej}}, \end{equation} \tag{ 17 }$$

where the mass of debris ejected in a collision is

$$\begin{equation} m_{\rm {ej}} = 0.5 (m_1 + m_2) \left(\frac{Q_c}{Q_d^*} \right)^{9/8}. \end{equation} \tag{ 18 }$$

This approach allows us to derive ejected masses for catastrophic collisions with Q_c ∼ Q*_d and for cratering collisions with Q_c ≪ Q*_d (see also Wetherill & Stewart 1993; Williams & Wetherill 1994; Tanaka et al. 1996; Stern & Colwell 1997; Kenyon & Luu 1999; O'Brien & Greenberg 2003; Kobayashi & Tanaka 2010). Consistent with numerical simulations of collision outcomes (e.g., Benz & Asphaug 1999; Leinhardt et al. 2008; Leinhardt & Stewart 2009), we set

$$\begin{equation} Q_d^* = Q_b r^{\beta _b} + Q_g \rho _g r^{\beta _g}, \end{equation} \tag{ 19 }$$

where $Q_b r^{\beta _b}$ is the bulk component of the binding energy, $Q_g \rho _g r^{\beta _g}$ is the gravity component of the binding energy, r is the radius of a planetesimal, and ρ_g is the mass density of a planetesimal.

To compute the evolution of the velocity distribution, we include collisional damping from inelastic collisions and gravitational interactions. For inelastic and elastic collisions, we adopt the statistical, Fokker–Planck approaches of Ohtsuki (1992) and Ohtsuki et al. (2002), which treat pairwise interactions (e.g., dynamical friction and viscous stirring) between all objects in all annuli (see also Stewart & Ida 2000). As in Kenyon & Bromley (2001), we add terms to treat the probability that objects in annulus i interact with objects in annulus j (see also Kenyon & Bromley 2004b, 2008).

2.2.1. Small Particles

In most of our previous calculations, we calculate the evolution of particles with radii larger than the "stopping radius," r_s≈ 0.5–2 m at 5–10 AU (Adachi et al. 1976; Weidenschilling 1977; Rafikov 2004). Although they are subject to gas drag, these particles are not well coupled to the gas. Here, we also consider the evolution of smaller particles entrained by the gas (see also Kenyon & Bromley 2009). Following Youdin & Lithwick (2007), we derive the Stokes number

$$\begin{equation} {\rm St} = \frac{r \rho _g \Omega }{\rho c_s}, \end{equation} \tag{ 20 }$$

where ρ = Σ/2H is the gas density. The vertical scale height of small particles is then

$$\begin{equation} H_s = \left\lbrace \begin{array}{@{}cll}H, & & {\rm St} \le \alpha \\ \sqrt{\frac{\alpha }{St}} H, & & {\rm St} > \alpha \\ \end{array}. \right. \end{equation} \tag{ 21 }$$

We assume that entrained small particles have vertical velocity v = H_sΩ and horizontal velocity h = 1.6v.

Estimating accretion rates of small particles by embedded protoplanets is a challenging problem. For particles with St ∼ 1, accretion rates depend on the strength of the planet's gravity relative to forces that couple particles to the gas. Recently, Ormel & Klahr (2010) began to address this issue with an analysis of interactions between protoplanets and particles (loosely) coupled to the gas. After identifying three classes of encounters, they derive growth rates as a function of m, St, and the local properties of the gas. For massive protoplanets with r≳ 1000 km, they derive short growth times, ∼10³ yr, for particles with 10⁻¹ < St < 10.

Here, we adopt a simple prescription for accretion of small particles with St≲ 1. Small protoplanets with r < 100 km do not have strong enough gravitational fields to wrest small particles from the gas. Thus, we assume small protoplanets cannot accrete particles with St < 1. The Ormel & Klahr (2010) results suggest large protoplanets accrete particles with St < 10⁻³ on long timescales, ≳1 Myr. Thus, we also assume protoplanets of any size cannot accrete these particles. For particles with St>10⁻³, we assume accretion proceeds as in our standard formalism with collision rates nσvf_g, where the cross-section σ includes the enhanced radius of the protoplanet discussed in Section 2.2.2. Our approach generally yields growth rates a factor of 1–10 smaller than those of Ormel & Klahr (2010). Thus, our formation times are longer than timescales derived from a more rigorous approach. For the calculations in this paper, however, this accuracy is sufficient to explore the impact of disk mass and viscosity on the formation of gas giant planets. In future papers, a more comprehensive treatment of small particle accretion and gas drag will yield a better understanding of the role of the local properties of the gas on planet formation.

Aside from enabling protoplanets to accrete smaller particles at large rates, including small particles with r≲ 1 m allows us to derive a more accurate calculation for the time evolution of debris from the collisional cascade. Compared to results in Kenyon & Bromley (2008, 2010), estimates for the abundance of very small grains with r∼ 1–10 μm now rely on extrapolation over three orders of magnitude in radius instead of six. Thus, this addition yields more accurate predictions for the variation of grain emission as a function of time (see, for example, Kenyon & Bromley 2008, 2010).

2.2.2. Planetary Atmospheres

As planets grow, they acquire gaseous atmospheres. Initially, the radius of the atmosphere R_a is well approximated by the smaller of the Bondi radius R_B and the Hill radius R_H:

$$\begin{equation} R_a = \left\lbrace \begin{array}{@{}cll}R_B = \frac{G M_p}{c_s^2}, & & R_B < R_H \\ R_H = \left(\frac{M_p}{3 M_\star } \right)^{1/3} a, & & R_B > R_H \\ \end{array}. \right. \end{equation} \tag{ 22 }$$

where M_p is the mass of the planet and a is the semimajor axis of the planet's orbit around the central star. For low-mass planets, R_H > R_B. Requiring R_B > R_P, where R_P is the radius of the planet, an icy planet with a mass density of 1 g cm⁻³ starts to form an atmosphere when M_p ≳ 3 × 10²⁵ g.

Once planets develop a small atmosphere, they accrete small particles more rapidly (Podolak et al. 1988; Kary et al. 1993; Inaba & Ikoma 2003; Tanigawa & Ohtsuki 2010). As small particles approach the planet, atmospheric drag reduces their velocities. Smaller velocities increase gravitational focusing factors, enabling very rapid accretion (see also Chambers 2008; Rafikov 2011). Because the extent of the atmosphere depends on the accretion luminosity, calculations need to find the right balance between the extent of the atmosphere and the accretion rate (and luminosity). Here, we follow Inaba & Ikoma (2003) and solve for the structure of an atmosphere in hydrostatic equilibrium around each planet with M_p ≳ 3 × 10²⁵ g. For each smaller planetesimal with mass m, we calculate the enhanced radius R_E(M_p, m), which replaces the physical radius in our derivation of cross-sections and collision rates.

This approach shortens growth times by factors of 3–10. In Chambers (2006), protoplanets with atmospheres grow 3–4 times faster at 5 AU than protoplanets without atmospheres (see also Inaba & Ikoma 2003). To confirm these results, we repeated the calculations of Kenyon & Bromley (2009) for protoplanets with atmospheres. In a power-law disk with x_m = 5, our previous calculations yielded 10 M_⊕ protoplanets on timescales of ∼3 Myr. Repeating these calculations for protoplanets with planetary atmospheres results in growth times shorter by factors of 5–10, ∼0.3–0.5 Myr. Several numerical tests suggest that different fragmentation laws produce factor of 1.5–2 changes in the growth time for protoplanets with atmospheres (see also Chambers 2006, 2008). Kenyon & Bromley (2009) derived only 10%–20% variations in growth times for a similar range in fragmentation laws. Thus, growth times are more sensitive to fragmentation when protoplanets have atmospheres.

2.2.3. Gas Accretion

As planets continue to grow, they begin to accrete gas from the disk. Although the minimum mass required to accrete gas varies with the accretion rate of the planet, the properties of the disk near the planet, and the distance of the planet from the central star, a typical minimum mass is ∼1–10 M_⊕ (Mizuno 1980; Stevenson 1982; Ikoma et al. 2000; Hori & Ikoma 2010; Rafikov 2006, 2011). Once gas accretion begins, the planet mass grows rapidly until it removes most of the gas in its vicinity, either by depleting the disk entirely or by opening up a gap between the planet and the rest of the disk. The planet then grows more slowly and may start to migrate through the disk (e.g., Lin & Papaloizou 1986; Lin et al. 1996; Ward 1997; D'Angelo et al. 2010).

Here, we approximate the complicated physics of gas accretion and the evolution of the atmosphere with a simple formula (see also Veras & Armitage 2004; Ida & Lin 2005; Alexander & Armitage 2009). The growth rates reviewed in D'Angelo et al. (2010) suggest

$$\begin{equation} \dot{M}_g \approx \dot{M}_0 \frac{M_p}{M_0} \frac{\Sigma \Omega ^2}{c_s} e^{-((\mu - \mu _0)/\sigma _m)^2}, \end{equation} \tag{ 23 }$$

where $\dot{M}_0$ is a typical maximum accretion rate, M₀≈ 0.1 M_J (1 M_J is the mass of Jupiter), μ = log M_p, μ₀ = log M₀, and σ_m≈ 2/3. This functional form captures the rapid increase in the accretion rate from the minimum core mass of 1–10 M_⊕ to larger core masses, the peak accretion rate at masses of roughly 10% the mass in Jupiter, and the decline in the accretion rate once the planet opens up a gap in the disk. Detailed numerical simulations suggest maximum accretion rates of 1 Jupiter mass every 0.01–0.1 Myr (e.g., Lissauer & Stevenson 2007; D'Angelo et al. 2010; Machida et al. 2010, and references therein). For simplicity, we consider $\dot{M}_0$ a free parameter; we set a rate appropriate for a Jupiter-mass gas giant in a disk with Σ = 100 g cm⁻², Ω = 1.78 × 10⁻⁸ s⁻¹ (5 AU), and c_s = 0.67 $\rm km\, s^{-1}$ and use Equation (23) to scale this rate throughout the disk.

To remove accreted gas from the disk in our numerical simulations of disk evolution, we set the change in surface density as

$$\begin{equation} \left(\frac{\partial \Sigma }{\partial t}\right)_{{\rm gas}} = -\frac{\dot{M}_g}{2 \pi R \Delta R}, \end{equation} \tag{ 24 }$$

where R is the semimajor axis of the planet and ΔR is the width of the annulus. Combined with photoevaporation from the central star, this expression yields

$$\begin{equation} \left(\frac{\partial \Sigma }{\partial t}\right)_{{\rm ext}} = \left\lbrace \begin{array}{@{}cll}(\frac{\partial \Sigma }{\partial t})_{{\rm diffuse}} + (\frac{\partial \Sigma }{\partial t})_{{\rm gas}}, & & R_{\rm in} \approx R_\star \\[3pt] (\frac{\partial \Sigma }{\partial t})_{{\rm direct}} + (\frac{\partial \Sigma }{\partial t})_{{\rm gas}}, & & R_{\rm in} \gg R_\star \\ \end{array}, \right. \end{equation} \tag{ 25 }$$

where R_in is the inner edge of the optically thick portion of the disk.

Although we use Equation (23) to derive gas accretion rates for the analytic disk model, we do not remove the accreted mass from the disk. To place some limit on the amount of mass accreted by gas giants, we halt gas accretion when the total mass in gas giants exceeds the total remaining mass in the disk.

2.2.4. Evolution of R and L of Planets

Throughout this phase, the radius of the planet R_p depends on the accretion rate and the properties of the atmosphere (e.g., Hartmann et al. 1997; Papaloizou & Nelson 2005; Chabrier et al. 2007; Lissauer & Stevenson 2007; Baraffe et al. 2009; D'Angelo et al. 2010). Before the planet opens a gap in the disk, accretion is roughly spherical; after the gap forms, the planet accretes material from a thin disk-like structure within the planet's Hill sphere (e.g., Nelson et al. 2000; D'Angelo & Lubow 2008). For planets with masses much larger than Neptune, disk accretion provides most of the gas. Thus, to derive a simple estimate for the evolution of R_p, we solve the energy equation for a polytrope of index n = 1.5 accreting material from a gaseous disk (e.g., Hartmann et al. 1997):

$$\begin{equation} \frac{\dot{R_p}}{R_p} = \left(\frac{7}{3} \eta - \frac{1}{3} \right) \frac{\dot{M_p}}{M_p} - \frac{7}{3} \left(\frac{L}{G M_p^2/R_p} \right). \end{equation} \tag{ 26 }$$

Here, L is the planet's photospheric luminosity in erg s⁻¹; 1 − η, where η ⩽ 1, measures the fraction of the accretion energy radiated away before material hits the planet's photosphere.³

Solving Equation (26) requires a relation for L. When $\dot{M}$ = 0, planets resemble n = 0.5–1 polytropes, with R_p ∝ M^m_p and −1/8 ⩽ m ⩽ 1/10 (e.g., Saumon et al. 1996; Fortney et al. 2010). For this phase, we derive a simple expression from published evolutionary tracks (e.g., Burrows et al. 1997; Spiegel et al. 2010; Baraffe et al. 2003, 2008):

$$\begin{equation} L \approx 10^{-9} \,{L_{\odot }}\left(\frac{M_p}{1\,\ M_J} \right)^{0.25} \left(\frac{R_p}{1\, R_J} \right)^{15}, \end{equation} \tag{ 27 }$$

where R_J is the radius of Jupiter. Accretion energy tends to expand planets considerably (e.g., Pollack et al. 1996; Marley et al. 2007; D'Angelo et al. 2010). To treat this phase, we derive a simple expression from published model atmosphere calculations applied to an n = 1.5 polytrope (Saumon et al. 1996; Fortney et al. 2010, and references therein):

$$\begin{equation} L \approx 10^{-6} \,\ {L_{\odot }}\left(\frac{M_p}{1\, M_J} \right) \left(\frac{R_p}{1\,\ R_J} \right)^{3}. \end{equation} \tag{ 28 }$$

When the planet has a radius R_t ≈ 1.78(M_p/1 M_J)^0.06R_J, these relations yield the same luminosity. Thus, we use Equation (27) for R(M_p)>R_t and $\dot{M}_g >$ 0 and Equation (28) otherwise.

Despite its simplicity, this approach yields reasonable results for L(t) and R_p(t). In the example of Figures 4 and 5, the calculation starts with a 10 M_⊕, 1 R_J planet accreting at a rate of 10⁻⁵ M_J yr^-1 until it reaches a mass of 1 M_J. Before the planet reaches its final mass, the luminosity increases exponentially with time. Once accretion stops, the luminosity declines. Peak luminosity and radius depend on η; planets that accrete hotter material (larger η) expand more than planets that accrete colder material (smaller η). For η≈ 0.3–0.4, our peak luminosities of ∼10⁻⁴ L_☉ are similar to those of more detailed calculations (Pollack et al. 1996; Marley et al. 2007; D'Angelo et al. 2010). At late times, evolution for all η converges on a single track for L(t) and R_p(t). This track yields a radius of 1.03 R_J at 4.5 Gyr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although this approach is not precise, it serves our purposes. The model yields planetary radii with sufficient accuracy (±10%) to derive the merger rate of growing planets from dynamical calculations with our N-body code. The estimated luminosity is also accurate enough to serve as a reasonable starting point for more detailed evolutionary calculations, as in Burrows et al. (1997) or Baraffe et al. (2003). This formalism is also flexible: it is simple to adopt better prescriptions for the luminosity or an energy equation with a different polytropic index.

Bromley & Kenyon (2006) provide a description of the N-body component of our hybrid code. To solve the equations of motion for a set of interacting particles, we use an adaptive algorithm with sixth-order time steps, based either on Richardson extrapolation (Bromley & Kenyon 2006) or a symplectic method (Yoshida 1990; Kinoshita et al. 1991; Wisdom & Holman 1991; Saha & Tremaine 1992). The code calculates gravitational forces by direct summation and evolves particles accordingly. In our earlier version, we performed integrations in terms of phase-space variables defined relative to individual Keplerian frames, following the work of Encke (Encke 1852; Vasilev 1982; Wisdom & Holman 1991; Ida & Makino 1992; Fukushima 1996; Shefer 2002); our latest version can track center-of-mass frame coordinates instead. This feature may be useful in instances where stellar recoil is important. As before, the code evolves close encounters—including mergers— between pairs by solution of Kepler's equations in the pairs' center-of-mass frame. The algorithm also includes changes in eccentricity, inclination, and semimajor axis, derived from the Fokker–Planck and gas drag algorithms in the coagulation code.

The 2006 version and the new version of N-body code can track the force from a dusty or gaseous disk. The older code handles only toy models for the gaseous disk, consisting of a power-law surface density profile and an exponential decay in time. The new code includes the more realistic disk potential derived in Section 2.1. The code represents a disk in annular bins with time-varying surface density specified by the physical model. Although the current code does not treat the lack of gravitational force from gaps in the disk produced by gas giant planets, this extra force is small compared to the forces from the rest of the disk and gas giant planets. Once this feature is included, we can then self-consistently track the dynamics of gas giants in an evolving gaseous disk.

Appendix A provides some background and details of our method to treat the gravity of the gaseous disk. To summarize, we approximate the disk as axisymmetric with a scale height H, which is small compared to the orbital distance (H/a ≲ 0.1). The acceleration that the disk produces on a planet near the disk midplane is then fast to calculate. When the disk has a power-law surface density, a simple analytical expression suffices. If the disk has more complicated radial structure (Section 2.1), then we split it up into constant-density annuli and calculate the contribution of each annulus to the acceleration at some point near the disk midplane. We repeat this calculation for a set of O(10⁴) radial points at the beginning of a simulation, and interpolate with a cubic spline thereafter to find the disk acceleration at the location of the N objects. As the surface density of the disk evolves, the code updates the date for the spline fit.

In addition to gas, the new version of the N-body code includes the gravity of non-physical objects representing dust or planetesimals. As in the SyMBA code of Duncan et al. (1998), we use massless "tracer" particles, which simply respond to the gravitational potential produced by the central star and planets. We use the evolution of the tracers to inform the coagulation code of the de/dt, di/dt, and da/dt induced by massive planets in the N-body code. We also have a second population of massive "swarm" particles, which have mutual interactions with planets but do not interact with either tracers or with each other. The swarm thus consists of super-particles which represent the ensemble of small objects that interact with each other statistically in the coagulation code and which also interact dynamically with planets in the N-body code. Independently of the tracers, we can use the swarm particles to inform the coagulation code of the de/dt, di/dt, and da/dt induced by massive planets in the N-body code. More importantly, we use the mutual interactions of swarm and N-body particles to calculate the response of planets to the surface density and velocity distributions of the swarm. In addition to the tests in Section 3, Bromley & Kenyon (2011) describe some applications of the interaction between swarm and N-body particles.

In terms of comparative computational cost, planets are expensive, tracers are cheap, and swarm particles are in between. Even without mutual interactions, the swarm particles can be a considerable burden in long-term integrations. To lighten this computational load, we evolve the star and planets independently of the swarm during a single coarse-grain time step, allowing the code to take as many substeps as necessary. Then, we interpolate the resulting planet trajectories to calculate the gravitational forces needed to update the swarm. When we use massive swarm particles, we allow the planets to deviate from their interpolated path in response to the smaller objects. This approach is valid only when the net forces on the planets from the swarm vary slowly in time, and when the our simple third-order interpolation of the planet's trajectory over a coarse time step is realistic. With this approach, there is a risk of lower accuracy for the very rare events when some members of the swarm interact with a pair of closely interacting N-body particles. Our third-order interpolation then does not yield a robust representation of the trajectories of these swarm particles. Detailed tests show that reduced accuracy only occurs for the very few swarm particles within 1–2 Hill radii of the pair of N-bodies, has a negligible impact on the trajectory of the N-body particles, and does not influence the outcome (merger or scattering event) of the interaction between the two N-bodies.

We put the new version of the N-body code through the same diagnostics as in Bromley & Kenyon (2006). To assess accuracy and dynamic range, we perform long-term orbit integrations of the major planets, as well as integrations of the planets with their masses scaled by a factor of 50 (Duncan et al. 1998; Bromley & Kenyon 2006). We also track the motion of tight binaries in orbit around a central star. In one limit, we evolve a close pair of Jupiter-mass objects (Duncan et al. 1998); in another, we simulate a surface-skimming satellite above the Earth as it orbits the Sun. We also set up pairs of planets on closely spaced circular orbits and evolve them to see whether we can resolve the critical separation that determines whether their orbits cross⁴ (Gladman 1993). With these tests we verify that the new code can replicate the results illustrated in Figures 1–4 of Bromley & Kenyon (2006).

To save computation time, the new code evolves tracers and swarm particles with lower temporal resolution than with the planets. To verify that the resulting orbits are reasonable, we follow Kokubo & Ida (1995) in simulating gravitational stirring of a planetesimal swarm by a pair of planets. Two 2 × 10²⁶ g planets are separated by 10 R_H and are embedded in the middle of an annulus of 800 2 × 10²⁴ g objects, which is 35 R_H wide and centered at a distance of 1 AU from a 1 M_☉ star. All objects are initially on circular orbits. We evolve this system treating all particles as mutually interacting massive bodies, then with the disk composed of massive swarm particles that interact with the planets only, and finally with a massless tracer disk. Figure 6 shows that all three methods yield output that is in good qualitative agreement with previous calculations (Kokubo & Ida 1995; Weidenschilling et al. 1997; Kenyon & Bromley 2001; Bromley & Kenyon 2006).

Zoom In Zoom Out Reset image size

The final set of tests for the N-body code concerns its ability to identify mergers, in particular between planets and swarm particles or tracers. We first confirm that minor changes to the code maintain its accuracy in high-speed mergers between small massive objects (the "grapefruit test" in Bromley & Kenyon 2006). We then test mergers between planets and tracers or swarm particles explicitly with the more sophisticated method proposed by Greenzweig & Lissauer (1990). We set up (1) a 10⁻⁶ M_☉ planet on a circular orbit at 1 AU and (2) test particles on orbits with eccentricity e = 0.007, inclination i = 0.2, and semimajor axes distributed in two rings (0.977–0.991 AU and 1.009–1.023 AU). We evolve the system through a single close encounter between the planet and each test particle to measure the fraction of test particles accreted by the planet. With the planet's physical radius set to 5300 km, we derive a merger fraction of 0.008 ± 0.001, consistent with previous estimates (Greenzweig & Lissauer 1990; Duncan et al. 1998; Bromley & Kenyon 2006).

As protoplanetary disks evolve, planetesimals become important to the large-scale orbital dynamics of the growing planets. Repeated weak interactions with planetesimals can push a planet radially inward or outward. Lin & Papaloizou (1979) and Goldreich & Tremaine (1980) first quantified aspects of this effect analytically; since then, others have explored it numerically (e.g., Malhotra 1993; Ida et al. 2000; Kirsh et al. 2009). Here, we show that our code can perform similar calculations.

Our first illustration of planet migration follows Kirsh et al. (2009). A 2 M_⊕ planet is embedded in a disk of planetesimals, each having a mass 1/600th of the planet's mass. The disk extends from 14.5 AU to 35.5 AU and has a surface density Σ = 1.2(a/a₀)⁻¹ g cm⁻² where a₀ = 25 AU. Initially, the planet has e = 0; the planetesimals have e_rms = 0.001 and i_rms = 0.5 e_rms, where "rms" is root-mean-squared. The planet stirs nearby planetesimals to large e and ejects some planetesimals from the system (Figure 8). As a reaction to the stirring by the planetesimals, the planet migrates from 25 AU at 10⁴ yr to 22 AU at 6 × 10⁴ yr, leaving excited planetesimals in its wake. Because the planetesimals have finite masses and encounter the planet sporadically, inward migration is not smooth (Figure 9). Our results agree with previous simulations (see Figures 3 and 4 of Kirsh et al. 2009).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our second demonstration of migration follows the work of Hahn & Malhotra (1999). Their disk has Σ ∝ a⁻¹ and total mass of 100 M_⊕ between 10 AU and 50 AU. The disk is composed of 1000 planetesimals with e_rms = 2i_rms = 0.01. With the solar system's four major planets initially packed between 5 AU and 23 AU, the system evolves dramatically in time. Dynamical interactions between the planets and planetesimals spread the orbits significantly. Jupiter migrates inward a small amount; the outer planets migrate outward in semimajor axis. After ∼30 Myr, Neptune reaches a ∼ 40 AU. Figure 10 shows the specifics of the migration with our code. The results compare well with Hahn & Malhotra (1999, lower left panel of their Figure 4).

Zoom In Zoom Out Reset image size

We begin with a discussion of planet formation in a disk composed of gas and Pluto-mass planetesimals. In these calculations, disk evolution scales with α. Disks with α = 10⁻² evolve rapidly. For M_d,0 = 0.1 M_☉, the disk mass declines to 0.07 M_☉ in 0.1 Myr, 0.03 M_☉ in 1 Myr, and 0.008 M_☉ in 10 Myr. Disks with α = 10⁻⁵ evolve slowly. For all initial masses, the disk mass declines by less than 1% in 10 Myr.

In disks with ensembles of Plutos, solid objects grow very slowly. In a disk with M_d,0 = 0.1 M_☉, it takes only ∼50–100 yr for the first Pluto–Pluto collision at ∼3 AU. Due to mutual viscous stirring, this first collision occurs at a modest velocity, e ≈ 10⁻³. Thus, this collision yields an object with a mass of roughly twice the mass of Pluto and some smaller collision fragments. Despite this promising start, it takes ∼0.1 Myr for this object to accrete another 10 Plutos and to reach a mass of roughly 10²⁶ g. Although these large objects accrete the collision fragments efficiently, the fragments have a very small fraction, ∼0.1%–1%, of the initial mass. Thus, accretion proceeds solely by the very slow collisions of large objects.

The dynamical evolution of the growing Plutos is also very slow. At ∼0.5 Myr, the first objects reach masses of ∼3 × 10²⁶ g and are promoted into the N-body code. With typical orbital separations of ≳0.1 AU, these objects are spaced at intervals of ≳10 Hill radii. Thus, they do not interact strongly. After ∼2–3 Myr, the largest objects have masses of 0.1–0.5 M_⊕ and typical orbital separations of 0.1–0.2 AU. A few of these are close enough to interact gravitationally. However, the interactions are slow. With masses much less than an Earth mass, these objects sometimes reach masses of 0.5–1 M_⊕ after 10–20 Myr. These objects will never accrete gas and become gas giant planets.

Figure 11 shows the mass distributions for these calculations at 3 Myr. Independent of α, the maximum planet mass scales with the initial disk mass as

$$\begin{equation} M_{{\rm max}} \approx 0.5 \, {M_{\oplus }}\ \left(\frac{M_{d,0}}{0.1\, M_\odot } \right)^{-1.5}. \end{equation} \tag{ 29 }$$

Because the most massive planets result from several rare collisions, the median mass is much smaller, ∼0.05 M_⊕ for disks with M_d,0≈ 0.02–0.1 M_☉. Disks with smaller initial masses cannot produce planets more massive than ∼0.01 M_⊕.

Zoom In Zoom Out Reset image size

For massive disks, the mass distribution is roughly a power law, with a cumulative probability p ∝ m^−β. Our results suggest that more massive disks have shallower power laws, with β≈ 2 for M_d,0= 0.1 M_☉ and β≈ 4 for M_d,0= 0.04 M_☉.

These calculations demonstrate that disks composed of Pluto-mass objects can never become gas giant planets. This result is not surprising. For a disk with surface density Σ, the accretion time is roughly t ∝ RP/Σ, where P is the orbital period (e.g., Lissauer 1987; Goldreich et al. 2004). Thus, large planetesimals grow more slowly than small planetesimals (see also Kenyon & Bromley 2008). In the limit of large objects growing without gravitational focusing, Equation (56) of Goldreich et al. (2004) suggests that an ensemble of Plutos can produce 10 M_⊕ objects in 200–400 Myr at 3 AU. Our numerical simulations confirm this long growth time.

Disks composed of a few Plutos and many pebbles evolve rapidly. For a model with M_d,0 = 0.1 M_☉, R_d,0 = 30 AU, and α = 10⁻³, it takes ∼1.5–1.8 Myr to promote the first object into the N-body code. Within 0.5 Myr, another 10–15 objects reach the promotion mass of 3 × 10²⁶ g. These objects have small atmospheres; they rapidly accrete the remaining pebbles in their feeding zones. It typically takes ∼0.2 Myr for protoplanets to reach masses of 1–2 M_⊕ and another 0.2–0.5 Myr to begin to accrete gas from the disk. Within a total evolution time of 2.5–3 Myr, some protoplanets grow into gas giant planets.

The disk viscosity parameter α sets the timescale for the early portions of this evolution. In calculations with α = 10⁻², the timescale to produce the first N-body is ∼2 Myr, much longer than the 0.2–0.4 Myr timescale for calculations with α = 10⁻⁴–10⁻⁵. In these models, α sets the scale height of the pebbles (Equation (21)). Gravitational focusing factors grow with smaller scale heights; thus, the Pluto-mass seeds accrete pebbles more rapidly when α is small. Our results suggest that planets grow much more rapidly in disks with α ≲ 10⁻⁴.

Stochastic processes are another feature of the evolution. In calculations with α = 10⁻⁴, it takes 2–5 × 10⁴ yr to produce an ensemble of massive oligarchs (Kokubo & Ida 1998) with M_p ≳ 3 × 10²⁶ g (Figure 12). For a large set of calculations with identical starting conditions, the random timing of oligarch formation and variations in accretion rates for each oligarch lead to a broad range in $N_{o,\rm max}$, the maximum number of oligarchs. Often a rapidly accreting oligarch prevents a neighboring oligarch from accreting pebbles and reaching the promotion mass. Thus, calculations with small $N_{o,\rm max}$ tend to have more massive oligarchs than calculations with large $N_{o,\rm max}$.

Zoom In Zoom Out Reset image size

Mergers of protoplanets are the final important feature of the growth of gas giant planets (see also Thommes et al. 2002; Kenyon & Bromley 2006; Ida & Lin 2010; Li et al. 2010). As protoplanets grow to masses of 5–10 M_⊕, their orbital separations are typically 10–20 Hill radii. Thus, mergers are rare. Once protoplanets start to accrete gas, their separations rapidly become smaller than 5–10 Hill radii. Strong dynamical interactions among pairs of protoplanets begin. Some interactions scatter protoplanets into the inner disk at 1–3 AU; others scatter protoplanets into the outer disk at ≳30–100 AU (see also Marzari & Weidenschilling 2002; Veras et al. 2009; Raymond et al. 2009a, 2009b, 2010). Throughout this chaotic growth phase, mergers rapidly deplete the number of growing protoplanets (Figure 12). Once the number of protoplanets reaches ∼5–10, planets have typical separations of ≳10 Hill radii. Chaotic growth ends.

The timescale for chaotic growth varies from one calculation to the next (Figure 12). Sometimes, protoplanets are widely separated, grow slowly, and reach chaotic growth at late times (indigo and blue curves in Figure 12). When several closely packed protoplanets begin to accrete gas at roughly the same time, their growth initiates an early phase of chaotic growth (magenta and turquoise curves in Figure 12). In nearly all cases, chaotic growth ceases on timescales of 3–10 Myr.

With growth rates sensitive to α, the mass distributions of planets also depend on α (Figure 13). In calculations with M_d,0 = 0.1 M_☉ and R_d,0 = 30 AU, disks with small α produce a variety of gas giants, with masses ranging from 10 M_⊕ to 10–20 M_J. Disks with larger and larger α yield smaller and smaller gas giants. For all α, the probability of a given mass is roughly a power law, p ∝ M⁻ⁿ_P, with n = 0.25–1. Thus, these calculations produce more Earth-mass planets than Jupiter-mass planets (see also Ida & Lin 2005; Mordasini et al. 2009). The exponent in the power law depends on α: models with smaller α produce broader probability distributions with smaller n.

Zoom In Zoom Out Reset image size

For all α, gas giants have a broad range of core masses.⁵ Although gas accretion begins when the core mass is 10 M_⊕, planets continue to accrete solid pebbles. During chaotic growth, protoplanets wander through a large range in orbital semimajor axis a, sweeping up (or scattering) pebbles and smaller oligarchs along their orbits. Mergers of two gas giants also augment the core mass. For disks with M_d,0 = 0.1 M_☉, core masses range from 10 M_⊕ to 100 M_⊕.

Our calculations suggest that Jupiter-mass or larger planets form throughout the disk (Figure 14). For α ≲ 10⁻⁴, Jupiters achieve fairly stable orbits at semimajor axes a_J≈ 2–30 AU. Although most 10 M_⊕ cores form at a∼ 4–10 AU, dynamical interactions scatter protoplanets to a≈ 1–3 AU and to a≳ 15–20 AU (see also Tsiganis et al. 2005). After the scattering events, these planets continue to accrete material and can grow to gas giant planets in their new orbits.

Zoom In Zoom Out Reset image size

Planets with smaller masses have broader distributions of semimajor axes than massive gas giants. For convenience, we divide lower mass planets into "Saturns" with masses of 15 M_⊕ to 1 M_J and "super-Earths" with masses of 1–15 M_⊕. For large α∼10⁻³–10⁻², Saturn-mass planets are often the most massive planet in the system and have orbits with a_S≈ 3–30 AU (Figure 15). When α is small, ≈10⁻⁴–10⁻⁵, Saturns are scattered into orbits with a_S≳ 100 AU, well outside those of Jupiter-mass planets. A few Saturn-mass planets are ejected from the system. Super-Earths have even broader distributions of semimajor axis. All super-Earths avoid regions close to gas giant planets (Figure 16). Although a few super-Earths are scattered into the inner disk, most are scattered into the outer disk. Many super-Earths are ejected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Lower mass disks produce correspondingly lower mass planets. Figure 17 shows predicted mass distributions for planets in a disk with M_d,0 = 0.04 M_☉. Disks with α ≳ 10⁻³ fail to produce Jupiter-mass planets; disks with smaller α rarely produce 3–20 M_J planets. The semimajor axes of these planets are much more concentrated toward the central star. The rare Jupiter-mass planets lie within 10 AU; Saturn-mass planets occupy 3–30 AU. Super-Earths occupy the largest range in a, with a_SE≈ 1–100 AU.

Zoom In Zoom Out Reset image size

Our results suggest that lower mass disks fail to produce gas giants for any α. In disks with M_d,0 = 0.02 M_☉ (Figure 18), disks with small α sometimes produce ice giants with M_p≈ 15–30 M_⊕. However, most planets in these simulations are super-Earths with M_p≈ 1–15 M_⊕. Nearly all of these planets lie at small semimajor axes, a_SE≈ 3–10 AU.

Zoom In Zoom Out Reset image size

Finally, these calculations produce a diverse set of multi-planet systems (Table 1). For any M_d,0, disks with smaller α yield a larger frequency of systems with at least two planets. Defining $M_{p,\rm max}$ as the maximum planet mass, planetary systems with smaller $M_{p,\rm max}$ tend to have more planets. Systems with 4 Jupiter-mass planets are rare: 5%–10% among 0.1 M_☉ disks and less than 1% among lower mass disks. Systems with 4–10 super-Earths are relatively common, with frequencies of 70% or more among 0.01–0.02 M_☉ disks.