The Concentration and Growth of Solids in Fragmenting Circumstellar Disks

Gas-giant planets are predominantly gaseous objects, but likely have solid central cores and overall metallicities closer to stars than terrestrial planets (Helled & Bodenheimer 2010; Kreidberg et al. 2014). Recent observations of the gravitational moments of Jupiter by the Juno mission indicate that it has a substantial solid or partially dissolved core of >10 M_⊕ (Bolton et al. 2017). Such a massive core is taken to be a strong indicator of formation by core accretion, as the accretion of a significant gas envelope requires a planetary core of around 10 M_⊕ (Pollack et al. 1996). Formation through gravitational instabilities (GI), on the other hand, involves the direct collapse of dense, gaseous regions into fragments. These massive, cool disks are most likely when a disk is young and devoid of considerable particle accumulations, neither a core or metallicity over the stellar value is expected through GIs; however, recent simulations suggest that fragments may have atmospheres that are distinct from the disk and envelope from which they form (Boley & Durisen 2010).

Observations of giant exoplanets and their atmospheres have led to the discovery of a number of trends in atmospheric chemical compositions and the development of theories that link these compositions to planet formation. Whether or not the solid material is sequestered in the core can affect the chemical balance of key molecules (Ilee et al. 2017). Atmospheric carbon to oxygen ratios (C/O) might be an indicator of formation location or formation mechanism (Öberg et al. 2011; Brewer et al. 2017) and can be inferred through spectroscopic measurements of transiting planets (Charbonneau et al. 2002) and directly imaged planets (Lee et al. 2013). While transiting planets are typically located at orbits too close to be formed in situ by GI, they still provide useful constraints on expected metallicity values of exoplanet atmospheres (Thorngren et al. 2016; Espinoza et al. 2017). In the case of the HR8799 system, the elemental abundances abundances are most consistent with superstellar C/O ratios in favor of core accretion, but the discrepancy between inner and outer planets suggests different accretion histories and formation pathways are possible (Lavie et al. 2017).

Gravitationally unstable disks may still be capable of producing large gaseous objects with significant solid cores by considering that even in low dust-to-gas environments, dust can collect in fragments (Helled et al. 2006). While spiral arms produced by GI are well-documented locations of particle concentration (Gibbons et al. 2012; Boss 2015), it is still unknown just how enriched a fragment formed therein may become enriched with solids during formation and immediately after, when temperatures are not high enough to evaporate solids. Boley & Durisen (2010) investigated particle concentration in self-gravitating disks with global simulations that included two particle sizes, and found considerable solid concentrations (i.e., cores and embryos), interior and exterior to fragments. With high-resolution shearing-box simulations of fragmentation, we aim to explore this process further.

In the linear analysis of the gravitational collapse of a gaseous disk, stability is defined according to the Toomre parameter (Toomre 1964; Goldreich & Lynden-Bell 1965)

$$\begin{eqnarray}&&Q=\displaystyle \frac{{c}_{{\rm{s}}}{\rm{\Omega }}}{\pi G{\rm{\Sigma }}},\end{eqnarray} \tag{ 1 }$$

where values above Q = 1 represent marginal stability to fragmentation, but non-axisymmetric structures may form up to Q ∼ 1.5 (Durisen et al. 2007). Here, c_s is the isothermal sound speed and Ω = (GM/R³)^1/2 is the orbital frequency corresponding to the stabilizing effects of thermal pressure and shear, respectively, which balance the gravitational collapse of a significantly overdense mass surface density Σ, with constants π and gravitational constant G. The formation of fragments in nonlinear gravitoturbulent evolution requires the fulfillment of the Gammie cooling criterion, whereby the disk cools on sufficiently short timescales such that thermal pressure support is removed on approximately the orbital timescale t_c = βΩ⁻¹, with a critical value around β_crit = 3 (Gammie 2001; Baehr et al. 2017; Deng et al. 2017).

A number of studies have looked into the behavior of solid particles in marginally gravitationally unstable disks, but have mostly focused on particle evolution in gravitoturbulent disks. These range from looking into the collection of particles in vortices (Gibbons et al. 2015) and fragments (Boley & Durisen 2010), shock fronts (Gibbons et al. 2012), and other relations with the GI turbulence (Shi et al. 2016). These find that the particles of typically moderate sizes will collect and form high solid concentrations through aerodynamic drag and gas-dust coupling. Shi & Chiang (2013) considered a concentration of particles well above the typical local dust-to-gas ratios for gas disks so that the particles will naturally fragment even if the gas disk is not gravitationally unstable (i.e., Q_gas ≥ 10). We are interested in the regime where the gas disk is gravitationally unstable (Q_gas ≤ 1), but the dust is not, therefore the gas will fragment but initial dust overdensities are not formed through the direct gravitational collapse of the dust, but by the hydrodynamic interactions between the gas and the dust.

