The impact of numerical oversteepening on the fragmentation boundary in self-gravitating disks

TL;DR

This study shows that the choice of numerical limiters in simulations of self-gravitating disks significantly affects the predicted fragmentation boundary, with oversteepening leading to overestimated surface densities and shifted fragmentation criteria.

Contribution

It demonstrates how numerical oversteepening influences fragmentation boundaries in self-gravitating disk simulations, highlighting the importance of limiter choice in accurate modeling.

Findings

Inappropriate limiters cause overestimation of surface density.

Oversteepening effect amplifies over dynamical timescales.

Fragmentation boundary shifts to larger cooling timescales.

Abstract

Context. It is still an open issue whether a self-gravitating accretion disk fragments. There are many different physical and numerical explanations for fragmentation, but simulations often show a non-convergent behavior for ever better resolution. Aims. We investigate the influence of different numerical limiters in Godunov type schemes on the fragmentation boundary in self- gravitating disks. Methods. We compare the linear and non-linear outcome in two-dimensional shearingsheet simulations using the VANLEER and the SUPERBEE limiter. Results. We show that choosing inappropriate limiting functions to handle shock-capturing in Godunov type schemes can lead to an overestimation of the surface density in regions with shallow density gradients. The effect amplifies itself on timescales comparable to the dynamical timescale even at high resolutions. This is exactly the environment, where clumps are expected to form. The effect is present without, but scaled up by, self-gravity and also does not depend on cooling. Moreover it can be backtracked to a well known effect called oversteepening. If the effect is also observed in the linear case, the fragmentation limit is shifted to larger values of the critical cooling timescale.

Equations

$$Q=\frac{c_s\kappa }{\pi G\Sigma }$$

$$\beta ={\tau }_{\mathrm{c}}\Omega \lesssim 3,$$

$$\frac{\partial \Sigma }{\partial t}+\nabla \cdot \left(\Sigma \mathbf{v}\right)=$$

$$\frac{\partial \Sigma \mathbf{v}}{\partial t}+\nabla \cdot \left(\Sigma \mathbf{v}\otimes \mathbf{v}+P\mathbb{I}\right)=$$

$$-\Sigma \nabla {\Phi }_{\mathrm{sg}}$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left(\left(E+P\right)\mathbf{v}\right)=$$

$$-\Sigma \nabla {\Phi }_{\mathrm{sg}}\cdot \mathbf{v}+\Sigma Q_{\mathrm{cool}}$$

$$\Delta {\Phi }_{\mathrm{sg}}=4\pi G\Sigma \delta (z),$$

$$Q_{\mathrm{cool}}=\frac{p\Omega }{(\gamma -1)\beta }$$

$$p=(\gamma -1)\Sigma ϵ,$$

$$E=\frac{1}{2}\Sigma {\mathbf{v}}^2+\Sigma ϵ.$$

$$\frac{\partial \Sigma }{\partial t}+\nabla \cdot \left(\Sigma {\mathbf{v}}^{\prime }\right)-q\Omega x\frac{\partial \Sigma }{\partial y}=0$$

$$\frac{\partial \Sigma {\mathbf{v}}^{\prime }}{\partial t}+\nabla \cdot \left(\Sigma {\mathbf{v}}^{\prime }\otimes {\mathbf{v}}^{\prime }+p\mathbb{I}\right)-q\Omega x\frac{\partial \Sigma {\mathbf{v}}^{\prime }}{\partial y}=$$

$$\Omega \Sigma \text{pmatrix}2v_y^{\prime }$$

$$E^{\prime }=\frac{1}{2}\Sigma {{\mathbf{v}}^{\prime }}^2+\Sigma ϵ.$$

$$f\left(x,y\right)=f\left(x,y+L_y\right)$$

$$f\left(x,y\right)=f\left(x+L_x,y-q\Omega L_{x}t\right).$$

$${\Phi }_{\mathrm{sg}}(\mathbf{k}(t))=2\pi G\frac{\Sigma (\mathbf{k}(t))}{|\mathbf{k}(t)|}.$$

$$\alpha =\frac{1}{q^2\gamma (\gamma -1)\beta }.$$

$$⟨H_{xy}⟩=⟨\Sigma v_x^{\prime }{v}_y^{\prime }⟩,$$

$$⟨G_{xy}⟩=\sum _{\mathbf{k}}\frac{\pi G{k}_x{k}_y}{|\mathbf{k}{|}^3}{\left|\Sigma (\mathbf{k})\right|}^2.$$

$$\alpha =\frac{⟨H_{xy}⟩+⟨G_{xy}⟩}{q⟨\Sigma c_{\mathrm{s}}^2⟩}.$$

