Convergence of the Critical Cooling Rate for Protoplanetary Disk Fragmentation Achieved: The Key Role of Numerical Dissipation of Angular Momentum

Massive rotating gas disks can become unstable as a result of their own self-gravity. Using linear perturbation theory applied to an axisymmetric, razor-thin disk, Toomre (1964) showed that the exponential growth of local density perturbations is controlled by the parameter

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

where c_s is the sound speed in the disk, κ is the epicyclic frequency, which equals the angular frequency Ω for a Keplerian disk, and Σ is the mass surface density (G is the gravitational constant). When $Q\lt {Q}_{\mathrm{crit}}\approx 1$, gravitational instability occurs locally, and sections of the disk will collapse into axisymmetric rings as gravity overcomes both pressure gradients and shear induced by rotation.

In realistic finite-thickness disks subject to generic perturbations the character of the instability is global, whereby the disk first undergoes a phase of nonaxisymmetric instability, namely, develops a prominent spiral pattern, before becoming locally unstable to collapse (Durisen et al. 2007). Inside and across spiral arms fluid elements are subject to torques that lead to angular momentum transfer (Cossins et al. 2009). Angular momentum transfer is indeed a key feature of gravitational instability and ultimately determines the evolution of the density field in other overdense regions inside spiral arms (Mayer et al. 2004; Boley & Durisen 2010). The other central aspect is that the spiral pattern causes shocks in the flow, releasing gravitational energy into heat, allowing the disk to evolve toward a marginally stable, gravito-turbulent state. However, if the cooling is too fast to be balanced by heating due to turbulence dissipation, then the disk fragments. Gammie (2001) proposed to parameterize the local cooling timescale t_cool as

$$\begin{eqnarray}&&\beta ={t}_{\mathrm{cool}}{\rm{\Omega }},\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=u{\left(\displaystyle \frac{{du}}{{dt}}\right)}^{-1},\end{eqnarray} \tag{ 3 }$$

with u being the specific internal energy. He showed that when $\beta \lt {\beta }_{\mathrm{crit}}\approx 3$ fragmentation occurs in local shearing sheet simulations with a ratio of specific heats $\gamma =2$. Physically this means that cooling has to occur on a timescale comparable to the local orbital time. Using 3D smoothed particle hydrodynamic (SPH) simulations, Rice et al. (2005) showed ${\beta }_{\mathrm{crit}}\approx 6\mbox{--}7$ for disks with $\gamma =5/3$. A similar value of ${\beta }_{\mathrm{crit}}$ was also found by Mayer et al. (2005) using a different 3D SPH code while comparing the evolution of isolated and binary protoplanetary disks. However, Meru & Bate (2011b) carried out similar SPH simulations to those in Rice et al. (2005) and found nonconvergence of the critical cooling rate when the resolution of the simulations was increased. They found that disks can fragment at a higher ${\beta }_{\mathrm{crit}}$ value in higher-resolution simulations. Meru & Bate (2012) suggested that the nonconvergent behavior was caused by extra heating due to artificial viscosity occurring at low resolution. At higher resolution, less heating is generated through artificial viscosity and a lower critical cooling rate (i.e., a larger β) is capable of balancing the combined heating by gravitational instability and artificial viscosity. Earlier on, Paardekooper et al. (2011) had run 2D simulations with the grid-based code FARGO, showing that the nonconvergence problem is not specific to SPH. By running a large set of simulations with both 3D SPH and the 2D grid-based code FARGO, Meru & Bate (2012) concluded that ${\beta }_{\mathrm{crit}}$ should converge to an asymptotic value $\approx 20\mbox{--}30$. This has important implications since it suggests that fragmentation can occur even with relatively long cooling times, corresponding to almost an order of magnitude longer than the local dynamical time. If clumps produced by fragmentation can later evolve into gas giant planets (see, e.g., Helled et al. 2014), it follows that giant planet formation via disk instability is more common than previously thought.

Paardekooper et al. (2011) suggested that the nonconvergence is at least partly due to the emergence of special locations in the disk, at the boundary between the turbulent and the laminar region, which appear when the disk is still adjusting toward a well-developed gravito-turbulent state. In order to avoid this edge effect, Paardekooper (2012) carried out high-resolution 2D local shearing sheet simulations, but found a similar increase of ${\beta }_{\mathrm{crit}}$ as resolution increased. Fragmentation also exhibited a stochastic nature in his simulations. However, recently Klee et al. (2017) performed a numerical study arguing that stochastic fragmentation is probably caused by oversteepening when performing slope limiting in grid-based codes. Moreover, Young & Clarke (2015) showed that the fragmentation boundary is strongly related to the gravitational softening in 2D simulations. When the softening is comparable to the resolution scale, ${\beta }_{\mathrm{crit}}$ increases with resolution, while when the softening is comparable to the disk scale height, ${\beta }_{\mathrm{crit}}$ varies little. Young & Clarke (2015) also showed that clumps form from overdensities with length scales $\sim H$. This is consistent with the earlier result of Gammie (2001), who found no power at scales much smaller than H, which were resolved in the simulations. These results combined suggest that whatever is the nature of nonconvergence, it should involve physics and/or numerics acting on a scale of order H.

In this paper, we attempt to explore again the issue of nonconvergence by comparing different Lagrangian hydro methods that have by design different numerical dissipation. We are driven by the notion that nonconvergence might be caused by numerical errors, perhaps associated with the boundary regions identified by Paardekooper et al. (2011). Our aim is to identify the exact source of nonconvergence, as well as a way to overcome the problem, in order to place self-gravitating disk simulations on a firm ground.

We run simulations using exactly the same initial condition of Meru & Bate (2012) with three different Lagrangian hydro methods. Specifically, we use the vanilla SPH implementation in the GIZMO code (Hopkins 2015), which is equivalent to the original GADGET3 code (see GADGET2 in Springel 2005) with standard Monaghan artificial viscosity and Balsara switch (Balsara 1995), a modified version of the same code with the artificial viscosity scheme of Cullen & Dehnen (2010) (inviscid SPH, hereafter ISPH), and the new Lagrangian meshless finite mass (MFM) method in GIZMO (Hopkins 2015), which does not need any artificial viscosity. Finally, a smaller set of additional simulations (not listed in Table 1) exploring recent variants of the SPH method are presented at the end of the paper.