For the study of fragment formation including particles, Boley & Durisen (2010) used global radiation hydrodynamic simulations, investigating the large-scale structure of gravitationally unstable disks with two particle sizes. In these simulations, particles are first concentrated into the spiral arms, reaching surface densities that start contributing to the self-gravity in these regions, warranting the inclusion of particle self-gravity. At this point, the particles and gas then collapse together, and the particles rapidly concentrate at the center and potentially constitute a new mode of planet formation (Nayakshin et al. 2014). Boley & Durisen (2010) ultimately found that there are significant concentrations of small grains at the center of fragments; however, their resolution is insufficient to capture the effects of small-scale turbulence caused by parametric instabilities near the midplane (Riols & Latter 2017).

Because young gravitationally unstable disks are largely dominated by their gas content, the distribution and evolution of solids can often be overlooked or ignored for convenience. A handful of studies have looked into the interaction between the gas and dust in gravitoturbulent disks (Gibbons et al. 2012, 2014, 2015; Shi & Chiang 2013; Shi et al. 2016), and have found varying degrees of solid concentration within disk features such as spirals arms, shocks, and vortices. However, none of these considered the effect of a fragmenting disk on particle distribution and concentration. Thus, this study will focus on the evolution of solids in a disk in which the gas is gravitationally unstable, but the dust, being at a realistic metallicity of the interstellar medium (ISM), is not initially susceptible to direct gravitational collapse. This may help determine the initial metallicities of gas-giant planets formed though disk instability as the size and mass of any significant solid clumps. In this paper we will use

$$\begin{eqnarray}&&Z\equiv {M}_{\mathrm{solid}}/{M}_{\mathrm{gas}},\end{eqnarray} \tag{ 2 }$$

to denote the mass ratio of solid to gas, within the entire simulation domain, unless otherwise stated. For the local dust-to-gas density ratio within a single grid cell we use

$$\begin{eqnarray}&&\epsilon \equiv {\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}.\end{eqnarray} \tag{ 3 }$$

The vorticity of the gas is a measure of the local rotation of a fluid and is a simple diagnostic of turbulent motion, defined as the curl of the velocity field ${\boldsymbol{u}}$: $\omega ={\rm{\nabla }}\times {\boldsymbol{u}}$. In the case of our 2D midplane figures, the value of the vertical component is plotted, which also serves to show the rotation in and around a fragment

$$\begin{eqnarray}&&{\omega }_{z}=\left(\displaystyle \frac{\partial {u}_{y}}{\partial x}-\displaystyle \frac{\partial {u}_{x}}{\partial y}\right)\hat{{\boldsymbol{z}}}.\end{eqnarray} \tag{ 4 }$$

Post-formation enrichment was studied by Helled & Schubert (2008) and Helled & Bodenheimer (2010), assuming that fragments formed initially as pure gas overdensities which subsequently accreted nearby material. While they found enrichment of the atmosphere is possible through the accretion of planetesimals, core formation is excluded due to (1) the high temperatures at the center, which evaporate all solids that do accrete, and (2) convection, which prevents dust settling. However, sublimation temperatures at the fragment core are only reached after tens of thousands of years, a large enough interval such that particles accumulated during formation may settle and form a core by the time sublimation occurs (Nayakshin 2010). Convection is less of a hindrance in the distant protoplanetary disk (∼100 au) where densities and optical depths are lower (Nayakshin 2010). This is the primary focus of the present investigation: attempting to quantify how much material is present immediately after fragment formation that could potentially survive sublimation at the center. However, additional attention will also be paid to the solid content in the atmospheres and its consequences regarding directly observed planets and their formation.

We continue by introducing the important hydrodynamic equations and processes necessary to describe a self-gravitating particle-gas mixture in Section 2, followed by a brief description of the numerical implementation of particles and initial conditions in Section 3. In Section 4 we present the results of our simulations of 3D gravitationally unstable disks with particles of various sizes, including estimates of core masses and fragment metallicities. We conclude by discussing the implications on planet formation theory and observations in Sections 5 and 6.

For this study we conduct 3D hydrodynamic shearing-box simulations of a self-gravitating disk with Lagrangian "swarm" particles embedded in the Eulerian mesh using the Pencil² code. Both gas and particles are treated as self-gravitating and particle back reaction is calculated on the gas by assuming particles are mapped to the grid with triangular-shaped clouds. Local simulations allow the Toomre wavelength ∼2πH to be well-resolved in the radial and azimuthal coordinates (x and y in the linearized coordinates, respectively) and keep the boundary conditions periodic. This Toomre wavelength is not resolved in the vertical direction z, but this has little consequence on fragmentation, which is dominated by radial and azimuthal collapse. Shearing-box simulations use hydrodynamic equations that are linearized and transformed into corotating Cartesian coordinates, where q = lnΩ/lnR = 3/2 is the shear parameter:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial t}-q{\rm{\Omega }}x\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial y}+{\rm{\nabla }}\cdot ({\rho }_{{\rm{g}}}{\boldsymbol{u}})={f}_{D}({\rho }_{{\rm{g}}}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}-q{\rm{\Omega }}x\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial y}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{\boldsymbol{u}}=-\displaystyle \frac{{\rm{\nabla }}P}{{\rho }_{{\rm{g}}}}+q{\rm{\Omega }}{v}_{x}\hat{{\boldsymbol{y}}}\\ \quad -\,2{\rm{\Omega }}\times {\boldsymbol{u}}-{\rm{\nabla }}{\rm{\Phi }}+\displaystyle \frac{\epsilon }{{\tau }_{s}}({\boldsymbol{u}}-{\boldsymbol{w}})+{f}_{\nu }({\boldsymbol{u}}),\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial s}{\partial t}-q{\rm{\Omega }}x\displaystyle \frac{\partial s}{\partial y}+({\boldsymbol{u}}\cdot {\rm{\nabla }})s\\ \quad =\ \displaystyle \frac{1}{{\rho }_{{\rm{g}}}T}\left(2{\rho }_{{\rm{g}}}\nu {{\bf{S}}}^{2}-{\rm{\Lambda }}+{f}_{\chi }(s)\right).\end{array}\end{eqnarray} \tag{ 7 }$$

