ON THE POSSIBILITY OF ENRICHMENT AND DIFFERENTIATION IN GAS GIANTS DURING BIRTH BY DISK INSTABILITY

We use CHYMERA (Boley 2007), with the particle modifications described below, to run a series of hydrodynamic simulations of gravitationally unstable disks with a significant population of rocks. CHYMERA is an Eulerian code that solves the equations of hydrodynamics, with self-gravity, on a regular, cylindrical grid. The star's motion is solved self-consistently, as described in B2009. The internal energy of the gas accounts for the vibrorotational states of H₂ as described in Boley et al. (2007b). Because we are studying extended disks that have very long orbital timescales, we assume that the ortho–para ratio of H₂ remains in equilibrium for the local temperature. For radiative cooling, CHYMERA can use flux-limited diffusion (Boley et al. 2006), a hybrid ray-tracing+flux-limited diffusion scheme (Boley et al. 2007a), or a local radiative cooling approximation (B2009). The hybrid scheme is the most accurate, but because it solves for radiative transport explicitly, it can become unstable whenever the radiative time step becomes smaller than the hydrodynamics step. Under most conditions, the algorithm is stabilized by limiting the amount of energy a cell can change for a given step. This sacrifices accuracy for speed, but appears to be reasonable for many situations (see Boley 2007). To avoid using limiters, a subcycling version of the hybrid scheme has been developed, but it is still in its testing phase and was not ready for use when this study began. Because we are interested in how the solids distribution behaves in a disk that does fragment, we choose to use the fast and stable local radiative cooling approximation.

In order to suppress artificial fragmentation, we ensure that the Truelove et al. (1997) criterion (see also Nelson 2006) is obeyed at all times by adding a cold pressure floor to our equation of state, i.e., P = max(P_ideal, P_cold). We have chosen to set the minimum Jeans wavelength to be at least four radial grid cells, which sets P_cold = (4Δr)²Gρ²/π. As will be discussed in more detail below, the cold pressure only becomes necessary at the very center of clumps.

For modeling solids, we have augmented CHYMERA to include particles, where each particle represents an ensemble of solids with some radius s. The Epstein and Stokes drags are calculated separately, and an interpolated solution is then found. First, consider the Epstein limit, which is appropriate whenever the Knudsen number Kn = λ/2s ≫ 1. Here, λ = μm_p(πR²ρ)⁻¹ is the mean-free path of a gas particle, where ρ is the gas density with mean molecular weight μ, m_p is the proton mass, and R is the typical gas atom/molecule radius. The internuclear spacing of H₂ is r₀ ∼ 0.74 Å, and the typical diameter for a molecule is ∼3r₀ (Allen 1976), so we take R = 10⁻⁸ cm. This size scale is slightly smaller than what has been used in other work (cf. Rice et al. 2004), but gives a cross section that is within a factor of a few to expected cross sections for atomic and molecular gases (∼10⁻¹⁵ cm²; Allen 1976).

We set the velocity difference between the gas and solids to be δv = v_g − v_p. Following Paardekooper (2007), we write the change of this difference in the Epstein limit as

$$\begin{equation} \frac{d \delta v}{dt}\bigg \vert _{\rm Epstein}=-\left(\frac{8}{\pi }\right)^{1/2}\frac{\rho c_a}{\rho _s s} f_e\delta v, \end{equation} \tag{ 1 }$$

for adiabatic sound speed c_a and internal particle density ρ_s. The term f_e is used to interpolate smoothly between low- and high-Mach speed flows, where the Mach speed here is defined relative to the gas such that $\mathcal {M}=\delta v/c_a$ and f_e = (1 + 9πδv²/(128c²_a))^1/2 (Kwok 1975). It is also convenient to define the stopping time for a particle moving through the gas as (e.g., Weidenschilling 1977)

$$\begin{equation} t_s=\frac{\rho _s s}{\rho c_a}. \end{equation} \tag{ 2 }$$

Assuming that the gas density and sound speed are constant over a time step Δt, Equation (1) can be integrated analytically. Carrying out the integration, we find

$$\begin{equation} \delta v \vert _{\rm Epstein} = 2 \delta v_0 \beta \mathcal {Q} \exp \left(-\alpha \Delta t\right) \left(\mathcal {Q}^2-\delta v_0^2\exp \left(-2\alpha \Delta t\right)\right)^{-1}, \end{equation} \tag{ 3 }$$

where δv₀ is the initial velocity difference, $\alpha =(\frac{8}{\pi })^{1/2} \frac{\rho c_a}{\rho _s s}$, β = (128c²_a/(9π))^1/2, and $\mathcal {Q}=\beta +(\beta ^2+\delta v_0^2)^{1/2}$ (see also Paardekooper & Mellema 2006 for an equivalent expression).

Now consider the Stokes limit, which is appropriate for small Kn. The change in the velocity difference can be written as

$$\begin{equation} \frac{d\delta v}{dt}\bigg \vert _{\rm Stokes}=-3\alpha {\rm Kn} k_d \delta v. \end{equation} \tag{ 4 }$$

The factor k_d depends on the Reynolds number ${\rm Re}=3\left(\frac{\pi }{8}\right)^{1/2}\frac{\mathcal {M}}{\rm Kn}$, where (e.g., Paardekooper 2007)