Table 1.  Simulations in Which Fragments Form and Survive for More Than One Orbital Period Are Marked with an F

Simulation	$\beta =$
	3	3.5	4	5	6	8	Boundary
TSPH-LR				F	F	NF	6–8
ISPH-LR				F	NF		5–6
MFM-LR	F	F	NF	NF			3.5–4
TSPH-HR						F	>8
ISPH-HR					F	NF	6–8
MFM-HR	F	F	F	F	NF	NF	5–6
TSPH-LR-TIC		F	NF				3.5–4
TSPH-HR-TIC			F	NF			4–5
MFM-LR-TIC	F	NF					3–3.5
MFM-HR-TIC	F	NF					3–3.5
MFM-UHR-TIC	F	NF					3–3.5

Note. Otherwise, they are marked with an NF. The boundary between fragmentation and nonfragmentation is denoted by the green background color.

Download table as:  ASCIITypeset image

In Section 2, we describe the setup of the simulations. In Section 3, we present the main results. We discuss the interpretation and implications of our results in Section 4, and finally we draw our conclusions in Section 5. Results of ancillary numerical tests are also presented in Appendices A and B.

In order to compare with previous work, we use the same disk model as in Meru & Bate (2012). A $0.1\,{M}_{\odot }$ disk spans a radial range, $0.25\,\mathrm{au}\lt R\lt 25\,\mathrm{au}$, surrounding a solar-mass star. The initial surface mass density and temperature profiles are ${\rm{\Sigma }}\propto {R}^{-1}$ and $T\propto {R}^{-1/2}$, respectively. The temperature is normalized so that at the outer edge of the disk the Toomre Q parameter equals 2. We use an adiabatic gas equation of state with $\gamma =5/3$. We note that the problem is essentially scale free, and we use 1 M_⊙, 1 au, and 1 yr as the mass, length, and time unit, respectively, in our calculation. The simulations are run with different cooling rates, as well as three different resolutions. The low-resolution model (LR) comprises 250,000 particles, while 2 million and 16 million particles are used, respectively, in the high-resolution (HR) and ultra-high-resolution models (UHR). The central star is modeled with a sink particle with sink radius of 0.25 au. Gas particles reaching the sink radius are deleted and their mass and momentum are added to the star to ensure mass and momentum conservation. We use adaptive gravitational softening (Price & Monaghan 2007).

The three disk models are run with the three different hydro methods, namely, standard traditional SPH as in Meru & Bate (2012; hereafter TSPH; see Wadsley et al. 2017); ISPH, which uses the implementation of artificial viscosity by Cullen & Dehnen (Cullen & Dehnen 2010); and the MFM method. The first two methods differ only in the way artificial viscosity is implemented. In vanilla SPH we solve the fluid equations using the standard density-entropy formulation (Springel 2005) and standard Monaghan artificial viscosity with Balsara switch (Balsara 1995), designed to reduce viscosity in flows with high vorticity. In ISPH we use the same formulation of the hydro equations, but with the Cullen & Dehnen viscosity, which uses the time derivative of the velocity divergence as a shock indicator in order to eliminate numerical viscous dissipation away from shocks. In this formulation, particles have individual adjustable viscosity coefficients between an ${\alpha }_{\max }$ and an ${\alpha }_{\min }$ value and are further assigned a coefficient l that sets the decay time of viscosity away from shocks. Therefore, the latter formulation requires that we set three parameters instead of the conventional two coefficients ${\alpha }_{\mathrm{sph}}$ and ${\beta }_{\mathrm{sph}}$ in the standard Monaghan viscosity. While previous artificial viscosity schemes had appeared that used switches acting on individual particles (Morris & Monaghan 1997), including with similar parameterization (Rosswog et al. 2000), the ISPH method is a further improvement toward a truly inviscid method owing to its shock tracking procedure. The MFM method is instead an entirely different method to solve the hydrodynamical equations, but shares with SPH the Lagrangian approach and the use of particles as tracers of the fluid. Most importantly, it does not need any artificial viscosity in order to generate stable solutions, even in high Mach number flows; hence, it is by nature inviscid. While we defer the reader to Hopkins (2015) for a detailed description of the code, here we recall that MFM uses the SPH-like kernel of a particle distribution to construct a volume partition of the domain. The fluid equations are then solved on the unstructured mesh thus generated by the partition using a Riemann solver (HLLC in this paper). As the volume elements are constructed from the particles, the smoothing length h of the assigned kernel function still defines the fundamental resolution length as in SPH. The method conserves mass, momentum, and angular momentum by design as in an SPH code, but this is even more strictly true for angular momentum relative to SPH since there is no added artificial viscosity, as pointed out and demonstrated with a series of tests in Hopkins (2015). Being an integral formulation, in this sense a hybrid with finite-volume grid-based methods, it should also not suffer from the errors associated with poor estimates of gradients proper of SPH. Gravity is solved using the same tree code in all three methods, inherited from the progenitor GADGET3 code.

For all codes we employ two different techniques to start the simulations. In the first one we start directly from the initial disk model without any relaxation phase, while in the second one we first relax the system by running the disk for four outer rotation periods (ORPs) with relatively long cooling time, using $\beta =12$. By relaxing the initial condition in the second case, we allow the development of a gravoturbulent state, with sustained spiral structure but no fragmentation, and then switch to the desired β value. We do not ramp down β to the desired value as in Clarke et al. (2007) because the 16M particle simulation is computationally expensive. The gravito-turbulent state at $\beta =12$, with mild spirals, is able to avoid the transition from a smooth disk to a turbulent disk when β is switched to a smaller value. Fragments in TIC runs form through collapse of spirals and are not linked to any transitional structures. We will show that this transition can lead to numerical fragmentation. The latter setup, which we indicate with TIC (turbulent initial condition), resembles that adopted by Paardekooper et al. (2011). For both initial setups, once the desired cooling timescale is turned on, simulations are run either for at least six ORPs or until the disk fragments.

