The Concentration and Growth of Solids in Fragmenting Circumstellar Disks

TL;DR

This study uses high-resolution simulations to explore how solid cores form within gas fragments in circumstellar disks, highlighting the influence of metallicity and particle size on core mass and composition.

Contribution

It demonstrates that solid cores can develop within gas fragments through aerodynamic concentration, with core mass and metallicity effects detailed.

Findings

Solid cores up to tens of Earth masses can form within fragments.

Higher metallicity leads to larger fragments and more massive cores.

Enrichment of atmospheres depends on core formation presence.

Abstract

Due to the gas rich environments of early circumstellar disks, the gravitational collapse of cool, dense regions of the disk form fragments largely composed of gas. During formation, disk fragments may attain increased metallicities as they interact with the surrounding disk material, whether through particle migration to pressure maxima or through mutual gravitational interaction. In this paper, we investigate the ability of fragments to collect and retain a significant solid component through gas-particle interactions in high-resolution 3D self-gravitating shearing box simulations. The formation of axissymmetric perturbations associated with gravitational instabilities allows particles of intermediate sizes to concentrate through aerodynamic drag forces. By the onset of fragmentation, the mass of local particle concentrations within the fragment are comparable to that of the gas component and the sebsequent gravitational collapse results in the formation of a solid core. We find that these cores can be up to several tens of Earth masses, depending on grain size, before the fragment center reaches temperatures which would sublimate solids. The solid fraction and total mass of the fragment also depend on the metallicity of the young parent protoplanetary disk, with higher initial metallicities resulting in larger fragments and larger solid cores. Additionally, the extended atmospheres of these soon-to-be gas giants or brown dwarfs are occasionally enriched above the initial metallicity, provided no solid core forms in the center and are otherwise lacking in heavier elements when a core does form.

Tables

Table 1. Table 1 Performed simulations and their various resolutions, particle sizes St, initial metallicities Z 0 subscript 𝑍 0 Z_{0} , Fragment metallicity Z / Z 0 𝑍 subscript 𝑍 0 Z/Z_{0} , fragment and core masses

Simulation	Grid Cells	Particle number	Stokes #	$Z_0$	$Z/Z_0$	Fragment mass $M_{Jup}$	Core mass $M_{\oplus }$
G512t2	${512}^2\times 64$	$-$	$-$	$-$	$-$	14.8	$-$
P512t2pss	${512}^2\times 64$	${10}^6$	$0.01$	${10}^{-2}$	3.4	5.24	57.9
P512t2ps	${512}^2\times 64$	${10}^6$	$0.1$	${10}^{-2}$	14.6	1.32	62.3
P512t2p	${512}^2\times 64$	${10}^6$	$1$	${10}^{-2}$	1.8	16.2	89.1
P512t2pl	${512}^2\times 64$	${10}^6$	$10$	${10}^{-2}$	1.7	2.3	$-$
P1024t2pss	${1024}^2\times 128$	${10}^6$	$0.01$	${10}^{-2}$	1.7	$-$	$-$
P1024t2ps	${1024}^2\times 128$	${10}^6$	$0.1$	${10}^{-2}$	5.1	4.8	69.0
P1024t2p	${1024}^2\times 128$	${10}^6$	$1$	${10}^{-2}$	3.8	5.7	70.7
P1024t2pl	${1024}^2\times 128$	${10}^6$	$10$	${10}^{-2}$	1.5	7.19	$-$
P1024t2psl	${1024}^2\times 128$	${10}^6$	$0.1$	$8\times {10}^{-3}$	5.0	3.83	49.3
P1024t2psm	${1024}^2\times 128$	${10}^6$	$0.1$	$1.33\times {10}^{-2}$	4.8	8.01	147.6

Equations

$$Q=\frac{c_{\text{s}}\Omega }{\pi G\Sigma },$$

$$Z\equiv M_{\mathrm{solid}}/M_{\mathrm{gas}},$$

$$ϵ\equiv {\rho }_{\mathrm{d}}/{\rho }_{\mathrm{g}}.$$

$${\omega }_z=\left(\frac{\partial u_y}{\partial x}-\frac{\partial u_x}{\partial y}\right)\hat{\mathbf{z}}.$$

$$\frac{\partial {\rho }_{\mathrm{g}}}{\partial t}$$

$$\frac{\partial \mathbf{u}}{\partial t}$$

$$-2\Omega \times \mathbf{u}-\nabla \Phi +\frac{ϵ}{{\tau }_s}(\mathbf{u}-\mathbf{w})+f_{\nu }(\mathbf{u})$$

