The Fragmentation Criteria in Local Vertically Stratified Self-gravitating Disk Simulations

Accretion driven by turbulence during the early self-gravitating phase of a disk is expected to be strong, and according to the α-formalism (Shakura & Sunyaev 1973), typically approaches unity α ≈ 0.1–1, a few orders of magnitude higher than that of magnetorotational instability (Balbus & Hawley 1991). To measure the α-stress in our simulations, we will need to consider the contribution of the gravitational and Reynolds stresses:

$$\begin{eqnarray}&&\alpha =-{\left(\displaystyle \frac{d\mathrm{ln}{\rm{\Omega }}}{d\mathrm{ln}r}\right)}^{-1}\displaystyle \frac{\langle {G}_{{xy}}\rangle +\langle {R}_{{xy}}\rangle }{\rho {c}_{s}^{2}},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&\langle {G}_{{xy}}\rangle ={\int }_{-\infty }^{\infty }\displaystyle \frac{{g}_{x}{g}_{y}}{4\pi G}{dz},\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&\langle {R}_{{xy}}\rangle =\langle \rho {u}_{x}{u}_{y}\rangle .\end{eqnarray} \tag{ 11 }$$

An analytic expectation of α in a gravitoturbulent thin disk is given by Gammie (2001):

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{4}{9}\displaystyle \frac{1}{\gamma (\gamma -1){t}_{c}{\rm{\Omega }}}.\end{eqnarray} \tag{ 12 }$$

Shi & Chiang (2014) compare the standard formulation of α (Equation (9)) with one based on density-weighted values of the pressure and stress, and find that while the density-weighted α' is a factor of 1.25–2 higher, the usual formulation more closely follows the theoretical expectation (12). From Equations (9) and (12), one can compare theoretical expectations of viscous stresses in the gravitoturbulent state with the stresses generated by our simulations. It is important to note that as a local description of viscous stress, α does not necessarily model gravitationally unstable disks, as gravity inherently acts on non-local scales (Balbus & Papaloizou 1999).

Proper modeling of disk cooling is complex even without using sophisticated ray-tracing and flux-limited diffusion schemes, depending significantly on the local temperature and density of the disk from midplane to atmosphere (Hubeny 1990; Rafikov 2005):

$$\begin{eqnarray}&&{t}_{{\rm{c}}}=\displaystyle \frac{4}{9\gamma (\gamma -1)}\displaystyle \frac{{\rm{\Sigma }}{c}_{{\rm{s}}}^{2}}{\sigma ({T}^{4}-{T}_{{\rm{irr}}}^{4})}f(\tau ),\end{eqnarray} \tag{ 13 }$$

where σ is the Stefan–Boltzmann constant, f(τ) = τ + 1/τ for optical depth τ, T_irr is the reference temperature and adiabatic index γ.

Contraction to a fragment requires that thermal pressure support does not counter the self-gravity of overdense perturbations. Heat generated through viscous heating and shock dissipation should be released on a cooling timescale short enough to keep pressure support from building up, or shorter than a few dynamical timescales (see Equation (2)). Cooling with a β-prescription is implemented in Pencil using a Newtonian cooling scheme,

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

where ${c}_{{\rm{s,irr}}}^{2}$ is proportional to a non-zero background temperature, corresponding to some constant irradiation, which as demonstrated by Lin & Kratter (2016) adds additional stability against small scale density perturbations, and t_c is a β-cooling timescale (2). The sound speed relates to entropy as

$$\begin{eqnarray}&&{c}_{{\rm{s}}}^{2}=\gamma \displaystyle \frac{P}{\rho }={c}_{{\rm{s,0}}}^{2}{\rm{\exp }}(\gamma s/{c}_{{\rm{p}}}+(\gamma -1)\mathrm{ln}(\rho /{\rho }_{0})),\end{eqnarray} \tag{ 15 }$$

where c_p = 1, initial sound speed ${c}_{{\rm{s,0}}}=\pi$, density ρ₀ = 0.293 and temperatures are derived from the sound speed using $T={c}_{{\rm{s}}}^{2}/{c}_{{\rm{p}}}(\gamma -1)$. This is a simple way to model disk cooling which is not as computationally expensive as more realistic methods, particularly at the resolutions used here.