In order to facilitate comparisons, we adopt standard values of artificial viscosity parameters in both TSPH and ISPH runs. For TSPH runs we use the Monaghan viscosity with ${\alpha }_{\mathrm{sph}}=1$ and ${\beta }_{\mathrm{sph}}=2$, which is the default choice for flows in which high Mach numbers can arise (Hernquist & Katz 1989) and has been advocated also by Mayer et al. (2004) and Meru & Bate (2012). For ISPH we set ${\alpha }_{\max }=2,{\alpha }_{\min }=0$, and l = 0.05 (Cullen & Dehnen 2010).

We declare a disk to undergo fragmentation (labeled "F" in Table 1) when clumps form that are at least two orders of magnitude denser than their surroundings and can survive at least one ORP without been sheared apart. Transient overdensities can form in the runs that we label "NF," but no robust clumps form in the way just defined. The results of simulations are summarized in Table 1, where we indicate whether or not the disk fragments, and for which value of β.

3.1. Overview

3.2. Numerical Dissipation in MFM Runs

When comparing TSPH and ISPH with MFM, differences in the hydro implementation are more subtle. Indeed, MFM does not have explicit numerical viscosity; rather, it solves the hydro equations analogously to a finite-volume method using a Riemann solver on volume elements mapped from the particle distribution. There could be some residual numerical dissipation associated with the nonanalytical Riemann problem solution, but at the same time MFM does not advect fluid elements through a grid, a common source of numerical dissipation in finite-volume grid-based methods. In the Keplerian ring test shown in Hopkins (2015) MFM conserves angular momentum better than all the other methods it was compared to, including SPH with a modified hydro force (pressure-entropy, hereafter PSPH), which helps to remove unwanted numerical effects, such as artificial tension, at fluid interfaces. This would suggest that MFM is inherently less viscous than all variants of SPH, irrespective of the artificial viscosity scheme adopted. However, in the latter test the disk is 2D, adiabatic, and non-self-gravitating, which is very different from the configuration we are studying here. Therefore, in order to reassess that MFM is indeed less viscous in a relevant albeit simple configuration, we run two different test problems. First, we consider the rotating uniform isothermal sphere collapse test in Boss & Bodenheimer (1979). The results are shown in Appendix A. We find that angular momentum is transported slightly faster outward in TSPH runs as expected from stronger numerical viscous transport, albeit differences are small and overall the angular momentum profiles are comparable between the two codes at late times. Furthermore, in Appendix B we model an accreting planet in a shearing sheet configuration. Here the flow responds more abruptly to the strong gravitational pull of the planet, and differences among the three methods emerge more clearly, establishing that MFM has better angular momentum conservation in flows with large velocity gradients. Although the Cullen & Dehnen switch suppresses the artificial viscosity ${\alpha }_{\mathrm{sph}}$ coefficient, the mass accreted onto the planet is almost identical to TSPH. Instead, MFM leads to consistently smaller density increase, and the difference remains over time.

In addition to angular momentum dissipation, the evolution of the internal energy is another important proxy for how strong the numerical dissipation by artificial viscosity is, especially in spiral shocks. Figure 2 shows the internal energy measured at the disk midplane after 1.07 ORPs in a nonfragmenting high-resolution run (for which we used a large enough $\beta =8$ to avoid fragmentation). A nonfragmenting run is a conservative choice for our purpose since spiral shocks are weaker. The figure highlights how the MFM run exhibits the cooler profile everywhere because, at an equivalent cooling rate, it suffers less from artificial viscous heating, even in the region near the inner boundary where SPH is notoriously problematic in handling the rapidly increasing shear (Mayer et al. 2004; Kaufmann et al. 2007). The TSPH and ISPH simulations are almost indistinguishable. The similar internal energy peak occurs slightly above 20 au in all runs because that is where a strong ring-like overdensity develops (similar to those in the top panel of Figure 3), suggesting that where the flow is intrinsically more compressive in all codes the dominant heating comes from pressure-volume work rather than numerical viscosity.

Having now assessed that MFM is the least viscous among the methods considered, we can turn again to the interpretation of the results in Figure 1. If, as suggested in Meru & Bate (2012), reduced numerical viscous heating as the resolution is increased is the source of nonconvergence, one would expect that MFM runs should be more prone rather than less prone to fragmentation relative to TSPH and ISPH runs. This is because gas can cool more effectively at a given β, as we have just seen by comparing the internal energy evolution. Instead, the opposite trend is seen—MFM is less prone to fragmentation and its behavior is nearly convergent with resolution (exactly convergent in TIC runs). Furthermore, ISPH and TSPH appear to be equivalent when we inspect the internal energy evolution in Figure 2, yet ISPH is less prone to fragmentation. All this suggests that, in order to understand the nature of our results, among the effects associated with numerical viscous dissipation, enhanced angular momentum transport, as opposed to enhanced heating, is key. We discuss this crucial point, as well as dependence on resolution, in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. Numerical Dissipation and the Resolution-dependent Flow Properties

In standard SPH methods artificial viscosity is needed for the stability of the flow near discontinuities. Without that, particle interpenetration occurs and shocks cannot be properly resolved. The conventional form of artificial viscosity employs two terms (Monaghan & Gingold 1983; Springel 2005): one term linear in the divergence of the flow, controlled by the ${\alpha }_{\mathrm{sph}}$ coefficient, and one term with quadratic dependence on the flow divergence, controlled by ${\beta }_{\mathrm{sph}}$. The quadratic term is dominant in shocks. When ${\beta }_{\mathrm{sph}}=0$, the linear term can be recast in a shear viscosity $\eta =(1/2){\alpha }_{\mathrm{sph}}\kappa {{hc}}_{s}\rho$ and a bulk viscosity $\zeta =(5/3)\eta$, where h is the SPH smoothing length, which controls the resolution, and κ is a constant that should be chosen according to the kernel (Cullen & Dehnen 2010). The shear component of the viscosity will always be present, namely, also in nonshocking regions, and can generate artificial angular momentum transport in shearing flows such as those inherent to astrophysical disks, even in the presence of damping factors dependent on the local vorticity such as the Balsara switch (Kaufmann et al. 2007). Cullen & Dehnen (2010), building on previous work (Morris & Monaghan 1997; Rosswog et al. 2000), proposed a new artificial viscosity switch with a high-order shock indicator that results in a much more effective damping of shear viscosity away from shocks. It is thus expected that, at fixed resolution, ISPH runs, which employ the latter viscosity scheme, should be less affected by numerical dissipation of angular momentum. At the same time, we caution that one expects mild shocks to arise as the spiral pattern grows in amplitude and the disk becomes marginally unstable, which could then result in significant viscous dissipation even in ISPH.