In Equations (5)–(7), ${\boldsymbol{u}}={({v}_{{\rm{x}}},{v}_{{\rm{y}}}+q{\rm{\Omega }}x,{v}_{{\rm{z}}})}^{T}$ is the gas flow plus shear velocity in the local box, ${\boldsymbol{w}}$ is the particle velocity that imparts a back reaction onto the gas proportional to the local dust-to-gas ratio , ρ_g is the gas density, and the thermodynamic variable is the gas entropy s. Viscous heat is generated by ${ \mathcal H }=2{\rho }_{{\rm{g}}}\nu {{\boldsymbol{S}}}^{2}$, with rate-of-strain tensor ${\boldsymbol{S}}$. Hyperdissipation is applied with the terms f_D(ρ_g), ${f}_{\nu }({\boldsymbol{u}})$, f_χ(s) which for each has the form

$$\begin{eqnarray}&&f(\xi )=\nu ({{\rm{\nabla }}}^{6}\xi ),\end{eqnarray} \tag{ 8 }$$

with constant ν = 2.5 H⁶Ω (Yang & Krumholz 2012). Heat is lost via simple β-cooling prescription

$$\begin{eqnarray}&&{\rm{\Lambda }}=\displaystyle \frac{\rho ({c}_{{\rm{s}}}^{2}-{c}_{{\rm{s,irr}}}^{2})}{(\gamma -1){t}_{{\rm{c}}}}\end{eqnarray} \tag{ 9 }$$

with t_c being given by t_c = βΩ⁻¹ and background irradiation term ${c}_{{\rm{s,irr}}}^{2}$, which ensures fragmentation is a result of growing mass perturbations rather than a transient region of effectively no thermal pressure support. This cooling prescription has no dependence on the optical depth and thus all regions cool with the same efficiency. In reality, the opacity will be dominated by the dust grains and an increase of the particle density will increase the local cooling timescale.

Self-gravity is solved in Fourier space by transforming the surface density to find the potential at wavenumber k and transforming the solution back into real space. The solution to the Poisson equation in Fourier space at wavenumber k is

$$\begin{eqnarray}&&{\rm{\Phi }}({\boldsymbol{k}},t)=-\displaystyle \frac{2\pi G\rho ({\boldsymbol{k}},t)}{{{\boldsymbol{k}}}^{2}},\end{eqnarray} \tag{ 10 }$$

where ${\rm{\Phi }}={{\rm{\Phi }}}_{{\rm{g}}}+{{\rm{\Phi }}}_{{\rm{d}}}$ and $\rho ={\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}$ are the potential and density of the gas plus dust particles combined.

Finally, we use an ideal equation of state, with internal energy ε, and specific heat ratio γ

$$\begin{eqnarray}&&P=(\gamma -1)\rho \varepsilon .\end{eqnarray} \tag{ 11 }$$

Particles are added such that the initial distribution maintains a physically motivated metallicity of 1:100, roughly that of the ISM, although this may be an overestimate of the gas content of the disk (Ansdell et al. 2016; Miotello et al. 2017). Stars may form from clouds of a wide range of metallicities and this may have an impact on planet formation, so deviations from this value will be investigated later. Particles evolve as a collection of solids, also known as a superparticle, i, with position ${{\boldsymbol{x}}}^{(i)}$ and velocity ${{\boldsymbol{w}}}^{(i)}$ as in Youdin & Johansen (2007):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{{\boldsymbol{w}}}^{(i)}}{{dt}} & = & 2{\rm{\Omega }}{w}_{y}^{(i)}\hat{{\boldsymbol{x}}}-\displaystyle \frac{1}{2}{\rm{\Omega }}{w}_{x}^{(i)}\hat{{\boldsymbol{y}}}-{\rm{\nabla }}{\rm{\Phi }}\\ & & +\,\displaystyle \frac{1}{{\tau }_{s}}\left({{\boldsymbol{w}}}^{(i)}-{\boldsymbol{u}}({{\boldsymbol{x}}}^{(i)})\right),\end{array}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{x}}}^{(i)}}{{dt}}={{\boldsymbol{w}}}^{(i)}-\displaystyle \frac{3}{2}{\rm{\Omega }}{x}^{(i)}\hat{{\boldsymbol{y}}},\end{eqnarray} \tag{ 13 }$$