$$\bar{\Delta }=\varphi (s)\Delta ,$$

$${\varphi }_{\mathrm{vl}}(r)=\frac{r+|r|}{1+|r|},$$

$${\varphi }_{\mathrm{sb}}(r)=$$

$$v_x(t)$$

$$v_y(t)$$

$$E_{\mathrm{epi}}={\Sigma }^2\left(v_x^2+\frac{2}{2+q}{v}_y^2\right)$$

$$\frac{{\mathrm{d}}^2\delta \Sigma }{{\mathrm{dt}}^2}-\frac{\mathrm{d}\delta \Sigma }{\mathrm{d}t}2q$$

$$\delta \Sigma$$

$$-2\Omega$$

$$\delta \xi =\frac{1}{{\Sigma }_0}\left(\frac{\partial v_y^{\prime }}{\partial x}-\frac{\partial v_x^{\prime }}{\partial y}\right)-\frac{(2-q)\Omega \delta \Sigma }{{\Sigma }_0^2}$$

Full text

*Aims. *We investigate the influence of different numerical limiters in Godunov type schemes on the fragmentation boundary in self-gravitating disks.

*Methods. *We compare the linear and non-linear outcome in two-dimensional shearingsheet simulations using the VANLEER and the SUPERBEE limiter.

*Results. *We show that choosing inappropriate limiting functions to handle shock-capturing in Godunov type schemes can lead to an overestimation of the surface density in regions with shallow density gradients. The effect amplifies itself on timescales comparable to the dynamical timescale even at high resolutions. This is exactly the environment, where clumps are expected to form. The effect is present without, but scaled up by, self-gravity and also does not depend on cooling. Moreover it can be backtracked to a well known effect called oversteepening. If the effect is also observed in the linear case, the fragmentation limit is shifted to larger values of the critical cooling timescale.

1 Introduction

If a disk becomes massive enough, compared to its enclosed central object, self-gravity plays a major role in the further evolution of the disk. This is of relevance in protoplanetary and AGN-disks, at least in the early phase. The resulting gravitational instabilities can either drive angular momentum transport by settling into a gravito-turbulent state (Paczynski, 1978), or lead to fragmentation and thus yield a possible mechanism for planet and star formation (Boss, 1998; Levin & Beloborodov, 2003). A disk becomes gravitationally unstable if the Toomre-parameter

[TABLE]

fulfills $Q\lesssim 1$ (Toomre, 1964). Here $c_{\mathrm{s}}$ is the speed of sound, $\kappa$ the epicyclic frequency, $\Sigma$ the surface density and $G$ the gravitational constant. Throughout this paper we will assume that the angular velocity $\Omega$ is close to its Keplerian value. In that case $\kappa=\Omega$ holds.

The exact conditions which lead to fragmentation is subject of debate. Gammie (2001) showed that if the cooling timescale $\tau_{\mathrm{c}}$ falls short of the dynamical timescale $\Omega^{-1}$, i.e., the dimensionless relative cooling time scale

[TABLE]

fragmentation occurs. More recent publications, however, show that with better resolutions, fragmentation can also occur for larger values of $\beta_{\mathrm{crit}}$, i.e., slower cooling (Meru & Bate, 2011, 2012). This non-convergence of numerical simulations remains an ongoing problem. Rice et al. (2005) even stress that $\beta_{\mathrm{crit}}$ should be used at all for investigating the onset of fragmentation. Instead, the maximum stress that a disk can carry before it is fragmenting should be used. This is represented by a maximum value for the dimensionless parameter $\alpha$ (Shakura & Sunyaev, 1973). Furthermore, it is not clear whether numerical or physical effects are responsible for this fragmentation. Cossins et al. (2009) check for a wide range of parameters in order to characterize self-gravitating disks in a more encompassing manner. The current state of the debate was reviewed recently by Kratter & Lodato (2016) and Rice (2016).

