---
abstract: 'We carry out simulations of gravitationally unstable discs using a Smoothed Particle Hydrodynamics (SPH) code and a grid-based hydrodynamics code, [fargo]{}, to understand the previous non-convergent results reported by [@Meru_Bate_resolution]. We obtain evidence that convergence with increasing resolution occurs with both SPH and [fargo]{} and in both cases we find that the critical cooling timescale is larger than previously thought. We show that SPH has a first-order convergence rate while [fargo]{} converges with a second-order rate. We show that the convergence of the critical cooling timescale for fragmentation depends largely on the numerical viscosity employed in both SPH and [fargo]{}. With SPH, particle velocity dispersion may also play a role. We show that reducing the dissipation from the numerical viscosity leads to larger values of the critical cooling time at a given resolution. For SPH, we find that the effect of the dissipation due to the numerical viscosity is somewhat larger than had previously been appreciated. In particular, we show that using a quadratic term in the SPH artificial viscosity (${\beta_{\rm SPH}}$) that is too low appears to lead to excess dissipation in gravitationally unstable discs, which may affect any results that sensitively depend on the thermodynamics, such as disc fragmentation. We show that the two codes converge to values of the critical cooling timescale, ${\beta_{\rm crit}}> 20$ (for a ratio of specific heats of $\gamma=5/3$), and perhaps even as large as ${\beta_{\rm crit}}\approx 30$. These are approximately $3-5$ times larger than has been found by most previous studies. This is equivalent to a maximum gravitational stress that a disc can withstand without fragmenting of $\alpha_{\rm GI,crit}\ \approx 0.013-0.02$, which is much smaller than the values typically used in the literature. It is therefore easier for self-gravitating discs to fragment than has been concluded from most past studies.'
author:
- |
    Farzana Meru$^{1,2,3}$[^1] and Matthew R. Bate$^3$\
    $^1$Institut für Astronomie, ETH Zürich, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland\
    $^2$Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany\
    $^3$School of Physics, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK\
bibliography:
- 'allpapers.bib'
title: 'On the convergence of the critical cooling timescale for the fragmentation of self-gravitating discs'
---

accretion, accretion discs - protoplanetary discs - planets and satellites: formation - gravitation - instabilities - hydrodynamics

Introduction
============

Historically, there have been two key quantities that have been used to determine whether a self-gravitating disc is likely to fragment. The first is the stability parameter [@Toomre_stability1964],

$$\label{eq:Toomre}
  Q=\frac{c_{\rm s}\kappa_{\rm ep}}{\pi\Sigma G},$$

where $c_{\rm s}$ is the sound speed in the disc, $\kappa_{\rm ep}$ is the epicyclic frequency, which for Keplerian discs is approximately equal to the angular frequency, $\Omega$, $\Sigma$ is the surface mass density and $G$ is the gravitational constant. [@Toomre_stability1964] showed that for an infinitesimally thin disc to fragment, the stability parameter must be less than a critical value, $Q_{\rm crit} \approx 1$.

[@Gammie_betacool] showed that in addition to the stability criterion above, the disc must cool at a fast enough rate. Using shearing sheet simulations, he showed that if the cooling timescale can be parametrized as

$$\label{eq:beta}
  \beta = t_{\rm cool}\Omega,$$

where

$$\label{eq:tcool}
  t_{\rm cool} = u \left ( \frac{{\rm d}u_{\rm cool}}{{\rm d}t} \right )^{-1},$$

$u$ is the specific internal energy and ${\rm d}u_{\rm cool}/{\rm d}t$ is the total specific cooling rate, then for fragmentation we require $\beta \lesssim \beta_{\rm crit}$. According to [@Gammie_betacool], $\beta_{\rm crit} \approx 3$ for a ratio of specific heats, $\gamma = 2$. [-@Rice_beta_condition] carried out three-dimensional simulations using a Smoothed Particle Hydrodynamics (SPH) code and showed that this cooling parameter is dependent on the equation of state. They showed that $\beta_{\rm crit} \approx 6-7$ for discs with $\gamma = 5/3$ and $\beta_{\rm crit} \approx 12-13$ for discs with $\gamma = 7/5$.

[@Gammie_betacool] and [@Rice_beta_condition] also showed that in a steady state disc where the dominant form of heating is that due to gravitational instabilities, since the gravitational stress in a disc can be linked to the cooling timescale in the disc using $$\label{eq:stress}
\alpha_{\rm GI} = \frac{4}{9} \frac{1}{\gamma (\gamma -1)} \frac{1}{\beta},$$ the rapid cooling required for fragmentation, ${\beta_{\rm crit}}$, can also be interpreted as a maximum gravitational stress that a disc can support without fragmenting, which they showed to be $\alpha_{\rm GI,crit} \approx 0.06$.

Recently, [@Meru_Bate_resolution] showed using SPH calculations of gravitationally unstable discs similar to those that have been performed by [@Rice_beta_condition] that the previous results on the critical cooling timescale had not converged. In particular, they found that the critical value of the cooling timescale, $\beta_{\rm crit}$, below which a disc would fragment increased linearly with increasing spatial resolution. This implied that the critical cooling rate might be much greater than that found from past studies (which would, for example, have implications for where in a real disc planets may form by the gravitational instability method). It also opened the question of whether or not a critical cooling rate indeed exists. Instead, a self-gravitating disc that is subject to a fixed cooling rate might fragment regardless of the value, given sufficient resolution (i.e. a disc may never be able to settle into a self-regulated state).

[@Lodato_Clarke_resolution] suggested that the non-convergent results may be an artefact of SPH artificially smoothing the density enhancements, or may be due to artificial viscosity if its effect was much larger than expected. [@Rice_cooling_convergence] suggested that the implementation of cooling in SPH may be to blame for the lack of convergence. However, [-@Paardekooper_convergence] showed using the two-dimensional grid-based hydrodynamics code, [fargo]{}, that the non-convergent problem was not specific to SPH. The source of non-convergence therefore cannot be constrained to SPH or to three-dimensional codes. [@Paardekooper_convergence] suggested that the boundary between the turbulent inner disc region and the smooth outer disc region (a consequence of starting the simulations with smooth initial conditions) may be the cause of the non-convergent results presented by [@Meru_Bate_resolution].

[@Bate_disc_frag_res] carried out SPH simulations of the collapse of molecular clouds to form protostars and discs. For particular initial conditions that lead to disc fragmentation, he noted that higher resolution simulations resulted in more fragments. Unlike the simulations of gravitationally unstable isolated protoplanetary discs discussed above, these simulations were of very early stage discs that formed prior to stellar core formation and were subject to rapid accretion from the surrounding molecular envelope. However, the interesting aspect here is that the fragmentation is more prevalent in higher resolution simulations of discs modelled using both isolated discs as initial conditions and using a parametrized cooling function [@Meru_Bate_resolution] as well as discs formed in molecular cloud collapse simulations using radiative transfer where no such smooth initial conditions are involved [@Bate_disc_frag_res].

[@Meru_Bate_resolution] expressed a concern about a lack of convergence with numerical resolution. However, even if convergence with numerical resolution is achieved, convergence between different numerical models is also important, i.e. the result is not believable if two different codes that can, in principle, model the same physical processes, produce physically different results.

In this paper, we present additional SPH results to those presented by [@Meru_Bate_resolution]. Rather than confining our investigations to a single hydrodynamics code, we also carry out a code comparison by performing further calculations using the grid-based Eulerian hydrodynamics code, [fargo]{}. We particularly focus on the dependence of the critical cooling timescale on the artificial viscosities employed in both codes.

In Section \[sec:numerics\] we describe the numerical methods adopted and discuss how numerical viscosity may affect the critical cooling timescale in discs in Section \[sec:num\_visc\]. We describe the simulations performed and present our results in Sections \[sec:sim\] and \[sec:results\], respectively. We discuss and make conclusions in Sections \[sec:disc\] and \[sec:conc\], respectively.