$$\begin{equation} k_d=\left\lbrace \begin{array}{@{}l@{\quad \quad }l} 1+0.15\ {\rm Re}^{0.687} & {\rm for\enspace\enspace Re\le 500,} \\[3pt] 3.96\times 10^{-6}{\rm Re}^{2.4} &{\rm for\enspace\enspace Re \le 1500,}\\[3pt] 0.11\ {\rm Re} &{\rm otherwise.}\end{array}\right. \end{equation} \tag{ 5 }$$

Each Reynolds regime can be integrated separately, yielding the following solutions:

$$\begin{eqnarray} &&\delta v \vert _{\rm Stokes}\nonumber\\ &&=\left\lbrace \begin{array}{@{}l@{\quad \enspace }l} \delta v_0 e^{ -3\alpha {\rm Kn} \Delta t}(\Theta \vert \delta v_0 \vert ^{0.687} & {\rm for\enspace\enspace Re\le 500,} \\[5pt] \quad(1-e^{-2.061 \alpha {\rm Kn}\Delta t})+1)^{-1.4556}\\[5pt] \delta v_0 (1+7.2 \vert \delta v_0 \vert ^{2.4} {\rm Kn} \alpha \Upsilon \Delta t)^{-0.4167} &{\rm for\enspace\enspace Re \le 1500,}\\[5pt] \delta v_0 \big(1+0.99 \big(\frac{\pi }{8}\big)^{1/2} \frac{\vert \delta v_0 \vert \alpha }{c_a}\Delta t\big)^{-1} &{\rm otherwise.}\end{array}\right.\nonumber\\ \end{eqnarray} \tag{ 6 }$$

For convenience, we have defined $\Theta = 0.15 ((\frac{\pi }{8})^{1/2}\frac{3}{c_a \rm Kn})^{0.687}$ and $\Upsilon =3.96\times 10^{-6}((\frac{\pi }{8})^{1/2}\frac{3}{c_a \rm Kn})^{2.4}$, which are constant over a time step.

Finally, the Epstein and Stokes solutions can be coupled to find an interpolated solution (modified from Woitke & Helling 2003)³

$$\begin{equation} \delta v \vert _{\rm total}=\frac{(3{\rm Kn})^2}{(3{\rm Kn})^2+1} \delta v \vert _{\rm Epstein}+\frac{1}{(3{\rm Kn})^2+1}\delta v \vert _{\rm Stokes}. \end{equation} \tag{ 7 }$$

Using conservation of momentum, this new velocity difference can be used to find the new momentum and velocity for the gas and solids. We do include the back reaction of the particles on the gas in this algorithm.

We should point out that a limitation of our approach is that the analytic solutions for the velocity difference consider only the drag force. However, differences between forces on the gas compared with those on the particles due to pressure require that either these forces be included in the solution of the velocity difference or that the hydrodynamics time step remains smaller than the stopping time of the particles. If these requirements are not met, then the asymptotic velocity drift of the particles through the gas will be modeled poorly. For the particle sizes and density regimes we consider here, this should not be a major difficulty. We present a drift test in the Appendix that addresses this issue.

We use a crude particle-in-cell method for evolving the particles. At the start of the simulation, particles are assigned to cells randomly, weighted by the gas mass in a given cell. Inside the cell, each particle is given a random position, an orbital frequency that is Keplerian, and a random vertical and radial velocity as explained in Section 3.2. Whenever the density of solids is needed, it is calculated by summing the mass of all particles in a given cell and dividing by that cell's volume. As discussed below, this density is included when calculating the disk's self-gravity. For future studies, we would like to include density estimates using a smoothing kernel over the nearest neighbors, but that is not implemented here.

We integrate the particles using a leap-frog method in Cartesian coordinates. During a given hydrodynamic step Δt, the particles are integrated after the second hydrodynamics sourcing step, which occurs after the potential update (see Figure 2.2 in Boley 2007), using the following scheme:

$$\begin{eqnarray} &\displaystyle v_0=v_{-1/2}+{\textstyle\frac{1}{2}}a_{\rm grav} \Delta t_{-1}\nonumber\\ &\displaystyle v_{0}^{\prime }=G_{\rm drag}(\delta v, v_0; \Delta t_{-1})\nonumber\\ &\displaystyle v_{1/2}=v_{0}^{\prime }+{\textstyle\frac{1}{2}}a_{\rm grav} \Delta t_{0}\\ &\displaystyle v_{1/2}^{\prime }=G_{\rm drag}(\delta v, v_{1/2}; \Delta t_{0})\nonumber\\ &\displaystyle x_{1}=x_0+v_{1/2}^{\prime }\Delta t_0.\nonumber \end{eqnarray} \tag{ 8 }$$

Here, v and x are the velocity and position values in Cartesian coordinates. The subscript 0 represents the values at the beginning of the current step. Using the previous time step with the current acceleration for the first update allows for better accuracy during variable time stepping. The velocity is first updated using the acceleration due to gravity. Next, using the updated velocity, a new velocity v' is found by applying the drag force using Equations (3), (6), and (7) and using conservation of momentum (see below), all represented in Equation (8) as G_drag. After the particle velocities are updated, the positions are advanced using v'. This algorithm permits 10⁵ particles to be evolved with a modest performance penalty.