On the numerical side, there are several effects which can be present in different kind of codes. Lodato & Clarke (2011) take artificial viscosity into account as a possible explanation for fragmentation. This is investigated further by Meru & Bate (2012). They tune their parameters defining the artificial viscosity in SPH and grid-based simulations (Fargo, see Masset 2000) in a way to maximize $\beta_{\mathrm{crit}}$, thus minimizing artificial viscosity, which goes along with a reduction of numerical heating. They also state that the tuned parameters are not in the ideal regime against particle penetration. Their results lead to large critical relative cooling timescales $\beta_{\mathrm{crit}}\lesssim 30$, at which fragmentation occurs, which differ from results for grid-based and other SPH simulations. Rice et al. (2012, 2014) argue that this procedure causes unforeseeable numerical errors, and state that the problem is more about the cooling function, which is not properly introduced within the SPH codes, leading to more likely fragmentation at higher resolutions. With this they reach a fragmentation boundary at around $6\lesssim\beta_{\mathrm{crit}}\lesssim 8$. Müller et al. (2012) stress that in two-dimensional simulations gravity needs to be smoothed in order to avoid singular behavior at higher resolutions. Simulations by Young & Clarke (2015) show indeed that their results change fundamentally when smoothing out the gravity to Keplerian vertical pressure scale height $H$, because of the neglected vertical structure. They state that within their framework, where two possible fragmentation modes are present, the second one, driven by quasi-static collapse, shows up only if the third dimension is taken into account. This approach needs to resolve gravitational forces on scales below the vertical scale-height $H$, which is not possible in 2D, because of the above mentioned limitations. In order to investigate non-quasi-static fragmentation a two-dimensional model is, however, still valid and has been used, for instance, by the same authors later (Young & Clarke, 2016). Another numerical effect may arise by additional edge effects, emerging from the boundaries in global simulations (Paardekooper et al., 2011) and leading to fragmentation at higher values of $\beta$. They solve this problem by adjusting $\beta$ from a very large value slowly to the actual one.

On the physical side Hopkins & Christiansen (2013) show that a turbulent disk is never stable. They use a stochastic approach and find that there is always a probability of a turbulent, self-gravitating disk to fragment. Rice (2016), however, claims that this is due to the isothermality of the disks and to the feedback cycle of gravitational instabilities is not being treated appropriately. Similar assumptions let Lin & Kratter (2016) revise the Toomre criterion itself, which actually applies only to a laminar, inviscid disk in which there are no net thermal losses. They conclude that adding cooling and viscosity leads to secondary instabilities which can drive fragmentation and that there formally exists no region for $\beta$ or the maximum viscous stress $\alpha$ where the disk is stable. That there are two modes of fragmentation, one for $\beta_{\mathrm{crit}}\lesssim 3$ at free fall timescale and one for $\beta_{\mathrm{crit}}\lesssim 12$ with a quasi-static collapse, is mentioned by Young & Clarke (2015). For the second mode, $\beta_{\mathrm{crit}}$ is derived by comparing cooling timescales with the timescales of encountering a spiral wave assuming the tightly-wound approximation. Young & Clarke (2015) conclude that two-dimensional calculations are not suitable to simulate the quasi-static part of the fragmentation, because these require gravitational smoothing, which suppresses a quasi-static collapse at the same time.

Moreover Paardekooper (2012) shows in the shearingsheet that fragmentation has a stochastic nature. He runs simulations with the same setup multiple times, which differ only in the subsonic random velocity field. Young & Clarke (2016) quantify this effect and come to the conclusion that stochastic fragmentation can, in general, not change the radii of fragmentation significantly.

Outline: In section 2 we will explain the underlying numerics. In section 3 we will show tests of the code, making extensive use of the linear theory test used by Gammie (2001) and Paardekooper (2012). The nonlinear outcome will be presented in section 4 where the value of $\beta_{\mathrm{crit}}$ and the stochastic nature is examined. In section 5 further numerical effects and their impacts on the results are discussed as well as the expansion to other kinds of codes. In section 6 we summarize and discuss our results.

2 Methods

We use the hydrodynamic simulation suite Fosite111http://www.astrophysik.uni-kiel.de/fosite/ (Illenseer & Duschl, 2009) which is based on an enhanced method of a Godunov type scheme (Kurganov & Tadmor, 2000). Generally, it can calculate on arbitrary two-dimensional curvilinear, orthogonal grids and is easily expandable, because of its modular design.

The shearingsheet-approximation is implemented within this framework, i.e., modules for fictitious forces, sheared boundary conditions, a Poisson-solver and the timescale-parameterized cooling by Gammie (2001) are implemented.

2.1 System of Equations

Similarly to Hawley et al. (1995) and Gammie (2001) we consider a small comoving region within a flat disk in plane-polar geometry at a fixed radius $r=r_{0}$ and the angle $\varphi=\varphi_{0}+\Omega(r_{0})t$, where $t$ is the time and $\Omega(r_{0})=\Omega_{0}$ the angular velocity at $r_{0}$. For abbreviation we drop the index of $\Omega_{0}$ in the following. We introduce local coordinates $x=r-r_{0}$ and $y=r_{0}(\varphi-\Omega t)$ and expand the equations of motion to first order in $|x|/r_{0}$ (Goldreich & Lynden-Bell, 1965). This yields the following system of equations:

[TABLE]