where τ_f (also referred to as the Stokes number) is the particle stopping time normalized by the dynamical time Ω⁻¹

$$\begin{eqnarray}&&\mathrm{St}={\tau }_{f}={\tau }_{{\rm{s}}}{\rm{\Omega }}.\end{eqnarray} \tag{ 14 }$$

This gives a hydrodynamic sense of the particle size in the Epstein regime, where the stopping time is defined as (Weidenschilling 1977)

$$\begin{eqnarray}&&{\tau }_{{\rm{s}}}=\displaystyle \frac{a{\rho }_{\cdot }}{{c}_{{\rm{s}}}{\rho }_{{\rm{g}}}},\end{eqnarray} \tag{ 15 }$$

where a is the particle diameter and ρ_· is the material density of a superparticle. Larger particles are less coupled to small-scale gas motions, retain their initial perturbations for longer, and thus have higher Stokes numbers. As a corollary, smaller particles are well-coupled and will quickly take the form of the gas motions in its vicinity and thus have low Stokes numbers. Particle mass is calculated from the gas through the initial condition of the metallicity.

We include no radial pressure gradient, so the radial migration of particles is not included. Because we include particles at the ISM dust-to-gas ratio  = 10⁻², the contribution of particles to the gravitational potential is initially miniscule, thus the potential is dominated by the gas distribution. All simulations include a million particles that are initially uniformly distributed throughout the domain, but rapidly settle to the disk midplane and concentrate with the gas overdensities. Because the gas is stratified and the particles are not, the initial dust-to-gas values near the vertical boundaries will be higher than the overall metallicity of 1/100 and the midplane values will be lower, although the overall metallicity is still at the desired value.

All simulations were run using the high-order finite-difference code Pencil (Brandenburg 2003) which is well-suited for turbulent flows. The code is stabilized through a six-order hyperdissipation scheme (Yang & Krumholz 2012), ensuring power is preserved at large scales, but numerical noise is damped at small scales. The 3D gas-only simulation shown in Figure 1 and included here for reference is based on those in Baehr et al. (2017) but altered to match the dimensions of the other runs presented here; more detailed descriptions of the numerical methods and initial conditions can be found therein.

Zoom In Zoom Out Reset image size

3.1. Particles in Pencil

3.2. Initial Conditions

Initiated with Q₀ = 1, all simulations quickly proceed to collapse, with particles settling to the midplane at different rates depending on their Stokes number. Smaller particles, i.e., St = 0.01, will take longer as minute gas disturbances will deflect their motion toward the midplane, while larger particles, St > 10, will be able to fall largely unimpeded. Figure 2 shows the evolution of the maximum gas and particle densities with time for particle sizes (from top to bottom) St = [0.1,1,10]. Gas densities above ∼100ρ₀ are considered to be the condition for fragmentation. One can see a clear trend toward higher particle concentrations with smaller particle size.

Zoom In Zoom Out Reset image size

The rapid concentration of particles within the collapsing axisymmetric and non-axisymmetric structures can be seen in the different panels of Figure 3. Each panel shows the collapse of the gas and particles for a different particle size, with background color showing the gas surface density with particle surface density overplotted in semi-transparent grayscale. Even though the particles are intended to be uniformly distributed, in practice particles are first evenly distributed among the processors and then within the grid cells of that processor. This introduces artificial gaps between groups of particles, depending on which processor they were originally allocated. These gaps become exaggerated by the large discrepancy between the number of total grid cells and particle number and the shear motions during the linear collapse, apparent in the top left panel of Figure 3. The particles are not unstable to gravitational collapse when these features are present and thus are not expected to cause any premature fragmentation.

Zoom In Zoom Out Reset image size

As one might expect, the particles that are most coupled to the gas show the most concentration in the gas overdensities. In the case of our largest particles (St = 10), the inability to couple strongly with the gas results in a steadily decreasing dust-to-gas ratio, unlike in the case of the smaller particles. This is similar to the simulations with 1 km sized rocks from Boley & Durisen (2010), which showed the delayed response of the larger particles to the gas drag resulted in wide particle arms that were in different positions than the gas spiral arms. This could lead to scenarios where solids can concentrate outside of the gas overdensities, enriching very different regions of the disk.

Our low-resolution runs showed particle concentration within fragments for particle sizes St = [0.01,0.1,1], but in the case of St = 0.01, the core is delayed but only because fragmentation of the gas, roughly defined when the maximum gas density is 100ρ₀, is also delayed. This may indicate that there is some dependency of particle fragmentation on the gravitational collapse of the gas, as might be expected from the initial conditions of the gas and dust (gravitationally unstable and stable, respectively).

In our high-resolution simulations, we again find that small to intermediate particles eventually collect the best within the collapsing fragment and eventually result in at least moderate cores. This is apparent in the radial profiles of the gas fragments, which form in our simulations and are shown in Figure 4. The gas (teal line) shows a roughly Bonnor–Ebert density profile, with constant density within a radius of ∼0.1 H and a power-law decline thereafter. This is consistent with the isothermal temperature structure assumed for Bonnor–Ebert spheres, which here is shown with the dark blue line.