For calculating a_grav, we first apply the exact force on the particle due to the star using Cartesian coordinates. Next, the force due to self-gravity from the disk is calculated on the cell faces of the cylindrical grid. These forces are then interpolated linearly, in the same way as the velocities (see below), to the particle's position in cylindrical coordinates. These interpolated values are then transformed to Cartesian coordinates and added to the star's force. The drag force is calculated by transforming the ith particle's Cartesian velocities to cylindrical components v_r,i, v_z,i, and v_ϕ,i. The initial velocity differences used in Equations (3), (6), and (7) are then calculated by linearly interpolating the gas velocities to the position of the particle. These interpolations, for now, are only one dimensional. For example, regardless of a particle's ϕ or z position within a given cell, the face-centered radial velocity component is used for interpolating to the r position for finding the gas v_r. Likewise, regardless of a particle's ϕ or r position in a cell, the face-centered vertical velocity component is used for interpolating to the z position for finding the gas v_z. When calculating the local gas orbital frequency Ω within the Jth radial cell, we use Ω = Ω_J + (Ω_J+1 − Ω_J−1)(r_i − r_J)/(2Δr), regardless of the particle's ϕ or z position within the given cell, for particle radial position r_i and cell center position r_J. This interpolation scheme is coarse but sufficient for this study (see the Appendix). Using the time step Δt/2, the new δv is calculated for each component. So far, we only have the new velocity differences. We now use conservation of momentum to find the new velocity for the particle and for the gas, as affected by this particle. For example, in the radial direction, we have δv_r,0 = v_r,g,0 − v_r,p,i,0 and initial total radial momentum density S = ρ_gv_r,g,0 + ρ_p,iv_r,p,i,0, where v_r,g,0 is the current gas radial velocity at the interpolated position, ρ_g is the gas density for the cell (cell center value), and v_r,p,i,0 is the ith particle's radial velocity. The density for the ith particle, ρ_p,i, is found by dividing the particle's mass by the volume of the cell that contains the particle. The naughts represent values before the drag from the ith particle is applied. The updated δv_r is then found by using δv_r,0 in Equations (3), (6), and (7). The new velocity in the radial direction for the ith particle is then v_r,p,i = (S − ρ_gδv_r)/(ρ_g + ρ_p,i). The new gas momentum at the location of the particle is then calculated. In the radial and vertical directions, the momentum densities are updated on cell faces by extrapolating the change in the gas momentum density to the cell faces using the same interpolation weights that were used to calculate the gas velocity at the particle's position. This is done to maintain the staggered grid formalism used in the code. Because the angular momentum density is cell centered, the change in the gas's angular momentum is applied directly to the cell in which the particle is contained. Finally, we should make it clear that every time the drag force is applied, the velocities for both the gas and the ith particle are updated. If the ith+1 particle resides within the same cell as the ith particle, then the drag on the ith+1 particle will be affected by the drag update for the ith particle. Waiting to update the gas velocities until all particles have had their drag forces taken into account leads to numerical instabilities when the particle density approaches the gas density.

A principal objective of this study is to address whether clumps that form by disk instability can be enriched at birth. To proceed, we explore two disks, each of which is forced to become highly unstable in order to effect fragmentation in the system. The first initial condition (IC1) is a 0.4 M_☉ disk around a 1.5 M_☉ star. We use a softening for the star of $5\, \rm AU$ as described in B2009. The ICs were created using the analytic approximations described in Boley & Durisen (2008). Unfortunately, this approximation uses Keplerian rotation, assumes a polytropic vertical temperature profile, and has sharp inner and outer boundaries. As a result, the analytic disk was evolved at low azimuthal resolution for about two orbital periods at r ∼ 200 AU to allow the disk to adjust. During this evolution, the disk was irradiated using a temperature profile of $T_{\rm irr}=T_0 (r /\rm {AU})^{-1/2} + 10$ K, with T₀ = 600 K. Figure 1 shows the corresponding Toomre (1964) Q and surface density Σ profiles. Recall that a disk is unstable to gravitational instabilities (GIs) whenever Q = c_aκ/(πGΣ) ≲ 1.7, where κ is the epicyclic frequency and c_a is the adiabatic sound speed (e.g., Durisen et al. 2007). In general, instability does not mean fragmentation, as GIs are capable of reaching a state of self-regulation (e.g., Gammie 2001; Lodato & Rice 2004; Mejía et al. 2005), even when radiation transport is included (Nelson et al. 2000; Boley et al. 2006, 2007a; Stamatellos & Whitworth 2008). However, at least for isothermal disks, fragmentation becomes very likely whenever Q < 1.4 (e.g., Nelson et al. 1998; Mayer et al. 2004). Figure 1 shows that the high irradiation temperature used for the initial model stabilizes the disk against fragmentation, but not against the development of GIs. During the low-resolution evolution, the surface density profile has relaxed to an exponential that roughly follows Σ ∼ 410exp(−r/45AU). This is not a detailed fit, as we are looking at a small range of disk radii, but it does give the basic sense of the surface density profile. Hereafter, simulations that use IC1 are labeled SIM1.

Zoom In Zoom Out Reset image size