The cooling described in Equation (14) is balanced by two sources of heat generation: viscous heating and shock dissipation. Viscous heating is calculated as

$$\begin{eqnarray}&&{ \mathcal H }=2\rho \nu {{\boldsymbol{S}}}^{2},\end{eqnarray} \tag{ 16 }$$

with kinematic viscosity ν = 2.5 H⁶Ω for all simulations, as seen in Equation (5), but ultimately the heating of the simulation will be dominated by the dissipation of shocks produced by the collapse of self-gravitating filaments and clumps. The heat generated by the dissipation of shocks is determined by the strength of the constant of hyperdissipation, so in addition to our standard collection of simulations spanning resolution and cooling timescale, we include a small sample of simulations with various dissipation values to note the effect on fragmentation.

In practice, this may be overestimating the heating resulting from strong shocks generated by self-gravity. It is well understood that in outer disk regions, disk heating is dominated by stellar irradiation instead of viscous accretion (D'Alessio et al. 1999). Additionally, recent work by Rafikov (2016) suggests that shocks are not a dominant source of heat in passive disks such as those considered here.

A proper implementation of a 3D shearing box simulation requires vertical stratification of the gas density. If this is not correctly initialized, the simulation will adjust on its own, causing undesired rapid heating and cooling as the disk swells up or contracts. Thus, we assume the disk begins in hydrostatic equilibrium with self-gravity

$$\begin{eqnarray}&&\displaystyle \frac{1}{\rho }\displaystyle \frac{{dP}}{{dz}}=-{{\rm{\Omega }}}^{2}z-4\pi G{\int }_{0}^{z}\rho ({z}^{{\prime} }){{dz}}^{{\prime} }.\end{eqnarray} \tag{ 17 }$$

The initial vertical density distribution is determined by rewriting the above as a dimensionless, differential equation, following Shi & Chiang (2014):

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}\tilde{\rho }}{d{\tilde{z}}^{2}}+\displaystyle \frac{\gamma -2}{\tilde{\rho }}{\left(\displaystyle \frac{d\tilde{\rho }}{d\tilde{z}}\right)}^{2}+{Q}_{0}^{2}{\tilde{\rho }}^{2-\gamma }+\displaystyle \frac{2}{h}{\tilde{\rho }}^{3-\gamma }=0.\end{eqnarray} \tag{ 18 }$$

We use an isothermal approximation to the polytropic equation of state for which Equation (18) was derived ($\gamma =1,K={c}_{{\rm{s}}}^{2}$), and assume ${{\rm{\Sigma }}}_{0}=2H\rho$. Solving for vertical hydrostatic equilibrium including self-gravity we find

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}\tilde{\rho }}{d{\tilde{z}}^{2}}-\displaystyle \frac{1}{\tilde{\rho }}{\left(\displaystyle \frac{d\tilde{\rho }}{d\tilde{z}}\right)}^{2}+{Q}_{0}^{2}\tilde{\rho }+\displaystyle \frac{2}{h}{\tilde{\rho }}^{2}=0,\end{eqnarray} \tag{ 19 }$$

where $\tilde{\rho }=\rho /{\rho }_{0}$, $\tilde{z}=z/({c}_{{\rm{s,0}}}^{2}/\pi G{{\rm{\Sigma }}}_{0})$, ${Q}_{0}={c}_{{\rm{s,0}}}{\rm{\Omega }}/\pi G{{\rm{\Sigma }}}_{0}$, and $h=H/({c}_{{\rm{s,0}}}^{2}/\pi G{{\rm{\Sigma }}}_{0})$. The approximation of the vertically stratified initial density creates some small oscillations in as the simulation equilibriates to the true solution. This results in a small flux of mass away from the regions between $1\lt | z| \lt 2$ and toward the outer disk atmosphere, but these result in density deviations on the order of 10⁻⁴ and do not have a significant impact on the outcome of the simulation. Due to the large empty regions above and below the disk midplane, it is necessary to impose a minimum density value to avoid extreme density differences which would slow down the code immensely.