Zoom In Zoom Out Reset image size

Regardless of the particle size, we still find that the particles will concentrate compared to the background distribution, coinciding with the location of the gas fragment. This is apparent in the top left and bottom right panels of Figure 4, which correspond to our largest and smallest particle species, and are least likely to concentrate strongly. While no core is formed, there is still increasing particle density toward the center of the fragment.

In the case of St = 0.1, the core is especially pronounced from the surrounding protoplanetary envelope, and is able to form a distinct particle clump within the fragment itself. This constitutes the most significant core we form in our simulations and as we detail in the following chapter, is tens of Earth masses in size. The fate of core formation in the simulation P1024t2pss, with the highest resolution and smallest particles, does not fragment before the time step becomes prohibitively small, but considering that the lower-resolution simulation formed a small but late core, it is reasonable to expect a core is possible in the high-resolution case as well.

The simulation P512t2ps represents an unusual case where the disk fragments and the particles proceed to collapse into a solid clump, but the gas dissipates shortly thereafter, resulting in the peculiarly high metallicity of the clump in that case. Something similar was seen in Boley & Durisen (2010), where fragments are tidally stripped, leaving a bare core that is proposed as an alternative planet formation pathway, in which tidally stripped cores are free to slowly accrete gas (Nayakshin 2010). This was, however, an isolated incident and none of the other simulations were observed to behave in this way.

4.1. Masses and Solid Concentrations

Importantly, from these simulations we want to determine the size and mass of the solid clumps and how enriched the atmospheres of gas giants or brown dwarfs that form are. First, we need to define our initial parameters in physical units consistent with a shearing box located 100 au from a 1 M_⊙ star. The orbital frequency at 100 au is 1.98 × 10⁻¹⁰ s⁻¹ and the sound speed of a 10 Kelvin monatomic gas (with γ = 5/3) is 262 ${\rm{m}}\,{{\rm{s}}}^{-1}$. The Toomre surface density Σ_T, which is the local surface density required to result in gravitational collapse for a given unstable Q value, can be used with our initial conditions to establish the initial surface density in the simulation

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{T}},0}=\displaystyle \frac{{c}_{{\rm{s}},0}{\rm{\Omega }}}{\pi {{GQ}}_{0}}=37\,{\rm{g}}\,{\mathrm{cm}}^{-2}.\end{eqnarray} \tag{ 16 }$$

This is roughly an order of magnitude greater surface density than the minimum mass solar nebula density at this radius.

We calculate the masses by determining the unit density at the desired orbital radius and at the desired orbital radius and (multiply this value by the total mass in code units of the gas or solids.) From Σ_T we calculate our unit density from the Gaussian profile of the initial vertical gas mass distribution

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{T}}}\equiv {\int }_{-H}^{H}{\rho }_{0}{e}^{-{z}^{2}/2{H}^{2}}{dz},\end{eqnarray} \tag{ 17 }$$

where ${\rho }_{0}=0.29\hat{\rho }$. Solving for our unit volume density $\hat{\rho }$ gives a value of 5 × 10⁻¹² ${\rm{g}}\,{\mathrm{cm}}^{-3}$, with a disk aspect ratio of H/R = 0.1. The mass is added up in code units and averaged at all radii around the densest point before being integrated over radial shells and multiplied by $\hat{\rho }$ to calculate the total mass.

The results of these mass calculations, taken at times shortly after the fragmentation, are shown in Table 1. The profiles of the fragments in the table are shown in Figure 5, as the cumulative mass is shown as a function of radius of the fragment r, beginning from the outside and going inward. Since the density of the mass of particles that constitute the solid core fall well short of the density of a solid body of similar mass, around 5 ${\rm{g}}\,{\mathrm{cm}}^{-3}$, these are not true cores at a solid material density, but a cloud of particles at the size of the grid resolution ∼0.02 H. Beyond this point, any further concentration of the particles is limited by the resolution and can be modeled by other means (see Nayakshin 2018).

Zoom In Zoom Out Reset image size

Table 1.  Performed Simulations and Their Various Resolutions, Particle Sizes St, Initial Metallicities Z₀, Fragment Metallicity Z/Z₀, and Fragment and Core Masses

Simulation	Grid Cells	Particle Number	Stokes #	Z0	Z/Z0	Fragment Mass (MJup)	Core Mass (M⊕)
G512t2	5122 × 64	⋯	⋯	⋯	⋯	14.8	⋯
P512t2pss	5122 × 64	106	0.01	10−2	3.4	5.24	57.9
P512t2ps	5122 × 64	106	0.1	10−2	14.6	1.32	62.3
P512t2p	5122 × 64	106	1	10−2	1.8	16.2	89.1
P512t2pl	5122 × 64	106	10	10−2	1.7	2.3	⋯
P1024t2pss	10242 × 128	106	0.01	10−2	1.7	⋯	⋯
P1024t2ps	10242 × 128	106	0.1	10−2	5.1	4.8	69.0
P1024t2p	10242 × 128	106	1	10−2	3.8	5.7	70.7
P1024t2pl	10242 × 128	106	10	10−2	1.5	7.19	⋯
P1024t2psl	10242 × 128	106	0.1	8 × 10−3	5.0	3.83	49.3
P1024t2psm	10242 × 128	106	0.1	1.33 × 10−2	4.8	8.01	147.6