Having shown in more than one way, in the previous section, that MFM is a less viscous method, and knowing that ISPH is by construction less viscous than TSPH, we have inferred that the different susceptibility to fragmentation of the three types of runs must be somehow caused by the different degree of angular momentum dissipation. The next step is to address how angular momentum dissipation affects fragmentation. Moreover, there is a related puzzling fact emerging from Table 1. This is that one would expect the results of hydro methods whose main difference is the viscous dissipation of angular momentum to eventually converge as the resolution is increased because, owing to the dependence of shear viscosity on the smoothing length h just highlighted, viscous dissipation should decrease as resolution increases. Instead, the opposite is seen: ${\beta }_{\mathrm{crit}}$ keeps increasing for TSPH as the resolution is increased, hence departing more and more from the ${\beta }_{\mathrm{crit}}$ value determined with MFM. However, we caution already at this point that this conventional way of reasoning is correct only if no other property of the flow changes with increasing resolution. If, for some reason, velocity gradients become intrinsically stronger as resolution is increased, one could imagine that convergence might be hard to achieve irrespective of resolution.

Here we argue that the nature of the flow in gravitationally unstable disks does change as resolution is increased. It develops regions of stronger shear and vorticity as sharper flow features become resolved, and such regions then become sites of higher numerical viscous dissipation. The resolution-dependent nature of the flow also explains why the way the initial setup is prepared matters for the fragmentation outcome at varying resolution, as we illustrate below. The changing nature of the flow with resolution is well illustrated in Figure 3. The simulations using unrelaxed initial conditions offer the best tool to understand what happens. In these the spiral pattern originates in the inner disk and propagates outward (Meru & Bate 2011a). We refer to the stage before spirals reach the outer edge of the disk as the transition phase. The duration of such a transition phase depends on the cooling rate, but in general lasts about two ORPs. By visual inspection it is clear that at higher resolution nonaxisymmetric modes are better resolved. This is the result of an improved gravitational force resolution, a point already emphasized in Mayer et al. (2004) and Mayer & Gawryszczak (2008) for both SPH and adaptive mesh refinement (AMR) codes.

Zoom In Zoom Out Reset image size

We refer to the boundary zones between the turbulent, highly nonaxisymmetric flow component generated by gravitational instability and the background smooth, axisymmetric disk as the interface region. As seen in the top panels of Figure 3, these are regions roughly a few code units in size, and they occur during the transition phase (see the rectangular region in Figure 3). Afterward the entire disk becomes nonaxisymmetric and gravito-turbulent. These interface regions are relevant because they are the location of the largest local velocity gradients, associated with both strong shear and vorticity.

We will use the analysis of vorticity as a proxy for high-velocity derivatives in the flow, as vorticity also carries information on the local angular momentum, which we posit plays a central role. The azimuthally averaged vorticity profile in the midplane is plotted in Figure 4 for MFM runs at different resolution. In all of them it is evident that vorticity can be an order of magnitude higher than the Keplerian flow vorticity. These vorticity peaks occur near interface regions and are more pronounced and more frequent as resolution increases. The latter is an effect of increased force resolution as higher overdensities appear that induce higher local vorticity. Indeed, if numerical dissipation of angular momentum does not affect the properties of the flow, one would expect density and vorticity to remain correlated. Let us consider an approximately spherical cusp of the flow located at the interface, namely, an overdense region of radius R_c and enclosed M_c, so that its mean density is $\rho \sim {M}_{c}/{{R}_{c}}^{3}$. When the cusp density is large enough, a fluid element at its boundary will move with a (circular) velocity V_c such that ${V}_{c}/{R}_{c}\propto \sqrt{G\rho }$. In other words, the cusp dominates the local gravitational field, even if collapse has not ensued, so that locally the vorticity of the flow should be correlated with local density, since, dimensionally, the vorticity ${\boldsymbol{\nabla }}\times {\boldsymbol{v}}\sim \tfrac{{V}_{c}}{{R}_{c}}\propto \sqrt{G\rho }$ (${\boldsymbol{v}}$ is the velocity vector of the gas flow). We can now turn back to our simulations and check whether or not vorticity and density in the cuspy regions grow at the relative rates implied by the proportionality just highlighted. First of all, we notice that the peak density at interface regions in MFM runs in Figure 5 is almost an order of magnitude larger in the UHR simulation relative to the LR simulation. As expected, vorticity is correspondingly higher. Note that TSPH runs share the same trend. ISPH runs yield a similar vorticity map to TSPH because at interface regions the flow has high-velocity derivatives, and hence artificial viscosity is not suppressed by the Cullen & Dehnen switch.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We cannot compare the vorticity in the MFM and SPH simulations in Figure 5 directly because the disk is much cooler in MFM (see Figure 2), so that the spiral pattern evolves faster than in the SPH simulations. But we can still compare the evolution of maximum density and maximum vorticity in a relative sense. We focus on the transition phase of a set of fast cooling simulations with $\beta =4$. We plot the peak volume density and peak vorticity in the $r\gt 8$ region in Figure 6. At t = 2 ORPs the spiral pattern reaches the outer edge of the disk in all simulations. In SPH simulations very compact, dense clumps have formed at this point, while in MFM only loosely bound overdensities are present. The peak density in MFM is indeed lower than in SPH simulations at comparable resolution. Also, SPH simulations show an exponential growth of the peak density, which is absent in the MFM simulation. The figure also shows that, in the early stage of the transition phase, the MFM run experiences a faster growth of density perturbations, reflecting the lower temperature (see Figure 2). Yet the key result emerges when we compare density and vorticity evolution in the different runs. It is noteworthy that in MFM runs vorticity appears to grow at the rate expected based on the relation ${\boldsymbol{\nabla }}\times {\boldsymbol{v}}\propto \sqrt{G\rho }$. Indeed, at 2 code units it has grown by about an order of magnitude while density has grown by nearly two orders of magnitude. Conversely, a comparable increase in vorticity occurs in TSPH and ISPH runs at this time, but there density increases much more, by about four orders of magnitude. We argue that this remarkable mismatch reflects strong numerical dissipation of angular momentum near cusps, which damps vorticity and promotes collapse.