$$\frac{\partial s}{\partial t}$$

$$\frac{1}{{\rho }_{\mathrm{g}}T}\left(2{\rho }_{\mathrm{g}}\nu {\mathbf{S}}^2-\Lambda +f_{\chi }(s)\right).$$

$$f(\xi )=\nu ({\nabla }^6\xi ),$$

$$\Lambda =\frac{\rho (c_{\text{s}}^2-c_{\text{s,irr}}^2)}{(\gamma -1)t_{\text{c}}}$$

$$\Phi (\mathbf{k},t)=-\frac{2\pi G\rho (\mathbf{k},t)}{{\mathbf{k}}^2},$$

$$P=(\gamma -1)\rho \epsilon .$$

$$\frac{d{\mathbf{w}}^{(i)}}{dt}$$

$$+\frac{1}{{\tau }_s}\left({\mathbf{w}}^{(i)}-\mathbf{u}({\mathbf{x}}^{(i)})\right)$$

$$\frac{d{\mathbf{x}}^{(i)}}{dt}$$

$$\mathrm{St}={\tau }_f={\tau }_{\mathrm{s}}\Omega .$$

$${\tau }_{\mathrm{s}}=\frac{a{\rho }_{\text{vbox}\text{scalebox}0.5\bullet }}{c_{\mathrm{s}}{\rho }_{\mathrm{g}}}$$

$${\Sigma }_{\mathrm{T},0}=\frac{c_{\mathrm{s},0}\Omega }{\pi G{Q}_0}=37{\mathrm{g}\mathrm{cm}}^{-2}.$$

$${\Sigma }_{\mathrm{T}}\equiv {\int }_{-H}^H{\rho }_0{e}^{-z^2/2H^2}dz$$

$$M\approx {\lambda }_{\mathrm{T}}^2{\Sigma }_{\mathrm{T},0}=17.2M_{\mathrm{Jup}},$$

$${\rho }_{\mathrm{T}}=\frac{{\Omega }^2}{\sqrt{2\pi }\pi G{Q}_0},$$

$$w_{\mathrm{sed}}=w-u=\frac{4\pi }{3}G{\rho }_{g}R{\tau }_s.$$

$${\tau }_{\mathrm{sed}}=\frac{{\lambda }_{\mathrm{T}}}{2(w_{\mathrm{sed}})}=\frac{15\sqrt{2\pi }}{4}\frac{Q_0}{\mathrm{St}}{\Omega }^{-1}.$$

Full text

1. Introduction