Here, $\mathbf{v}$ is the velocity and $E$ the total energy density. Additionally, we have $P$ the vertically integrated pressure, $\mathbf{\Omega}=\Omega\mathbf{\hat{e}}_{z}$ the oriented angular velocity and $q=-\left.\frac{\mathrm{d}\ln{\Omega}}{\mathrm{d}\ln{r}}\right|_{r_{0}}$ the shearing parameter with $q=3/2$ in the keplerian case. $q=3/2$ is used throughout this paper. $\Phi_{\mathrm{sg}}$ is the potential induced by self-gravity in the midplane of the disk. It connects to the system via Poisson’s equation

[TABLE]

where $\delta(z)$ is the vertical density structure in the $z$-direction. The parametrized cooling is

[TABLE]

with $\gamma$ the two-dimensional heat capacity ratio.

In order to close the system of equations, we assume an ideal gas equation of state

[TABLE]

with the pressure $p$ and the specific internal energy $\epsilon$. The latter is related to the total energy density according to

[TABLE]

In order to allow for larger timesteps we split the constant velocity offset by introducing the fiducial velocities $\mathbf{v}^{\prime}=\mathbf{v}+q\Omega x\mathbf{\hat{e}}_{y}$, with $v_{x}^{\prime}=v_{x}$ and $v_{y}^{\prime}=v_{y}+q\Omega x\mathbf{\hat{e}}_{y}$ (Masset, 2000; Mignone et al., 2012).This yields the new system of equations:

[TABLE]

with the reduced energy density

[TABLE]

The Fargo scheme (Masset, 2000) was already implemented in Fosite (Jung, 2016) and can be enabled for $\mathbf{v}=\mathbf{v^{\prime}}+\mathbf{w}$, where $\mathbf{w}$ is an arbitrary solenoidal velocity (Mignone et al., 2012).

2.2 Fictitious Forces

The source terms for fictitious forces are implemented in two different ways according to the right hand side of eqs. 3b, 3c and eqs. 8, 8, respectively. The implementation differs depending on whether advection splitting is enabled or not.

2.3 Boundary Conditions

We chose periodic boundary conditions in $y$-direction and sheared periodic boundary conditions in $x$-direction (Hawley et al., 1995). With $f\in\{\Sigma,v_{x},v_{y}^{\prime},p\}$ and $L_{x}$, $L_{y}$ the size of the field in both directions we have

[TABLE]

in $y$-direction and

[TABLE]

in $x$-direction. If necessary, the data is linearly interpolated between adjacent cells (Hawley et al., 1995).

2.4 Self-gravity

We consider a razor thin disk, which means that $\delta(z)$ in eq. 4 is the $\delta$-distribution, and solve this equation via common practice FFT routines222We use the following libraries: FFTW (Frigo & Johnson, 2012) on scalar, fft (from MathKeisan package, see www.mathkeisan.com) on NEC-SX systems. similar to the way Gammie (2001) did. We shift the surface density field to the next periodic point, transform the result and calculate components of the potential via

[TABLE]

Again transforming and shifting back gives us the solution for $\Phi_{\mathrm{sg}}(x,y)$. The gravitational acceleration is obtained computing second order finite differences. We thereby also cut of all wavenumbers $\mathbf{|k|}=\pi/\left(\sqrt{2}\Delta x\right)$ (Gammie, 2001). Young & Clarke (2015) argue that this mimics the effect of a gravitational softening with length scale $\sim$0.3\text{,}\mathrm{H}$$, where $\mathrm{H}$ is the scale height of the disk.

2.5 Viscosity

In order to measure the viscous stresses in a disk, the $\alpha$ prescription can be used (Shakura & Sunyaev, 1973). Taking into account the cooling function in eq. 5 an analytical expression for the dimensionless parameter $\alpha$ can be found (Gammie, 2001):

[TABLE]

Simultaneously the $\alpha$-parameter can be calculated from the simulations assuming a Jacobi-like integral (Gammie, 2001), by splitting the stress-tensor $T_{xy}$ into a hydrodynamical part

[TABLE]

with the residual velocities $v_{x}^{\prime}$, $v_{y}^{\prime}$ and a gravitational part

[TABLE]

Together these yield

[TABLE]

Thus we can compare the numerical solution eq. 8p with the analytical prediction in eq. 8m.

2.6 Numerical Viscsosity