Zoom In Zoom Out Reset image size

We can attempt to understand the nature of the problem by considering resolution requirements in self-gravitating disks. In order to properly follow the fragmentation, the most unstable Toomre wavelength ${\lambda }_{T}$ must be resolved (Nelson 2006). This is

$$\begin{eqnarray}&&{\lambda }_{T}=\displaystyle \frac{2{c}_{s}^{2}}{G{\rm{\Sigma }}},\end{eqnarray} \tag{ 4 }$$

where ${\rm{\Sigma }}\sim \rho H$. In a marginally unstable disk $Q\approx 1,{\lambda }_{T}\sim Q\tfrac{{c}_{s}}{{\rm{\Omega }}}\sim H$ (Cossins et al. 2009). We substitute Σ with $\rho {\lambda }_{T}$ in Equation (4). We can then write ${\lambda }_{T}\sim \tfrac{{c}_{s}}{\sqrt{G\rho }}$. We recall that the Toomre wavelength is borne out of linear perturbation theory for axisymmetric local perturbations. It is a conservative resolution marker for our purpose since it has been shown that the nonlinear stage of fragmentation in spiral arms occurs on a characteristic wavelength that is 5–6 times smaller than the Toomre wavelength (Boley 2009; Tamburello et al. 2015). In the interface regions the spiral pattern dominates the flow; hence, the same considerations would apply. Recalling again that, near cusps, we can write ${\boldsymbol{\nabla }}\times {\boldsymbol{v}}\,\sim \tfrac{{V}_{c}}{{R}_{c}}\,\propto \sqrt{G\rho }$, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{\lambda }_{T}}{h}\sim \displaystyle \frac{{c}_{s}}{h| {\boldsymbol{\nabla }}\times {\boldsymbol{v}}| }.\end{eqnarray} \tag{ 5 }$$

We define

$$\begin{eqnarray}&&q\equiv \displaystyle \frac{{c}_{s}}{h| {\boldsymbol{\nabla }}\times {\boldsymbol{v}}| },\end{eqnarray} \tag{ 6 }$$

where h is the kernel smoothing length. Since the Toomre wavelength should be resolved by at least one kernel (Nelson 2006), the resolution will not be adequate if $q\lt 1$, bearing in mind that the Toomre wavelength is a generous estimate for the characteristic wavelength of fragmentation. Likewise, the equation above highlights how, as vorticity increases, the resolution has to increase proportionally, namely, the smoothing length h has to decrease, in order to keep q large and thus follow a self-gravitating fluid appropriately. The top panel of Figure 7 shows that, as we expected, $q\lt 1$ in the interface regions and decreases as resolution is increased, reflecting the dominant effect of increasing vorticity relative to the smoothing length decrease.

Zoom In Zoom Out Reset image size

Conversely, outside these regions q increases as resolution is increased, as the dominant effect is now the decrease of the smoothing length. Note that Figure 7 uses the MFM runs. The bottom panel of Figure 7 also shows that q is always high at all resolutions when the initial conditions are relaxed (TIC runs), which reflects the fact that, even for MFM, only in this case there is exact convergence for the critical cooling time (see Table 1).

Indeed, in TIC runs the velocity field of the flow varies more smoothly across the disk, with no transition phase and thus no cusps (Figure 3). This confirms the claim of Paardekooper et al. (2011) on the importance of the initial setup for convergence studies. This is because artificial angular momentum transport during the transition from a smooth disk to a fully developed spiral pattern is avoided. Yet Table 1 also shows that only MFM shows exact convergence among TIC runs (${\beta }_{\mathrm{crit}}=3$). The SPH runs still have excessive viscous dissipation of angular momentum even in the presence of more moderate velocity gradients in shearing flows. While we did not run ISPH with TIC conditions, we have already shown that the viscosity switch alone becomes rather ineffective at high resolution, resulting in a behavior similar to TSPH (Figure 6).

Fragmentation in self-gravitating disks is known to be a process highly sensitive to the numerical implementation of the hydrodynamical equations and resolution. Earlier work carried out with both grid-based and SPH codes highlighted the importance of numerical approaches, including the effect of artificial viscosity, for the simple case of locally isothermal disks (Pickett et al. 2001; Mayer et al. 2004; Durisen et al. 2007; Mayer & Gawryszczak 2008). Strong dependence on numerics has also been seen in cloud fragmentation simulations. For example, Bate (2011) carried out SPH simulations of molecular cloud collapse to form protostars and disks and found more fragmentation with increasing resolution. Hu et al. (2014) ran galaxy simulations with the SPH code GADGET (Springel 2005) and showed that different artificial viscosity switches can lead to very different fragmentation results.

In this paper we have shown that simply increasing the resolution does not yield a more correct answer for disk fragmentation; rather, the nature of viscous dissipation and how it couples with the properties of the flow is the crucial aspect. Since the properties of the flow change with resolution, studying convergence as a function of resolution alone is insufficient. Instead, comparing different hydro methods, and in particular different methods with varying numerical dissipation, reveals the real nature of the problem. Angular momentum dissipation in strong shearing regions, which we study using vorticity, is the main reason behind nonconvergence in SPH methods. ISPH, which suppresses artificial viscosity in pure shearing flows, withstands faster cooling without fragmentation than TSPH. MFM is even more resistant to fragmentation.