The initial temperature distribution is not stratified but instead isothermal throughout the entire vertical extent, equal to the background irradiation level. Since distant regions of the disk will be dominated by disk irradiation rather than accretion heating this is a reasonable assumption. As the disk settles to hydrostatic equilibrium, shocks are generated from self-gravitational collapse and generate significant heat, but cooling returns the vertical temperature profile to its original flat distribution.

Introducing vertical structure affects the stability and dynamics of a simulation, and should be accounted for in comparison to a thin disk approximation. Mass density should be stratified so that it is as close to hydrostatic equilibrium as possible, otherwise the simulation will adjust itself and if the difference is significant, the simulation will be adversely affected. In a simulation on the border between fragmentation and stability, too much heating initial heating may suppress fragmentation and too much cooling may result in premature fragmentation.

The stratification of mass density also results in the weakening of self-gravity, meaning a simulation established with Q₀ = 1 will not become unstable. In the thin disk approximation, the gas density is uniformly distributed over two scale heights, but in our 3D simulations, it is instead concentrated in at the midplane but distributed over twelve vertical scale heights. With more mass further away from the midplane, a disk needs to be more massive to become unstable. Theoretical expectations indicate a critical Toomre criterion closer to Q₀ ≈ 0.7 (Behrendt et al. 2015; Kratter & Lodato 2016) and our simulations reflect that.

An important question regarding the stability of self-gravity is the role of artificial viscosity on the onset of fragmentation. We use a hyperdiffusion method, implemented in Pencil according to Appendix B of Yang & Krumholz (2012):

$$\begin{eqnarray}&&{f}_{D}({\rm{\Sigma }})={\zeta }_{D}({{\rm{\nabla }}}^{6}{\rm{\Sigma }}),\end{eqnarray} \tag{ 20 }$$

which smooths perturbations on small scales without affecting the power on larger scales. It has been suggested that such viscosities are the reason for the non-convergence of the 2D simulations (Rice et al. 2014). The choice of a stronger dissipation scheme, such as hyperviscosity, is justified because of the effects of numerical oversteepening (Klee et al. 2017), which may lead to the overestimation of overdense perturbations and contribute to non-convergent behavior.

The shearing box is set up with periodic y and z boundaries, corresponding to the azimuthal and radial physical extents, respectively. This means that all quantities that exit one boundary are reproduced at the other side of the box, including the self-gravitational potential in the vertical direction.

$$\begin{eqnarray}&&f(x,y,z,t)=f(x,y+L,z,t)\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&f(x,y,z,t)=f(x,y,z+L,t)\end{eqnarray} \tag{ 22 }$$

The boundaries x = 0 and x = L require a different boundary condition on account of the shear velocity ${u}_{{\rm{y}}}={v}_{{\rm{y}}}+\tfrac{3}{2}{\rm{\Omega }}x$ (Hawley et al. 1995).

$$\begin{eqnarray}&&f(x,y,z,t)=f(x+L,y-\displaystyle \frac{3}{2}{\rm{\Omega }}{Lt},z,t)\end{eqnarray} \tag{ 23 }$$

When calculating this over a shear periodic x-boundary, the displacement due to the shear is taken into account by shifting the entire y-direction to make the x-direction periodic before proceeding with the transform in the x-direction. After the calculation of the potential in Fourier space the process is reversed to get back to real space.

Initial values of density and sound speed are set by establishing an overall Toomre Q (Equation (1)) that is borderline unstable. To this end, we use Ω = G = 1 and must consider how we calculate the sound speed and surface density. By setting the sound speed to ${c}_{{\rm{s,0}}}=\pi$, we are left with the surface density Σ, which while normally one would ensure ${{\rm{\Sigma }}}_{0}=\int {\rho }_{0}(z){dz}=1$ to achieve Q₀ = 1, we find it necessary to also consider the case where Q₀ ≈ 0.676 and initialize the surface density accordingly. To account for the gravitational field of the central star, we include a sinusoidal acceleration profile, such that mass settles toward the midplane but acceleration approaches 0 at the vertical boundaries. Such a profile is not as realistic as a linear profile, but it improves stability.