At the start of the simulation, SIM1's ICs are loaded onto a high-resolution grid, with r, z, ϕ = 256, 64, 512 zones. The initial radial and vertical cell widths remain unchanged from the relaxation process, with Δr = Δz = 2 AU. The disk is given a 10% random density perturbation and the prefactor T₀ in the irradiation temperature profile is reduced to 130 K. Dropping the irradiation temperature by such a large factor causes a sudden loss of disk stabilization and allows Q to drop below unity temporarily. The decision to change T₀ dramatically is mainly to drive SIM1 to instability violently; however, it does have the physical interpretation of sudden and prolonged shielding of the outer disk by, for example, the inner disk.

The second set of ICs (IC2) uses the B2009 simulation called SIMD. We start our realization of the simulation about half of an orbit before fragmentation (Figure 2). To recapitulate, the disk is 0.33 M_☉ and surrounds a 1 M_☉ star. We do not use softening for the interactions between the star and gas in this simulation. The grid has the same dimensions and physical resolution as the SIM1 grid. The disk is grown by allowing mass accretion from an envisaged envelope at $\dot{M}\sim 10^{-4} M_{\odot }$ yr⁻¹ until 0.33 M_☉ is reached, after which the accretion is abruptly halted. This mass accretion rate is high and is chosen mainly as a computational time consideration (see B2009 for more details). The radial density profile for the infalling mass is set to ρ ∼ r^−0.5 (note that r^−1.5 profiles were also explored in B2009). This shallow density profile causes mass to build up at large radii, which creates a surface density enhancement in the disk's profile (Figure 3) and causes Q to plummet. The surface density profile reflects that, for this particular case, disk formation occurs faster than GIs can activate and smooth out the density structure. During the formation of the disk and its subsequent evolution, an irradiation temperature of 30 K is incident on the disk everywhere. We include this irradiation for the simulation presented here, which we call SIM2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We assume that half of the solid mass is contained in large particles that do not contribute to the opacity and that half goes into a perfectly gas-entrained dust distribution that is solely responsible for the opacity. We use the D'Alessio et al. (2001) opacity tables and only consider the Rosseland mean for simplicity. The grain size distribution is assumed to be a power law, with dn/ds ∝ s^−3.5 and a minimum grain size of 0.005 μm. The maximum grain size in the distribution can be varied, which can have a strong effect on the cooling (e.g., Cai et al. 2006). SIM1 is run with a maximum grain size of 1 mm and with 1 μm, while SIM2 is only run with 1 μm. Hereafter, when we use the term "1 mm opacities," we are referring to a dust grain distribution with a maximum s = 1 mm. The total solids-to-gas mass ratio is assumed to be 0.02 (including ices), with 0.01 in small grains, which sets the opacity, and 0.01 in rock-size solids composed of silicates and ices. We refer to the solids in this manuscript simply as "rocks" because most simulations are run with 10 cm size particles, but rocks can also refer to the solids in the 1 km size run. We choose this simple division between large solids and grains for multiple reasons, which are detailed at the end of this subsection. Note that SIMD in B2009 was assumed to have solar metallicity, so our assumption here reduces the original dust opacity for SIM2 by a factor of two.

Particles are distributed according to a density-weighted Monte Carlo sampling of the disk. The initial particle surface densities are plotted in Figures 1 and 3 for SIM1 and SIM2, respectively, and show good agreement with the gas surface density. We assume that each particle in the simulation represents an ensemble of solids. The particles in any given run are assumed to have a single size population of either 10 cm or 1 km in radius and an internal density of 3 g cm⁻³. For the initial velocity distribution, the particles are placed on Keplerian orbits for the star. The disk potential in the initial orbital frequency was not included because we intended to focus on just 10 cm size particles, which respond strongly to gas drag. The 1 km size run was originally intended to be presented in a complementary follow-up study, but for providing a comparison with the low-drag limit, it was later decided to run SIM2 with km-size particles. To make sure that any difference we find can be attributed to the change in solid size only, we also put the 1 km size particles on Keplerian orbits. In SIM1, the radial and vertical velocities were set to zero and given a random perturbation with a magnitude no greater than 0.01 of the Keplerian speed. The perturbation was included in an attempt to avoid all particles at a given disk radius from passing through the midplane simultaneously. In SIM1, we found the 0.01 perturbation to be too weak, so the perturbation was increased to 0.1 of the Keplerian speed for SIM2. The larger perturbation kept strong caustics from forming during particle midplane passages. However, the fast midplane settling times (∼1 orbit) marginalize this detail for the 10 cm size solids.

For this study, we do not change the fraction of solids that go into dust grains. Although our 50–50 choice is arbitrary, the effects of decreasing or increasing the amount of dust in the system on GI amplitudes have already been addressed in previous research (see Cai et al. 2006, 2008; Boley & Durisen 2008). As more solids remain in a dusty phase, the opacity of the disk will increase, which will tend to stabilize the disk against fragmentation. Likewise, as the opacity is lowered due to more dust going into larger solids, the amplitude of GIs can increase due to shorter cooling times, as long as the emissivity does not approach zero. Another consideration is whether there is a qualitative change in the disk's behavior for different large solids-to-gas ratios, which could occur for different disk metallicities or for different grains/large solids fractions. The 50–50 choice seems to be a reasonable starting point, as it is not at either extreme, and, as we will show, does demonstrate enrichment and differentiation at birth. To begin to understand how the large solids affect the large-scale structure of the gaseous disk, we run one simulation without large solids, while keeping the opacity the same.