MFM is the only method for which we were able to prove convergence of the critical cooling timescale among the methods explored so far, but even in this case exact convergence requires relaxation of the initial conditions in order to avoid transients that lead to sharp flow velocity gradients and high local vorticity at interface regions. The important role of the initial conditions is supported also by a recent study of Backus & Quinn (2016) with locally isothermal disks using a modern formulation of the SPH hydro force, but still standard artificial viscosity, using the ChaNGa code. It is also in line with previous findings by Paardekooper et al. (2011) with the FARGO code.

Clearly we have only considered a small set of hydro methods in this paper. Even in the context of SPH many different revised formulations of the hydro force have appeared in the literature over past few years that would be worth exploring, such as SPHS (Hayfield et al. 2011), GDSPH (Keller et al. 2014; Tamburello et al. 2015), and the pressure-entropy formulation (PSPH) in Hopkins (2013). Thermal diffusion alone has been shown to be effective at removing artificial surface tension (Price 2008; Shen et al. 2010). It may play a role as it does in all problems where sharp fluid interfaces appear (Agertz et al. 2007). These modifications of other aspects of the SPH formulation could also affect ${\beta }_{\mathrm{crit}}$. The smoothed cooling approach of Rice et al. (2014) is also another example of how the fragmentation problem is sensitive to the details of the numerical implementation within the same category of methods. In order to at least partially address this, we rerun the HR simulation from relaxed initial conditions (TIC) using PSPH and thermal diffusion in the GIZMO code (Hopkins 2015) and with the Cullen & Dehnen viscosity formulation. We find that the disk fragments for β = 4–5. Therefore, using PSPH does not lead to an improvement relative to SPH, confirming the central role of residual viscous dissipation as opposed to other aspects of the SPH scheme, for example, the SPH "E0" zeroth-order errors (Read et al. 2010).

Dehnen & Aly (2012) showed that the kernel function does affect numerical dissipation in shear flow simulations. The Wendland C4 kernel shows the best conservation properties in their numerical tests. Although MFM solves the fluid equations on the volume elements constructed from the particle distribution, the volume elements themselves are constructed using a specified kernel function; hence, it is relevant to ask what effect the choice of the kernel function has on the flow solution. Therefore, we rerun the MFM-HR-Beta3-TIC simulation with the Wendland C4 kernel, using 200 neighbors. We found that the simulation still fragments at ${\beta }_{\mathrm{crit}}=3$. This is reassuring, as it strengthens the conclusion that the results have truly converged, namely, we cannot push ${\beta }_{\mathrm{crit}}$ even further by reducing numerical dissipation using a more conservative kernel function. However, a word of caution must be said. With so many neighbors the effective gravitational force resolution is half the force resolution in the standard HR simulation, being equal to that in the LR simulations, and the force resolution at later stages is also not directly comparable owing to its full adaptive nature.

Meru & Bate (2012) could not prove exact convergence for a set of FARGO simulations, for which they also considered varying the coefficient of viscosity. Indeed, despite showing a trend with resolution that was pointing toward convergence, they did not find ${\beta }_{\mathrm{crit}}$ to maintain the exact same value when they increased their resolution by a factor of 2. The asymptotic value of ${\beta }_{\mathrm{crit}}$ suggested by their results was also higher than the corresponding value found in their SPH runs. In light of our findings that show how convergence occurs to a very small vaue of ${\beta }_{\mathrm{crit}}=3$, we argue that the published FARGO simulations might suffer from even stronger numerical dissipation than SPH, resulting not only from numerical viscosity but also from advection errors occurring in cuspy regions that are not at all axisymmetric and hence cannot be captured properly by the cylindrical geometry of the grid. For the same reason, significant improvements can be expected using relaxed initial conditions, as deviations from axisymmetry are less severe owing to the absence of the transition phase. This is once again in line with the findings of Paardekooper et al. (2011). Advection errors in the interface regions could be even more severe in Cartesian grids, unless aggressive AMR is used to better capture the curvature of the flow. As usual in numerical convergence studies, one would like to show that the converged answer is also the physically correct answer. As with highly nonlinear phenomena such as gravitational instability, we do not have analytical solutions with which to compare. One strategy might be to impose the same explicit viscosity to all codes, hence the same Reynolds number. If viscosity is the dominant effect in the nonconverge problem, as we argue here, the SPH variants, MFM, and possibly grid-based codes, all aligned in the latter way, should converge. The asymptotic behaviour for decreasing viscosity should then define the correct physical answer. Such a complimentary code comparison project would be a valuable follow-up to the work presented in this paper.

Finally, at fixed smoothing length the density increase in self-gravitating flow is strongly affected by gravitational softening. Adaptive softening is expected to play a key role, as it increases the force resolution as particles converge toward a cusp or sharp interface. We run additional TSPH simulations with fixed softening to quantify the importance of the choice of softening for the problem under study.

Figure 8 indicates that fixed softening simulations also experience an exponential peak density growth. In the TSPH-LR run, when we choose a gravitational softening of 0.1 code length units, which equals half the smoothing length in the disk midplane at r = 9 ( see the curve of the TSPH-LR-Beta4-S0.1 run), the peak volume density evolution is similar to that in the the adaptive softening run. If we decrease the softening to 0.05 (TSPH-LR-Beta-S0.05), the peak volume density increases significantly. The TSPH-HR-Beta4-S0.05 has even higher peak density than TSPH-LR-Beta4-S0.05, and also even higher than the peak density in the adaptive softening run, resulting from the denser cusps forming at higher resolution. As a result, runaway fragmentation in the presence of vorticity damping can be even more problematic than in the reference adaptive softening run.

Zoom In Zoom Out Reset image size