In order to resolve shocks, numerical schemes need to introduce non-linear terms. In SPH or finite difference schemes this is, for example, achieved by artificial viscosity. Fosite is based on finite volume methods which introduce limiters in order to handle the shock regions appropriately. Otherwise artificial numerical oscillations will occur in the vicinity of discontinuities. According to Godunov’s theorem (Godunov, 1959) this is a fundamental limitation of linear difference schemes with higher than first order spatial accuracy. Essentially, the limiter is switching to first-order accuracy in the shock regions. In smooth regions an ideal limiter should not affect the numerical scheme. In the real world, however, any limiter effects the solutions even though the slopes are moderate. Slope-limiters, like in our case, are applied to the slopes of the physical quantities itself whereas flux-limiters, which were used in Paardekooper (2012), are applied to the fluxes. With a slope $\Delta$, we introduce

[TABLE]

where $\phi(s)$ is the limiter which depends on the ratio of neighboring slopes $s$.

Within this text, we use two different limiters, VANLEER (van Leer, 1974) and SUPERBEE (Roe & Baines, 1982). VANLEER tends to lead to more numerical dissipation and resolves shocks less good than SUPERBEE. The latter one, however, can oversteepen the slopes in shallow regions. A more thoroughly discussion on this issue can be found in Toro (2009), who gives an overview of some limiters applied to linear advection. This goes conform with the guidelines of Hirsch (2007, ch. 8.3.4), who additionally states that ’limiters with continuous slopes generally favor convergence’. Whereas VANLEER has a continuous slope, limiters like SUPERBEE or MINMOD (Roe, 1986) show discontinuities in their slopes. The VANLEER limiter is described by

[TABLE]

and SUPERBEE by

[TABLE] &r≤0 max(0,min(2|r|,1),min(|r|,2))r>0

.

3 Tests & Verification

We test our code with epicyclic motion (Stone & Gardiner, 2010; Paardekooper, 2012) and linear theory (Gammie, 2001; Paardekooper, 2012). Thereby the linear theory test is being used extensively in order to obtain insights into the deviations between linear and nonlinear simulations.

3.1 Epicyclic Motion

The epicyclic motion test is adapted from Paardekooper (2012). We use exactly the same initial conditions with a resolution of $N_{x}=N_{y}=128$, a field size of $L_{x}=L_{y}=1$, an initial velocity perturbation $v_{x,0}=0.1c_{\mathrm{s}}$ and $v_{y,0}=0$, with a speed of sound $c_{\mathrm{s}}=0.01$. Then, we have oscillating solutions in the absence of any gradients

[TABLE]

which holds for the keplerian case.

Stone & Gardiner (2010) stress that the epicyclic energy

[TABLE]

should be conserved, which is achieved by a symmetric time integrator. Since Fosite is a purely explicit solver, we cannot conserve the epicyclic energy to round-off error. However, with the Dormand-Prince method (Dormand & Prince, 1980), we run the test over $1000\text{\,}\Omega^{-1}$ and get a relative error of order $10^{-5}$ for $E_{\mathrm{epi}}$. Additionally Dormand-Prince is of fifth order accuracy, compared to the Crank-Nicholson method, which is symmetric but only of second order. In fig. 1 the velocity evolution $v_{x}$ is shown over the first $100\text{\,}\Omega^{-1}$.

3.2 Linear Theory

In order to test especially the gravitational solver (but also the boundaries and the fictitious forces) the equations of motion are linearized in $\Sigma=\Sigma_{0}+\delta\Sigma\exp{\left(i\mathbf{k}(t)\cdot\mathbf{x}\right)}$ (similarly for other variables) and the evolution of the amplitude of a small density wave is examined. For the isothermal case this is derived in Gammie (1996). Without magnetic fields and viscosity an ordinary differential equation

[TABLE]

is the result, where

[TABLE]

is the initial perturbation in the vorticity $\xi$. See also Paardekooper (2012) for a thorough derivation. In Gammie (2001) the test is done including the energy equation, however we stay with the isothermal version to be comparable with Paardekooper (2012).

Similar to Paardekooper (2012) we use eq. 3.2 with an initial perturbation of $\delta\Sigma=$5\text{\times}{10}^{-4}\text{,}\Sigma_{\mathrm{0}}$$ and $\Sigma_{0}=1$. The wavenumbers are $k_{x,0}=-2(2\pi)$ and $k_{y,0}=2\pi$.

In fig. 2 we show convergence for the isothermal case with the VANLEER limiter. Typically, the first maximum is approximated well with all resolutions and limiters. At later stages the numerical viscosity plays a larger role with decreasing resolutions.

Fig. 3 shows the same setup, but using the SUPERBEE limiter. There, we do not see convergence, but an overestimation of the analytical theory at higher resolution. The error is largest for a resolution of $N=128$ and then declines for higher resolutions. For $N=64$ however, the behavior is similar to the one in fig. 2. Thus it seems that we have two opposing effects, one has a declining, viscous nature, which one should expect, since numerical dissipation usually declines with increasing resolutions. The other is strongly changing with the limiting function and it is only present or noticeable, if numerical dissipation is low enough, but leads to an excess in the density amplitude.