Download table as:  ASCIITypeset image

To get a rough estimate of the maximum initial mass expected for these simulations, one assumes that all mass within a critical wavelength λ_T = 2πH

$$\begin{eqnarray}&&M\approx {\lambda }_{{\rm{T}}}^{2}{{\rm{\Sigma }}}_{{\rm{T}},0}=17.2{M}_{\mathrm{Jup}},\end{eqnarray} \tag{ 18 }$$

which is the maximum mass our fiducial particleless simulation ever reaches. Because not all material will collapse into the fragment, a lower limit can be approximated by a factor of 4 decrease (Kratter & Lodato 2016). The majority of the fragments form with initial masses within this 4.3–17.2 M_Jup range and are generally closer to the lower-mass end of this range when including particles with back reaction.

For the purposes of determining the local gravitational instability of the clumps in Figure 4, we also define a midplane Toomre volume density ρ_T, from the 3D Toomre parameter Q_3D of Mamatsashvili & Rice (2010)

$$\begin{eqnarray}&&{\rho }_{{\rm{T}}}=\displaystyle \frac{{{\rm{\Omega }}}^{2}}{\sqrt{2\pi }\pi {{GQ}}_{0}},\end{eqnarray} \tag{ 19 }$$

where we have used ${\rm{\Sigma }}=\rho H\sqrt{2\pi }$. This sets a threshold midplane mass above which the disk is expected to collapse locally.

For grains in a spherical cloud, sedimentation timescales are inversely proportional to the particle size expressed in terms of the Stokes number. Considering our particles in the Epstein regime, the sedimentation velocity of particles w affected by gravity and drag in a gas of velocity u is (Nayakshin 2010).

$$\begin{eqnarray}&&{w}_{\mathrm{sed}}=w-u=\displaystyle \frac{4\pi }{3}G{\rho }_{g}R{\tau }_{s}.\end{eqnarray} \tag{ 20 }$$

Assuming the radius of collapse R is equal to the half the Toomre wavelength and the central gas density ρ_g to be roughly an order of magnitude greater than the Toomre density, we approximate the sedimentation timescale as

$$\begin{eqnarray}&&{\tau }_{\mathrm{sed}}=\displaystyle \frac{{\lambda }_{{\rm{T}}}}{2({w}_{\mathrm{sed}})}=\displaystyle \frac{15\sqrt{2\pi }}{4}\displaystyle \frac{{Q}_{0}}{\mathrm{St}}{{\rm{\Omega }}}^{-1}.\end{eqnarray} \tag{ 21 }$$

This is at odds with the evolution of the three particle species in Figure 2, which show that after the gas has fragmented (at around t = 9Ω⁻¹) the particles concentrate to saturated levels considerably faster for the smaller St = 0.1 particles, slower for St = 1, and not at all for the largest size. This indicates that gravity and aerodynamics are not a sufficient description of the particle collapse and the initial angular momentum and transfer throughout the collapse affect the sedimentation with size. In this case the drift of particles inward becomes comparable to the circumstellar disks where drift is fastest for particles around 0.1 < St < 1 (Weidenschilling 1977).

The solid mass ratios of the clumps formed are greater than the initial Z₀ = 0.01 and in the case of St ≤ 1, often more comparable to the ratios observed in Jupiter or Saturn. Figure 5 shows the total mass in gas (teal) and solids (yellow) as a function of radius for the fragment profiles in Figure 4. From this one can see that we consistently form gas fragments of a few Jupiter masses with varying degrees of particle concentration. Particularly, in the case of St = 0.01 and St = 10, we can see that although there are several Earth masses worth of material in the entire fragment, it is distributed in a similar fashion to the gas (i.e., remains constant throughout the inner regions of the fragment), and thus no significant core is formed. In the cases of St = 0.1 and St = 1, keeps rising as when approaching the fragment center and the formation of a core begins to take shape.

4.2. Metallicity

In simulations that do not form a core, the fragment envelopes show radial dust-to-gas ratios consistently larger than their initial ₀ = 0.01 level, although the scatter is sometimes significant, particularly for values below that level. Clumps that do form cores have atmospheres depleted in heavier grains, but may have significant peaks in the local dust-to-gas ratio that are the locations of smaller secondary solid bodies that have not immediately accreted to the fragment center. As later accretion is expected to further alter the C/O ratio by introducing material from chemically distinct regions of the disk, the amount of solid material available at formation is important to understand the observations of gas-giant atmospheres (Espinoza et al. 2017). Whether or not a core forms could have a significant impact on the atmospheric metallicity and C/O ratio.

The satellites are perhaps due to the transfer of angular momentum of the particles with the spinning gas fragment. As evident in Figure 1, without particles, a fragment forms with significant rotation, mostly cyclonic, but also a weak region of anticyclonic rotation. Falling in toward the center of the fragment, the individual particles do not feel the pressure gradient of the gas, but will interact with the rotating gas velocity field. This is borne out not only in the satellites, but also in spiral arms of particles around the fragment center visible in Figure 6. These bodies have masses on the order of an Earth mass and remain in an orbit around the primary for the duration of the short simulation runtime.