We explore two particle sizes, but not a distribution of particle sizes. There are several reasons for this choice. First, we seek to run as few supercomputer simulations as necessary to demonstrate the enrichment of fragments at birth. We feel this can be done well by picking a solid size that is likely to maximize enrichment and differentiation at birth, which is represented by 10 cm rocks (see the Appendix). Particles that are much smaller will remain entrained with the gas, while particles much larger will be unaffected by gas drag. The former does not require additional simulations, as this has been the assumption in most GI hydrodynamics simulations to date. The latter is demonstrated using the 1 km solid size simulation, but it is not expected to lead to significant enrichment (e.g., Helled & Bodenheimer 2010). Instead of exploring a single population, one might argue that a distribution of particle sizes would be more realistic. While this is correct qualitatively, any distribution that we could have chosen would be as equally arbitrary as a single solid size inasmuch as the actual distribution of solid sizes is unknown, particularly during the early stages of disk evolution explored here. Even during later stages, planetesimal growth is a matter of ongoing research (e.g., Johansen et al. 2007; Morbidelli et al. 2009). Consider, for example, the results of the Brauer et al. (2008) study. They find that the particle size of solids can become strongly peaked between 1 and 10 cm at 100 AU for the conditions they explore, similar to what is explored in this manuscript, but only if they ignore collisional destruction as well as inward drift due to gas drag. If they include destruction and drift, they find very little dust growth, which is inconsistent with observations. For example, Greaves et al. (2008) and Ricci et al. (2010) both find evidence for particles with sizes of a few centimeters or larger in the extended regions of disks. Even if we use these observations to infer a power-law distribution for particle sizes up to 10 cm (e.g., Draine 2006), we are still left without knowing the distribution of solids larger than 10 cm, due to observational limits, or whether the distribution appropriate for disks with ages of a few Myr is also appropriate for very young disks. Finally, there is a danger in assuming some distribution of solids without having a model for grain growth and destruction that can be evolved self-consistently. As solid sizes differentiate due to disparate drag timescales, the local population of solids will also evolve. The link between the spatial and size evolution of the solids will alter the differentiation process itself, as well as local opacities. This problem is by no means removed by considering a single particle size, as a local concentration of rocks will likely change the local dust opacity, but a single size does allow as clean of a numerical experiment as possible. Complexity can be added in later studies after the simple case is understood, which we present in detail here.

Overall, four variations of the simulations are run. As is described in the next section, the SIM1 disk is evolved with 1 mm opacities, with 1 μm opacities, and the same 1 μm opacities, but without rocks. We refer to these simulations as SIM1mm, SIM1mu, and SIM1muNoRocks, respectively. SIM1mu and SIM1muNoRocks are both branched from the SIM1mm disk after the first 620 yr, as described below. When referring generally to characteristics of the disk, we only write SIM1. Two realizations of SIM2 are presented in this study, one using 10 cm size particles (simply SIM2) and 1 km size particles (SIM2km), both of which are run with 1 μm opacities.

Due to the short cooling times, the sudden change from T₀ = 600 K to 130 K causes the disk to contract quickly, and after about 620 yr of evolution, the SIM1 disk has contracted enough that it becomes possible to evolve the same disk on a grid with half the initial physical size in r and z. Using direct injection, i.e., no interpolation, the disk is loaded onto a higher resolution grid with Δr = Δz = 1 AU, while keeping the number of radial, azimuthal, and vertical grid cells the same. Figure 4 shows the Q and surface density profiles shortly after being loaded onto the high physical resolution grid. The surface density profile has been altered significantly after reducing the irradiation temperature. The choppy Q profile is due to direct injection, which adds noise to the epicyclic frequency. According to Figure 4, the disk is strongly unstable.

Zoom In Zoom Out Reset image size

SIM1mm is only run for 1260 yr. The surface density evolution is shown in Figure 5. Although very flocculent structure develops, including knot formation, clumps do not form. If, however, the maximum grain size is reduced to 1 μm at 620 yr, which marks the start of SIM1mu, the disk fragments into many clumps. At the time of the opacity switch, dense spiral arms have not yet fully developed. Figure 6 shows the evolution of SIM1mu from 620 yr to 1870 yr. The effect of opacity on the spiral structure is illustrated by the ρ–T phase diagram in Figure 7. For the same temperature, SIM1mu reaches much higher densities than SIM1mm. Figure 8 shows the D'Alessio Rosseland mean opacity for a maximum grain size of 1  μm and 1 mm. For T ≲ 100 K, the 1 mm opacity is considerably more opaque, which does not allow the spiral arms in SIM1mm to cool fast enough for fragmentation to proceed. One might worry, instead, that the cold pressure (see Section 2.1) employed in these simulations could be affecting fragmentation. The thick line in Figure 7 shows where the cold pressure becomes important, assuming an equivalent thermodynamic temperature for the pressure. This threshold is only met in clump cores, so this cold pressure cannot be responsible for the lack of fragmentation in the 1 mm case. It is possible that, after significantly longer evolution, SIM1mm would fragment, too. We chose not to investigate this possibility for two reasons: (1) we wanted to focus our computing resources on disks that we know will fragment; and (2) even this short simulation demonstrates that the grain size in the outer disk can strongly affect clump formation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, before focusing on the evolution of the rocks in these simulations, we compare disk structures between SIM1mu and SIM1muNoRocks (Figure 9). Although the disks are qualitatively similar, differences have developed. This will be discussed in more detail in the next section.