Gas giant planets are predominantly gaseous objects, but likely with 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 $>10M_{\oplus}$ (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 $10M_{\oplus}$ (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 gravitational instabilities; however recent simulations suggest that fragments may have atmospheres which 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 which 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 a useful constraints on expected metallicity values of exoplanet atmospheres (Espinoza et al., 2017; Thorngren et al., 2016). 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)

[TABLE]

where values above $Q=1$ represent marginal stability to fragmentation, but non-axisymmetric structures may form up to $Q\sim 1.5$ (Durisen et al., 2007). Here, $c_{\textnormal{s}}$ is the isothermal sound speed and $\Omega=(GM/R^{3})^{1/2}$ is the orbital frequency corresponding to the stabilizing effects of thermal pressure and tidal shear respectively, which balance the gravitational collapse of a significantly overdense mass surface density $\Sigma$, with constants $\pi$ 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_{\mathrm{c}}=\beta\Omega^{-1}$, with a critical value around $\beta_{\mathrm{crit}}=3$ (Gammie, 2001; Deng et al., 2017; Baehr 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 voritices (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) consider 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_{\rm{gas}}\geq 10$). We are interested in the regime where the gas disk is gravitationally unstable ($Q_{\rm gas}\leq 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 which 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 find 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).

Since 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, 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

[TABLE]

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

[TABLE]

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 $\bm{u}$: $\omega=\nabla\times\bm{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

[TABLE]

Post-formation enrichment was studied by Helled & Schubert (2008) and Helled & Bodenheimer (2010), assuming that fragments are formed initially as pure gas overdensities which subsequently accrete nearby material. While they find enrichment of the atmosphere is possible through the accretion of planetesimals, core formation is excluded due to the high temperatures at the center which evaporate all solids which do accrete and 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 ($\sim 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 which 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.

2. Theory

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 111http://pencil-code.nordita.org/ 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 $\sim 2\pi 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 which are linearized and transformed into co-rotating Cartesian coordinates, where $q=\mathrm{ln}\Omega/\mathrm{ln}R=3/2$ is the shear parameter:

[TABLE]

In equations (5) - (7), $\mathbf{u}=(v_{\textnormal{x}},v_{\textnormal{y}}+q\Omega x,v_{\textnormal{z}})^{T}$ is the gas flow plus shear velocity in the local box, $\bm{w}$ is the particle velocity which imparts a back reaction on to the gas proportional to the local dust-to-gas ratio $\epsilon$, $\rho_{\mathrm{g}}$ is the gas density, and the thermodynamic variable is the gas entropy $s$. Viscous heat is generated by $\mathcal{H}=2\rho_{\mathrm{g}}\nu\mathbf{S}^{2}$, with rate-of-strain tensor $\mathbf{S}$. Hyperdissipation is applied with the terms $f_{D}(\rho_{\mathrm{g}})$, $f_{\nu}(\bm{u})$, $f_{\chi}(s)$ which for each has the form

[TABLE]

with constant $\nu=2.5\,H^{6}\Omega$ (Yang & Krumholz, 2012). Heat is lost via simple $\beta$-cooling prescription

[TABLE]

with $t_{\textnormal{c}}$ given by $t_{\mathrm{c}}=\beta\Omega^{-1}$ and background irradiation term $c_{\textnormal{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

[TABLE]

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

Finally, we use an ideal equation of state, with internal energy $\varepsilon$, and specific heat ratio $\gamma$

[TABLE]

Particles are added such that the initial distribution maintains a physically motivated metallicity of 1:100, roughly that of the interstellar medium (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 value will be investigated later. Particles evolve as a collection of solids, also known as a superparticle, $i$, with position $\bm{x}^{(i)}$ and velocity $\bm{w}^{(i)}$ as in Youdin & Johansen (2007)

[TABLE]

where $\tau_{f}$ (also referred to as the Stokes number) is the particle stopping time normalized by the dynamical time $\Omega^{-1}$

[TABLE]

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

[TABLE]

where $a$ is the particle diameter and \rho_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}} 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 $\epsilon=10^{-2}$ the contribution of particles to the gravitational potential is initially miniscule and thus the potential is dominated by the gas distribution. All simulations include a million particles which 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.

3. Numerical Model

All simulations were run using the high-order finite-difference code Pencil (Brandenburg, 2003) 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 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.

3.1. Particles in Pencil

Particles are included in the Pencil code using Lagrangian ’swarm’ particles which represent collections of similar-sized particles all moving together (Johansen & Youdin, 2007; Youdin & Johansen, 2007; Schreiber & Klahr, 2018). These particles feel the effect of self-gravity and feel drag forces from the gas as well as feedback of the particles onto the gas but do not collide, merge or otherwise interact with the gas or each other. We include particles of various sizes, i.e. with a range of Stokes numbers $\mathrm{St}=[0.01,0.1,1,10]$ that also quantify how quickly the particle couples to the gas, to evaluate the ability of these various sizes to be captured by fragments.

The evolution of the particles is governed by equations (12) and (13), and are mapped to the grid using a triangle-shaped-cloud method, which distributes particle mass with a weight function $W_{\mathrm{I}}(\bm{x}^{(i)}-\bm{x})$ which puts most of the mass in the nearest cell and the rest in the 26 surrounding cells (Youdin & Johansen, 2007). This keeps the particles from becoming too discretely distributed and for the correct calculation of the back-reaction of the gas onto the dust, which uses the same distribution function to weight the velocity contributions towards the superparticle from each of the 27 cells (Yang & Johansen, 2016).

3.2. Initial Conditions

Our simulations cover $18H$ in radial and azimuthal directions and $L_{\rm z}=2H$ in the vertical, with a resolution of either 30 or 60 grid points per scale height in all spatial dimensions. This fulfills the resolution requirements to avoid truncation errors which can become unstable (Truelove et al., 1997; Nelson, 2006) as well as the requirements for individual particles, which require that the grid spacing $\Delta x<c_{\textnormal{s}}\tau_{s}$ (Bai & Stone, 2010). Additionally, the number of particles must be enough such that at the midplane layer there is 1 particle per cell, which is marginally met in our high grid resolution simulations for a $1024\times 1024$ midplane with $10^{6}$ total particles. Each simulation contains a single particle species of constant size defined by its Stokes number.

Density and temperature are established such that $Q_{\mathrm{0}}$ is on the edge of stability. While Baehr et al. (2017) found that instability requires $Q_{\mathrm{0}}\approx 0.676$, this was likely due to the small radial domain, which excluded the collapse of large unstable wavelengths. Thus for the larger radial-azimuthal domains considered here $Q_{\mathrm{0}}\approx 1$ remains the critical Toomre value, for these simulations meaning that for $c_{\textnormal{s}}=\pi$ and $G=\Omega=1$ the vertically integrated volume density is $\Sigma=1$.

Figure 1 shows the particle-less fiducial simulation which demonstrates the fragmentation of the gas which is expected in all simulations. Particularly interesting is the presence of a strong coherent anticyclonic vortex at the fragment center as well as a second rotating region separated by a largely stationary gap. If particles concentrate at the center of fragments the existence of a vortex could aid the process.

4. Results

Initiated with $Q_{0}=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. $\mathrm{St}=0.01$ will take longer as minute gas disturbances will deflect their motion towards the midplane, while larger particles, $\mathrm{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) $\mathrm{St}=[0.1,1,10]$. Gas densities above $\sim 100\rho_{0}$ are considered to be the condition for fragmentation. One can see a clear trend towards higher particle concentrations with smaller particle 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 greyscale. Even though the particles are intended to be uniformly distributed, in practice particles are first even 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.

As one might expect, the particles which are most coupled to the gas show the most concentration in the gas overdensities. In the case of our largest particles ($\mathrm{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 which 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 $\mathrm{St}=[0.01,0.1,1]$, but in the case of $\mathrm{St}=0.01$, the core is delayed but only because fragmentation of the gas, roughly defined when the maximum gas density is $100\rho_{0}$, 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 $\sim 0.1H$ and a power law decline thereafter. This is consistent with the isothermal temperature structure assumed for Bonnor-Ebert spheres, here shown with the dark blue line.

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 towards the center of the fragment.

In the case of $\mathrm{St}=0.1$, the core is especially pronounced from the surrounding protoplanetary envelope, 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 timestep 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 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 and is speculated 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_{\odot}$ star. The orbital frequency at 100 au is $1.98\times 10^{-10}$ $\mathrm{s}^{-1}$ and the sound speed of a 10 Kelvin monatomic gas (with $\gamma=5/3$) is $262$ $\mathrm{m\,s}^{-1}$. The Toomre surface density $\Sigma_{\mathrm{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

[TABLE]

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 multiplying this value by the total mass in code units of the gas or solids. From $\Sigma_{\mathrm{T}}$, we calculate our unit density from the Gaussian profile of the initial vertical gas mass distribution

[TABLE]

where $\rho_{\mathrm{0}}=0.29\hat{\rho}$. Solving for our unit volume density $\hat{\rho}$ gives a value of $5\times 10^{-12}$ $\mathrm{g\,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 profile 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 inwards. Since the density of the mass of particles which constitute the solid core fall well short of the density of a solid body of similar mass, around 5 $\mathrm{g\,cm}^{-3}$, these are not true cores at a solid material density, but a cloud of particles at the size of the grid resolution $\sim 0.02H$. Beyond this point, any further concentration of the particles is limited by the resolution and can be modeled by other means (c.f. Nayakshin, 2018).

To get a rough estimate of the maximum initial mass expected for these simulations, one assumes that all mass within a critical wavelength $\lambda_{\rm T}=2\pi H$

[TABLE]

which is the maximum mass our fiducial particleless simulation ever reaches. Since 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_{\rm Jup}$ range and 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 $\rho_{T}$, from the 3D Toomre parameter $Q_{\rm 3D}$ of Mamatsashvili & Rice (2010)

[TABLE]

where we have used $\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).

[TABLE]

Assuming the radius of collapse $R$ is equal to the half the Toomre wavelength and the central gas density $\rho_{g}$ to be roughly an order of magnitude greater than the Toomre density, we approximate the sedimentation timescale as

[TABLE]

This is at odds to the evolution of the three particle species in Figure 2, which show that after the gas has fragmented (at around $t=9\Omega^{-1}$) the particles concentrate to saturated levels considerably faster for the smaller $\mathrm{St}=0.1$ particles, slower for $\mathrm{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 affects the sedimentation with size. In this case the drift of particle inwards becomes comparable to the circumstellar disks where drift is fastest for particles around $0.1<\mathrm{St}<1$(Weidenschilling, 1977).

Solid mass ratios of the clumps formed are greater than the initial $Z_{0}=0.01$ and in the case of $\mathrm{St}\leq 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 $\mathrm{St}=0.01$ and $\mathrm{St}=10$, we can see that although there is several Earth masses worth of material in the entire fragment, it is distributed in a similar fashion to the gas (i.e.: $\epsilon$ remains constant throughout the inner regions of the fragment), and thus no significant core is formed. In the cases of $\mathrm{St}=0.1$ and $\mathrm{St}=1$, $\epsilon$ keeps rising as one approaches the fragment center and the formation of a core begins to take shape.

4.2. Metallicity

In simulations which do not form a core, the fragment envelopes show radial dust-to-gas ratios consistently larger than their initial $\epsilon_{0}=0.01$ level, although the scatter is sometimes significant, particularly for values below. Clumps which do form cores have atmospheres depleted in heavier grains, but may have significant peaks in the local dust-to-gas ratio which are the locations of smaller secondary solid bodies which have not immediately accreted to the fragment center. Since 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 anti-cyclonic rotation. Falling in towards 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.

There is a noted correlation of giant planet occurrence with the stellar metallicity, with observed planets more common around metal-rich stars, a trend which typically supports core accretion planet formation (Fischer & Valenti, 2005; Zhu, 2019). We added a pair of simulations which deviated from the initial dust abundance to $Z=0.008$ and $Z=0.013$ for our simulations which formed the most prominent core, those with particles of size $\mathrm{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_{0}=[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.

5. Discussion

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 particles sizes $\mathrm{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 which could lead to the high dust-to-gas ratios which enable the formation of planetesimals (Raettig et al., 2015). Even in the early, gravitationally unstable disk, 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, one can see that our particle-less fiducial simulations fragment and a strong cyclonic vortex spins at its center. However, cyclonic vortices are not good dust traps because since they rotate in the same direction as the disk spins, they will increase the angular momentum of solid material which enter into the vortex, preventing settling.

Theoretical and observational work is coming to the agreement that planetary masses are increasingly rare if formed by gravitational instabilities. 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) into the relation between giant planets masses and the metallicity of their host stars finds that companions below a mass of $4\,M_{\rm Jup}$ are predominantly found around metal-rich stars while those larger than $10\,M_{\rm Jup}$ do not preferentially form around stars of any metallicity. Thus they find that there are two distinct populations of objects which are separated into planets and brown dwarfs/stars at roughly $10\,M_{\mathrm{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 notice 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 the case 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 a another fragment, contribute to the eventual formation of a terrestrial planet or may be quickly accreted into the fragment and 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, 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 it’s eventual fragmentation (Sengupta et al., 2019).

Since we use a particle size defined though 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 would mean that particles which are originally on the order of centimeters in size would instead have an effective size of meters once inside a gas overdensity one hundred 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 which maintains their aerodynamic properties (Birnstiel et al., 2011).

Since we do not model the further collapse of a fragment which may significantly increase central temperatures over 1400 K and evaporate the material which 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 $\beta$-cooling used here lacking. One possible solution which remains computationally inexpensive is to use a modified cooling prescription which 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 which takes this into account, especially for 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 since 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 $\mathrm{St}=10$ particles do not collect particularly well within fragments, even $\mathrm{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 $\mathrm{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 $\mathrm{St}=1$ would migrate quickly away from the gravitationally unstable regions of interest. The spiral arms induced by gravitational collapse are however 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 since 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 which cannot be captured in such a local model and will affect the overall metallicity (Mercer & Stamatellos, 2017; Ilee et al., 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 which may be gas 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.

6. Conclusions

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 planetery cores or enrich the metallicity of gas giant planets. We summarize our findings in the following points:

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. 2.

Once collected within gas overdensities through aerodynamic drag, solids can concentrate to the point where their contribution to the total self-gravitational potential is non-negligible and thus form considerable core within proto-gas giants. Similar to the results of Boley & Durisen (2010), we find that for gas fragments up to several Jupiter masses, with central solid cores typically of a few tens of $M_{\oplus}$. 3. 3.

We find that intermediate particle sizes $\mathrm{St=0.1\--1}$, fragments form early with already sizable cores before they have core temperatures in excess of $\sim$1400 K which would sublimate solid material. This would mean that sedimentation after formation may not be necessary to explain core in planets formed by gravitational instability. 4. 4.

We find fragments in total to have solid fractions are 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 formation mechanism. 5. 5.

Our fragments have atmospheres with metallicities above the initial value $Z_{0}=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 which improved the quality of the paper. HB would like to thank 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.