Zoom In Zoom Out Reset image size

There is a noted correlation of giant planet occurrence with the stellar metallicity, with observed planets more common around metal-rich stars, a trend that typically supports core accretion planet formation (Fischer & Valenti 2005; Zhu 2019). We added a pair of simulations that deviated from the initial dust abundance to Z = 0.008 and Z = 0.013 for our simulations that formed the most prominent core, those with particles of size St = 0.1. As seen in Figure 7, the result is a fragment with a core, regardless if the initial metallicity is (from top to bottom) Z₀ = [1/125, 1/100, 1/75]. The difference lies in the fragment and core masses, which are significantly higher in both gas and dust components when more solid material is included. Still, the overall metallicity enhancement remains constant for the clumps in all three simulations, with an increase by a factor of around 5.

Zoom In Zoom Out Reset image size

The concentration of solids in non-fragmenting disks is well-documented (Gibbons et al. 2012; Shi et al. 2016), and our simulations show many similar features. For example, intermediate particle sizes St = 0.1,1 concentrate particularly well, while smaller and larger particles are either too well-coupled to the gas' smallest gas motions or barely influenced by the gas at all, respectively. This is apparent in the non-axisymmetric features of both 2D and 3D simulations, as both cases show most particles concentrating in these structures (Boley et al. 2010; Humphries & Nayakshin 2018).

Vortices are known as potential particle traps and can lead to the high dust-to-gas ratios that enable the formation of planetesimals (Raettig et al. 2015). Even in early, gravitationally unstable disks, vortices have shown the potential to efficiently collect material (Gibbons et al. 2015); however, vortices may not persist long enough in vertically stratified simulations (Lin & Pierens 2018), affecting the ability of particles to form planetesimals or cores. In Figure 1, our particleless fiducial simulations fragment and a strong cyclonic vortex spins at its center. However, cyclonic vortices are not good dust traps because they rotate in the same direction that the disk spins, so they increase the angular momentum of solid material that enters the vortex, preventing settling.

Theoretical and observational work is coming to the agreement that planetary masses are increasingly rare if formed by GIs. Observations have found candidates to be rare occurrences (Vigan et al. 2017) and simulations and analytic work have frequently had difficulty forming planetary masses at the right positions and masses and keeping them from migrating away (Baruteau et al. 2011; Kratter & Murray-Clay 2011). Our simulations here show that initial fragments are already quite large, and assuming that they will further accrete material, their masses could very likely increase into the brown dwarf regime. If they do not acquire much further mass the masses here are mostly within the range of many directly observed planets of a few Jupiter masses.

The analysis of Schlaufman (2018) of the relation between giant planets masses and the metallicities of their host stars finds that companions below a mass of 4 M_Jup are predominantly found around metal-rich stars, while those larger than 10 M_Jup do not preferentially form around stars of any metallicity. Thus, they find that there are two distinct populations of objects that are separated into planets and brown dwarfs/stars at roughly 10 M_Jup. While metallicity does not affect the ability of our simulations to fragment, it does appear to have an effect on the total mass, with higher initial metallicities producing fragments with more massive gaseous components. In this respect we find no evidence of a separation between two populations from the initial fragment properties, but that excludes all later dynamics and evolution of the fragment.

In addition to the formation of a core and some satellites, we also note the formation of a number of smaller solid bodies, particularly in the spiral arms, where particles tend to concentrate when not within the fragment. These objects become much more numerous in cases where fragmentation does not occur but instead gravitoturbulence prevails and fragmentation of the dense particle layer only results in many smaller fragmentation events rather than a single blob. Planetesimals such as these may be the seed for another fragment, and can contribute to the eventual formation of a terrestrial planet or be quickly accreted into the fragment. This will be the subject of subsequent investigations.

5.1. Limitations

Many gas–dust interactions have been neglected or simplified in these simulations, including grain heating and cooling, and collisional destruction and aggregation, among other effects (Nayakshin 2010). Collision speeds are expected to be high and may prevent the growth of larger solids, but may be low enough within the spiral arms to allow for aggregation to occur. Additionally, grain growth alone may contribute to the overall Toomre instability of the disk and aid in its eventual fragmentation (Sengupta et al. 2019).

Because we use a particle size defined through a Stokes number, the effective size of a particle changes with the properties of the surrounding fluid, which we illustrate in Figure 8. A particles' size in the Epstein regime scales linearly with the gas density (see Equation (15)), so as a particle moves into a high gas density region, its size will increase, which may not necessarily reflect the true nature of the system. This means that particles that are originally on the order of centimeters in size would instead have an effective size of meters once inside a gas overdensity 100 times the initial value. However, as they move through fluid regions of varying density, particles of a particular Stokes number will tend to grow to a mass that maintains their aerodynamic properties (Birnstiel et al. 2011).

Zoom In Zoom Out Reset image size