Zoom In Zoom Out Reset image size

By the time spiral arms in the gas form, most rocks have settled into a very thin subdisk, contained within one vertical grid cell. This is consistent with expectations: as discussed in the Appendix, the stopping time for 10 cm rocks is ∼10–100 yr for the gas densities in these disks. Most particles will also make their first midplane passage after 1/4 of an orbit due to the ICs, i.e., about 200 yr at 100 AU. As seen in Figure 6, the spiral arms remain weak until about 780 yr, giving enough time for particles to settle before spiral waves fully develop.

Figure 10 shows the face-on evolution of the gas and particles in SIM1mu. As rocks encounter gas overdensities in spiral arms, they become trapped and concentrated, which is consistent with the results reported by Rice et al. (2004). By the end of the simulation, most rocks are associated with spiral arms or clumps. To determine whether there is an enhancement of solids in these structures, we plot in Figure 11 the ratio of rock surface density to gas surface density for four snapshots. The data are shown using the cylindrical coordinates of the grid, where the abscissa is the r coordinate and the ordinate is the ϕ coordinate, to marginalize the chance that a high peak in the ratio is due to interpolation error. In about one dynamical time, the solids become filtered and caught in the spiral arms, with some regions exceeding a tenth of the gas surface density. The right panel in Figure 11 shows the rock-to-gas ratio for gas and particles in only the midplane cells, which accounts for most of the rock mass. Thin, wispy rock arms dominate the midplane density in some regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to concentrating the rocks into arms, the gas–rock coupling in this highly perturbed disk stops rapid inward drift. Rocks still move inward, but they are locked with the progression of spiral arms. Figure 12 (left panel) shows the azimuthally averaged surface density evolution for the rocks. For the first ∼600 yr of evolution, particles outside r ∼ 70 AU move inward relative to the gas, while particles inside 70 AU move outward and collect in an annular surface density maximum. After ∼600 yr, spiral arms develop in the 70 AU surface density enhancement, and undermine simple inward drift due to nontrivial gas dynamics and the evolving pressure maxima. Note that rocks move inward and outward over tens of AU. The spikes in the rock surface density correspond to gaseous spiral arms and/or clumps. As shown in Figure 12 (right panel), the surface density enhancements in rocks coincide with maxima in the azimuthally averaged gas surface density.

Zoom In Zoom Out Reset image size

Occasionally, features that do not have an obvious gas surface density counterpart form. The most obvious type of feature is a rock bridge between gaseous spiral arms. These structures are created when a gas spiral arm dissipates, but the thin, dense rock arm does not. As the spiral arm dissipates, it does so by diffusing away from the pressure maximum. This would allow for particles that are highly concentrated and not perfectly entrained with the gas to remain concentrated. Bridges mostly have irregular shapes, but rare, straight bridges do form after interacting arms create a temporary leading rock arm that is then sheared into a trailing arm. Knotty structure in the rocks can also form due to buckling of the rock arm after incomplete dissolution of the gaseous arm (see x, y ≈ 50, 50 AU in the 1410 yr panel of Figure 10).

Finally, the strong clumping and density enhancements in the rock distribution seem to be causing noticeable changes in the evolution of the gas. In particular, as seen in Figure 9, SIM1mu shows additional clumping and fragmentation compared with SIM1muNoRocks. Unfortunately, it is difficult to distinguish whether the differences between the simulations are brought about by the chaotic nature of GIs in the nonlinear regime, or whether they are directly related to the solids. This is particularly difficult because variations of a percent could be amplified over a dynamical time for either reason. For example, by removing the solids from the simulation, we instantaneously decrease the mass of the disk by 1%. On the other hand, significant concentration in spiral arms could create small, but important, gravitational potential enhancements, aiding fragmentation. Additional effects, such as the momentum exchange between the gas and the rocks, may have some effect as well (e.g., Johansen et al. 2007). Moreover, a greater fraction of large solids will likely exacerbate the differences, if they can indeed be attributed to solid–gas interactions. This result is a topic by itself and needs to be addressed in a separate study.

SIM2 allows us to explore fragmentation and concentration of solids in a disk that is driven toward instability by mass loading. It is not as violently unstable as SIM1 and has much more ordered spiral arms. In addition, visible, weak spiral structure is already present at the start of the evolution presented here. Because the particles are distributed according to gas mass, the rocks have not fully settled to the midplane by the time dense spiral arms develop, unlike SIM1. This allows for the possibility that the spiral waves might significantly hinder midplane rock settling.