Since we expect that methods introduced to resolve shocks only add dissipation to the system, we would assume a lower amplitude over time. However, our test case inhibits no shocks, but shallow gradients. We therefore have strong indications that oversteepening caused by the SUPERBEE limiter is the reason for this effect (see e.g. Leveque, 2002; Hirsch, 2007, ch. 6; ch. 8).

It is always present at shallow regions but comes only into effect when resolution is high enough, because the inherent discretization error of the numerical scheme which leads to numerical dissipation decreases. Shallow, shock free regions are the expected environment where clumps form. To substantiate our thesis, we have a closer look at our calculation with gravity from fig. 2 (VANLEER) and fig. 3 (SUPERBEE) at a resolution of $N=64$. In fig. 4

a cut through the backshifted density field (by $\Delta y=q\Omega x\mathbf{\hat{e}}_{y}$) along the $x$-axis at $y=0$ and time $t=520\Omega^{-1}$ is shown. This is shortly after the first transition through $\delta\Sigma=0$. Comparing the solutions, we see first that the SUPERBEE simulation is overestimating the theoretical prediction, but this is mainly because the numerical solution is forerunning the analytical one at low resolutions (see fig. 3). Second and most importantly the wave shows a clear plateau at the top of the waves. This is the common behavior and a clear indicator333Together with steeper gradients, but the analytical solution gives us only the Fourier-component. for an oversteepened function, which leads to considerable, non-physical deviations from the analytical solution.

We want to highlight under which circumstances the problem occurs in order to make predictions for a non-linear run. Therefore in fig. 5

the test is shown twice as long for the highest resolution of $512\times 512$ with different limiters. Over these timescales it can be observed that the error amplifies increasingly and even if the first $10\text{\,}\Omega^{-1}$ seem to be converged, at later stages the numerical solution diverges. Without the gravitational module in fig. 6

the behavior remains the same, so gravitational forces are not the cause, but scale up initial errors.

We also tested the influence of a different time integrator (strong-stability-preserving Runge-Kutta method of same order, Macdonald 2003), a different Riemann solver (HLLC, Toro et al. 1994; Jung et al. 2015), with and without the advection splitting. We also carried out calculations with the energy equation. There are minor differences for different values of $\gamma$, but in all tests we performed, the limiting function plays the dominant role for the diverging nature of the test. We cannot exclude that the effect does not appear for the VANLEER limiter at higher resolutions and over longer timescales. Lastly, we also crosschecked the simulations with additional limiters, namely MINMOD (Roe, 1986) and OSPRE (Waterson & Deconinck, 1995), which also show convergence up to the resolution of $N=512$ (very slowly for MINMOD, however).

There is an important difference to mention between our code and for example the one used by Paardekooper (2012). Our code is limiting the slopes in order to achieve a TVD-system (Illenseer & Duschl, 2009) whereas Paardekooper (2012) is using flux limiting functions for the same purpose. Having a look at comparing simulations of linear advection, it seems that slope-limiters are more prone to oversteepening than flux-limiters (see Toro, 2009, ch. 13.10). The effect exists however, in both cases.

4 Nonlinear Results

The non-linear results with the VANLEER and SUPERBEE limiting functions are tested and the fragmentation boundary is investigated. If they both only add dissipation to the system, they should converge at high enough resolutions to the same macroscopic behavior, since the artificial heating added to the system should be low enough to exclude this as the major driver of non-convergence, as, for instance, Meru & Bate (2012) stated.

4.1 Standard run

Similar to Gammie (2001) we set $L=$320\text{,}\mathrm{G}\Sigma\mathrm{/}\Omega^{2}$\sim$50\text{,}\mathrm{H}$$ with a resolution of $N=1024$ and $\beta=2$ or $\beta=10$, respectively, and in order to be comparable we also use $\gamma=2.0$. Until $100\text{\,}\Omega^{-1}$, we can reproduce previous results with fragmentation for $\beta=2$ and the gravito-turbulent state for $\beta=10$.

For the following calculations, we use an initial $\beta=20$ until $t=$50\text{,}\Omega^{-1}$$ and afterwards decrease $\beta$ for the next $50\text{\,}\Omega^{-1}$ linearly to the desired value in order to avoid numerical problems from initializations (Paardekooper et al., 2011).