Because we do not model the further collapse of a fragment, that may significantly increase central temperatures over 1400 K and evaporate the material that constitutes the core (Helled & Schubert 2008). Additionally, the grains within the clump would begin to dominate the opacity and likely reduce the effectiveness significantly if not completely, rendering the β-cooling used here lacking. One possible solution that remains computationally inexpensive is to use a modified cooling prescription that drastically reduces cooling efficiency in overdense regions (Baehr & Klahr 2015). While an already-present core would survive the high temperatures, further growth by dust and pebbles would likely be halted. The long-term fate of a core in our simulations will someday benefit from more extensive radiative transfer models, which take into account the grain opacities and heating at the fragment centers, which could lead to grain evaporation. It should be noted that our cooling prescription includes no metallicity dependence, which is particularly important due to the dominance of dust opacities. For this reason, results might change significantly for a treatment of cooling that takes this into account, especially for a the particle-dense central region.

These simulations also assume that the there is a uniform distribution of particle sizes, all with a dust-to-gas ratio of 1:100. This is most certainly not the case, and we expect the smaller particles to be more common because there has been little opportunity for collisional growth and aggregation (Mathis et al. 1977). There is, however, evidence of some growth in the molecular cloud phase, but with a limit at micron-sized particles (Steinacker et al. 2015), but early gravitoturbulence and particle growth could lead to larger grain sizes (Sengupta et al. 2019). While larger St = 10 particles do not collect particularly well within fragments, even St = 1 particles may be far less common than smaller species, resulting in less solid accumulation within a fragment.

We do not include any radial pressure gradient, which would result in migration velocities, particularly for St = 1 particles, and affect their ability to concentrate at any one location. While these simulations may appropriately describe the local concentration of particles, the absence of a pressure gradient and thus particle drift means that particle concentrations are likely overestimated, as particles around St = 1 would migrate quickly away from the gravitationally unstable regions of interest. The spiral arms induced by gravitational collapse, however, are significant pressure maxima, which decrease the particle velocity dispersions and facilitate growth (Booth & Clarke 2016). Shearing boxes of this size are large enough to have gradients corresponding to global disk properties, but because these simulations are not intended to explore the longer-term nature of the fragments that form within, we consider the shearing-box approximation suitable to model the local collapse and initial conditions of a fragment.

As fragments migrate they may accrete more material, whether it be gaseous or solid, after formation that cannot be captured in such a local model and will affect the overall metallicity (Ilee et al. 2017; Mercer & Stamatellos 2017). A process similar in principle to pebble accretion (Ormel & Klahr 2010) could result in later metallicity enhancements (Humphries & Nayakshin 2018). Fragments are also prone to rapid migration (Baruteau & Masset 2008), which may take them to new regions of the disk that may be gas-rich or dust-rich and accretion will alter their composition, which cannot be captured in this study. What is presented in this paper represents a look at the possible initial conditions of gas-giant planets formed by GI and their later evolution is a subject of further study.

We have performed 3D hydrodynamic simulations of gravitationally unstable gaseous disks using high-resolution finite-difference simulations with the Pencil code. Using Lagrangian superparticles within the Eulerian mesh, we modeled the evolution of solids of various sizes in simulations, which fragment into dense gas clumps. The overdense arms and fragments are gravitational and hydrodynamic traps for particles, which can then concentrate into planetary cores or enrich the metallicity of gas-giant planets. We summarize our findings in the following points:

• 1.  
  Fragmentation of the gas disk concentrates particles of intermediate sizes initially through aerodynamics, but once concentrated enough they may be comparable gravitationally to the gas fragments themselves. Thus, even though the dust is not initially gravitationally unstable, there is still potential for significant solid accumulations to form in fragmenting gravitationally unstable disks, potentially resulting in the growth of solid bodies.
• 2.  
  Once collected within gas overdensities through aerodynamic drag, solids can concentrate to the point that their contribution to the total self-gravitational potential is non-negligible and thus form considerable cores within proto-gas giants. Similar to the results of Boley & Durisen (2010), we find that gas fragments up to several Jupiter masses, with central solid cores typically of a few tens of M_⊕.
• 3.  
  We find that for intermediate particle sizes St = 0.1– 1, fragments form early with already sizable cores before they have core temperatures in excess of ∼1400 K, which would sublimate solid material. This would mean that sedimentation after formation may not be necessary to explain cores in planets formed by gravitational instability.
• 4.  
  We find fragments in total to have solid fractions around 1.5%–5%, similar to the amounts suggested by core accretion and by observations of Jupiter and Saturn. This suggests that core formation and enhanced metallicity are not unique to core accretion and more information is required of a planet to determine its formation mechanism.
• 5.  
  Our fragments have atmospheres with metallicities above the initial value Z₀ = 1/100 in the case when a core does not form, but below initial metallicities when a core is formed. Whether or not a core forms may have an impact on the atmospheric metallicity of the resulting gas-giant planet.

The authors thank the anonymous referee for thorough comments that improved the quality of the paper. H.B. thanks Andreas Schreiber, Chao-Chin Yang, Dominic Batzler, and Wlad Lyra for useful discussions. All simulations made use of the Isaac cluster at the Max-Planck Center for Data and Computing in Garching.