A reference temperature proportional to ${c}_{{\rm{s,irr}}}^{2}={c}_{{\rm{s,0}}}^{2}={\pi }^{2}$ is used in Equation (14) because the sound speed is initialized as ${c}_{{\rm{s,0}}}=\pi$. The shearing box needs to be large enough to contain a critical Toomre wavelength, λ = 2π/k_cr = 2πH, while also resolving this critical wavelength with at least four grid cells and avoid artificial fragmentation (Nelson 2006). We meet these conditions using ${L}_{{\rm{x}}}={L}_{{\rm{y}}}={L}_{{\rm{z}}}=(40/\pi )H\approx 12.7H;$ see Figure 1 for an example of the box scale. These boxes are smaller in radial and azimuthal extent than other comparable simulations (Hirose & Shi 2017; Riols et al. 2017), allowing for higher resolution per scale height comparatively.

Zoom In Zoom Out Reset image size

Artificial fragmentation due to poor resolution is a significant concern with self-gravitating simulations, as insufficient resolution can lead to non-negligible errors which form spurious fragments (Truelove et al. 1997). Nelson (2006) found that for gravitationally unstable disk simulations, the critical Toomre wavelength needs to be resolved by at least four resolution elements to avoid such fragmentation, described by a maximum Toomre number of

$$\begin{eqnarray}&&T=\displaystyle \frac{\delta x}{{\lambda }_{T}}=0.25.\end{eqnarray} \tag{ 24 }$$

Our simulations contain almost exactly one critical Toomre wavelength and with a resolution of N³ = 256³ cells in each direction, our minimum resolution easily meets this resolution criterion with T = 0.008, decreasing by a factor of one-half when increasing the resolution to N³ = 512³ and N³ = 1024³. The simulations of Boss (2017) only consider fragments at or near this resolution criterion, and are thus potentially susceptible to fragmentation due to numerical errors. Additionally, the Toomre length only roughly describes the feeding zone of an unstable clump which may or may not fragment further (Kratter & Murray-Clay 2011), and is not necessarily a length indicative of a fragment.

In the 2D case, fragmentation can be easily established with a Roche surface density threshold (See Section 2.4 of Baehr & Klahr 2015), but this threshold does not translate to 3D simulations and a volume density threshold. Instead, we find the pressure-modified Hill condition of Kratter & Murray-Clay (2011) to result in a suitable threshold density. This threshold considers that pressure support stabilizes a fragment against collapse and leaves it more susceptible to tidal shear. Starting with the condition that a fragment has a radius that is a fraction of it's Hill radius,

$$\begin{eqnarray}&&{r}_{\mathrm{frag}}={\left(\displaystyle \frac{1}{f}\right)}^{1/3}{\left(\displaystyle \frac{{M}_{\mathrm{frag}}}{3{M}_{* }}\right)}^{1/3}R,\end{eqnarray} \tag{ 25 }$$

where M_* is the mass of the central star, R is the location away from the center, and f is factor determined by the choice of specific heat ratio γ. Defining the fragment mass as the integration of density over radial shells,

$$\begin{eqnarray}&&{M}_{\mathrm{frag}}=4\pi {\int }_{{h}_{50}}^{{r}_{\mathrm{frag}}}{\rho }_{{\rm{h}}}{\left(\displaystyle \frac{{r}_{{\rm{h}}}}{r}\right)}^{k}{r}^{2}{dr},\end{eqnarray} \tag{ 26 }$$

with the density at the center of the fragment ρ_h at a size one-fiftieth the scale height h₅₀, a fragment outer radius of r_frag and a power-law index of the fragment density distribution k. Then, inserting the fragment mass into Equation (25) and assuming ${h}_{50}\ll {r}_{\mathrm{frag}}$, we then solve for the central density ρ_h:

$$\begin{eqnarray}&&{\rho }_{{\rm{h}}}=\displaystyle \frac{3f(3-k){M}_{* }}{4\pi {R}^{3}}{\left(\displaystyle \frac{{r}_{\mathrm{frag}}}{{h}_{50}}\right)}^{k}.\end{eqnarray} \tag{ 27 }$$

We use γ = 5/3 which makes f = 4, whereas a softer choice of γ = 7/5 would mean f = 16, resulting in a slightly higher central density. Equation (27) will serve as our threshold density for our simulations, provided that the maximum density is above this threshold for several timescales, thus separating fragments from overdensities that have not yet collapsed enough.

Finally, we wish to apply our fragment threshold, Equation (27), to our simulations, from which we need to extract the density profile of a fragment. We find that a fragment scales with k = 0.7 to a fragment radius of r_frag = 2/3 H. With this, we find that for an overdensity in our simulations to withstand tidal disruption, as is shown in Figure 7, it must satisfy

$$\begin{eqnarray}&&{\rho }_{{\rm{h}}}\geqslant 25.57\displaystyle \frac{{M}_{* }}{{R}^{3}},\end{eqnarray} \tag{ 28 }$$

which, for a protoplanetary disk around a solar-mass star, results in a densities of $1.2\times {10}^{-10}{\rm{g}}\,{\mathrm{cm}}^{-3}$ at 50 au and 1.5 × 10⁻¹¹ g cm⁻³ at 100 au. Comparing this with our simulations and with the Roche density in Figure 7, we see that this density threshold more adequately delineates between fragments and non-fragments. A disk with the temperatures as in these simulations will have an aspect ratio H/R ≈ 0.11 at 50 au and H/R ≈ 0.16 at 100 au. With these densities and aspect ratios, we approximate threshold fragment masses of 5 M_J and 15 M_J at 50 and 100 au, respectively.

Zoom In Zoom Out Reset image size

Applying the density threshold to our low-resolution simulations, which failed to form sustained fragments as can be seen in Figure 8, simulations with very short cooling times (β = 1, 0.5) surpass this threshold. If the maximum density of a fragment can surpass this critical density for more than an orbit, it can be considered a fragment even if it is subsequently disrupted. Thus, this threshold may provide a useful sink cell formation criterion for simulations with resolution too low to fragment.

Zoom In Zoom Out Reset image size

Following the results of Young & Clarke (2015), which demonstrated the need for a smoothing length to achieve convergent disk fragmentation in thin disk simulations, we expected our 3D simulations would converge, as a smoothing length is no longer required for a disk with finite thickness. There might be additional effects which affect fragmentation at different resolutions including, but not limited to, cooling, viscosity and numerical effects, but investigating the extension of local self-gravitating disks to include a vertical component is the most natural extension of that work and our results support their findings in this regard.

Recent work has focused on numerical effects and their influence on non-convergence. Klee et al. (2017) found that numerical oversteepening could be responsible for the overestimation of overdense peaks resulting in fragmentation. Through our use of sixth-order hyperdiffusion, we do not expect oversteepening to be a significant factor. Deng et al. (2017) find that dissipation of angular momentum in SPH methods removes shear support from regions with high flow velocities and as this depends on resolution, fragmentation becomes easier with increasing resolution.

While our results strongly suggest convergence, at our lowest resolution fragmentation is suppressed due to insufficient resolution for both 2D and 3D models, similar to the results of Turk et al. (2012), which found that turbulent velocity fluctuations suppressed the formation of small scale magnetic dynamos when the Jeans length is not resolved by at least 64 grid cells. At higher resolutions, the results follow the standard fragmentation scenario, with clear fragmentation at β = 2 and a lack of fragmentation evident at β = 5, 10, 20. The β cooling prescription is not a true substitute for realistic cooling, only a simple parametrization of the underlying physics. The use of a fixed background term may have an effect on convergence (Rice et al. 2012), thus one must be careful that convergence is not lost if this term is removed (Lin & Kratter 2016). However, as simulations using β cooling and background irradiation converge, it supports our understanding of the physics.