With the SUPERBEE limiter and $\beta=10$ fragmentation occurs at $t\simeq$300\text{,}\Omega^{-1}. We also can see the stochastic nature in this case. Therefore we ran the simulation with $\beta=15$ twice until $t=$1100\text{\,}\Omega^{-1}, one simulation shows fragmentation, the other not. However, there is a difference to Paardekooper (2012): Over the whole simulation time, fragmentation occurs most probably up to $\beta\lesssim 12$. At the same time, with VANLEER limiter we see no fragmentation at all for $\beta>2$ at a resolution of $N=1024$.

4.2 The non-fragmenting case

We compare the non-linear outcome for both limiters and a large value of $\beta$, where the gravito-turbulent state is present throughout the simulation time. In fig. 7 and fig. 8 the macroscopic behavior can be seen at the end of the run with $t=$1100\text{,}\Omega^{-1}$$, for the VANLEER and SUPERBEE limiter, respectively.

The simulation with SUPERBEE inhibits more structure and details. It shows small density enhancements, which did not lead to full fragmentation. These features are absent in the simulation with VANLEER limiter. Instead the gas is dominated by shocks running through the shearingsheet. Comparing our results with those of Gammie (2001, fig. 4) and Baehr & Klahr (2015, fig. 1 - left), (with different cooling function) look much more like fig. 7 at least in a qualitative way, whereas fig. 8 can be associated with results presented in Paardekooper (2012, fig. 6 - top).

In fig. 9

we show the maximum density $\Sigma_{\mathrm{max}}$ together with the measured value of $\alpha$. The simulation with the SUPERBEE limiter has generally a larger value for $\Sigma_{\mathrm{max}}$ and a larger amplitude of variations in $\Sigma$. If $\beta$ were smaller these excesses may lead to fragmentation. At the same time, $\alpha$ yields similar values, with slightly better results with SUPERBEE limiter. The VANLEER limiter shows a larger variance around a value slightly above the theoretical value of $\alpha=0.011$ calculated by eq. 8m. This is probably the case, because the VANLEER limiter is more dissipative, leading to additional heating. It can, however, also be seen that this value is not very large and shows that numerical heating has a minor, but not fully neglectable impact on the simulations.

In fig. 10

we see an overview of the fragmenting behavior. It shows fragmentation for $\beta<3$ throughout. For a resolution of $N=512$ both limiters have their fragmentation boundary directly at that value, very similar to Gammie (2001). At a larger resolution of $N=1024$ this behavior changes dramatically. The VANLEER limiter has its boundary value still at the same value for $\beta$ wheres the SUPERBEE limiter leads to fragmentation for at least up to $\beta=12$. At $N=2048$ SUPERBEE behaves very similar. For the VANLEER limiter the largest value with which we see fragmentation is $\beta=5$ at $N=2048$. It should be mentioned that we do not get the slow way of fragmentation in that case, but local collapse directly at $t=$100\text{,}\Omega^{-1}$$, where the final value of $\beta$ is reached. Interestingly, Paardekooper (2012) obtains the same maximum result in $\beta_{\mathrm{crit}}$ for MINMOD, which generally is more dissipative. However, we cannot see any stochastic fragmentation for VANLEER limiter. It may still exists but then in a much tighter region around $\beta_{\mathrm{crit}}\sim 3-5$. Generally the simulations at a resolution of $N=2048$ where not performed over the whole $1100\text{\,}\Omega^{-1}$ but only up to $300\text{\,}\Omega^{-1}$, which is most probably the reason that we do not get fragmentation for $\beta=12$ and the SUPERBEE limiter.

In order to check our results coherently in comparison to other studies we also made simulations with MINMOD limiter and with $\gamma=1.6$ ourself. MINMOD shows also fragmentation for $\beta<3$ at a resolution of $N=1024$ and no fragmentation for larger values. When choosing the smaller value for $\gamma$ we also see fragmentation for VANLEER limiter and $\beta=5$ at $N=1024$. This is consistent with Rice et al. (2005), who state that lower a $\gamma$ leads to larger critical $\beta_{\mathrm{crit}}$.

5 Discussion

We showed that the choice of the most appropriate limiting function is crucial to the non-linear outcome of the simulations and cannot be seen just as a numerical helper adding more or less dissipation to the system. We showed this for the special case of Goduvov type schemes which make use of slope limiters in order to achieve total-variation-diminishing.

Rice et al. (2014) showed that the cooling function needs to be adjusted in SPH simulations in order to converge. Young & Clarke (2015) however state that the gravitational smoothing leads to very different results in 2D simulations. Both publications change thereby physical quantities on small scales that are crucial for the fragmentation itself. Gravitational instabilities inherently depend on the cooling function and of course also on the gravitation itself. Thus, changing these quantities will always change the fragmentation behavior and could suppress any other physical or numerical effect on these scales. Still, we cannot exclude that these effects are the cause of fragmentation. Young & Clarke (2015) argue that smoothed 2D simulations cannot resolve the quasi-static collapse within their description. In other contexts, however, like stochastic fragmentation, this poses no problem as 2D calculations by Young & Clarke (2016) demonstrate. The question of how strong this effect evolves in three dimensions and together with gravity is an open one.