Within about 1/2 orbit at r ∼ 100 AU, an arm becomes very dense and fragments into a clump. The simulation is evolved for an additional 1/2 orbit after fragmentation. At the time the simulation is stopped, a second arm is showing signs of clump formation, but we do not investigate it further here. The surface density evolution of the disk is given in Figure 13. As done for SIM1, the same surface density snapshots are shown in Figure 14, but with the particles superimposed. The rocks respond quickly to the density perturbations, and the spiral arms become enhanced despite the short duration of the simulation, as highlighted in Figure 15. Such rapid concentration indicates that enhancement of solids will occur as soon as spiral arms form. Even with large-scale spiral arms, significant settling to the midplane takes place. Within one orbit, the midplane has a very high solids concentration, with >95% of rocks mass contained in the midplane cells. In some spiral arms, the rock density dominates over the gas density.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 16, we show the azimuthally averaged rock surface density profile for several snapshots. Particles move inward, but are thrown back out to larger radii as the spiral arms increase in amplitude. As seen in SIM1, rocks are concentrated in azimuthally averaged density maxima. In contrast to the 10 cm size particles, the solids in SIM2km do not collect in the spiral gaseous spiral arms (Figure 17). Spiral arms do form in the solids, but because the 1 km size solids respond to gas drag slowly, the particles do not collect in the high-density gas regions. This, in turn, allows for super nebular solids-to-gas ratios outside of the gaseous spiral arms (Figure 18). The enhancement is partly due to the much broader solid spiral arms, causing enhancement in the low-density regions immediately before and after a gaseous arm, and partly due to misalignments between the gaseous and solid arms.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

GIs form dense spiral structures that are very efficient at concentrating solids with sizes ∼10 cm, with surface density enhancements of solids 30 times larger than the average nebula value, in agreement with the Rice et al. (2004) study. Near the midplane, the mass density can become dominated by 10 cm rocks. Super rock-to-gas ratios even occur in SIM2km due to the spiral arms of solids being offset and broader than the gaseous arms. We also find that solids do not simply drift inward; instead, they follow the evolution of the gas spiral structure, including outward migration. Simulations should be extended over much longer periods than studied here so that the asymptotic behavior of GIs and solids can be addressed, but these results do show that solids need not be lost by drag and that large particles could be kept at large radii while GIs operate. In addition, the results of SIM1muNoRocks suggest that solids can affect the gas evolution. As discussed earlier, this result requires further study, as numerical chaos could be contributing to some of the differences that we see between SIM1muNoRocks and SIM1mu. However, SIM2km also showed differences in the evolutions of clumps SIM2:C1 and SIM2km:C1, indicating that more than just numerical chaos is at work.

Clump formation takes place exactly where rock-size solids are concentrated, i.e., the spiral arms. Although the midplane regions of these arms can be highly enriched, disk instability only leads to modest enrichment within the fragments, with a total enrichment of tens of percent over the nebula's average value. This value may change for different rock distributions, but we argue that these results provide an upper limit on the amount of enrichment in fragments, as spiral arms are very efficient at sweeping up 10 cm rocks for the conditions explored here. This is emphasized by SIM2km, which actually shows a depletion in solids relative to the gas mass. However, Figure 18 also shows that there are supernebular solids near the clump for the 1 km conditions. If solid evolution is considered, then some of the depletion seen in SIM2km may be marginalized. Using the limits that the 10 cm and 1 km simulations provide, we can provide an upper limit to the amount of enrichment at birth, assuming that rock concentration does not become even more efficient for larger metallicities or for different opacity/large solids fractions. If we assume that most of the solids grow to 10 cm sizes and that the rock enrichment seen in these simulations is doubled as a consequence, we would expect to have a typical f_RtoG ∼ 0.03, giving a total enrichment f_enr ∼ 2. This suggests that a disk instability planet, although enriched, will have at birth no more than twice the amount of metals compared with its host star, taking our results at face value. If a substantial amount of the protoplanet's envelope is tidally stripped, the final planet could become more enriched than estimated from the initial clump enrichment (see Section 6.4).

Enrichment at birth provides a base level for other enrichment mechanisms, such as planetesimal capture. However, Helled & Bodenheimer (2010) found that enrichment by planetesimals with sizes ⩾1 km is negligible at radii of r ∼ 70 AU for the conditions they studied. Because fragmentation is expected to be much more likely at large radii (e.g., Stamatellos et al. 2007; Stamatellos & Whitworth 2008; Boley 2009; Rafikov 2009; Clarke 2009), we expect enrichment at birth to be the dominant enrichment mechanism. Nevertheless, if rapid migration occurs before dissociative collapse or if fragmentation occurs, however rarely, at 10 s of AU, then planetesimal capture could also enhance disk instability objects.