The results found in this paper suggest that disks need to cool really fast to fragment, on a timescale shorter than the local orbital time. This reinforces the notion that it is in the outer disk regions, at $R\gt 30\,\mathrm{au}$ (Clarke & Lodato 2009; Rafikov 2009), that the conditions are favorable for fragmentation. However, our study focused on idealized isolated disk models, neglecting the effect of both external triggers, such as mass accretion onto the disk from the envelope in the early stages of disk evolution (Boley & Durisen 2010; Vorobyov & Basu 2010), and internal triggers, such as perturbations by a previously formed fragment (Meru 2015). Mass loading has already been shown to bring disks to fragmentation even when radiative cooling is slow (Boley 2009), or even formally absent as in the collapse of polytropic turbulent cloud cores (Hayfield et al. 2011). The correct treatment of radiative cooling and heating originating from accretion ultimately requires radiative transfer to be employed (see, e.g., Meru & Bate 2011a; Rogers & Wadsley 2011; Mayer et al. 2016). However, attention should be paid to the initial condition setup and accretion boundary condition. Complex flow structures and overdensities arising in turbulent collapse may lead to strong numerical dissipation and cause spurious fragmentation for the same reasons described in this paper. Therefore, a weakly dissipative method such as MFM holds promise to be most suitable for radiative disk formation studies as well.

The complexity inherent to model disk fragmentation emerging from our study suggests that it is too premature to use existing simulations to compare with exoplanet detection via imaging surveys for wide-orbit gas giants (see Vigan et al. 2017). Indeed, once the formation stage of clumps is robustly modeled, their evolution via further collapse (Galvagni et al. 2012; Szulágyi et al. 2017) and migration in the disk (Baruteau et al. 2011; Malik et al. 2015) is also not completely understood, and not yet firm from the numerical side. Fast inward migration of clumps is routinely found, which alone would invalidate any direct comparison with wide-orbit gas giants as a way to infer the role of disk instability, but the migration rate in self-gravitating disks might be affected directly or indirectly by numerical dissipation effects. A comparison of hydro methods is warranted in this case as well.

In the gravitotubulent state, the disk is in thermal equilibrium (Gammie 2001; Rice et al. 2005; Cossins et al. 2009) and the gravitational stress is linked to the cooling timescale by

$$\begin{eqnarray*}&&{\alpha }_{{GI}}=\displaystyle \frac{4}{9}\displaystyle \frac{1}{\gamma (\gamma -1)}\displaystyle \frac{1}{\beta }.\end{eqnarray*}$$

The ${\beta }_{\mathrm{crit}}$ corresponds to the maximum gravitational stress the gravito-turbulent state can support. With $\gamma =7/5$, Rice et al. (2005) found a larger ${\beta }_{\mathrm{crit}}$ than that for $\gamma =5/3$. We also expect a larger ${\beta }_{\mathrm{crit}}$ for $\gamma =7/5$ disks with MFM using relaxed initial conditions.

Finally, there is an interesting upshot of our study. This is that actual physical sources of viscosity in disks (molecular viscosity is negligible compared to ${\alpha }_{{GI}}$), for example, the magnetorotational instability (MRI; Balbus & Hawley 1991), might promote fragmentation if they can dissipate angular momentum efficiently in overdense regions with strong shear. However, MRI turbulence can heat up the disk and alleviate gravitational instability. How MRI affects fragmentation is hard to predict. This will be explored in our future work.

We performed 3D hydrodynamic simulations using SPH with different artificial viscosity implementations, as well as with a new Lagrangian method, MFM, which employs a Riemann solver on a volume partition generated by particles, thus not requiring any explicit numerical viscosity.

Our TSPH simulations agree well with previous studies by Meru & Bate (2011b, 2012). They confirm the nonconvergence of the critical cooling rate for disk fragmentation in a self-gravitating disk. MFM instead attains convergence, and to a significantly lower value of the critical cooling time, ${\beta }_{\mathrm{crit}}=3$. ISPH simulations, which adopt the more conservative artificial viscosity implementation by Cullen & Dehnen (2010), resulting formally in no shear viscosity, exhibit an intermediate behavior. Finally, exact convergence, even in MFM, occurs only with relaxed initial conditions. We find a coherent explanation for our results, showing that the driving effect is numerical dissipation of angular momentum in strong shearing flows. This is exacerbated as resolution is increased, despite shear viscosity in SPH formally decreasing with resolution, because the flow develops stronger cusps with large velocity derivatives where vorticity is artificially damped. In those regions the characteristic wavelength associated with fragmentation, conservatively measured with the Toomre wavelength, remains poorly resolved even when the number of particles is increased by nearly two orders of magnitude, hence the high sensitivity on the numerical dissipation. Relaxed initial conditions exhibit a more convergent behavior because they avoid a transient phase in which overdense cusps and associated high-velocity derivatives occur, thus suffering less from the impact of artificial angular momentum dissipation by numerical viscosity.

We stress that the small value ${\beta }_{\mathrm{crit}}=3$ in the converged MFM simulations is significantly smaller than any value previously found in the literature when trying to assess convergence with increasing resolution. This indicates that the trends suggestive of asymptotic convergence reported in the literature for both SPH codes and FARGO (Meru & Bate 2012; Rice et al. 2014) were just reflecting saturation of numerical angular momentum dissipation rather than a decreasing numerical dissipation with increasing resolution.

We have further assessed the validity of our statements by running additional simulations that explore other aspects of SPH implementations, such as adaptive softening, kernel function, and form of the hydro force. These tests confirm our claim that only MFM attains exact convergence.

In summary, our results indicate that the use of hydro methods with minimal numerical viscosity is the first necessary step toward predictive simulations of disk instability because the issues we highlighted would still play a role in complex simulations with radiative transfer and envelope accretion.

We thank Jim Stone, Frederic Masset, Tom Quinn, and Charles Gammie for useful discussions and comments. We are grateful to Philip Hopkins for helping run the code. We than the anonymous referee for helpful comments. Pynbody (Pontzen et al. 2013) is used for the visualization. We acknowledge support from the Swiss National Science Foundation via the NCCR PlanetS. F.M. acknowledges support from The Leverhulme Trust and the Isaac Newton Trust.