Additionally we showed that oversteepening is also present without gravity and generally it is not affected by the number of dimensions, because of its pure numerical nature. The discussed numerical effect is neither dependent on gravitation, nor on the cooling function. We can also exclude that fragmentation for larger $\beta\sim 10$ is connected to numerical heating, since our results yield reasonable values for $\alpha$. It takes place in exactly the region where clumps form and the error occurs on timescales which are comparable to the orbital timescales. Additionally it can be backtracked to a well known effect called oversteepening and is testable with the linear theory test, at least in shearingsheet simulations. SPH codes introduce artificial stresses, which are always present, not only in shock regions, and some of them are especially introduced in order to avoid numerical fragmentation, because of, for instance, the testile instability (Monaghan, 2000). This is why it cannot be excluded that SPH simulations also suffer from similar effects leading to an enhancement of density fluctuations which are in smooth regions scaled up by gravity.

This leads to an attractive property of oversteepening, as it may very well explain, why codes converge to different results or do not converge when $\beta>3$. The artificial viscosities or limiters, which are introduced in order to resolve shocks well, are handled completely differently by different kinds of codes and probably also handle smooth gradients completely different.

Here, we showed that the effect exists for Godunov-type schemes. This is for example also the case in Paardekooper (2012). In that case, the limiting functions are applied to other physical quantities (fluxes instead of slopes), which may change but not remove the effect. Our simulations run in concordance with the estimated behavior, because the effect tends to be stronger for slope limiters in the literature (Hirsch, 2007), which we also see for our fragmentation behavior compared to Paardekooper (2012), both in the linear theory test and in the non-linear outcome.

With the SUPERBEE limiter we can also see stochastic fragmentation like Paardekooper (2012). Running at a resolution $N=1024$ with $\beta=15$ over $t_{\mathrm{sim}}=$1100\text{,}\Omega^{-1}$$ twice, we observe fragmentation in one case only, but we did not investigate this issue in greater detail. With the other choice of the limiting function, we were not able to reproduce the stochastic features from Paardekooper (2012). They may indeed exist, but then in a much lower range then previously expected (around $\beta_{\mathrm{crit}}\lesssim 5$). This is also consistent with the results of Paardekooper (2012) himself, since he got $\beta_{\mathrm{crit}}\lesssim 5$ for the very dissipative limiter MINMOD, where most probably no oversteepening is present. With the VANLEER limiter, which is less dissipative, but also does not show oversteepening within the test, we get the same critical value for $\beta_{\mathrm{crit}}$. Furthermore, independent SPH simulations by Young & Clarke (2016) show that the effect is not as strong as expected before.

Another point important to mention is that the statement of Young & Clarke (2015) that there is no dependency on $\gamma$ in two-dimensional simulations is not coherent with our results. We can for example see fragmentation for $\beta=5$ and the VANLEER limiter, when choosing $\gamma=1.6$ and a resolution of $N=1024$. This is in concurrency with the results of Rice et al. (2005), where it was stated that smaller $\gamma$ lead to larger critical $\beta_{\mathrm{crit}}$. Additionally, heat capacity ratios are often compared between two and three dimensions equivalently and although they might behave similarly, they are not the same (Gammie, 2001).

6 Conclusions

In this paper, we investigated the impact of the limiting function on the fragmentation boundary in self-gravitating shearing-sheet simulations. By comparing the numerical outcome with linear theory, we showed that there is strong evidence that a unfortunate choice of the limiter can lead to an overestimation in the surface density in regions with shallow density gradients, which makes fragmentation easier for larger $\beta$. We accomplish this by using two-dimensional shearingsheet simulations and looking at deviations from to analytical results for linear theory. The behavior can be ascribed to a well known numerical effect called oversteepening, which occurs in regions with shallow density gradients. It is only observable if the resolution is high and the corresponding numerical diffusion low enough. Additionally we can exclude gravity as the driving force of the effect, but it is scaling the small initial errors up. For the limiting functions where oversteepening sets in place in the test, we see very strong differences in the non-linear outcome. In these calculations we need much smaller critical $\beta_{\mathrm{crit}}$, because of abundant density excesses in the field. We also cannot reproduce stochastic fragmentation with a well chosen limiter, however this might take place in a much smaller range for values of $\beta$.

Acknowledgements.

We thank an anonymous referee for very helpful comments on an earlier version of this paper.