Grain growth alters the opacity of the disk (D'Alessio et al. 2001) and changes the cooling rates, which in turn affects the behavior of GIs (Cai et al. 2006, 2008). Our results are consistent with this picture, where fragmentation is suppressed in SIM1mm due to the higher opacities. The problem is, unfortunately, complicated by the competing effects of grain growth, which suppresses fragmentation at large radii, and settling, which would reduce the opacity for a large volume of gas. This emphasizes that, to address this problem properly, grain growth models need to be included with the dynamical evolution of solids in order to model fragmentation of gaseous disks correctly. Concerning grain growth, we have assumed for this study that 10 cm or 1 km particles are present by the time the disk fragments. Whether this can be achieved is unclear, but the observations of Greaves et al. (2008) do suggest that 10 cm particles are present in massive, extended disks.

We have assumed that the clump masses calculated here are comparable to the final gas giant/brown dwarf masses. Based on the arguments presented in Boley et al. (2010), we do not expect rapid mass growth after a clump undergoes its second collapse due to dissociation of H₂ (the fragmentation process being the first collapse), which should occur after a few thousand years for these clump masses (see also Helled & Schubert 2008). Because the centrifugal radii of gas fragments are expected to be a few tenths to an AU at disk radii ∼100 AU, a substantial disk will likely form as the fragment collapses. Subsequent mass growth will be mediated by subdisk evolution, which has recently been demonstrated in simulations by Vorobyov & Basu (2010). Moreover, the clump masses calculated here represent objects that will become either gas giants or brown dwarfs with surrounding disks, where the disks themselves may be massive. Global instabilities in the clumps, e.g., the bar instability (Durisen et al. 1986), could also redistribute angular momentum before dissociative collapse, ejecting tens of percent of the clump's mass. If convection is unable to redistribute the centrally concentrated rocks before collapse, then such instabilities would preferentially eject the solid-poor material. However, even this scenario would only increase the enrichment of the clump by tens of percent. Only if most of the clump's initial mass could be stripped, e.g., due to tides, would the enrichment of the clump match that of Jupiter's (f_enr up to 6; see, e.g., Baraffe et al. 2010).

Tidal disruption of clumps with dense rock concentrations is seen in these simulations. For example, SIM1mu:C1 is an 8 M_J clump with 38 M_⊕ of solids at its center and becomes sheared away within half an orbit. Unfortunately, we do not have the resolution to follow core formation, so investigating whether these solid concentrations can separate from the gas (see, e.g., Rice et al. 2006) and form large cores will need to be addressed in a future study. However, it is tempting to speculate whether a third formation channel for gas giants, which we call "core assist plus gas capture," becomes available if cores do form in clumps that eventually become disrupted. In this mechanism, gas fragments behave as factories for the formation of rocky/icy cores. After fragments become tidally disrupted, their cores could be liberated into the disk and begin to capture gas slowly, possibly even in the later and less active phase of disk evolution. This is similar to other hybrid gas giant formation mechanisms where hydrodynamical instabilities accelerate core formation by concentrating solids (e.g., Durisen et al. 2005, Klahr & Bodenheimer 2006). We also acknowledge that core formation through differentiation at birth may be subject to complicated equation of state effects as the pressure and density increase as the clump contracts.

As we have shown here with the 1 mm opacity simulation, fragmentation is very sensitive to opacity, and dense knots or short-lived clumps could be more common than long-lived fragments (see, e.g., Vorobyov & Basu 2010). In this case, core assist plus gas capture may very well represent another channel for in situ formation of gas giants on wide orbits, if the gravitational torques and gas interactions keep the cores at large radii. If, instead, the cores migrate into the inner nebula as they capture gas, they will produce inner gas giants. Interestingly, we note that HR149026b has presented a puzzle for planet formation due to its extremely large core estimate (∼70 M_⊕; Sato et al. 2005). If SIM1mu:C3 were to be sheared away, its concentration of 70 M_⊕ of solids in its central regions could give rise to the formation of such a planet by core assist plus gas capture.

Our results suggest that there may be three formation channels for gas giants planets (core accretion plus gas capture, direct formation by disk instability, and core assist plus gas capture), which could give rise to three populations of gas giants. If this is indeed the case, it is crucial to understand how the formation channels could be distinguished observationally. For this discussion, let us refer to a Population I gas giant (Pop I GG) as one that forms via traditional core accretion plus gas capture, Population II GGs as disk instability gas giants, and Pop III GGs as core assist plus gas capture.

If Pop I GGs are typically as enriched as Jupiter and Saturn, then a metallicity difference between Pop II and Pop I GGs may indeed be observable. Unfortunately, the typical outcome of core accretion has not yet been established, as core accretion simulations follow planet formation up until the runaway gas accretion phase is established, but not the end of accretion and subsequent evolution (e.g., see recent work by Movshovitz et al. 2010). As a result, the final metallicity of Pop I GGs is unclear. To illustrate this point, take the enrichment of heavy elements in Jupiter to be about five times the solar value. If Jupiter had continued to accrete gas with nearly solar composition, it would need to grow to about 4 M_J before its enrichment would drop to twice solar and becoming indistinguishable, using metallicity, from Pop II GGs. Looking for trends in the frequency of GGs as functions of semimajor axis, host star metallicity, and time, as discussed in the introduction, may still be the best way to distinguish Pop I and II GGs.

Pop III GGs suffer the same uncertainty in final metallicity as do Pop I GGs. In the simulations presented here, we unfortunately cannot compute the formation of a solid core, let alone the subsequent runaway gas capture onto that core, and cannot state what the final enrichment of such an object will be. We do note however that, with the large concentration of solids seen in gas clumps, Pop III GGs may be very similar to Pop I GGs, but have very large cores for their mass (e.g., HR149026b). Whether Pop II or III objects can form in a given disk will likely depend on the protostellar core's angular momentum. The fraction of protostars for which Pop II and III GG formation will be important is uncertain, but we remind the reader, as also noted by Stamatellos et al. 2007, that the distribution of angular momentum within dark cores probably does permit the formation of extended disks (Goodman et al. 1993).