Numerical methods {#sec:numerics}
=================

Our SPH simulations are carried out using the exact same code as that used by [@Meru_Bate_resolution], originally developed by [@Benz1990], further developed by [-@Bate_Bonnell_Price_sink_ptcls] and [@Price_Bate_MHD_h] and parallelised using both OpenMP and MPI (see [@Meru_Bate_fragmentation] for details). Our simulations with a grid-based code are carried out using the Fast Advection in Rotating Gaseous Objects ([fargo]{}) two-dimensional fixed polar hydrodynamics code [@Masset_FARGO; @Baruteau_Masset_FARGO1; @Baruteau_Masset_FARGO2].

We include the heating effects in the disc due to work done on the gas and artificial viscosity. The cooling in the disc is taken into account using the cooling parameter, $\beta$ (equation \[eq:beta\]), first proposed by [@Gammie_betacool] which cools the gas on a timescale given by equation \[eq:tcool\]. For the SPH simulations carried out in this paper, we ensure that the timestepping is limited by the following timestep criterion (in addition to the Courant condition, the force condition and the viscous timestep condition; see [@SPH_Monaghan]):

$$\Delta t \leq C \frac{\beta}{\Omega},
  \label{eq:tcool_timestep}$$

where $C = 0.3$. [@Meru_Bate_fragmentation] show that this condition is adequate to ensure that the fragmentation results are not affected by the timestepping imposed. For the simulations performed using [fargo]{}, the timestep constraint (using $C = 1$) is also included for all simulations involving $\beta \lesssim 6$. This constraint appears as an additional term in the denominator of equation 15 of [@Masset_FARGO]. We have verified that for larger values of $\beta$ the effect of including this timestepping constraint is negligible.

To model the shocks in the discs, both codes use artificial viscosity. The SPH code uses the artificial viscosity method described by [@Chow_Monaghan_newAV] and [@Price_Monaghan_v_dot_r_AV], the implementation of which is summarised in equations \[eq:mom\_av\], \[eq:Pi\] and \[eq:mu\_2\] (see Appendix \[appendixA\] for details). The artificial viscosity is controlled by the parameters $\alpha_{\rm SPH}$ and $\beta_{\rm SPH}$. [fargo]{} uses the [@Vonneumann_Richtmyer_av] artificial viscous pressure with parameter $q$. We use the default values for the SPH and [fargo]{} artificial viscosity parameters of $(\alpha_{\rm SPH}, \beta_{\rm SPH}) = (0.1, 0.2)$ and $q = 1.41$, respectively. Where specified, we also vary the amount of SPH and [fargo]{} artificial viscosities to investigate their effects on the fragmentation boundary.

[fargo]{}’s grid is set up to use linear spacing in the azimuthal direction and logarithmic spacing in the radial direction. We use open boundary conditions at both the inner and outer radial boundaries and use a fixed gravitational softening length of $3 \times 10^{-4} H$, where $H$ is the vertical scaleheight of the disc, in all the [fargo]{} simulations.

The effect of numerical viscosity {#sec:num_visc}
=================================

In hydrodynamics codes, artificial viscosity is frequently applied to correctly capture shocks and to avoid post-shock oscillations. In trying to understand the evolution of gravitationally unstable discs, equation \[eq:stress\] seeks to link the dissipation rate in the disc to the cooling rate of the disc. However, equation \[eq:stress\] is derived for a steady-state disc which assumes that gravitational instability is the only heating source for the disc. In reality, in numerical simulations, there will be additional heating due to numerical dissipation.

In SPH, the artificial viscosity typically includes both a linear term controlled by $\alpha_{\rm SPH}$ and a quadratic term controlled by $\beta_{\rm SPH}$. The linear term provides a bulk viscosity which dissipates kinetic energy as particles approach each other to reduce particle oscillations following a shock, while the quadratic term is primarily important to stop particle interpenetration. [fargo]{}, on the other hand, only contains a quadratic term controlled by the artificial viscosity parameter, $q$, and provides a bulk viscosity.

The dissipation from the bulk viscosity in shocks (such as those generated by gravitational instability) is physical. In a gravitationally unstable disc, this provides the $\alpha_{\rm GI}$ term. However, in these discs the artificial viscosities will also provide some shear viscosity, the heating effects of which are not accounted for in equation \[eq:stress\].

SPH artificial viscosity {#sec:SPH_av}
------------------------

[@Monaghan1985] showed that in the continuum limit, the $\alpha_{\rm SPH}$ component mimics a Navier-Stokes viscosity with bulk and shear coefficients proportional to the resolution length (see [@Meglicki_av_continuum] and Appendix \[appendixA\]). This has been confirmed numerically [e.g. @Artymowicz_Lubow_av; @Lodato_Price_betaSPH]. The shear viscosity contributions from the linear and quadratic SPH terms due to the artificial viscosity can also be compared to the [@SS_viscosity] viscosity of the form $\nu = \alpha_{\rm SS} c_{\rm s} H$. For the SPH calculations discussed in this paper, where the viscosity is only applied between approaching particles, assuming a Keplerian disc it can be shown that (see Appendix \[sec:recent\_AV\], \[sec:new\_alpha\] and \[appendixD\])

$$\label{eqss}
\alpha_{\rm SS,lin} = \frac{31}{525} {\alpha_{\rm SPH}}\frac{h}{H}$$

and

$$\label{eqssquad}
\alpha_{\rm SS,quad} = \frac{9}{70 \pi} {\beta_{\rm SPH}}\left ( \frac{h}{H} \right )^2$$

where $\alpha_{\rm SS,lin}$ and $\alpha_{\rm SS,quad}$ are the contributions from the linear and quadratic terms, respectively, and $h$ is the particle smoothing length. Note that the coefficients in this SPH code are marginally different to some other SPH implementations. In SPH codes employing the older [@Monaghan_Gingold_art_vis] formalism, the coefficients would be $1/20$ and $3/(35 \pi)$ for the linear and quadratic terms, respectively (see Appendix \[sec:original\_AV\], \[sec:evaluating\] and \[appendixD\] for details).

For a [@SS_viscosity] disc model, the dissipation rate per unit mass is given by $\frac{9}{4} \alpha_{\rm SS} c_{\rm s}^2 \Omega$. In a purely gravitationally unstable disc the dissipation rate may be parametrized $\frac{9}{4} \alpha_{\rm GI} c_{\rm s}^2 \Omega$. However, using SPH we expect an additional heating due to numerical dissipation given by $\frac{9}{4} (\alpha_{\rm SS, lin} + \alpha_{\rm SS, quad}) c_{\rm s}^2 \Omega$. Thus, the combined heating rate per unit mass is expected to be

$$\frac{9}{4} \left ( \alpha_{\rm GI} + \frac{31}{525} {\alpha_{\rm SPH}}\frac{h}{H} + \frac{9}{70 \pi} {\beta_{\rm SPH}}\left ( \frac{h}{H} \right )^2 \right ) c_{\rm s}^2 \Omega.
\label{eq:av_dissipation}$$

Note that while the dissipation due to the quadratic term is often ignored when comparing the viscosity in SPH simulations of $\alpha$-discs, we show in Sections \[sec:diss\_nonGI\] and Appendix \[appendixD\] that its contribution can be non-negligible. In a numerical simulation, it is this heating rate that must be balanced by the imposed cooling (equation \[eq:tcool\]) for the disc to settle into a quasi-steady state. Thus, equation \[eq:stress\] can be rewritten as

$$\beta = \frac{4}{9} \frac{1}{\gamma (\gamma -1)} \frac{1}{(\alpha_{\rm GI} + \alpha_{\rm SS,lin} + \alpha_{\rm SS,quad})}.$$

For the disc to fragment, the combined heating must be insufficient to balance the cooling so that

$$\beta_{\rm crit} = \frac{4}{9} \frac{1}{\gamma (\gamma -1)} \frac{1}{(\alpha_{\rm GI, crit} + \alpha_{\rm SS, lin, crit} + \alpha_{\rm SS, quad, crit})},$$

where $\alpha_{\rm GI,crit}$ is the *true* value of the gravitational stress that a disc can support before fragmenting and $\alpha_{\rm SS, lin, crit}$ and $\alpha_{\rm SS, quad, crit}$ are the contributions to the heating from the artificial viscosity that allows a disc to fragment once $\beta \lesssim {\beta_{\rm crit}}$ for any one particular resolution.

Now, for a particular cooling self-gravitating disc calculation, let us suppose that there is some maximum gravitational stress that can be produced by the disc beyond which it will fragment. In this case, $\alpha_{\rm GI, crit}$ will be a constant, but $\alpha_{\rm SS, lin, crit}$ and $\alpha_{\rm SS, quad, crit}$ will decrease with increasing resolution. If $\alpha_{\rm SS, lin}$ and $\alpha_{\rm SS, quad}$ obey equations \[eqss\] and \[eqssquad\], then for a set of simulations with increasing resolution

$$\beta_{\rm crit} = \frac{4}{9} \frac{1}{\gamma (\gamma -1)} \left( \alpha_{\rm GI, crit} + \eta \frac{31}{525} \alpha_{\rm SPH} \frac{h}{H} + \zeta \frac{9}{70 \pi} {\beta_{\rm SPH}}\left ( \frac{h}{H} \right )^2 \right)^{-1},
\label{eqnumerical}$$

where $\eta$ and $\zeta$ are constants which we expect to equal 1.

In order to compare equation \[eqnumerical\] with the results of SPH simulations, we need to determine $h/H$ just before the fragmentation sets in. Assuming that the disc fragments when $Q \approx 1$, using equation \[eq:Toomre\] and noting $H=c_{\rm s}/\Omega$ and $\Omega^2=G M_*/R^3$, we have

$$H \approx \frac{\pi \Sigma R^3}{M_*},
\label{eq:H}$$

where $R$ is the radius in the disc and $M_*$ is the mass of the central object. The smoothing length in an SPH simulation for a disc that is resolved (i.e. $h<H$, which is true for all the simulations presented here at the radius of fragmentation) can be estimated, using equation 3 of [@Price_Bate_MHD_h] with density, $\rho \approx \frac{\Sigma}{2 H}$, as

$$\label{eqh}
h \approx 1.2 \left( \frac{2 H m_{\rm p}}{\Sigma} \right)^{1/3},$$

where $m_{\rm p}$ is the mass of an SPH particle. We use constant mass SPH particles and so the mass of an SPH particle is the disc mass divided by the number of SPH particles, $m_{\rm p} = M_{\rm d}/N_{\rm part}$. The ratio of the smoothing length to disc scaleheight can then be approximated to be

$$\frac{h}{H} \approx \frac{1.2}{\Sigma R^2} \left( \frac{2 M_*^2 M_{\rm d}}{\pi^2 N_{\rm part}} \right)^{1/3}.
\label{eq:h_H}$$

In Section \[sec:h\_H\] we verify that this is indeed the case for steady-state marginally stable discs that have a Toomre parameter, $Q \approx 1$. We therefore expect that

$$\begin{aligned}
\label{eq:betacrit_hH}
\beta_{\rm crit} &= \frac{4}{9} \frac{1}{\gamma (\gamma -1)} \left( \alpha_{\rm GI, crit} + \eta \frac{31}{525} \frac{1.2}{\Sigma R^2} \left( \frac{2 M_*^2 M_{\rm d}}{\pi^2 N_{\rm part}} \right)^{1/3} \alpha_{\rm SPH} \right. \nonumber \\
&\qquad \left. + \zeta \frac{9}{70 \pi} \frac{(1.2)^2}{\Sigma^2 R^4} \left ( \frac{2 M_*^2 M_{\rm d}}{\pi^2 N_{\rm part}} \right )^{\frac{2}{3}} {\beta_{\rm SPH}}\right )^{-1},\end{aligned}$$

i.e. we have three unknowns: $\alpha_{\rm GI,crit}$, $\eta$ and $\zeta$, since we can determine ${\beta_{\rm crit}}$ for any one resolution. If $\eta$ and $\zeta$ are unity, then the second and third terms in the denominator in equation \[eq:betacrit\_hH\] each have values of $\approx O(10^{-3})$ for the discs studied by [@Rice_beta_condition] and [@Meru_Bate_resolution] simulated with 250,000 particles (using the parameters described in Section \[sec:sim\], and given that the fragmentation occurs in the outer parts of the disc – see [@Meru_Bate_fragmentation]). The contribution from these terms are approximately a factor of $O(10)$ smaller than the original estimate of $\alpha_{\rm GI,crit} \approx 0.06$ [@Gammie_betacool; @Rice_beta_condition]. Thus it was assumed in earlier SPH studies that the heating due to artificial viscosity would be negligible compared to the dissipation due to gravitational instabilities (see Appendix A of [@Lodato_Rice_original]). Note, however, that if the SPH artificial viscosity plays a significant role then the ${\alpha_{\rm SPH}}$ term scales linearly with the smoothing length, $h$, such that the convergence of $\beta_{\rm crit}$ towards the true value is expected to be first order in $h$ (i.e. very slow as the numerical resolution is increased). On the other hand, if the dominant term is the ${\beta_{\rm SPH}}$ term, the convergence will be faster since the SPH artificial viscosity scales quadratically with the smoothing length.

[fargo]{} artificial viscosity {#sec:FARGO_av}
------------------------------

Most grid-based hydrodynamical codes are second order and do not have a linear viscosity. Therefore, one would expect that their rate of convergence towards the true value of $\beta_{\rm crit}$ will be second order in spatial resolution and thus possibly faster than SPH codes. Indeed, [fargo]{} uses the [@Vonneumann_Richtmyer_av] artificial viscous pressure given by [see also @Numerical_Methods]

$$\label{eq:visc_press}
q^2 \rho (\Delta x)^2 \left | \frac{\partial v}{\partial x} \right |^2$$

where $\Delta x$ is the cell size, $\rho$ is the density and $q = l/\Delta x$ is a constant which indicates the number of grid cells over which the shock is spread and whose value is dependent on the numerical scheme and is usually $0.05 \le q \le 2$ ($l$ indicates the strength of the artificial viscosity). This is a bulk viscosity. In a cylindrical code, if the gas travels in circles, there should be no shear viscosity at all. However, in a gravitationally unstable disc, this will not be the case and there will be some shear viscosity (that arises from the bulk viscosity) which we expect will scale in roughly the same way, i.e. proportional to the square of the size of the grid cell. Assuming the shear rate, $\left | \partial v / \partial x \right |$, is approximately Keplerian, then

$$\left | \frac{\partial v}{\partial x} \right | \sim \left | R \frac{{\rm d} \Omega}{{\rm d} R} \right |.
\label{eq:shear}$$

Using equation \[eq:shear\] and equating equation \[eq:visc\_press\] to the shearing force per unit area as defined in equation \[eq:FoverA\] yields a kinematic viscosity due to the artificial viscosity given by

$$\nu_{\rm av, {\sc fargo}} = q^2 (\Delta x)^2 \left | R \frac{{\rm d} \Omega}{{\rm d} R} \right | \approx \frac{3}{2} q^2 (\Delta x)^2 \Omega.
\label{eq:fargo_av}$$

where the final approximation assumes a Keplerian flow. This gives a [@SS_viscosity] type viscosity of the form $\nu = \alpha_{\rm SS} c_{\rm s} H$ with

$$\label{eqss_fargo}
\alpha_{\rm SS, {\sc fargo}} = \frac{3}{2} q^2 \left ( \frac{\Delta x}{H} \right )^2.$$

Analogous to the derivation of equation \[eqnumerical\] for SPH artificial viscosity this yields

$$\label{eqgrid}
\beta_{\rm crit} = \frac{4}{9} \frac{1}{\gamma (\gamma -1)} \frac{1}{ \left ( \alpha_{\rm GI, crit} + \xi \frac{3}{2} q^2 \left(\frac{\Delta x}{H}\right)^2 \right ) },$$

where we expect that $\xi$ is unity. Substituting for the scaleheight using equation \[eq:H\] yields a formula for ${\beta_{\rm crit}}$ that is equivalent to equation \[eq:betacrit\_hH\] for a grid-based code:

$$\beta_{\rm crit} = \frac{4}{9} \frac{1}{\gamma (\gamma -1)} \left( \alpha_{\rm GI, crit} + \xi \frac{3}{2} q^2 \left(\frac{\Delta x M_*}{\pi \Sigma R^3}\right)^2 \right )^{-1}.
\label{eq:betacrit_fargo}$$

Since the [fargo]{} simulations use a logarithmic grid we take the cell size at a radius of $R = 22$ au, i.e. close to the edge of the disc where we would expect fragmentation to occur when the cooling timescale is close to the critical value for any one particular resolution. We find that the cell size at this radius scales as

$$\Delta x \approx \frac{125}{N_{\rm cells}^{\frac{1}{2}}} \rm ~au.$$

where $N_{\rm cells}$ is the total number of cells used in the simulation. We therefore have two unknowns: $\alpha_{\rm GI, crit}$ and $\xi$, since we can determine ${\beta_{\rm crit}}$ at any one resolution.

Equation \[eqgrid\] shows that if the artificial viscosity plays a significant role in the dissipation in the disc, the convergence is expected to be second order in spatial resolution, i.e. potentially faster than with SPH.

Artificial viscosity effects with resolution {#sec:av_res}
--------------------------------------------

In an SPH code when the number of particles increases, $h/H$ is reduced and therefore as the resolution increases, $\alpha_{\rm SS, lin}$ and $\alpha_{\rm SS, quad}$, and thus the heating due to artificial viscosity, tend to zero. Similarly, as the resolution is increased in a grid-based code, the cell size decreases for any one problem, and the contribution to the dissipation from the artificial viscosity decreases. In the limit of infinite resolution, equations \[eqnumerical\] and \[eqgrid\] will return to equation \[eq:stress\]. But for a finite resolution, the value of $\beta_{\rm crit}$ [*obtained from a numerical simulation should always be smaller than the true value*]{}.

It is also important to note that simply reducing the value of ${\alpha_{\rm SPH}}$, ${\beta_{\rm SPH}}$ or $q$ is not necessarily a sufficient way in which to decrease the numerical dissipation and obtain the *true* value of ${\beta_{\rm crit}}$. By inspection of equations \[eqnumerical\] and \[eqgrid\] we might naively assume this to be the case. However, reducing these values may mean that shocks are treated inaccurately. For example, the shocks may not be spread over a large enough lengthscale to model them numerically and/or there may be post-shock oscillations that are eventually damped, resulting in dissipation. [@Lodato_Price_betaSPH] showed using SPH that setting ${\alpha_{\rm SPH}}= 0$ counter-intuitively led to a larger amount of dissipation. [@Price_Federrath_betaSPH] found that if an adequate value of ${\beta_{\rm SPH}}$ is not used, particle interpenetration may occur. In their case, they stated that this makes very little difference to their dissipation rate since their linear term dominates almost everywhere. However, their simulations explored a different regime to that being explored here and with a different artificial viscosity scheme. Their simulations employed the [@MM_viscosity] artificial viscosity switch where the value of ${\alpha_{\rm SPH}}$ ranges between 0.05 and 1.0 (the higher value being implemented close to shocks). Their simulations were of high Mach number shocks (${\mathscr M} = 10$) and so a large part of their simulations would require the use of the higher value of ${\alpha_{\rm SPH}}$ and thus this would dominate the dissipation. The simulations performed by [@Rice_beta_condition] and [@Meru_Bate_resolution] used a fixed value of ${\alpha_{\rm SPH}}= 0.1$ - such a low value opens up the possibility of the quadratic term being important and therefore decreasing the value of ${\beta_{\rm SPH}}$ may affect the overall dissipation rate and thus the fragmentation outcome.

In [fargo]{}, if $q$ is set to zero, there will be no controlled numerical dissipation (e.g. to capture shocks or other disturbances at the grid scale). However, the code will still have some level of numerical diffusion and dissipation which is not controllable other than that it too should decrease with increasing resolution.

Simulations {#sec:sim}
===========

The disc and star properties used to carry out the simulations in this paper are exactly the same as those used by [@Rice_beta_condition] and [@Meru_Bate_resolution]: a $0.1 {{\rm M_{\odot}}}$ disc surrounding a $1 {{\rm M_{\odot}}}$ star. The SPH simulations span a radial range, $0.25 \le R \le 25$ au, while the [fargo]{} simulations span a radial range, $1 \le R \le 25$ au (only marginally different to the SPH simulations for numerical reasons).

The initial surface mass density and temperature profiles are $\Sigma \propto R^{-1}$ and $T \propto R^{-1/2}$, respectively, and the temperature is normalised so that the minimum initial Toomre stability value at the outer edge of the disc, ${Q_{\rm min}}=2$. The discs are modelled with a ratio of specific heats, $\gamma = 5/3$.

Table \[tab:res\_sim\] shows a summary of the initial SPH simulations and the key fragmenting results carried out by [@Meru_Bate_resolution] (obtained from their Table 1) as well as those in this paper (bold text). We supplement the [@Meru_Bate_resolution] results by carrying out an additional simulation using 2 million particles with $\beta = 9$ and three additional simulations using 16 million particles with $\beta = 12$, 15 and 20.

Table \[tab:fargo\_sim\] summarises the initial simulations carried out using [fargo]{} and the key fragmentation results. We perform simulations, using $q = 1.41$, at five different resolutions and determine the critical value of $\beta$ at each of these resolutions. The lowest resolution simulations are carried out using 768 and 256 grid cells in the azimuthal and radial directions, respectively. We then increase the linear resolution by factors of 2, 4, 8 and 16 in both the azimuthal and radial directions.

The simulations were run either for at least 6 outer rotation periods (ORPs) or until the discs fragmented. Fragments are defined as regions whose surface mass densities are at least two orders of magnitude denser than their surroundings. In addition, we ensure that the fragments survive for at least one rotation to verify that they do not shear apart.

To investigate the effects of the different components of artificial viscosity in SPH, we carry out a number of simulations where we vary the values of ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$ (see Table \[tab:alpha\_beta\_inves\] for details). We perform a suite of simulations using 250,000 particles. Firstly, we set ${\alpha_{\rm SPH}}= 0.1$ and vary the value of ${\beta_{\rm SPH}}$ (Table \[tab:alpha\_beta\_inves\], top section). Secondly, we set ${\beta_{\rm SPH}}= 2.0$ and vary the value of ${\alpha_{\rm SPH}}$ (Table \[tab:alpha\_beta\_inves\], middle section). We then carry out simulations with $({\alpha_{\rm SPH}}, {\beta_{\rm SPH}}) = (0.1, 2.0)$ using 31,250, 250,000 and 2 million particles to determine the effect that changing the value of ${\beta_{\rm SPH}}$ has on the fragmentation boundary at each of these resolutions (Table \[tab:alpha\_beta\_inves\], bottom section).

In addition, we perform a number of SPH simulations without self-gravity (see Table \[tab:non-GI\] for details) using various values of ${\alpha_{\rm SPH}}$, ${\beta_{\rm SPH}}$ to compare the measured dissipation to the analytically expected values in equation \[eq:av\_dissipation\] (also see Appendix \[appendixD\]).

Finally, we investigate the effects that artificial viscosity in [fargo]{} has on the critical cooling timescale by varying the value of $q$ in equation \[eq:visc\_press\] between 0 and 2.5 (see Table \[tab:fargo\_av\]) for discs modelled using 786,432 grid cells (512 and 1536 cells in the radial and azimuthal directions, respectively). We then carry out simulations using an artificial viscosity parameter, $q = 0.5$, at all but the lowest resolutions considered in this paper to determine the effect that this has on the fragmentation boundary.

Results {#sec:results}
=======

The convergence rate of ${\beta_{\rm crit}}$ with SPH {#sec:nonconverge}
-----------------------------------------------------

![Graph of $\beta$ against resolution of the non-fragmenting (open squares) and fragmenting (solid triangles) SPH simulations. This figure contains the results presented by @Meru_Bate_resolution [their Figure 3] as well as the new simulations highlighted in Table \[tab:res\_sim\]. The solid line, obtained by fitting equation \[eq:fit\], shows a dividing line between the fragmenting and non-fragmenting cases and the grey region is where fragmentation can take place. The graph shows clear evidence of convergence of results with increased resolution. These simulations are carried out with (${\alpha_{\rm SPH}}$, ${\beta_{\rm SPH}}$) = (0.1, 0.2). The convergence rate is first order with spatial resolution. The dotted line (which coincides well with the solid line) is obtained by fitting equation \[eq:betacrit\_hH\].[]{data-label="fig:res_beta"}](Fig1.pdf){width="1.0\columnwidth"}

![image](Fig2.pdf){width="2.0\columnwidth"}

Table \[tab:res\_sim\] summarises the results of the SPH simulations, using ${\alpha_{\rm SPH}}$, ${\beta_{\rm SPH}}$) = (0.1, 0.2), carried out by [@Meru_Bate_resolution] and those performed for this paper. [@Meru_Bate_resolution] found *borderline* simulations which they defined to be discs which showed signs of fragmentation but the fragments sheared apart rapidly (within 1 ORP) and no further signs of fragmentation were seen. We find that *borderline* simulations can in fact range a span of $\beta$ values. However, since ultimately they are discs that do not end up fragmenting, we now simplify this terminology and refer to them as non-fragmenting simulations. Figure \[fig:res\_beta\] shows a summary of the SPH results. This figure is the same as Figure 3 of [@Meru_Bate_resolution] but with the the addition of the new SPH results presented in this paper. [@Meru_Bate_resolution] found no evidence for convergence, but the addition of the new high resolution calculations now provides evidence for a very slow convergence of $\beta_{\rm crit}$ with increasing resolution. Figure \[fig:image\_SPH\_beta0.2\] shows two of the highest resolution SPH simulations (16 million particles) carried out using $\beta = 10$ and $\beta = 12$. It can clearly be seen that at this resolution, fragmentation occurs for $\beta = 10$ but not for $\beta = 12$. To estimate the rate of convergence we fit a formula of the form

$$\beta = \frac{{\beta_{\rm crit}}}{1 + \lambda l^{\sigma}}
\label{eq:fit}$$

where $l$ is the linear spatial resolution, $\lambda$ is a constant and $\sigma$ is the convergence rate. For SPH, we simply take $l \propto N_{\rm part}^{-\frac{1}{3}}$. We fit this formula to the values of $\beta$ in Figure \[fig:res\_beta\] that lie half way between the lowest non-fragmenting value of $\beta$ and the highest fragmenting value of $\beta$ for each numerical resolution, i.e. the fragmentation boundary. We find that a good fit is obtained with ${\beta_{\rm crit}}= 15.6 \pm 1.0$ and $\sigma = 1.08 \pm 0.05$. The value of $\sigma$ shows that the rate of convergence is first order in spatial resolution. For the benefit of understanding the results presented in Section \[sec:frag\_bdry\_optimumSPH\] (which only have data points at the lowest three resolutions) we fit equation \[eq:fit\] to the data presented in Figure \[fig:res\_beta\] but *exclude* the highest resolution simulations. We find ${\beta_{\rm crit}}= 17.4$ and $\sigma = 1.03$, so excluding the last point does not alter the fit significantly due to the slow convergence rate of SPH.

In addition, we also fit equation \[eq:betacrit\_hH\] to this data. We find that $\alpha_{\rm GI, crit} = 0.024 \pm 0.001$, $\eta = 21.1 \pm 1.3$ and $\zeta = 1.7 \pm 1.5$. In the limit of infinite resolution, this value of the critical gravitational stress obtained is equivalent to a critical cooling timescale of ${\beta_{\rm crit}}\approx 17$. However, we point out that while $\zeta$ is reasonably close to unity, the value of $\eta$ is very large. This is suggestive of an additional source of dissipation present in the simulations over and above what we expect from artificial viscosity in a shear dominated disc.

The convergence rate of ${\beta_{\rm crit}}$ with [fargo]{} {#sec:FARGO}
-----------------------------------------------------------

![Graph of $\beta$ against resolution of the non-fragmenting (open squares) and fragmenting (solid triangles) [fargo]{} simulations carried out using $q = 1.41$. The solid line, obtained by fitting equation \[eq:fit\], shows a dividing line between the fragmenting and non-fragmenting cases and the grey region is where fragmentation can take place. The graph shows clear evidence of convergence of results with increased resolution. The convergence rate is second order with spatial resolution. The dotted line (which coincides well with the solid line) is obtained by fitting equation \[eq:betacrit\_fargo\].[]{data-label="fig:res_beta_fargo"}](Fig3.pdf){width="1.0\columnwidth"}

Table \[tab:fargo\_sim\] and Figure \[fig:res\_beta\_fargo\] summarise the results using the grid-based code, [fargo]{}. As with the SPH results, we also see that as the resolution increases numerical convergence does appear to take place. Figure \[fig:image\_fargo\_q1.41\] shows the surface mass density rendered images of the discs modelled using [fargo]{} at resolutions of 786,432, 3.1 million, 13 million and 50 million grid cells. The top two panels show that using a cooling timescale of $\beta = 10$, the discs do not fragment using 786,432 grid cells but when the resolution is increased to 3.1 million grid cells, fragmentation is seen. Similarly, using a cooling timescale of $\beta = 18$, the disc modelled with 13 million grid cells does not fragment whereas that modelled at the higher resolution of 50 million grid cells does fragment.

Again, we use equation \[eq:fit\] to estimate the rate of convergence. For [fargo]{}, we simply take $l \propto N_{\rm cells}^{-\frac{1}{2}}$, where $N_{\rm cells}$ is the number of grid cells. Note that the linear resolution is inversely proportional to the square root of the number of cells because the calculation is two dimensional. We fit this formula to the fragmentation boundary in Figure \[fig:res\_beta\_fargo\]) as done for the SPH results. We find that a good fit is obtained with ${\beta_{\rm crit}}= 22.3 \pm 2.3$ and $\sigma = 2.03 \pm 0.36$. The value of $\sigma$ shows that the rate of convergence is second order in spatial resolution. We then fit the data using equation \[eq:betacrit\_fargo\] and find that $\alpha_{\rm GI, crit} = 0.018 \pm 0.001$ and $\xi = 0.87 \pm 0.08$. In the limit of infinite resolution, this value of $\alpha_{\rm GI, crit}$ is equivalent to a critical cooling timescale, ${\beta_{\rm crit}}\approx 22$ (using equation \[eq:stress\]), similar to the value obtained using equation \[eq:fit\].

Thus the value of ${\beta_{\rm crit}}$ converges more rapidly using [fargo]{} than SPH. In Section \[sec:FARGO\_av\] we note that if artificial viscosity plays a significant role in the determination of ${\beta_{\rm crit}}$ then [fargo]{} might be expected to display second-order convergence since it only applies a quadratic artificial viscosity. On the other hand SPH includes both linear and quadratic artificial viscosities. If the linear term is dominant, this may lead to first-order convergence. This implies that artificial viscosity may be significant in determining ${\beta_{\rm crit}}$. We note that [fargo]{} appears to converge to a higher value of ${\beta_{\rm crit}}$ than SPH, but this result may also be caused by the different artificial viscosities. Therefore, in the following sections we investigate the dependence of ${\beta_{\rm crit}}$ on the artificial viscosities applied in both codes.

The effect of SPH artificial viscosity on convergence {#sec:SPH_av_effect}
-----------------------------------------------------

In Section \[sec:num\_visc\] we present analytical arguments that suggest that artificial viscosity may play a role in the numerically determined value of the critical cooling timescale. We show that the contribution to the dissipation due to the artificial viscosity is expected to decrease with increasing resolution (Appendix \[appendixD\] and equation \[eq:av\_dissipation\]). Therefore, if the slow convergence can be attributed to SPH artificial viscosity, this may be the reason why [@Meru_Bate_resolution] found that ${\beta_{\rm crit}}$ increases with increasing resolution and is a plausible explanation as to why the results presented in Section \[sec:nonconverge\] show a slow convergence. We test the role that artificial viscosity plays on the fragmentation of self-gravitating discs by varying the values of ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$ separately.

### The effect of ${\beta_{\rm SPH}}$ on the critical cooling timescale {#sec:betaSPH_test}

![Graph of $\beta$ against ${\beta_{\rm SPH}}$ of the non-fragmenting (open squares) and fragmenting (solid triangles) simulations carried out using 250,000 particles and ${\alpha_{\rm SPH}}= 0.1$. The solid line, included by eye, shows a dividing line between the fragmenting and non-fragmenting cases and the grey region is where fragmentation can take place. The graph shows an S-shaped curve: any additional particle oscillation that may exist is stopped using ${\beta_{\rm SPH}}\approx 2$ since the effect on the critical cooling timescale does not change above this value; at lower values of ${\beta_{\rm SPH}}$, the critical cooling timescale is smaller as the additional particle oscillation results in excess dissipation that needs to be overcome before fragmentation can take place. At very low values of ${\beta_{\rm SPH}}$ either so much particle oscillation occurs (either at the edge of the shock front or due to particle interpenetration) or that the ${\alpha_{\rm SPH}}$ term dominates, that the effects of lowering ${\beta_{\rm SPH}}$ does not result in more dissipation.[]{data-label="fig:S-shaped"}](Fig5.pdf){width="1.0\columnwidth"}

In Section \[sec:nonconverge\] we show that the first-order convergence seen for the SPH results suggests that the ${\alpha_{\rm SPH}}$ term may be responsible. However, as mentioned in Section \[sec:av\_res\], if the optimum value of ${\beta_{\rm SPH}}$ is not used (i.e. a value that minimises numerical dissipation) additional dissipation may occur and affect the fragmentation boundary. Therefore, while not immediately obvious from the results in Section \[sec:nonconverge\], the value of the ${\beta_{\rm SPH}}$ term may affect the fragmentation conclusions. Table \[tab:alpha\_beta\_inves\] (top panel) and Figure \[fig:S-shaped\] summarise the results of the simulations carried out to investigate what effect the value of ${\beta_{\rm SPH}}$ has on the critical cooling timescale using 250,000 particles and maintaining a fixed value of ${\alpha_{\rm SPH}}= 0.1$. It can be seen that the shape of the fragmenting/non-fragmenting boundary line appears to follow a somewhat S-shaped curve. At high values of ${\beta_{\rm SPH}}$, any potential particle interpenetration is appropriately dealt with as changing the value of ${\beta_{\rm SPH}}$ from 2 to 4 has no effect on the critical cooling timescale. As ${\beta_{\rm SPH}}$ is reduced, particle interpenetration and additional particle velocity dispersion can occur since the appropriate amount of the quadratic term of the artificial viscosity is not used. Eventually, these oscillations are damped down by the ${\alpha_{\rm SPH}}$ term resulting in dissipation. Consequently, a more rapid cooling is required to overcome the additional dissipation resulting in smaller critical cooling values. At very low values of ${\beta_{\rm SPH}}$ the critical cooling timescale remains the same. This may be due to one of two reasons: 1) a “saturation” of additional oscillations occurs such that reducing ${\beta_{\rm SPH}}$ further does not increase the dissipation - by this we mean that the cause of the oscillation (i.e. the incorrect modelling at the edge of the shock front and particle interpenetration) is so high that any reduction in ${\beta_{\rm SPH}}$ cannot cause more oscillation to occur; or 2) at such low values of ${\beta_{\rm SPH}}$, the linear artificial viscosity term dominates the dissipation such that any additional particle interpenetration does not increase the overall dissipation by much (an effect also noted by [@Price_Federrath_betaSPH]).

In any case, Figure \[fig:S-shaped\] shows that the amount of ${\beta_{\rm SPH}}$ that is required to deal with the particle oscillations in this problem is $\approx 2$. This value ensures that the dissipation resulting from artificial viscosity is as low as possible (since the value of ${\beta_{\rm crit}}$ that results is higher) and therefore is likely to give a result that is “closer to the real answer”. While these simulations are carried out at a single resolution (250,000 particles), [@Matthews_thesis] shows that ${\beta_{\rm SPH}}\approx 2$ is sufficient to stop particle interpenetration for Mach numbers across a shock of $\mathscr M \approx 3$. In our simulations we find that the Mach numbers across the shock are up to $\approx 3$. Note that the simulations presented by [@Rice_beta_condition] and [@Meru_Bate_resolution] were all carried out using ${\beta_{\rm SPH}}=0.2$. Consequently, the calculations in both papers will have been affected by this and thus the converged value of ${\beta_{\rm crit}}$ is expected to be even higher than that suggested by Figure \[fig:res\_beta\].

### The effect of ${\alpha_{\rm SPH}}$ on the critical cooling timescale {#sec:alphaSPH_test}

![Graph of $\beta$ against ${\alpha_{\rm SPH}}$ of the non-fragmenting (open squares) and fragmenting (solid triangles) simulations carried out using 250,000 particles and with ${\beta_{\rm SPH}}= 2.0$. The solid line, included by eye, shows a dividing line between the fragmenting and non-fragmenting cases and the grey region is where fragmentation can take place. At high viscosities, the dissipation is higher (Figure \[fig:diss\_alphaSPH\_betaSPH\], left panel), resulting in a faster cooling, i.e. a lower value of $\beta$, required to overcome the dissipation and cause fragmentation. As ${\alpha_{\rm SPH}}$ is decreased, the dissipation also decreases requiring a slower cooling for fragmentation. At very low values of ${\alpha_{\rm SPH}}$, additional dissipation occurs resulting in a lower value of ${\beta_{\rm crit}}$. This may be due to additional particle oscillation as the shocks are not modelled adequately with such a low value of ${\alpha_{\rm SPH}}$.[]{data-label="fig:beta_alphaSPH"}](Fig6.pdf){width="1.0\columnwidth"}

In Section \[sec:SPH\_av\] we show that the dissipation due to the linear term in the artificial viscosity may play a part in the fragmentation results, and hence the value of ${\beta_{\rm crit}}$ at any one resolution. In Section \[sec:nonconverge\] we show that this may indeed have been the case by considering the rate of convergence with increasing resolution. In this section we maintain a fixed resolution using 250,000 particles (i.e. keep the value of $h/H$ constant) and use a fixed value of ${\beta_{\rm SPH}}= 2.0$, but vary the value of ${\alpha_{\rm SPH}}$ to confirm that equation \[eqss\] does indeed play a part in determining the fragmentation boundary. Table \[tab:alpha\_beta\_inves\] (middle panel) and Figure \[fig:beta\_alphaSPH\] summarise the simulations performed and the key fragmentation results. At higher values of ${\alpha_{\rm SPH}}$ the dissipation due to the artificial viscosity is expected to increase. Consequently, the cooling required to overcome this additional dissipation is larger and as a result, the critical cooling timescale for fragmentation is lower. As the ${\alpha_{\rm SPH}}$ term is decreased, the amount of dissipation also decreases and thus the cooling does not have to be so rapid, resulting in a higher critical cooling timescale. At values below ${\alpha_{\rm SPH}}\approx 0.1$, however, the dissipation increases once again as there is not enough artificial viscosity to remove the oscillations at shock fronts. Examining the velocity dispersion of particles in the disc around their expected almost Keplerian values, we find that with very low viscosity, the velocity dispersion of the particles increases. The particles are ‘jostled’ by one another when the viscosity is lower and the relative motions grow larger (also see Section \[sec:diss\_nonGI\]). Therefore, even though the value of ${\alpha_{\rm SPH}}$ is decreased, the dissipation increases as the small amount of viscosity that is present tries to damp these larger velocities. This suggests that ${\alpha_{\rm SPH}}\approx 0.1$ is a happy medium whereby it minimises the dissipation and avoids large oscillations at shock fronts.

### Determining the fragmentation boundary using optimum values of ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$ {#sec:frag_bdry_optimumSPH}

![Graph of $\beta$ against resolution for non-fragmenting (open squares) and fragmenting (solid triangles) simulations carried out with (${\alpha_{\rm SPH}}$, ${\beta_{\rm SPH}}$) = (0.1, 2.0) using 31,250, 250,000 and 2 million particles. The solid line, obtained by fitting equation \[eq:betacrit\_hH\], shows a dividing line between the fragmenting and non-fragmenting cases and the grey region is where fragmentation can take place. The region to the right of $\approx 2$ million particles has not been shaded in as it is unclear from these results alone what the shape of the dividing line would lie. The dotted line shows the fit using equation \[eq:betacrit\_hH\] when $\zeta$ is set to the minimum value it can be, i.e. unity.[]{data-label="fig:factor1.5"}](Fig7.pdf){width="1.0\columnwidth"}

![Graph of ${\beta_{\rm crit}}$ against resolution of the SPH simulations carried out using ${\beta_{\rm SPH}}= 0.2$ (squares) and ${\beta_{\rm SPH}}= 2.0$ (triangles). The value of ${\alpha_{\rm SPH}}$ is 0.1. It can be seen that the effect of increasing ${\beta_{\rm SPH}}$ to 2.0 (i.e. to a value that minimises the additional dissipation) is to increase the critical cooling timescale.[]{data-label="fig:SPH_betaSPH0.2_2_comparison"}](Fig8.pdf){width="1.0\columnwidth"}

![image](Fig9.pdf){width="2.0\columnwidth"}

The SPH artificial viscosity parameters that appear to produce a minimum excess dissipation for this problem are $({\alpha_{\rm SPH}}, {\beta_{\rm SPH}}) \approx (0.1, 2.0)$. However, given that the previous simulations did not use these optimum values [@Rice_beta_condition; @Meru_Bate_resolution], it is important to correct for this. We therefore carry out a number of SPH simulations using 31,250, 250,000 and 2 million particles to determine what the critical cooling timescale is for the discs simulated using (${\alpha_{\rm SPH}}$, ${\beta_{\rm SPH}}$) = (0.1, 2.0). Table \[tab:alpha\_beta\_inves\] (bottom panel) and Figure \[fig:factor1.5\] summarise the results of these simulations. It can immediately be seen that the critical cooling timescale is higher than the equivalent simulations with ${\beta_{\rm SPH}}= 0.2$ (also see Figure \[fig:SPH\_betaSPH0.2\_2\_comparison\]). Figure  \[fig:disc\_betaSPH\] (left panel) shows an image of a fragmented disc modelled using 2 million SPH particles, (${\alpha_{\rm SPH}}$, ${\beta_{\rm SPH}}$) = (0.1, 0.2) and a cooling timescale of $\beta = 9$ which fails to fragment. Figure \[fig:disc\_betaSPH\] (right panel) shows the equivalent disc modelled using (${\alpha_{\rm SPH}}$, ${\beta_{\rm SPH}}$) = (0.1, 2.0) which fragments even though it is modelled using a slower cooling time of $\beta = 12$. The results (Figures \[fig:factor1.5\] and \[fig:SPH\_betaSPH0.2\_2\_comparison\]) still show that as the resolution increases, the critical cooling timescale increases, consistent with the results presented with a lower value of ${\beta_{\rm SPH}}$. However, since there are only three data points with ${\beta_{\rm SPH}}= 2.0$, it is firstly not clear whether convergence exists and secondly, what function should be used to fit this data. If we assume a functional form as given by equation \[eq:fit\], we find that ${\beta_{\rm crit}}= 36.6 \pm 6.9$, assuming a first-order convergence rate (i.e. $\sigma = 1.0$) as indicated in Section \[sec:nonconverge\]. Using the same convergence rate as found in Section \[sec:nonconverge\] (i.e. $\sigma = 1.08$) we find that ${\beta_{\rm crit}}= 29.2 \pm 2.4$. This implies that the true value may well be as high as $\approx 30$. We note that a fit assuming a second-order convergence rate gives a poor fit to the data. In Section \[sec:nonconverge\], omitting the 16 million particle data point makes no significant difference to the value of ${\beta_{\rm crit}}$ obtained due to the slow convergence rate. Therefore we do not expect the absence of the 16 million particle data point here to significantly affect the fit.

Furthermore, we attempt to fit the analytical formula from equation \[eq:betacrit\_hH\]. Allowing all three parameters to vary ($\alpha_{\rm GI, crit}$, $\eta$ and $\zeta$) we find that $\alpha_{\rm GI, crit} = 0.015$, $\eta = 14.0$ and $\zeta = 0.3$. There are several points to note here. Firstly, this value of the critical gravitational stress is equivalent to a critical cooling timescale in the limit of infinite resolution of ${\beta_{\rm crit}}\approx 27$. Secondly, the values of $\eta$ and $\zeta$ are much less than those obtained in Section \[sec:nonconverge\]. This implies that when increasing the quadratic artificial viscosity term to ${\beta_{\rm SPH}}= 2.0$, not only does the total dissipation decrease but the level of excess dissipation also decreases. However, we express caution here: the value of $\zeta$ obtained here is lower than unity which is not possible. This is likely to be an artefact of using three data points to fit three unknowns. We therefore refit the data by setting $\zeta$ to the minimum possible value it can be, i.e. unity. In this case, we find that $\alpha_{\rm GI, crit} = 0.024 \pm 0.005$ and $\eta = 6 \pm 3$ (see dotted line in Figure \[fig:factor1.5\]). Again, we emphasise that $\eta$ is smaller than previously obtained in Section \[sec:nonconverge\] suggesting that with ${\beta_{\rm SPH}}= 2.0$, the excess dissipation is significantly reduced.

We note that the analytical formula in equation \[eqnumerical\] is dependent on two aspects. Firstly it assumes that the ratio of the smoothing length to the disc scaleheight, $h/H$, is given by equation \[eq:h\_H\]. Secondly, it assumes that the dissipation due to the artificial viscosity is indeed given by equation \[eq:av\_dissipation\]. To check the analytical arguments presented in Section \[sec:SPH\_av\], it is important to test these two aspects.

### Testing the analytical formula for $h/H$ {#sec:h_H}

![image](Fig10a.pdf){width="0.66\columnwidth"} ![image](Fig10b.pdf){width="0.66\columnwidth"} ![image](Fig10c.pdf){width="0.66\columnwidth"}

Figure \[fig:h\_H\_profile\] shows the analytical estimate of the ratio of the smoothing length to the disc scaleheight, $h/H$, against the azimuthally averaged radial profile of $h/H$ for the non-fragmenting discs (i.e. marginally stable discs where $Q \approx 1$) carried out in Section \[sec:frag\_bdry\_optimumSPH\] (bottom section of Table \[tab:alpha\_beta\_inves\]). It is very clear from this graph that at all resolutions considered, the analytical estimate of $h/H$ is a good approximation in the outer parts of the disc where the fragments generally form. It is important to note, that the analytical formula assumes the initial surface mass density profile remains constant. However, the discrepancy in the inner regions is due to the change in surface mass density profile as the disc evolves into a state of mechanical equilibrium on a viscous timescale. The change in surface mass density profile thus changes the value of $h/H$ (equation \[eq:h\_H\]). Since the viscous timescale is shorter at small radii, the disc evolves more rapidly there and thus the discrepancy is larger. However, it is important to note that for a cooling timescale close to the critical one, fragmentation occurs in the outer parts of the discs for surface mass density profiles shallower than $\Sigma \propto R^{-2}$ ([@Meru_Bate_fragmentation]) since the resolution increases with radius (equation \[eq:h\_H\]). The outer parts are where the agreement is best between the analytically expected and azimuthally averaged values of $h/H$. We also note that the agreement between the analytical formula for $h/H$ and the simulation data is better with increasing resolution as the viscosity decreases (equations \[eqss\] and \[eqssquad\]) and so the effective viscous time is larger resulting in a slower evolution of the surface mass density profile. We therefore conclude that the reason why the value of $\eta$ in Section \[sec:frag\_bdry\_optimumSPH\] is not unity cannot therefore be put down to a mismatch between the analytical and actual values of $h/H$ in the region where fragmentation will occur.

### Testing the analytical formula for ${\beta_{\rm crit}}$ using non-self-gravitating discs {#sec:diss_nonGI}

![image](Fig11a.pdf){width="1.0\columnwidth"} ![image](Fig11b.pdf){width="1.0\columnwidth"}

![Setting up an initial disc model with purely Keplerian shear flow, we measure the instantaneous viscous dissipation averaged over $\approx 800$ particles in a thin radial extent. We then add increasing random particle velocities to the disc setup and measure the dissipation. In this figure, we plot the measured values of the dissipation divided by the analytic values expected for purely Keplerian flow due to the $\alpha_{\rm SPH}$ (solid line) and $\beta_{\rm SPH}$ (dashed line) terms separately versus the magnitude of the random velocity dispersion in units of the sound speed (i.e. we plot $\eta$ (solid line) and $\zeta$ (dashed line) as defined in equation \[eqnumerical\]). The excess dissipation due to small-scale particle velocity dispersion can be substantially higher than that produced by a pure shear flow.[]{data-label="fig:dissipation"}](Fig12.pdf){width="1.0\columnwidth"}

We carry out a number of simulations of discs with 250,000 particles with the same physical parameters as described in Section \[sec:sim\] but *without* self-gravity. Table \[tab:non-GI\] shows a summary of these simulations. The goal of this exercise is to see what effect the change in the SPH artificial viscosity parameters has on the dissipation in a *laminar* disc *which should only be due to shear* and whether this is as we would expect from the analytical formulae. To do this, we must start from exactly the same disc. However, the effects of the initial conditions must also be removed as this may affect the amount of dissipation. Therefore, we run a disc using a cooling time, $\beta = 20$, for 1.5 ORPs using the artificial viscosity parameters $({\alpha_{\rm SPH}}, {\beta_{\rm SPH}}) = (0.1, 2.0)$ (i.e. the values that we expect would minimise the additional dissipation). This is equivalent to $\approx 3$ orbital periods at 15 au (where this analysis is done). The cooling time (equation \[eq:beta\]) is also $\approx 3$ orbital periods at 15 au. Since the initial evolution time and the cooling time are approximately equal, the disc is then settled such that the heating matches the cooling. We then change the artificial viscosity parameters in the disc according to what is shown in Table \[tab:non-GI\] and run the simulations for a short period of time (0.1 ORPs or $\approx 0.2$ orbital periods at 15 au) and measure the total dissipation rate due to the artificial viscosity in the radial range $14.9 \le R \le 15.1$ au and compare these with the expected dissipation due to the artificial viscosity using equation \[eq:orig\_dissipation\] (using the actual values of the sound speed and smoothing length obtained from the simulation rather than the initial values). It is important to note that in order to make this comparison, we only calculate the dissipation over a short period of time as we are comparing the instantaneous expected dissipation rate with the instantaneous actual (azimuthally averaged) dissipation rate. If we allow the discs to run for a very long time before measuring the dissipation rate, the discs will evolve considerably and a like-for-like comparison is then not possible. This subsequent evolution takes place over a much smaller timescale than the orbital timescale. Therefore, there is no time for $h/H$ or the velocity field of the particles to change. Therefore any change in the disc’s dissipation *must* be due to the change in artificial viscosity parameters.

Figure \[fig:diss\_alphaSPH\_betaSPH\] (left panel) shows a graph of how the dissipation rate (measured and expected) changes with ${\alpha_{\rm SPH}}$ (using a fixed ${\beta_{\rm SPH}}= 2.0$). The expected dissipation rate is the sum of the dissipation due to the ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$ terms. It can immediately be seen that at low values of ${\alpha_{\rm SPH}}$, the expected dissipation due to the ${\beta_{\rm SPH}}$ term is very important (though its contribution is often thought to be negligible in comparison to the ${\alpha_{\rm SPH}}$ term). At high values of ${\alpha_{\rm SPH}}$ the measured dissipation matches the expected values very well. Most strikingly, the total dissipation is *higher* than the expected dissipation from the analytical formula at low values of ${\alpha_{\rm SPH}}$. The discrepancy is about a factor of two for ${\alpha_{\rm SPH}}\lesssim 0.1$.

Figure \[fig:diss\_alphaSPH\_betaSPH\] (right panel) shows the measured and expected dissipation rates against ${\beta_{\rm SPH}}$ (using a fixed ${\alpha_{\rm SPH}}= 0.1$). In this case the total dissipation is always approximately a factor of 2 higher than the expected values in these non-self-gravitating calculations. However, we note that there is no obvious additional dissipation at small values of ${\beta_{\rm SPH}}$ compared to large values, in contrast to the self-gravitating calculations in Figure \[fig:S-shaped\] which shows a definite difference in results between low and high ${\beta_{\rm SPH}}$ values. The ${\beta_{\rm SPH}}$ viscosity was originally introduced into SPH to stop particle interpenetration at shocks in supersonic flows. Shocks are not present in the non-self-gravitating calculations, so particle interpenetration is not an issue, but shocks play a significant role in the self-gravitating calculations. Thus, the apparent reduction of the dissipation in the self-gravitating calculations when the value of ${\beta_{\rm SPH}}$ is increased is likely to be because particle interpenetration is stopped more effectively with a higher value of ${\beta_{\rm SPH}}$. Regardless, both panels in Figure \[fig:diss\_alphaSPH\_betaSPH\] clearly show that the expected contribution to the dissipation from the quadratic artificial viscosity term can be larger than that from the linear term. In particular, for the simulation using $({\alpha_{\rm SPH}}, {\beta_{\rm SPH}}) = (0.1, 0.2)$, which were the values used by [@Rice_beta_condition] and [@Meru_Bate_resolution] as well as many others, the dissipation is more than three times larger than the analytically expected value from the ${\alpha_{\rm SPH}}$ viscosity alone.

The level of dissipation expected from the SPH artificial viscosity given in Appendix \[appendixD\] is lower than that measured from the actual simulations. What is the source of the excess dissipation that we find? The key is that the derivation assumes that the only contribution from the artificial viscosity to the thermal energy is due to shear flow in a purely Keplerian disc. Any other motions will add to this dissipation. In a gravitationally unstable disc, we also expect heating from the bulk component of the artificial viscosity due to shocks generated in the disc. Indeed, this is the assumed source of heating that is supposed to allow a gravitationally unstable disc to achieve a quasi-steady state when an imposed cooling timescale, $\beta$, is applied. However, it is exactly this maximum heating rate that the disc can provide without fragmenting that we are trying to measure when we try to determine $\beta_{\rm crit}$. The fact that the convergence rate of $\beta_{\rm crit}$ with increasing resolution is slow (first order) and that increasing $\beta$ from 0.2 to 2.0 increases the critical cooling timescale significantly implies that there is a third source of heating. Furthermore, as demonstrated earlier in this section, even when we measure the dissipation in a calculation without self-gravity, we still find some excess dissipation beyond what Appendix \[appendixD\] predicts. Figure A3 of [@Lodato_Rice_original] shows a calculation of the Reynolds stress in a non-self-gravitating disc with $({\alpha_{\rm SPH}}, {\beta_{\rm SPH}}) = (0.1, 0.2)$. They find that the $\alpha_{\rm SS}$ parameter due to this is a few $\times 10^{-3}$ in the range $0 \le R \le 25$ au. As a check to ensure we are consistent with previous results, we calculate the Reynolds stress in the same way as [@Lodato_Rice_original] (though we only average over 0.1ORPs) and also find the $\alpha_{\rm SS}$ value to be a few $\times 10^{-3}$ over the same radial range. Furthermore, our results are also consistent with [@Forgan_alpha] who also find their $\alpha_{\rm SS}$ parameter to be a few $\times 10^{-3}$ (in the inner parts of their self-gravitating discs where the effects of self-gravity are not very important; see their Figure 4).

In a non-self-gravitating calculation there are essentially only two possible contributions to the viscous heating. The first is from the Keplerian shear flow. The second is any additional particle motions. Without self-gravity, these can only come from ‘random’ particle motions. It is well known that in a typical SPH simulation the particles ‘jostle’ one another, resulting in a particle velocity dispersion. This velocity dispersion results from errors in the pressure gradients due to the finite number of particles within a smoothing kernel. In a compressible SPH simulation, such motions are typically at the level of some fraction of the sound speed. In order to determine the effect of these motions on the dissipation in a disc we perform a simple toy experiment whereby for illustrative purposes, we introduce different amoounts of particle velocity dispersion to see its effects on the dissipation in the disc. Setting up a purely Keplerian disc, we compute the instantaneous average values of the dissipation for particles in a small radial extent, due to the $\alpha_{\rm SPH}$ and $\beta_{\rm SPH}$ terms separately. Comparing these values to those expected from the ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$ terms respectively (Appendix \[appendixD\]), i.e.

$$D_{\alpha} = \frac{93}{700} \eta {\alpha_{\rm SPH}}c_s h \Omega^2
\label{eq:Dalpha}$$

and

$$D_{\beta} = \frac{81}{280} \zeta {\beta_{\rm SPH}}h^2 \Omega^3,
\label{eq:Dbeta}$$

respectively. We find that as expected, $\eta=\zeta=1$ to a high level of precision ($\approx 1-3$ per cent when averaging over $\approx 800$ particles). We then experiment with adding different levels of random velocities in addition to the underlying Keplerian motion. The results are displayed in Figure \[fig:dissipation\], where the magnitude of the particle velocity dispersion is given as a fraction of the local sound speed in the disc. We see that if random motions at the level of, e.g. 30 per cent of the sound speed are present, the dissipation increases by factors of $\eta=1.7$ and $\zeta=2.2$. This provides us with an explanation for the excess dissipation in the non-self-gravitating calculations. When the level of artificial viscosity is low, the velocity dispersion of the particles in the disc generates a non-negligible fraction of the dissipation. When the viscosity is high (in particular the linear $\alpha_{\rm SPH}$ term), this velocity dispersion is damped, and since the contribution to the dissipation from the shear flow is larger, no significant dissipation beyond that expected from the shear flow is found. We can see from Figure \[fig:diss\_alphaSPH\_betaSPH\] that $\alpha_{\rm SPH}=0.1$ is too low to effectively damp the particle velocity dispersion, resulting in dissipation rates that are approximately a factor of two larger than expected. This does not require a high level of particle velocity dispersion – from Figure \[fig:dissipation\] we see that a velocity dispersion of only $\approx 25$ per cent of the sound speed is enough to boost the dissipation due to the $\beta_{\rm SPH}$ viscosity by a factor of two. Thus, although we show in Sections \[sec:betaSPH\_test\] and \[sec:alphaSPH\_test\] that the minimum dissipation is obtained for $(\alpha_{\rm SPH}, \beta_{\rm SPH}) = (0.1, 2.0)$, this minimum dissipation is still larger than that expected from the analytic derivation.

In the self-gravitating disc calculations the situation is more complex. Here there are gravitational forces from the fluid, shocks and local pressure gradients in the disc which can stir up the particles. A low level of artificial viscosity (particularly $\beta_{\rm SPH}$) will allow particle penetration in shocks and a low value of $\alpha_{\rm SPH}$ will be ineffective at damping post-shock oscillations and other small-scale particle motions. If the random motions become a substantial fraction of the sound speed, which they may well do since we find the Mach numbers across a shock to be up to $\approx 3$, the factors can become very large ($\eta=3-9$ and $\zeta=7-24$ for random velocities of $60-100$ per cent of the sound speed). Thus, if random particle motions are indeed also playing a part in the self-gravitating calculations, it should be no surprise that we infer a level of dissipation that is well beyond that expected from equation \[eq:orig\_dissipation\], i.e. the counter-intuitive nature of ${\beta_{\rm SPH}}$ that leads to Figure \[fig:S-shaped\]. This leaves us with a problem with the SPH simulations. In order to obtain a level of dissipation that is close to that predicted by a purely Keplerian flow we can infer from Figure \[fig:diss\_alphaSPH\_betaSPH\] that we would need to use $\alpha_{\rm SPH} {\mathrel{\hbox{\rlap{\lower.55ex \hbox {$\sim$}}
                   \kern-.3em \raise.4ex \hbox{$>$}}}}1$ and from Figure \[fig:S-shaped\] that we require $\beta_{\rm SPH} \ge 2$. This should cut the particle penetration at shocks, post shock oscillations, and other particle velocity dispersion to low levels, thus making the analytic predictions of the dissipation accurate. On the other hand, the higher level of viscosity would increase the dissipation generated by the shear flow, thus reducing the measured value of $\beta_{\rm crit}$ at a given resolution (c.f. Figure \[fig:beta\_alphaSPH\]). Thus, although the simulations may be better behaved, an even higher numerical resolution would be needed to determine the converged value of the critical cooling timescale, $\beta_{\rm crit}$. We discuss other options in Section \[sec:disc\].

In summary, in a purely Keplerian disc, with no random particle velocity dispersion, the dissipation is as expected from the analytical values in equation \[eq:orig\_dissipation\]. However, in the simulations of non-self-gravitating discs, particularly with low values of ${\alpha_{\rm SPH}}$, the dissipation is somewhat higher than expected. We attribute this to random particle velocity dispersion, since there is no other source of heating in such discs over and above the viscous heating due to Keplerian shear flow. In self-gravitating discs additional particle dispersion will be present, which may well result in the counter-intuitive nature of the artificial viscosity that leads to Figure \[fig:S-shaped\], i.e. more dissipation with lower ${\beta_{\rm SPH}}$.

The effect of the [fargo]{} artificial viscosity on the critical cooling timescale {#sec:FARGO_av0.5}
----------------------------------------------------------------------------------

![Graph of $\beta$ against the [fargo]{} artificial viscosity parameter, $q$, of the non-fragmenting (open squares) and fragmenting (solid triangles) simulations carried out using 512 and 1536 cells in the radial and azimuthal directions, respectively. The solid line, included by eye, shows a dividing line between the fragmenting and non-fragmenting cases and the grey region is where fragmentation can take place. For low values of $q$, the dissipation due to artificial viscosity is low (and is most likely dominated by the intrinsic numerical diffusion) resulting in fragmentation occurring with high values of $\beta$. As the artificial viscosity is increased, a faster cooling is required to overcome the additional dissipation, resulting in lower values of ${\beta_{\rm crit}}$.[]{data-label="fig:beta_q_fargo"}](Fig13.pdf){width="1.0\columnwidth"}

In Section \[sec:FARGO\_av\] we show that the dissipation due to artificial viscosity present in [fargo]{} may play a part in the critical cooling timescale. Table \[tab:fargo\_av\] and Figure \[fig:beta\_q\_fargo\] summarise the results of the simulations carried out to investigate this. As the artificial viscosity parameter, $q$, is increased, it becomes harder for the disc to fragment due to the extra heating. Consequently, the critical value of $\beta$ required to overcome this and allow the disc to fragment decreases.

At lower values of $q$, the effect of artificial viscosity on the fragmentation boundary is much less obvious. Note from Table \[tab:fargo\_av\] that when the artificial viscosity parameter is set to zero, the fragmentation boundary decreases to a lower value of $\beta$ as with the SPH results in Figure \[fig:beta\_alphaSPH\] but the effect of reducing the viscosity is much less pronounced in [fargo]{} than in SPH. However, the reasoning is likely to be different because there is no dissipation associated with the artificial viscosity term since it is set to zero. We note from Figure \[fig:beta\_q\_fargo\] that ${\beta_{\rm crit}}$ increases rapidly as $q$ is decreased to $q \approx 0.5$ and then plateaus - this is most likely because at such low values of the artificial viscosity the dissipation is dominated by intrinsic dissipation in the code.

For fragmentation to occur, the dissipation associated with high values of the artificial viscosity parameter needs to be overcome with a faster cooling. As with the SPH results presented in Sections \[sec:betaSPH\_test\] and \[sec:alphaSPH\_test\], artificial viscosity clearly plays a part in whether these discs, modelled using a grid-based code, fragment or not. Figure \[fig:beta\_q\_fargo\] suggests that a value of $q \approx 0.5$ may be sufficient to avoid any excess dissipation.

### Determining the fragmentation boundary using the optimum value of the [fargo]{} artificial viscosity parameter

![Graph of $\beta$ against resolution of the non-fragmenting (open squares) and fragmenting (solid triangles) [fargo]{} simulations. These simulations are carried out with an artificial viscosity parameter, $q = 0.5$. The solid line, obtained by fitting equation \[eq:fit\], shows a dividing line between the fragmenting and non-fragmenting cases and the grey region is where fragmentation can take place. The graph shows clear evidence of convergence of results with increased resolution. The convergence rate is second order with spatial resolution. The dotted line (which coincides well with the solid line) is obtained by fitting equation \[eq:betacrit\_fargo\].[]{data-label="fig:res_beta_fargo_q0.5"}](Fig14.pdf){width="1.0\columnwidth"}

![Graph of ${\beta_{\rm crit}}$ against resolution of the [fargo]{} simulations carried out with an artificial viscosity parameter, $q = 1.41$ (squares) and $q = 0.5$ (triangles). It can be seen that the effect of reducing the artificial viscosity parameter to $q = 0.5$ (i.e. to a value that minimises the additional dissipation) is to increase the critical cooling timescale, with the effect being much greater at lower resolution.[]{data-label="fig:fargo_q1.41_q0.5_comparison"}](Fig15.pdf){width="1.0\columnwidth"}

In Section \[sec:FARGO\] we show that convergence appears to be reached at higher resolution with [fargo]{} and that the convergence is second order in spatial resolution. However, these simulations use a value of the artificial viscosity parameter, $q = 1.41$, which we show in Section \[sec:FARGO\_av0.5\] does not minimise the dissipation. We carry out simulations of self-gravitating discs using a value of $q = 0.5$ at various different resolutions. Table \[tab:fargo\_sim\_q0.5\] and Figure \[fig:res\_beta\_fargo\_q0.5\] summarise the simulations carried out to investigate this, and the key fragmenting results. It can be seen that the effect of using a lower value of $q$ is that ${\beta_{\rm crit}}$ is higher than obtained in Section \[sec:FARGO\]. Figure \[fig:disc\_fargo\_q\] shows a surface mass density rendered image of one of the highest resolution discs (modelled using 50 million grid cells) and shows clear fragmentation with a cooling time as high as $\beta = 26$. However, despite the critical cooling time being larger, we can see from Figure \[fig:res\_beta\_fargo\_q0.5\] that convergence is still being achieved. We firstly fit the data using equation \[eq:fit\] and find that ${\beta_{\rm crit}}= 28.0 \pm 0.2$ and $\sigma = 1.89 \pm 0.05$ showing that the convergence rate is second order with spatial resolution. We then fit the data using equation \[eq:betacrit\_fargo\] and find that $\alpha_{\rm GI, crit} = 0.0145 \pm 0.0001$ and $\xi = 3.16 \pm 0.04$. In the limit of infinite resolution, this value of $\alpha_{\rm GI, crit}$ is equivalent to a critical cooling timescale, ${\beta_{\rm crit}}\approx 28$ (using equation \[eq:stress\]).

Figure \[fig:fargo\_q1.41\_q0.5\_comparison\] shows the fragmentation boundary (with error bars) using $q = 1.41$ (as in Section \[sec:FARGO\]) and $q = 0.5$. As the resolution increases, the difference between the two sets of results decreases: since the convergence with [fargo]{} is fast, i.e. second-order, the effect of using different values of the artificial viscosity parameter (i.e. $q = 0.5$ versus $q = 1.41$) becomes negligible with 50 million grid cells compared to a lower resolution. This further corroborates that at a higher resolution, the artificial viscosity plays less of a part in the fragmentation boundary. As suggested by the analytics in Section \[sec:FARGO\_av\] it is expected that the artificial viscous dissipation should decrease both when the resolution is increased and when $q$ is decreased (until the numerical dissipation becomes dominated by intrinsic grid dissipation). Indeed, the higher value of $\xi$ obtained here in comparison to that in Section \[sec:FARGO\] suggests that the dissipation due to the artificial viscosity has been minimised and that the intrinsic grid dissipation is becoming more important, consistent with Figure \[fig:beta\_q\_fargo\]. We note that the value $q = 0.5$ effectively means that the shock is spread over approximately half a grid cell which will affect the treatment of shocks. We emphasise that we choose this value since it gives the lowest artificial heating rate that is possible with [fargo]{}, as done so with the SPH simulations. More importantly, we show that the choice of the value of $q$ has much less of an effect at higher resolution than at low resolution.

Discussion {#sec:disc}
==========

The non-convergence of results concerning the fragmentation of self-gravitating discs has opened up a number of questions concerning both the physics and the numerics involved in determining whether a disc will fragment into bound objects. Consequently, [@Meru_Bate_resolution] presented the possibility that either the critical cooling timescale was larger than originally thought, or the extreme possibility that the physics behind the fragmentation of discs needs to be reconsidered. In this paper we find that using both SPH and the [fargo]{} codes i.e. a three-dimensional particle-based Lagrangian code and a two-dimensional grid-based Eulerian code, respectively, the artificial viscosity that is used to accurately model shocks plays a significant part in the convergence rate. Not only do we show that convergence can occur, but we also show that the rate at which it occurs is as expected from analytical arguments involving artificial viscosity (first-order with linear resolution for SPH and second order for [fargo]{} i.e. a faster convergence with [fargo]{}). This affects the results on the fragmentation boundary in both SPH and grid-based calculations.

In particular, we conclude that oscillations at the shock front and particle interpenetration may not have been adequately accounted for in the previous SPH simulations (via the quadratic artificial viscosity term, ${\beta_{\rm SPH}}$). Those SPH simulations that used a value of ${\beta_{\rm SPH}}$ that was too low, may have counterintuitively resulted in more dissipation, causing fragmentation to have been *underestimated*. After minimising the additional dissipation associated with the artificial viscosity employed in both codes, we find that the critical cooling timescale is at least as high as 20 and perhaps even as much as $\approx 30$, for a ratio of specific heats, $\gamma = 5/3.$

Previous simulations that investigate the effects of artificial viscosity show that the strength of the gravitational instabilities are weakened [@Pickett_thermal_AV] and clump formation is reduced [@Boss_AV] when artificial viscosity is employed. These results are in the same sense as our results, i.e. excess effective viscosity reduces the propensity for fragmentation.

[@Mayer_AV] carried out a test on the fragmentation and disc evolution resulting from the inclusion of artificial viscosity in their three-dimensional SPH code. [@Pickett_Durisen_AV] carried out a similar test using a three-dimensional grid-based code. Both sets of authors perform their tests on isothermal simulations i.e. they only include the effects of artificial viscosity in the momentum equation and did not consider its heating effects in the energy equation. They suggest that artificial viscosity may reduce or even prevent clump formation from occurring. While our results are consistent with theirs with respect to preventing clump formation when artificial viscosity is increased, we stress that the dissipation associated with the artificial viscosity plays a key role in the fragmentation results.

Another possible numerical parameter that may affect the results is the gravitational softening used in the two-dimensional grid code. [@Muller_smoothing] show that an incorrect value of the gravitational softening length in two-dimensional disc simulations can significantly affect the fragmentation conclusions: a low value causes the gravitational forces on short distances to be over-estimated, resulting in the conclusion that fragmentation does occur, when the converse conclusion is reached for larger values of the softening parameter. Indeed, it is well known that three-dimensional discs are more stable than two-dimensional discs since the vertical component dilutes the effect of gravity [@Toomre_stability1964]. Thus, incorrectly taking into account the effects of the vertical direction in a two-dimensional simulation may cause the disc to be more prone to fragmentation than its three-dimensional equivalent. [@Muller_smoothing] show that a value of $\approx 0.6 H$ is required to model the gravitational forces correctly (though they do say that a comparison with 3D simulations is required). Since we are using a softening length of $3 \times 10^{-4} H$, our [fargo]{} simulations may overestimate fragmentation.

For the SPH simulations, although we see evidence for convergence of the critical cooling timescale, the convergence rate is only first order with increasing resolution. This is partly due to the larger dissipation than that predicted by the continuum limit of the SPH equations in a shear flow. We argue that the excess dissipation is due to small-scale particle velocity dispersion. In non-self-gravitating discs, this results from pressure gradient errors due to discretisation, but in the self-gravitating discs there are other potential sources: primarily particle penetration at shock fronts and post-shock oscillations (particularly when the levels of artificial viscosity are low), but perhaps also discretisation errors in the self-gravity. To achieve well-behaved dissipation (i.e. that which is close to that predicted by the continuum limit of the SPH equations) requires $\alpha_{\rm SPH}{\mathrel{\hbox{\rlap{\lower.55ex \hbox {$\sim$}}
                   \kern-.3em \raise.4ex \hbox{$>$}}}}1$ and $\beta_{\rm SPH}{\mathrel{\hbox{\rlap{\lower.55ex \hbox {$\sim$}}
                   \kern-.3em \raise.4ex \hbox{$>$}}}}2$. However, the dissipation from the shear flow is then relatively high meaning that even higher resolution would be necessary to obtain a converged value of the critical cooling timescale, $\beta_{\rm crit}$.

However, there are many possibilities that might improve the SPH performance. We have employed the most basic form of SPH artificial viscosity in this paper (i.e. constant values of $\alpha_{\rm SPH}$ and $\beta_{\rm SPH}$). An obvious aspect to investigate is whether a viscosity switch such as those proposed by [@MorrisMonaghan1997] or [@CullenDehnen2010] which increase $\alpha_{\rm SPH}$ and $\beta_{\rm SPH}$ in the presence of a shock and allow them to decay away from shocks can improve the convergence rate. This could potentially provide high viscosity to avoid particle penetration at shocks and post-shock oscillations (i.e. reducing small-scale particle velocity dispersion), but retain low viscosity (and thus low heating rates) in the bulk of the disc, thus minimising the heating due to the shear flow. Since some of the particle velocity dispersion originates from pressure gradient errors, another possibility is to try a more accurate kernel. The quintic spline kernel generally performs better than the cubic spline kernel and [@MorrisFoxZhu1997] reported that it significantly reduced velocity field noise in their calculations. Testing these variants is beyond the scope of this paper, but we expect that these and other SPH variations may be able to significantly improve the performance of SPH on this problem. We stress that these possible improvements to the SPH convergence rate will not decrease the value of the critical cooling timescale and thus the values obtained in this paper indicate lower limits for ${\beta_{\rm crit}}$.

Since [@Meru_Bate_resolution] published their results highlighting the convergence problem, a number of authors have attempted to explain the non-convergence. [@Lodato_Clarke_resolution] speculated that the cause may be the artificial smoothing of the density enhancements in SPH or a larger than expected level of artificial viscosity. Our results clearly show that artificial viscosity plays a major role in numerical determinations of the critical cooling rate.

[@Paardekooper_convergence] suggested that the boundary between the turbulent inner disc region and the laminar outer disc region (a natural consequence of starting with smooth initial conditions) may cause an edge in the disc that becomes more and more pronounced at higher resolutions, making it easier to fragment. They suggested that if the smooth initial conditions were removed, convergence could be achieved. In light of the new results presented in this paper, the effect of edges should be considered in more detail. If edge effects do play a part, it is unclear whether they should always continue to become sharper at higher resolution (and hence inconsistent with the results presented here), or whether they should “saturate” at higher resolution (and thus consistent with these results). It is important to note, however, that [@Bate_disc_frag_res] performed radiative transfer calculations of molecular cloud collapses and found that disc fragmentation is more prevalent in higher resolution calculations. These discs did not begin with smooth initial conditions, and yet a similar resolution dependence was seen.

More recently, [@Paardekooper_shearing_sheet] carried out shearing sheet simulations, similar to those performed by [@Gammie_betacool] where no such edge effects should play a part. He found that as the resolution was increased, the critical cooling timescale also increased, showing that the convergence issue is not restricted to global simulations, but also affects local simulations. He found fragmentation for at least as large as ${\beta_{\rm crit}}\approx 20$. However, this was for simulations carried out with a ratio of specific heats, $\gamma = 2$. This is equivalent to a maximum gravitational stress as least as small as $\alpha_{\rm GI, crit} \approx 0.011$, consistent with the value of the gravitational stress found using our global simulations.

It is important to note that many of the previous simulations that have attempted to explain the convergence problem highlighted by [@Meru_Bate_resolution] have tried to do so by carrying out simulations with resolutions in the *non-convergent* region of the resolution space shown here in this paper. It is therefore hard to interpret those results since they may have been affected by artificial viscosity. For those codes that do not use artificial viscosity, other sources of numerical diffusion relating to the specific implementation may be important and the effect of these need to be thoroughly understood. It would be interesting to try to understand the convergence problem with a Godunov scheme that does not implement an artificial viscosity, or to apply a fixed Navier-Stokes viscosity.

While this work focusses on the effects on fragmentation of self-gravitating discs, the key conclusion that artificial viscosity can play a significant role in the physical interpretation of simulations is more general. We emphasise that any simulations whose outcome is highly dependent on the thermodynamics of a problem should ensure that the effects of artificial viscosity in their code are well understood as well as highlighting the importance of convergence of results with both resolution and numerical method.

Implications for the fragmentation of real discs
------------------------------------------------

A critical cooling timescale of ${\beta_{\rm crit}}\approx 20$ or $\approx 30$ is equivalent to a maximum gravitational stress of $\alpha_{\rm GI, crit} \approx 0.02$ or $\approx 0.013$, respectively. [@Clarke2009_analytical] produced an analytical model for the structure of a gravitationally unstable disc which is subject to realistic cooling. She showed that for optically thick discs that are sufficiently low in temperature that they are dominated by ice grains,

$$\alpha_{\rm GI} = 0.4~\bigg(\frac{R}{100~\rm{au}}\bigg)^{\frac{9}{2}},$$

for a disc with interstellar opacities and surrounding a $1 {{\rm M_{\odot}}}$ star, where $R$ is the radius being considered. This relationship shows that for a maximum value of the gravitational stress, a critical radius, $R_{\rm crit}$, can be found outside of which fragmentation can occur (for a disc with a shallow surface mass density profile). While the previously accepted result of $\alpha_{\rm GI,max} \approx 0.06$ gives a critical radius of $R_{\rm crit} \approx 68$ au, the values of ${\beta_{\rm crit}}$ obtained here moves the critical radius inwards to $R_{\rm crit} \approx 47-51$ au (for a disc around a $1 {{\rm M_{\odot}}}$ star using interstellar opacities). The core accretion scenario is thought to occur out to $\approx 10$ au, while gravitational instability is historically thought to operate outside of $\approx 70-120$ au [@Rafikov_SI; @Clarke2009_analytical]. Therefore, an intermediate radial region exists where no one in situ formation method adequately seems to describe the formation of planets. Our results show that this gap can at least partly be bridged if the true critical cooling timescale is as much as ${\beta_{\rm crit}}\approx 20-30$.

We point out that equation \[eq:stress\] is derived by assuming that the dominant form of heating in a disc is that due to the gravitational instabilities. In a real disc, there may be a contribution to the stress from the magnetorotational instability (MRI), $\alpha_{\rm MRI}$. In this case, equation \[eq:stress\] may be written as

$$\beta = \frac{4}{9} \frac{1}{\gamma (\gamma -1)} \frac{1}{(\alpha_{\rm GI} + \alpha_{\rm MRI})}.$$

Therefore, while the critical cooling timescale for a purely gravitationally unstable disc is quite large, if the contribution to the gravitational stress from the MRI (or in fact other heating sources) becomes important, a faster cooling will be required to overcome this additional heating and allow the disc to fragment. Since we find that the critical stress may be as low as $\alpha_{\rm GI, crit} \approx 0.01$, the heating due to MRI will certainly be expected to be important if it provides an effective stress of approximately this level or higher. Even if $\alpha_{\rm MRI}$ is a factor of 10 smaller, this will still make a 10% difference to the heating which can change the critical cooling timescale required for fragmentation.

Conclusions {#sec:conc}
===========

We perform hydrodynamical simulations using a three-dimensional Smoothed Particle Hydrodynamics code and a two-dimensional Eulerian grid-based code of self-gravitating discs to investigate how the presence of artificial viscosity may affect fragmentation results. We present additional SPH results to those presented by [@Meru_Bate_resolution] as well as perform similar simulations using the grid-based hydrodynamics code, [fargo]{}, and show that convergence with resolution of the critical cooling timescale can be achieved with both codes. We show that the previous non-convergent results are largely due to the effects of artificial viscosity that play a more prominent role at lower resolution. We find that the convergence rate of the critical cooling timescale required for fragmentation is first order in spatial resolution using SPH and second order using [fargo]{}. Furthermore, we find that the dissipation from the artificial viscosities in SPH is exactly as we would expect in a purely laminar disc. However, if random particle motions are present, they can produce dissipation due to artificial viscosity that is larger than expected. In self-gravitating discs, using a value of the quadratic artificial viscosity term that is too low can result in counterintuitively high dissipation. This may be caused by additional random particle velocity dispersion due to the presence of shocks, causing the dissipation to significantly deviate from that expected from the SPH continuum limit equations. In addition, the dissipation due to the ${\beta_{\rm SPH}}$ term may not be as small as previously assumed and should not be ignored.

We show using analytical arguments and numerical simulations that as the resolution is increased, the artificial viscosity term becomes less important and the rate of convergence is as expected from the analytical arguments. With the particular setup adopted here, which is the same as that used by [@Rice_beta_condition] and [@Meru_Bate_resolution], we find that once the effects of artificial viscosity have been minimised, the critical cooling timescale converges with increasing resolution to a value at least as high as ${\beta_{\rm crit}}\approx 20$ and perhaps even as high as ${\beta_{\rm crit}}\approx 30$. However, a convergence between the two codes has not yet been achieved. We conclude that this is much more of a problem than had previously been supposed, in part due to the slow convergence rate of SPH and in part due to the enormous resolution required to obtain convergence. The critical cooling timescale is a factor of $\approx 3-5$ times larger than the value of ${\beta_{\rm crit}}\approx 6$ that has been used in the past, and is equivalent to a maximum gravitational stress of $\alpha_{\rm GI, crit} \approx 0.013-0.02$ that a disc can handle before it fragments (in contrast to the previously obtained value of $\alpha_{\rm GI, crit} \approx 0.06$).

We show that using values of the artificial viscosity that do not minimise the additional dissipation caused by it can significantly affect fragmentation results and we expect that any other results that sensitively depend on the thermodynamics of a problem, e.g. collapse of AGN discs and molecular clouds into stars may also be affected. This highlights the importance of ensuring that the artificial viscosity does not play a significant role when carrying out numerical simulations. We show that fragmentation of self-gravitating discs can be suppressed if the effects of artificial viscosity are not carefully considered. This suggests that fragmentation of discs into bound objects (e.g. for the formation of planets, binary companions and stars formation in galaxy simulations) is easier than previously thought.

Acknowledgments {#acknowledgments .unnumbered}
===============

We thank the referee, Cathie Clarke, for her thorough review and insightful comments. We thank Daniel Price, Jim Pringle, Neal Turner, Giuseppe Lodato, Ken Rice, Clément Baruteau, Sijme-Jan Paardekooper and Tom Quinn for helpful discussions. Most of the calculations reported here were performed using the University of Exeter’s SGI Altix ICE 8200 supercomputer. We gratefully thank the bwGRiD project for some of the computational resources. bwGRiD (http://www.bw-grid.de) is a member of the German D-Grid initiative, funded by the Ministry for Education and Research (Bundesministerium fuer Bildung und Forschung) and the Ministry for Science, Research and Arts Baden-Wuerttemberg (Ministerium fuer Wissenschaft, Forschung und Kunst Baden-Wuerttemberg). Some of the calculations reported here were performed using the [brutus]{} cluster at ETH Zürich. Some of the figures were produced using the publicly available [splash]{} visualisation software [@SPLASH]. MRB is grateful for the support of a EURYI Award which also funded FM. This work, conducted as part of the award “The formation of stars and planets: Radiation hydrodynamical and magnetohydrodynamical simulations" made under the European Heads of Research Councils and European Science Foundation EURYI (European Young Investigator) Awards scheme, was supported by funds from the Participating Organisations of EURYI and the EC Sixth Framework Programme. FM acknowledges the support of the German Research Foundation (DFG) through grant KL 650/8-2 within the Collaborative Research Group FOR 759: [*The formation of Planets: The Critical First Growth Phase*]{}. FM was also supported by the ETH Zurich Postdoctoral Fellowship Programme as well as by the Marie Curie Actions for People COFUND program.

Analytic Derivation of the Shear Viscosity Present in SPH {#appendixA}
=========================================================

In the following sections, for the sake of clarity, we make the simplifying assumptions that the SPH particle smoothing length $h$, the sound speed, $c_{\rm s}$, and the density, $\rho$, are all slowly varying (i.e. are constant).

The original SPH artificial viscosity {#sec:original_AV}
-------------------------------------

The standard SPH artificial viscosity method described by [@Monaghan_Gingold_art_vis], which is a time-independent fixed artificial viscosity, adds the following term to the momentum equation:

$$\frac{{\rm d} {\bf v}_i}{{\rm d} t} = - \sum_{j} {m_j \Pi_{ij} {\bf \nabla}_i W_{ij}}
\label{eq:mom_av}$$

where $$\Pi_{ij} =
\left\{
\begin{array}{l l}
\displaystyle    \frac{- \alpha_{\rm SPH}~c_{{\rm s}}~\mu_{ij}~+~\beta_{\rm SPH}~\mu_{ij}^2}{\rho} & {\bf v}_{ij} \cdot {\bf r}_{ij} < 0 \\[1.25ex]
   0 & {\bf v}_{ij} \cdot {\bf r}_{ij} > 0,
\end{array}
\right.
\label{eq:Pi}$$ $$\mu_{ij} = \frac{h {\bf v}_{ij} \cdot {\bf r}_{ij}}{{\bf r}_{ij}^2 + \eta^2},
\label{eq:mu}$$ $m_j$ is the mass of particle $j$, $W_{ij}$ is the smoothing kernel adopted, $h$ is the smoothing length and ${\bf v}_{ij}={\bf v}_i - {\bf v}_j$ is the velocity difference between particles $i$ and $j$. The quantity $\eta^2=0.01h^2$ is included to avoid divergence for small separations between neighbouring particles. The viscosity involves two terms, the strengths of which are controlled by the parameters ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$. Note that the artificial viscosity is only applied when particles approach each other and is turned off when they recede from each other.

It has been recognised for some time that in the continuum limit, the $\alpha_{\rm SPH}$ viscosity term applies both a bulk and a shear viscosity which has the form of a Navier-Stokes type viscosity [@Monaghan1985; @Pongracic1988; @Meglicki_av_continuum]. In particular, [@Meglicki_av_continuum] provide a clear derivation showing that the viscous acceleration due to the $\alpha_{\rm SPH}$ viscosity is given by $$\label{eq:viscous_force}
\frac{{\rm d}{\bf v}}{{\rm d}t} = \frac{\alpha_{\rm SPH} h \kappa}{2 \rho} \left[ \nabla \cdot (c_s \rho {\bf S})  + \nabla (c_s \rho \nabla \cdot {\bf v}) \right],$$ where $$S_{ij} = \frac{\partial v^i}{\partial x^j} +  \frac{\partial v^j}{\partial x^i},$$ is the deformation tensor. The first term in equation \[eq:viscous\_force\] is a shear viscosity, while the second term is a bulk viscosity. The constant, $\kappa$, depends on the number of spatial dimensions and the kernel used by the SPH code. If the $\alpha_{\rm SPH}$ viscosity is applied in three dimensions and between both approaching and receding particles (unlike equation \[eq:Pi\]) then (see Appendix \[sec:evaluating\]) $$\kappa = - \frac{4 \pi}{15} \int r^3 \frac{{\rm d}W}{{\rm d}r}~{\rm d}r.$$ The value of this integral depends on the kernel that is being used. We use the standard cubic spline kernel $$W(q,h) = \frac{\sigma}{h^d}
\left\{
\begin{array}{l l}
   1-\frac{3}{2}q^{2}+\frac{3}{4}q^{3} & {\rm for}  \ 0 \leq q <  1, \\[1.25ex]
   \frac{1}{4}(2-q)^{3} & {\rm for} \ 1 \leq q < 2, \\[1.25ex]
   0 & {\rm otherwise},
\end{array}
\right.
\label{eq:kernel_1}$$ where $d$ is the number of dimensions, $\sigma$ is the normalisation constant equal to 2/3, $10/(7\pi)$ and $1/\pi$ in one, two and three dimensions respectively, and $q = r/h$. In this case, the integral has the value $-3/(4\pi)$, such that $\kappa=1/5$. We find that in general the shear viscosity contribution to the momentum equation in two and three dimensions can be written $$\frac{{\rm d}{\bf v}}{{\rm d}t}  = \nu \nabla \cdot {\bf S} = \frac{1}{2(2+d)} \alpha_{\rm SPH} c_s h \nabla \cdot {\bf S},$$ where $\nu$ is the kinematic shear viscosity. Thus, for example, when simulating accretion discs, [@Artymowicz_Lubow_av] used $\nu = \frac{1}{8} \alpha_{\rm SPH} c_s h$ for their two dimensional SPH simulations, while [@Lodato_Price_betaSPH] used $\nu = \frac{1}{10} \alpha_{\rm SPH} c_s h$ for their three dimensional SPH simulations.

However, as expressed in equation \[eq:Pi\], it is usual in SPH simulations to only apply the artificial viscous force between approaching particles. Thus, in general the kinematic shear viscosity in SPH simulations is a factor of two smaller and, for three dimensional calculations, is $$\nu = \frac{1}{20} \alpha_{\rm SPH} c_s h.
\label{eq:orig_nu}$$

A more recent variation of SPH artificial viscosity {#sec:recent_AV}
---------------------------------------------------

Recently, a slightly different form of artificial viscosity has been applied in SPH codes [@Chow_Monaghan_newAV; @Price_Monaghan_v_dot_r_AV] where, essentially, equation \[eq:mu\] is replaced by $$\mu_{ij} = {\bf v}_{ij} \cdot {\bf \hat{r}}_{ij}.
\label{eq:mu_2}$$ This is the form of the artificial viscosity that we use and is constructed in analogy to dissipative terms in Riemann methods. This form also avoids the arbitrary quantity $\eta^2$ which was included to avoid numerical divergences. It also means that the magnitude of $\mu_{ij}$ differs by a factor of $r_{ij}/h$ from the original form of the viscosity. Note that different implementations can differ slightly in the way that the coefficients $\alpha_{\rm SPH}$ and $\beta_{\rm SPH}$ enter the equations. For example, in some implementations the coefficient of the quadratic viscosity term is actually given by the product of $\alpha$ and $\beta$, while in others the coefficient is fixed to be $2\alpha$. Some implementations may also differ in the value of $\Pi_{ij}$ by factors of two, so care needs to be taken when evaluating the continuum limit of the viscosity in any particular SPH code.

It can be shown (see Appendix \[appendixB\]), that when using equations \[eq:Pi\] and \[eq:mu\_2\] and the standard cubic spline kernel in three dimensions that the shear viscosity is actually 18% larger than for the original artificial viscosity such that $$\nu = \frac{31}{525} \alpha_{\rm SPH} c_s h.
\label{eq:new_nu}$$

The $\beta_{\rm SPH}$ viscosity
-------------------------------

Although many past studies have considered the continuum limit of the linear $\alpha_{\rm SPH}$ viscosity in SPH, to our knowledge, nobody has considered the contribution of the quadratic $\beta_{\rm SPH}$ viscosity. By inspection of equations \[eq:Pi\] and \[eq:mu\] we can determine how the shear viscosities due to both the $\alpha_{\rm SPH}$ and $\beta_{\rm SPH}$ should scale. For a pure shear flow, $\mu_{ij}$ provides an estimate of the shear rate of the fluid multiplied by $h$. Thus, since the kinematic viscosity is the ratio of the shear stress (given by $\rho \Pi_{\rm ij}$) to the shear rate, we expect the kinematic viscosity due to the $\alpha_{\rm SPH}$ term to scale as $$\nu_{\alpha} = q_{\alpha} {\alpha_{\rm SPH}}c_{\rm s} h
\label{eq:nu_alpha}$$ where $q_{\alpha}$ is a constant of proportionality. This is consistent with the analysis given above. An alternative method of arriving at this equation is to use the fact that from kinetic theory, kinematic viscosity is proportional to the characteristic speed of interchange between particles and to the characteristic distance over which interchange occurs. In this case, the characteristic speed of particle interchange is given by $c_s$ and the distance over which particles interact is a smoothing length, $h$.

The [@Vonneumann_Richtmyer_av] type viscosity (SPH $\beta$-viscosity) is a second order viscosity with the viscous forces depending on the square of the relative speed of the particles. Therefore, although the characteristic distance over which the viscosity acts is still a smoothing length, the characteristic speed is now the relative speed between particles over a smoothing length. In a flow directed in the $x$-direction which is sheared in the $y$-direction this is given by

$$h \left( \frac{{\rm d}v_x}{{\rm d}y}\right)$$

and thus the $\beta$-viscosity is expected to scale as

$$\nu_{\beta} = q_{\beta} {\beta_{\rm SPH}}h^2 \left( {{\rm d}v_x \over {\rm d}y} \right)
\label{eq:nu_beta}$$

where $q_{\beta}$ is a constant of proportionality.

SPH artificial viscosity in the continuum limit {#appendixB}
===============================================

One way to evaluate the constants of proportionality, $q_{\alpha}$ and $q_{\beta}$, is to use the defining equation for kinematic viscosity. This can be obtained by considering the shearing force produced by a viscous fluid on a plane running parallel to the direction of motion of the fluid. If the fluid flows in the $x$-direction and there is a velocity gradient across the flow in the $y$-direction, then the kinematic viscosity of the fluid is defined by

$${F \over A} = \nu \rho {{\rm d}v_x \over {\rm d}y}
\label{eq:FoverA}$$

which gives the force per unit area exerted on the plane surface as a function of the kinematic viscosity, $\nu$, the density of the fluid, $\rho$, and the shear rate of the fluid.

The force exerted on a volume element of fluid is determined by considering two planes parallel to the flow. If the two planes are parallel to the $x$-$z$ plane and they are separated by a distance $\delta y$, then the net force on the fluid element is

$$F_1 - F_2 = \left[{\nu_1 \rho_1 \left({{\rm d}v_x \over {\rm d}y}\right)_1 - \nu_2 \rho_2 \left({{\rm d}v_x \over {\rm d}y}\right)_2 \over \delta y}\right]  \delta x \delta y \delta z$$

so that the force per unit volume is given by

$${F \over V} = \nu \rho {{\rm d}^2v_x \over {\rm d}y^2}$$

in the limit that $\delta z \rightarrow 0$, assuming the density and kinematic viscosity are constant, which is a valid assumption to make if the two planes are sufficiently close to each other. The force per unit volume for a fluid element is simply the acceleration of the fluid element multiplied by its density. Thus, the equation

$${{\rm d}v_x \over {\rm d}t} = \nu {{\rm d}^2v_x \over {\rm d}y^2}$$

is obtained, which can be compared directly to the SPH momentum equation \[eq:mom\_av\]. Using equations \[eq:nu\_alpha\] and \[eq:nu\_beta\] we therefore find that the specific force can be expressed as

$$\frac{{\rm d}v_x}{{\rm d}t} \bigg |_{{\alpha_{\rm SPH}}} = q_{\alpha} {\alpha_{\rm SPH}}c_{\rm s} h \bigg ( \frac{{\rm d}^2 v_x}{{\rm d}y^2} \bigg )
\label{eq:mom_alphaSPH}$$

and

$$\frac{{\rm d}v_x}{{\rm d}t} \bigg |_{{\beta_{\rm SPH}}} = q_{\beta} {\beta_{\rm SPH}}h^2 \bigg ( \frac{{\rm d}v_x}{{\rm d}y} \bigg ) \bigg ( \frac{{\rm d}^2 v_x}{{\rm d}y^2} \bigg ),
\label{eq:mom_betaSPH}$$

for the ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$ terms, respectively.

Evaluating the constant of proportionality, $q_{\alpha}$, with the original form of the SPH artificial viscosity {#sec:evaluating}
----------------------------------------------------------------------------------------------------------------

To evaluate the constant of proportionality for the $\alpha_{\rm SPH}$ term in the artificial viscosity, we can derive the continuum limit and compare it to equation \[eq:mom\_alphaSPH\] to determine the magnitude of $q_\alpha$. Using equations \[eq:mom\_av\] - \[eq:mu\], the force per unit mass due to artificial viscosity on particle $i$ is

$$\begin{aligned}
\frac{{\rm d}{{\bf v}_i}}{{\rm d}t} &= -
\sum_{j} \frac{m_j}{\rho} \left [ - {\alpha_{\rm SPH}}h c_{\rm s} 
\frac{{\bf v}_{ij} \cdot {\bf r}_{ij}}{({\bf r}_{ij}^2 + \eta^2)} + {\beta_{\rm SPH}}h^2 \frac{({\bf v}_{ij} \cdot {\bf r}_{ij})^2}{({\bf r}_{ij}^2 + \eta^2)^2} \right ] \nabla_i W_{ij} \nonumber \\
&\approx \int_{\rm Kernel} \left [{\alpha_{\rm SPH}}h c_{\rm s} \frac{{\bf v}_{ij} \cdot {\bf r}_{ij}}{{\bf r}_{ij}^2} - {\beta_{\rm SPH}}h^2 \frac{({\bf v}_{ij} \cdot {\bf r}_{ij})^2}{({\bf r}_{ij}^2)^2} \right ] \nabla_i W_{ij} {\rm d}^3x.
\label{eq:F_full}\end{aligned}$$

The artificial viscosity in the continuum limit due to the linear term (${\alpha_{\rm SPH}}$) is, $$\frac{{\rm d}{{\bf v}_i}}{{\rm d}t} \approx \int_{\rm Kernel} {\alpha_{\rm SPH}}h c_{\rm s}  \frac{{\bf v}_{ij} \cdot {\bf r}_{ij}}{({\bf r}_{ij}^2 + \eta^2)} \nabla_i W_{ij} {\rm d}^3x.
\label{eq:F_alpha}$$ Following Appendix A of [@Meglicki_av_continuum], we expand ${\bf v}_j$ around ${\bf r}_i$ to give $${\bf v}_{ij} \cdot {\bf r}_{ij} = \left [ \Delta x^p \frac{\partial v_i}{\partial x^p} + \frac{\Delta x^p \Delta x^q}{2} \frac{\partial^2 v_i}{\partial x^p \partial x^q} + \frac{\Delta x^p \Delta x^q \Delta x^a}{6} \frac{\partial^3 v_i}{\partial x^p \partial x^q \partial x^a} \right ] \Delta x^r .$$ In addition, $$\nabla_i W_{ij} = -\frac{{\rm d}W_{ij}}{{\rm d}r} \frac{\Delta x^k}{r}
\label{eq:grad_W}$$ where $r = |{\bf r}_{ij}| = |r_i - r_j|$. Inserting equation \[eq:grad\_W\] into the first part of equation \[eq:F\_alpha\] and retaining terms of order $O(h^4)$, which are the lowest non-vanishing terms, gives $$\frac{{\rm d}{v^k_i}}{{\rm d}t} \approx - \int_{\rm Kernel} \frac{{\alpha_{\rm SPH}}h c_s}{r^3} \left[ \frac{\partial v_i^r}{\partial x^p} + \frac{1}{2} \Delta x^q \frac{\partial^2 v_i^r}{\partial x^p \partial x^q} \right] \Delta x^p \Delta x^r \Delta x^k \frac{{\rm d}W_{ij}}{{\rm d}r} {\rm d}^3 x + O(h^6).
\label{eq:Falpha_full}$$ We note that integrating a term with $(\Delta x)^t$, where $t$ is odd, over a symmetric kernel yields a zero result. Therefore, the result of the integration of the first term in equation \[eq:Falpha\_full\] is zero. Simplifying equation \[eq:Falpha\_full\] gives $$\frac{{\rm d}{v^k_i}}{{\rm d}t} \approx - \frac{{\alpha_{\rm SPH}}h c_{\rm s}}{2}  \frac{\partial^2 v_i^r}{\partial x^p \partial x^q}   \int_{\rm Kernel} \Delta x^p \Delta x^r \Delta x^k \Delta x^q \frac{1}{r^3} \frac{{\rm d}W_{ij}}{{\rm d}r} {\rm d}^3 x.
\label{eq:Falpha_simplified}$$ The integral is a fourth order symmetric isotropic tensor which can be written in the form $$- \int_{\rm Kernel} \Delta x^p \Delta x^r \Delta x^k \Delta x^q \frac{1}{r^3} \frac{{\rm d}W_{ij}}{{\rm d}r} {\rm d}^3 x = \kappa_{\alpha} ( \delta_{pq} \delta_{rk} + \delta_{pr} \delta_{qk} + \delta_{pk} \delta_{rq} )
\label{eq:fourth_order}$$ where $$\kappa_{\alpha} = A \int_{\rm Kernel} r^3 \frac{{\rm d}W_{ij}}{{\rm d}r} ~{\rm d}r.
\label{eq:kappa}$$ Equation \[eq:Falpha\_simplified\] can therefore be written $$\frac{{\rm d}{v^k_i}}{{\rm d}t} \approx \frac{\kappa_\alpha {\alpha_{\rm SPH}}h c_{\rm s}}{2} \frac{\partial^2 v_i^r}{\partial x^p \partial x^q}  ( \delta_{pq} \delta_{rk} + \delta_{pr} \delta_{qk} + \delta_{pk} \delta_{rq} ).$$ Contracting with the delta terms yields $$\frac{{\rm d}{v^k_i}}{{\rm d}t}  \approx \frac{\kappa_\alpha {\alpha_{\rm SPH}}h c_{\rm s}}{2} \left [
\frac{\partial^2 v_i^k}{\partial x^p \partial x^p} +
\frac{\partial^2 v_i^p}{\partial x^p \partial x^k} +
\frac{\partial^2 v_i^r}{\partial x^k \partial x^r}  \right ]
\label{eq:contraction}$$ If we now assume that locally there is a constant flow in one direction (e.g. the $x$-direction) with a velocity gradient in an orthogonal direction (e.g. the $y$-direction) so that $\partial v_i^k / \partial x^k = 0$ and $\partial / \partial x^k (\partial v_i^k / \partial x^p) = 0$, as we would expect for a shear flow, then the above equation can be simplified to give $$\frac{{\rm d}{v^k_i}}{{\rm d}t}  \approx \frac{\kappa_\alpha {\alpha_{\rm SPH}}c_{\rm s} h}{2} \left (
\frac{\partial^2 v_i^k}{\partial x^p \partial x^p}
\right )
\label{eq:Falpha_continuum}$$ Comparing this to equation \[eq:mom\_alphaSPH\] yields $\kappa_\alpha / 2 = q_{\alpha}$.

Equation \[eq:fourth\_order\] defines $\kappa_\alpha$. To calculate $\kappa_\alpha$ in three dimensions, we need to sum over all possible combinations of $r$, $p$ and $q$ in equation \[eq:fourth\_order\]. The simplest case that yields a non-zero value of the right hand side of equation \[eq:fourth\_order\] involves $k = r = p = q$ such that equation \[eq:fourth\_order\] gives

$$- \int_{\rm Kernel} (\Delta x^k)^4 \frac{1}{r^3} \frac{{\rm d}W_{ij}}{{\rm d}r} {\rm d}^3 x = 3 \kappa_{\alpha}
\label{eq:kappa_alpha1}.$$

There are three additional cases where $k$ is equivalent to one other letter while the remaining two letters are equal but in an orthogonal direction to $k$. Since in three dimensions there are two orthogonal directions to $k$, equation \[eq:fourth\_order\] gives $$- 6 \int_{\rm Kernel} (\Delta x^k)^2 (\Delta x^q)^2 \frac{1}{r^3} \frac{{\rm d}W_{ij}}{{\rm d}r} {\rm d}^3 x = 6 \kappa_{\alpha}
\label{eq:kappa_alpha2}.$$ Summing equations \[eq:kappa\_alpha1\] and \[eq:kappa\_alpha2\] together yields $$- \int_{\rm Kernel} \left [ (\Delta x^k)^4 + 6 (\Delta x^k)^2 (\Delta x^q)^2 ) \right ] \frac{1}{r^3} \frac{{\rm d}W_{ij}}{{\rm d}r} {\rm d}^3 x = 9 \kappa_{\alpha}.
\label{eq:kappa_alpha_sum}$$ Without loss of generality, we use $x = r \sin \theta \cos \phi = \Delta x^k$, $y = r \sin \theta \sin \phi = \Delta x^q$, and $d^3x = r^2 \sin \theta ~{\rm d}r ~{\rm d}\theta ~{\rm d}\phi$ and substitute into equation \[eq:kappa\_alpha\_sum\]. For the $\theta$-component we integrate over $\theta = [0, \pi]$. However, for the $\phi$-component, care must be taken to integrate over the correct range since we only consider particles that are approaching each other. Figure \[fig:schematic\] shows that in the frame of the particle being considered the ranges $\phi = [0, \pi/2]$ and $\phi = [\pi, (3\pi)/2]$ involve particles approaching each other whereas outside this range particles recede from each other. Integrating yields

$$\kappa_{\alpha} = - \frac{2 \pi}{15} \int_{\rm Kernel} r^3 \frac{{\rm d}W}{{\rm d}r}~{\rm d}r.
\label{eq:kappa_alpha_full}$$

Using the standard cubic spline kernel given above (for three dimensions) $$\int r^3 \frac{{\rm d}W(q,h)}{{\rm d} r} = -\frac{3}{4 \pi}.$$ Therefore, substituting into equation \[eq:kappa\_alpha\_full\] yields $\kappa_\alpha = 1/10$, and so equation \[eq:contraction\] becomes $$\frac{{\rm d}v_x}{{\rm d}t} \bigg |_{{\alpha_{\rm SPH}}} \approx  \frac{1}{20} {\alpha_{\rm SPH}}c_{\rm s} h \frac{{\rm d}^2 v_x}{{\rm d}y^2}.$$ We note that the same constant is achieved if the integration in the $\phi$-direction is done over all space and simply divided by two.

Evaluating the constant of proportionality, $q_{\alpha}$, with the recent form of the SPH artificial viscosity {#sec:new_alpha}
--------------------------------------------------------------------------------------------------------------

For the more recent form of SPH artificial viscosity (equation \[eq:mu\_2\]), the procedure is identical, but the expansion and, thus, the integral is slightly different. The required expansion is $${\bf v}_{ij} \cdot {\bf \hat{r}}_{ij} =  \frac{1}{r} \left [ \Delta x^p \frac{\partial v_i}{\partial x^p} + \frac{\Delta x^p \Delta x^q}{2} \frac{\partial^2 v_i}{\partial x^p \partial x^q} + \frac{\Delta x^p \Delta x^q \Delta x^a}{6} \frac{\partial^3 v_i}{\partial x^p \partial x^q \partial x^a} \right ] \Delta x^r .$$ so that equation \[eq:Falpha\_simplified\] becomes $$\frac{{\rm d}{v^k_i}}{{\rm d}t} \approx - \frac{{\alpha_{\rm SPH}}c_{\rm s}}{2}  \frac{\partial^2 v_i^r}{\partial x^p \partial x^q}   \int_{\rm Kernel} \Delta x^p \Delta x^r \Delta x^k \Delta x^q \frac{1}{r^2} \frac{{\rm d}W_{ij}}{{\rm d}r} {\rm d}^3 x.
\label{eq:Falpha_simplified_2}$$ Note that $h$ is missing from this equation and the integral differs by a factor of $r$. Thus, the equivalent of equation \[eq:kappa\_alpha\_full\] that needs to be solved is $$\kappa_{\alpha} = A \int_{\rm Kernel} r^4 \frac{{\rm d}W_{ij}}{{\rm d}r} ~{\rm d}r$$ where $A = (-2 \pi)/15$ (ensuring that we account for the fact that the viscosity is only applied between approaching particles), but this time the integral has a value of $-31 h/(35 \pi)$. Therefore, we obtain $$\frac{{\rm d}v_x}{{\rm d}t} \bigg |_{{\alpha_{\rm SPH}}} \approx \frac{31}{525} {\alpha_{\rm SPH}}c_{\rm s} h \frac{{\rm d}^2 v_x}{{\rm d}y^2},$$ which as noted above is approximately 18% larger than the value obtained for the original form of the artificial viscosity.

The dissipation associated with ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$ {#appendixD}
============================================================================

![Schematic diagram showing the fluid motion (denoted by solid arrows) surrounding an SPH particle (black dot) in its rest frame, as the gas disc orbits the central star. In the region of $\phi$-space where $\phi = [0, \pi/2]$ and $\phi = [\pi, (3\pi)/2]$ (shaded region) the fluid is approaching the SPH particle and in this region artificial viscosity is turned on while in all other areas the fluid is receding from the SPH particle and here the artificial viscosity is turned off (see equation \[eq:Pi\]). The origin is the location of the SPH particle.[]{data-label="fig:schematic"}](Fig17.pdf){width="0.3\columnwidth"}

For this present paper we are more interested in the magnitude of the thermal dissipation provided by the viscosity than the angular momentum transport as such. The dissipation rate per unit mass in a viscous accretion disc is given by $$\frac{{\rm d} u}{{\rm d} t}  = \frac{T_{R\phi}}{\Sigma} \left (R \frac{{\rm d}\Omega}{{\rm d} R} \right ) = \nu R \frac{{\rm d}\Omega}{{\rm d} R} \left ( R \frac{{\rm d}\Omega}{{\rm d} R} \right ) = \frac{9}{4} \nu \Omega^2,
\label{eq:Keplerian_dissipation}$$ where $T_{R \phi}$ is the stress tensor and the final equality assumes a Keplerian disc. In this case the shear rate is $$\frac{{\rm d}v_x}{{\rm d}y} = R \frac{{\rm d} \Omega}{{\rm d} R} = - \frac{3}{2} \Omega.$$

Rather than derive the dissipation rate from the continuum limit of the SPH momentum equation, we can derive the dissipation rate directly from the continuum limit of the SPH energy equation $$\frac{{\rm d} u_i}{{\rm d} t} = \frac{1}{2} \sum_{j} {m_j \Pi_{ij} {\bf v}_{ij} \cdot {\bf \nabla}_i W_{ij}}.$$ Taking the continuum limit of this equation and using the more recent form of the artificial viscosity gives $$\frac{{\rm d} u}{{\rm d}t} \approx \frac{1}{2} \int_{\rm Kernel} \left [ - {\alpha_{\rm SPH}}c_{\rm s} {\bf v} \cdot {\bf \hat{r}} + {\beta_{\rm SPH}}({\bf v} \cdot {\bf \hat{r}})^2 \right ] {\bf v} \cdot  \hat{\bf r}~ \frac{{\rm d}W}{{\rm d}r}~ {\rm d}^3x.
\label{eq:cont_diss}$$

The simplest way to obtain the dissipation rate due to the artificial viscosity in a Keplerian disc is to consider a small patch of the disc around a particle such the the local velocity field is given by $$v_x \approx - \frac{3 y \Omega_0}{2},$$ where $\Omega_0$ is the angular velocity of the fluid at the radius being considered and $y$ is the displacement in the inward radial direction (see Figure \[fig:schematic\]). Inserting this into equation \[eq:cont\_diss\] and taking $x=r\sin\theta\cos\phi$ and $y=r\sin\theta\sin\phi$, we obtain $$\frac{{\rm d} u}{{\rm d}t} \approx \frac{9}{4} \left [\frac{31}{525} \alpha_{\rm SPH} c_{\rm s} h \Omega_0^2 + \frac{9}{70 \pi} \beta_{\rm SPH} h^2  \Omega_0^3 \right],
\label{eq:orig_dissipation}$$ taking care to integrate only over the ranges $\phi=[0,\pi/2]$ and $\phi=[\pi,3\pi/2]$ (i.e. where the flow is approaching and not receding; see Figure \[fig:schematic\]).

For the original form of the artificial viscosity, the dissipation is $$\frac{{\rm d} u}{{\rm d}t} \approx \frac{9}{4} \left [\frac{1}{20} \alpha_{\rm SPH} c_{\rm s} h \Omega_0^2 + \frac{3}{35 \pi} \beta_{\rm SPH} h^2  \Omega_0^3 \right],
\label{eq:new_dissipation}$$ again only taking the integral over the regions where the flow is approaching.

Note that the coefficients preceding $\alpha_{\rm SPH}$ in equations \[eq:orig\_dissipation\] and \[eq:new\_dissipation\] are the same as the coefficients appearing in the the kinematic viscosity as given by equations \[eq:orig\_nu\] and \[eq:new\_nu\], respectively, as is expected from equation \[eq:Keplerian\_dissipation\].

Using $H = c_s/\Omega$, the ratio of the dissipation rates associated with the linear and quadratic artificial viscosity terms using the recent form of the SPH artificial viscosity is given by $$\frac{D_{\alpha}}{D_{\beta}} = \frac{62 \pi}{135} \frac{{\alpha_{\rm SPH}}}{{\beta_{\rm SPH}}} \frac{H}{h}.$$ Typical values of ${\alpha_{\rm SPH}}$ and ${\beta_{\rm SPH}}$ are frequently within a factor of 2 of each other (with ${\beta_{\rm SPH}}= 2 {\alpha_{\rm SPH}}$). Therefore, the important variable that determines the relative magnitude of the dissipation associated with the linear and quadratic SPH artificial viscosity terms is the ratio of the smoothing length to disc scaleheight, $h/H$. In a well resolved disc, $h/H \ll 1$ such that $D_{\alpha} \gg D_{\beta}$. However, if the disc is poorly resolved and/or $\beta_{\rm SPH} \gg \alpha_{\rm SPH}$ the dissipation associated with the quadratic viscosity may be significant. Finally, note that the above equations assume that the only viscous dissipation comes from the shear in a purely Keplerian disc.

[^1]: farzana.meru@phys.ethz.ch