A large part of this investigation has been about the differences between 2D and 3D simulations of self-gravitating disks, but it is also necessary to discuss how they are similar. Vertical structure and the accompanying gravitational dilution mean disk instability and fragmentation require even higher gas surface densities to collect locally before GI can set in. Once instabilities develop, the results of 2D simulations with smoothed self-gravity and 3D simulations show similar structures and features on similar scales (Figure 9). Gravitoturbulent simulations in 3D show non-axisymmetric structures of the same scale as in 2D, but do not show as strong a contrast as 2D simulations due to the weaker self-gravity on scales smaller than H.

Zoom In Zoom Out Reset image size

Ultimately, both 2D and 3D simulations converge to the same critical cooling timescale, but 3D simulations take significantly more resources. While 3D simulations are convergent and can capture effects at small scales, there is growing evidence that scales smaller than ∼H contribute little to the growth of gravitationally unstable modes (Young & Clarke 2015; Riols et al. 2017). Thus, for the modeling of self-gravitating disks, 2D simulations may be sufficient so long as an appropriate smoothing length is used to account for the vertical scale of the disk.

Planet formation by gravitational instability is hindered by the large disk masses required. With the requisite Q₀ for fragmentation decreasing in these simulations, it is be even more difficult to form planets by GI by requiring either more mass or a colder disk for fragmentation. Larger relative disk sizes required for fragmentation also suggest that the resulting fragments will draw from a larger mass reservoir to overcome stabilizing effects and form larger objects. Thus, it is more likely fragments already have enough mass to be considered brown dwarfs or low-mass stars, without taking into account post formation accretion and growth.

The density threshold of Equation (28) and the subsequent fragment mass suggest objects formed in these simulations are formed with at least $4.9\,{M}_{{\rm{J}}}$ at both 50 au and 100 au around a solar-like star. Since companions will generally form larger than these thresholds and subsequently accrete more mass, these fragments will likely form planets many Jupiter masses in size and brown dwarfs as well. Recent observations show fragments are indeed quite large, as is the case of L1448 IRS3B (Tobin et al. 2016) which is consistent with theoretical expectations (Stamatellos & Whitworth 2009).

Global simulations can model the growth, migration, and potential destruction of fragments formed in a disk, but with our shearing box simulations are able to resolve the small scales of GI structure and development. For a resolution and convergence study this makes local simulations ideal for the relatively small patches of the disk where fragmentation occurs.

The Fourier solver in Pencil is limited to domains with the same number of grid cells in the vertical direction as in the x-direction, meaning that a simulation with a flattened box would have significantly higher resolution in z, resulting in different values of hyperviscosity in each direction, possibly leading to unphysical effects. For that reason, we conducted our simulations spanning roughly 12.7 scale heights in the vertical direction, much of which was not crucial to the simulation of fragmentation. Accommodating these largely empty spaces require carefully establishing a minimum density low enough to prevent undesired additional mass from settling onto the disk and affecting fragmentation. A more efficient approach would be to simulate the vertical direction in four to six scale heights while still retaining the radial and azimuthal extents. Additionally, the periodic vertical boundary conditions for the self-gravitational potential rather than typical outflow conditions and may contribute to the somewhat erratic α stresses calculated in Figure 6.

While the simple cooling with irradiation used here might result in convergence of the cooling criteria, a better description of collapse and fragmentation will ultimately come from using radiative transfer calculations, whether in the form of approximations such as flux-limited diffusion (Boley et al. 2006; Mayer et al. 2007) or more sophisticated but intensive methods like ray tracing. These methods take into account the gas opacities and optical depths to more appropriately model variations of cooling from one region of the disk to the next.

One feature of protoplanetary disks that is typically ignored in self-gravitating contexts but might have an effect on fragmentation is the disk magnetic field. While ionization is expected to be low in the cold outer regions of circumstellar disks, the effect of magnetic fields on fragmentation is uncertain. Only a pair of studies have looked into the affects of magnetic fields on the fragmentation criterion (Riols & Latter 2016; Forgan et al. 2017), and came to opposite conclusions about the importance of magnetic fields on disk fragmentation.