Software: GIZMO code (Hopkins 2015), Pynbody (Pontzen et al. 2013).

To further assess the low numerical viscosity in MFM, we run the standard isothermal test case for the collapse of a rotating molecular cloud. This test, first described by Boss & Bodenheimer (1979), is a typical test case for numerical codes studying fragmentation (Bate et al. 1995). An initially isothermal, spherical, self-gravitating hydrogen molecular cloud starts to collapse from rest. Physical viscosity and the magnetic field are not included, so the angular momentum of the system should be conserved. The gas is assumed to be ideal gas, and the initial temperature is 10 K. The mass of the cloud is $1\,{{M}}_{\odot }$, and the radius of the cloud is $3.2\times {10}^{16}\,\mathrm{cm}$. The angular velocity is ${\rm{\Omega }}=1.6\times {10}^{-12}\,\mathrm{rad}\,{{\rm{s}}}^{-1}$. The ratio of thermal energy to the magnitude of gravitational energy and that of rotational energy to the magnitude of gravitational energy are $\alpha =0.25$ and $\beta =0.20$, respectively. The mean density is $\rho =1.44\times {10}^{-17}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. No density perturbation is imposed initially. The free-fall time for the initial cloud is ${t}_{\mathrm{ff}}=5.52\times {10}^{11}\,{\rm{s}}$.

We run this setup with 10,000, 100,000, and 1 million particles, using both MFM and TSPH. All the simulations share the same properties discussed below. In Figure 9, we plot the specific angular momentum profile of the sphere at different times.

Zoom In Zoom Out Reset image size

Initially the sphere is uniformly rotating, and the special angular momentum $j\propto {r}^{2}$, where r is the radius. The initial cloud radius is about 0.6 code length units. When material starts to flow inward, angular momentum must be transported outward in order to conserve the total angular momentum of the system. Particles at the outer part of the cloud drift further. We show the angular momentum profile out to only one code length unit because the outer part is noisy owing to the small number of particles.

The angular momentum profile evolves less at higher resolution in Figure 10, which is a sign of less numerical dissipation. The is because of the nice continuous flow property and because numerical viscosity scales with ${\text{}}h$. From Figure 11, we see very similar evolution of the angular momentum profile by MFM and TSPH. However, at fixed time, the TSPH angular momentum profile (dashed lines) has evolved a bit faster, i.e., more angular momentum has been transported outward, which is reflected by the fact that the MFM curve is above the TSPH curve inside 0.6 code length units and below outside.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We caution that the changing nature of the flow with increasing resolution, which we have shown to be a central aspect of unstable disk evolution, does not play a role here. This also means that there are no sites of strong shear and correspondingly strong numerical dissipation as in the disk simulations. Therefore, the fact that MFM appears to perform only marginally better than TSPH is not surprising, as this is an easy flow configuration to handle.

We are interested in assessing how numerical viscosity affects the angular moment transport, especially in cusps at interface regions where high-velocity derivatives occur. To avoid the complexity of resolving a full disk but carry out a test that can address a configuration with velocity derivatives, we use a 2D shearing sheet (Hawley et al. 1995) with an active accreting planet in the center to mimic the collapse under self-gravity. The simulations are similar to those in Ormel et al. (2015a), while the gravity of the planet is added following Ormel et al. (2015b). We choose the local sound speed and the disk scale height as the natural units of speed and length, respectively. We set the gravitational constant to 1. The planet's dimensionless mass $m\equiv {GMp}/({c}_{s}^{3}/{\rm{\Omega }})$ equals its Bondi radius in dimensionless units ${R}_{\mathrm{Bondi}}/H$. In these units the dimensionless Hill sphere ${r}_{\mathrm{Hill}}={R}_{\mathrm{Hill}}/H\,={(m/3)}^{1/3}$.

The gravity of the planet is added, $F(t)\,={F}_{0}\left\{1-\exp \left[-\tfrac{1}{2}{\left(\tfrac{t}{{t}_{\mathrm{inj}}}\right)}^{2}\right]\right\}$. The force increases gradually to mimic the collapse due to self-gravity. Here t_inj is the injection time. The gravity from the planet is smoothed to avoid numerical instability and increase the efficiency of calculation, ${F}_{0}=-\tfrac{m}{{r}^{2}}\exp \left[-A{\left(\tfrac{{r}_{\mathrm{in}}}{r}\right)}^{p}-B{\left(\tfrac{r}{{r}_{\mathrm{out}}}\right)}^{q}\right]$. We set A = 10, p = 8 and B = 1, q = 4. The force transits from a pure Newtonian force ${F}_{0}=-m/{r}^{2}$ for ${r}_{\mathrm{in}}\lt r\lt {r}_{\mathrm{out}}$ to zero continuously and quickly. We set ${r}_{\mathrm{in}}=0.1,{r}_{\mathrm{out}}=m$ to prevent boundary condition corruption by the gravity of the planet.

We run simulations with $m=0.2,{t}_{\mathrm{inj}}=2$ at a resolution of 64 × 64 particles (this means 64 particles per scale height, which is even higher resolution than the 16 million particles of UHR 3D disk runs used in this paper). The artificial viscosity ${\alpha }_{\mathrm{sph}}$ coefficient is shown in Figure 12 to compare directly TSPH and ISPH. ISPH, as expected, has a lower artificial viscosity coefficient. We use the time evolution of peak density to measure how much mass is accreted over time inside the planet radius r_in as shown in Figure 13. The faster the growth of peak density, the faster must be the outward transport of angular momentum resulting from any form of numerical dissipation, including, but not only, explicit artificial viscosity. Figure 13 highlights the lower dissipation of MFM. Indeed, MFM always has the lowest peak density, while ISPH and TSPH are almost identical despite the lower value of the ${\alpha }_{\mathrm{sph}}$ coefficient in ISPH. We argue that the faster accretion in SPH methods is due to higher numerical viscous dissipation. The Cullen & Dehnen switch is most effective away from shocks, and it is not able to effectively suppress numerical viscosity in a flow with a large velocity gradient as in this test.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size