Stable anisotropic heat conduction in smoothed particle hydrodynamics

TL;DR

This paper presents a stable method for simulating anisotropic heat conduction in Smoothed Particle Hydrodynamics, emphasizing entropy increase and proposing a two-first-derivative approach for improved stability and larger timesteps.

Contribution

It introduces a novel stable SPH method using two first derivatives with alternating operators for anisotropic heat conduction, avoiding instability of second derivative methods.

Findings

Stable method requires entropy increase.

Two first derivatives with alternating operators are stable.

Larger timesteps are possible with the proposed method.

Abstract

We investigate how to simulate anisotropic heat conduction in a stable manner in Smoothed Particle Hydrodynamics. We show that the requirement for stability is that entropy must increase. From this, we deduce that methods involving direct second derivatives in SPH are unstable, as found by previous authors. We show that the only stable method is to use two first derivatives with alternating differenced and symmetric SPH derivative operators, with the caveat, that one may need to apply smoothing or use an artificial conductivity term if the initial temperature jump is discontinuous. Furthermore, we find that with two first derivatives the stable timestep can be 3--8 times larger even for isotropic diffusion.

Tables

Table 1. Table 1: Summary of stability conditions for anisotropic diffusion in SPH using different methods.

Method	Isotropic stability	Anisotropic stability
Brookshaw/Espanol-Revenga	Always	If Eq. 24 satisfied
Direct second derivatives	If Eq. 27 satisfied	If Eq. 26 satisfied
Two first derivatives	Always	Always

Table 2. Table 2: The integration constants used in the numerical tests.

	$M_4$	$M_6$	$M_8$
Brookshaw/Espanol-Revenga	0.15	0.18	0.2
Direct second derivatives	0.2	0.35	0.4
Two first derivatives	0.6	1.2	1.6

Equations

$$\frac{\mathrm{d}u}{\mathrm{d}t}=\frac{1}{\rho }\nabla \cdot (\mathbf{\kappa }\nabla T),$$

$$A({\mathbf{r}}_a)\equiv A_a=\sum _b\frac{m_b}{{\rho }_b}{A}_{b}W(|{\mathbf{r}}_a-{\mathbf{r}}_b|,h),$$

$${\nabla }_a{A}_a=\sum _b\frac{m_b}{{\rho }_b}\left(A_b-A_a\right){\nabla }_a{W}_{ab}.$$

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\sum _b\frac{m_b}{{\rho }_a{\rho }_b}{\overline{\kappa }}_{ab}(T_b-T_a)\frac{-2F_{ab}}{|r_{ab}|}.$$

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\sum _b\frac{m_b}{{\rho }_a{\rho }_b}{\overline{\kappa }}_{ab}(T_b-T_a){\nabla }_a^2{W}_{ab}.$$

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\sum _b\frac{m_b}{{\rho }_a{\rho }_b}(T_b-T_a)\left({\overline{\kappa }}_{ij}^{ab}{\nabla }_a^i{\nabla }_a^j{W}_{ab}\right),$$

$${\nabla }_a^i{\nabla }_a^j{W}_{ab}=\frac{1}{C_{\nu }{h}^{\nu +2}}\left[\left(f^{\prime \prime }-f^{\prime }{q}^{-1}\right){\hat{r}}_{ab}^i{\hat{r}}_{ab}^j+f^{\prime }{q}^{-1}{\delta }^{ij}\right].$$

$${\nabla }^2{Y}_{ab}\equiv \frac{-2F_{ab}}{|r_{ab}|},$$

$$F_{ab}=\frac{C_{\nu }}{h^{\nu +1}}{f}^{\prime }(q),$$

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=-\sum _b\frac{m_b}{{\rho }_a{\rho }_b}({\kappa }_a+{\kappa }_b)(T_b-T_a)\frac{F_{ab}}{|r_{ab}|}.$$

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=-\sum _b\frac{m_b}{{\rho }_a{\rho }_b}\frac{4{\kappa }_a{\kappa }_b}{({\kappa }_a+{\kappa }_b)}(T_b-T_a)\frac{F_{ab}}{|r_{ab}|}.$$

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=-\sum _b\frac{m_b}{{\rho }_a{\rho }_b}{T}_{ba}{\overline{\kappa }}_{ij}^{ab}\left[(\nu +2){\hat{r}}_{ab}^i{\hat{r}}_{ab}^j-{\delta }^{ij}\right]\frac{F_{ab}}{|r_{ab}|}.$$

$$F_{ab}^{ij}=-\left[(\nu +2){\hat{r}}_{ab}^i{\hat{r}}_{ab}^j-{\delta }^{ij}\right]\frac{F_{ab}}{|r_{ab}|},$$

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\sum _b\frac{m_b}{{\rho }_a{\rho }_b}{T}_{ba}{\overline{\kappa }}_{ij}^{ab}{F}_{ab}^{ij}.$$

$${\tilde{F}}_{ab}^{ij}=\gamma {\kappa }_{ij}^{ab}{F}_{ab}^{ij}+\frac{1-\gamma }{3}{\delta }^{ij}.$$

$$F^j$$

$$\frac{\mathrm{d}u}{\mathrm{d}t}$$

$$F_a^j$$

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}$$

$${\Omega }_a\equiv 1+\frac{3h_a}{{\rho }_a}\sum _b{m}_b\frac{\partial W_{ab}(h_a)}{\partial h_a}.$$

$$\frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\sum _a{m}_a{u}_a=\sum _a{m}_a\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=0.$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}=\sum _a{m}_a\frac{\mathrm{d}{s}_a}{\mathrm{d}t}=$$

$$\sum _a\sum _b\frac{m_a{m}_b}{{\rho }_a{\rho }_b}\left(\frac{1}{T_a}-\frac{1}{T_b}\right)(T_b-T_a){\overline{\kappa }}_{ij}^{ab}{F}_{ab}^{ij}.$$

$${\overline{\kappa }}_{ii}^{ab}{F}_{ab}^{ii}=\kappa \frac{-2F_{ab}}{|r_{ab}|}\ge 0.$$

$$-{\kappa }_{xx}[(\nu +2){\hat{r}}^x{\hat{r}}^x-1)]-{\kappa }_{yy}[(\nu +2){\hat{r}}^y{\hat{r}}^y-1]\ge 0.$$

$$\frac{{\kappa }_{xx}}{{\kappa }_{yy}}\in \left[\frac{1}{\nu +1};\nu +1\right].$$

$$-{\kappa }^{xx}[(f^{\prime \prime }-f^{\prime }{q}^{-1}){\hat{r}}^x{\hat{r}}^x+f^{\prime }{q}^{-1}]$$

$$-{\kappa }^{yy}[(f^{\prime \prime }-f^{\prime }{q}^{-1}){\hat{r}}^y{\hat{r}}^y+f^{\prime }{q}^{-1}]\ge 0,$$

$$-\kappa [(f^{\prime \prime }-f^{\prime }{q}^{-1})+\nu f^{\prime }{q}^{-1}]\ge 0.$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}$$

$$\sum _a{m}_a\frac{1}{T_a}\sum _b{m}_b\left[\frac{{\kappa }_{ij}^a{F}_a^i{\nabla }_a^j{W}_{ab}(h_a)}{{\Omega }_a{\rho }_a^2}+\frac{{\kappa }_{ij}^b{F}_b^i{\nabla }_a^j{W}_{ab}(h_b)}{{\Omega }_b{\rho }_b^2}\right].$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}$$

$$\Delta t\le \frac{C{c}_v\rho h^2}{{\kappa }_{ij}},$$

Full text

††pubyear: 2018††pagerange: Stable anisotropic heat conduction in smoothed particle hydrodynamics–Stable anisotropic heat conduction in smoothed particle hydrodynamics

1 Introduction

In this paper, we aim to find stable discretisations of thermal diffusion equations of the form

[TABLE]

where $u=c_{v}T$ is the internal energy per unit mass, $\rho$ is the density, $T$ is the temperature, and $\boldsymbol{\kappa}\equiv\kappa_{ij}$ is a heat conduction tensor (with $\kappa_{ij}\equiv\kappa\delta_{ij}$ when the conduction is isotropic). The heat conduction tensor must be a symmetric tensor ($\kappa_{ij}=\kappa_{ji}$) with positive diagonal elements, to ensure that heat flows from hot to cold and not vice versa.

Such equations arise in astrophysics when considering heat (Parrish & Stone, 2005) or cosmic ray diffusion (Forman & Gleeson, 1975; Girichidis et al., 2016) in the presence of magnetic fields, and the anisotropic diffusion of radiation (Petkova & Springel, 2009).

Smoothed particle hydrodynamics (SPH; for reviews see Monaghan 1992, 2005) is a Lagrangian method for fluid dynamics. The method consists of two main steps. First, discretise matter onto a set of finite particles. For each particle $a$ we define properties such as mass $m_{a}$, density $\rho_{a}$, and position $\mathbf{r}_{a}$. Second, discretise the set of equations using the summation interpolant (Gingold & Monaghan, 1977; Lucy, 1977)

[TABLE]

where $A_{a}$ is the interpolated property, and $W(|\mathbf{r}_{a}-\mathbf{r}_{b}|,h)\equiv W_{ab}$ is the kernel function describing the weight assigned to neighbours. As the kernel $W_{ab}$ is the only part of (2) that depends on $\mathbf{r}_{a}$, derivatives of $A_{a}$ with respect to the position $\mathbf{r}_{a}$ simply involve the gradient of the kernel function, e.g.,

[TABLE]

Discretisation of the isotropic heat equation in SPH was considered in Lucy (1977), but most modern implementations follow the method presented by Brookshaw (1985) and Cleary & Monaghan (1999) whereby

[TABLE]

In the above, $F_{ab}$ is the scalar part of the kernel gradient, and $\overline{\kappa}$ is either an arithmetic or harmonic mean of the (scalar) conductivity between the particle pair e.g. $\overline{\kappa}_{ab}=(\kappa_{a}+\kappa_{b})/2$.

Other methods have been proposed to reduce errors in second derivatives computed with SPH by correcting the gradient operators (e.g. Korzilius, Schilders & Anthonissen, 2017). We do not consider such methods in this paper because they do not guarantee the conservation of energy. By contrast, all of the schemes we examine conserve energy to the precision of the timestepping algorithm.

More recently, Español & Revenga (2003) generalised the Brookshaw method to the anisotropic case (i.e. when $\boldsymbol{\kappa}$ is a tensor). However, they found that achieving accuracy to within a few percent required at least 50 neighbours in 2D, which would be equivalent to 250 neighbours in 3D. This is costly.

Petkova & Springel (2009) found a more significant problem: The method proposed by Español & Revenga (2003) is unstable when the diffusion is highly anisotropic. While Petkova & Springel (2009) proposed an ‘anisotropy limiter’ to fix the stability problem, Hopkins (2017) showed that this could lead to incorrect results.

The paper is organised as follows: In Section 2 we outline methods for isotropic and anisotropic heat conduction in SPH, and assess their stability. This shows why the Espanol & Revenga method is unstable — it does not guarantee positive entropy. Table 1 summarises the stability conditions. We test our ideas in Section 3, and show how to remove oscillations in the solution obtained with two first derivatives when the initial conditions are discontinuous. We summarise and draw conclusions in Section 4 and 5.

2 Heat conduction in SPH

2.1 Methods

2.1.1 Direct second derivatives

By taking second derivatives of (2), an SPH discretisation of (1) in the case where $\kappa$ is isotropic is given by

[TABLE]

For the anisotropic case with a tensor heat conduction coefficient $\kappa_{ij}$ we can generalise this to

[TABLE]

where we assume that repeated $i$ or $j$ indices imply summation. The second derivative of the kernel can be written in terms of the dimensionless kernel function according to (e.g. Price, 2010)

[TABLE]

We refer to this method as the ‘direct second derivative’ operator for heat conduction.

2.1.2 Brookshaw method for isotropic heat conduction

Brookshaw (1985) showed that the direct second derivative operator (5) could be very inaccurate, explaining the errors reported by Lucy (1977). Instead, Brookshaw proposed to use the first derivative of the kernel function to compute the second derivative. The method is equivalent to defining a new kernel function $Y$, such that (e.g. Price, 2012)

[TABLE]

where $F_{ab}$, in terms of the dimensionless kernel function $f(q)$, is

[TABLE]

with $q\equiv|r_{ab}|/h$, $C_{\nu}$ is the kernel normalisation constant, $h$ is the smoothing length and $\nu$ is the number of dimensions. The operator for isotropic diffusion in SPH is therefore

[TABLE]

Cleary & Monaghan (1999) proposed an alternative form using the harmonic mean for the case where $\kappa$ is discontinuous, given by

[TABLE]

We do not consider the harmonic mean in this paper since Price & Laibe (2015) found that it could produce incorrect results in dust diffusion problems.

2.1.3 Espanol & Revenga: Anisotropic heat conduction

Español & Revenga (2003) generalised the method further to be applicable to anisotropic diffusion. In this, more general case, they showed that the correct expression is

[TABLE]

where $T_{ba}\equiv T_{b}-T_{a}$. Their method is equivalent to defining

[TABLE]

such that the anisotropic discretisation of (1) is given by

[TABLE]

2.1.4 Petkova & Springel: Anisotropy-limited diffusion

Petkova & Springel (2009) showed that the method proposed by Español & Revenga (2003) is unstable when the diffusion is highly anisotropic. To stabilise the method, Petkova & Springel (2009) proposed to modify the operator according to

[TABLE]

In order to obtain a stable solution in 3D, one should use at least $\gamma=\frac{2}{5}$. This converts the fully anisotropic problem to $\frac{2}{5}$ of anisotropic diffusion plus $\frac{3}{5}$ of isotropic diffusion. However limiting the anisotropy of the diffusion in this way produces incorrect results, as discussed by Hopkins (2017).

2.1.5 Two first derivatives

A simple alternative approach is to compute the diffusive flux explicitly before taking the gradient, i.e.

[TABLE]

We show in Section 2.2.4 that this method is stable as long as we choose the derivative operators carefully. In particular, one should implement this ‘two first derivatives’ method using a combination of differenced and symmetric derivative operators (Cummins & Rudman, 1999; Tricco & Price, 2012). Considering the general case where the smoothing length on each particle is different, we discretise this using (e.g. Price, 2012)

[TABLE]

In the above, $\Omega$ is the term used in the evaluation of the smoothing length Monaghan (2002); Springel & Hernquist (2002)

[TABLE]

For completeness, in Section 3.1.3 we tested the derivative operators in the reverse order — differenced after symmetric. Although we find no significant difference in the results (see Figures 4 and 5), the order above is necessary to conserve total energy, $E$, since

[TABLE]

The last sum is zero when inserting (19) because of the antisymmetry of the kernel gradient producing a double summation which is antisymmetric in the particle index $a$ and $b$ and therefore equal to zero. The above version is also efficient because it does not require adding a third loop over the particles between the density and force evaluations — since evaluating the flux does not require prior knowledge of the density.

2.2 Stability

As stated by Monaghan (2005), the main advantage of the Brookshaw method over a direct second derivative of the SPH kernel is the guarantee of increasing entropy. For heat conduction problems, this means that heat always flows from hot to cold particles (which guarantees stability but says nothing about accuracy). To show that this is true, we need to prove that the rate of change of the total entropy $S$ is positive. For the Espanol-Revenga operator, we have

[TABLE]

As the terms involving mass, density and temperature are all positive, positive entropy requires $\overline{\kappa}_{ij}^{ab}F_{ab}^{ij}$ to be positive definite.

2.2.1 Brookshaw: Isotropic heat conduction

When heat conduction is isotropic ($\boldsymbol{\kappa}=\kappa\delta^{ij}$), the Español & Revenga (2003) operator reduces to the Brookshaw (1985) operator. The stability condition becomes

[TABLE]

Because $\kappa$ is always positive and $F_{ab}$ is negative (assuming SPH kernels where the weight decreases monotonically with radius), the kernel part (23) and the total entropy (22) are always positive, and therefore the method is stable.

2.2.2 Espanol & Revenga: Anisotropic heat conduction

As the following idea is the same for more than one dimension, we will discuss the simplest case of 2D anisotropic heat conduction. Knowing that, from the energy argument, the heat conduction tensor is a real symmetric tensor, it is always possible to find a principal axis and to represent the heat conduction tensor as a diagonal (orthotropic) tensor using spectral decomposition with respect to the principal axis. This means that for any anisotropic heat conduction problem we can find a new local coordinate axis in which $\boldsymbol{\kappa}$ is a diagonal tensor. Thus, to have positive total entropy (22) in 2D, we need to satisfy the inequality

[TABLE]

To satisfy the inequality (24), the ratio of components should be

[TABLE]

As in practice this ratio can be arbitrary, anti-diffusive behaviour is possible. To illustrate this, the grey shading in Figure 1 shows the regions inside the kernel compact support radius from which neighbouring particles contribute anti-diffusive terms. The instability grows from these regions resulting in negative temperatures.

2.2.3 Direct second derivatives

While we found this approach to be more stable than Espanol-Revenga method both practically and theoretically, after the same analysis as in Section 23 we arrive at a similar condition for stability as for the Espanol-Revenga operator, namely

[TABLE]

or in the isotropic case

[TABLE]

This means that for direct second derivatives with standard kernels there is always a region in which anti-diffusive behaviour occurs. These regions are shown for different kernels in grey in Figure 1 for the test case described in Section 3.2 where $\kappa^{xx}=1$ and $\kappa^{yy}=0$.

2.2.4 Two first derivatives

We compute the entropy evolution using

[TABLE]

Taking into account that $\nabla_{a}W_{ab}=-\nabla_{b}W_{ba}$ we can rearrange the double summation to give

[TABLE]

Hence, as long as the flux is calculated using (18), entropy increase is guaranteed. Importantly, stability is guaranteed independent of the choice of SPH kernel. This result is similar to the proof given in Price et al. (2017) for physical viscosity or in Price (2012) for magnetohydrodynamics.

Table 1 summarises the stability conditions for each of the methods discussed above.

2.3 Timestep constraints

For diffusion problems, the timestep requirement is given by (e.g. Cleary & Monaghan 1999)

[TABLE]

where $C$ is a constant. The above is a local constraint. Since we used the same timestep for all particles, we take the minimum over all particles.

Table 2 gives the values of $C$ we found to be stable in our 2D tests. We obtained these values empirically, although they could, in principle, be determined analytically by a stability analysis. The constants for Brookshaw and direct second derivative methods apply only to isotropic diffusion, as the methods are unstable for anisotropic diffusion regardless of the choice of timestep, which we confirmed down to $C=0.01$.

The change of the stability constant between kernels (left to right in Table 2) is not surprising. This simply follows the standard deviation of the kernel, as discussed by Dehnen & Aly (2012). Even so this effect helps to mitigate the computational cost of smoother spline kernels.

The surprising aspect is the factor of 3–8 increase in timestep possible with the two first derivatives method (bottom row of Table 2). We confirmed that this also holds in 3D. We attribute this to the extra stability provided by the double convolution of the kernel gradient inherent in the two first derivatives method. This gives us compelling reason to use this method even for isotropic diffusion.

3 Numerical tests

3.1 Diffusion in a slab with constant heat conduction tensor

To test our ideas in practice, we considered diffusion in a 3D slab: $x\in[-1,1]$ with Dirichlet boundary condition, and $y,z\in[-4\Delta r,4\Delta r]$, where $\Delta r$ is the particle spacing, with periodic boundary conditions. We set up the problem using $64\times 8\times 8$ particles in 3D. Initially, $T_{L}=1$ for $x<0$, and $T_{R}=2$ for $x>0$. In the isotropic case the exact solution is one dimensional and can be approximated for some time (until diffusion hits the boundary) by the exact solution for one dimensional heat conduction in an infinite slab

[TABLE]

This solution holds for anisotropic diffusion when $\kappa^{xx}$ is the only non-zero component

[TABLE]

3.1.1 Diffusion in the direction of the heat gradient

Figure 2 shows the results of this test at a resolution of 64 particles along the x-direction, placed on a uniform lattice. We compare the two first derivatives method (left) to direct second derivatives (right), finding that both two first derivatives and direct second derivative methods produce stable results. Using noisier kernels (M4 instead of M6) or higher resolution leads to numerical instability with the direct second derivative.

At this resolution and kernel (M6), a ‘carbuncle mode’ appears in the solution obtained with the two first derivatives method caused by the initially discontinuous temperature profile, while the solution obtained with direct 2nd derivative is less noisier. This error decreases with resolution, and can be eliminated by smoothing the initial conditions (see Section 3.1.4).

3.1.2 Diffusion perpendicular to the heat gradient

When the only non-zero component of heat conduction tensor is $\kappa^{yy}$, but heat conduction $T(x,t)$ depends only on the $x$ component, there should be no diffusion at all.

Figure 3 shows the results of the diffusion test with conduction allowed only in the direction perpendicular to the heat gradient. In contrast to the previous case, two first derivatives give the correct solution (left), while the direct second derivatives method is unstable (right).

3.1.3 Does the order of derivative operators matter?

Comparing left and right panels in Figures 4 and 5 shows that the order of operators (Symmetric after Differential in the left column or Differential after Symmetric in the right columns) does not affect the results. However, total energy is only exactly conserved when the symmetric operator is used in the thermal energy equation.

3.1.4 Eliminating the carbuncle mode

The reason for the oscillations in Section 3.1.1 is the discontinuity in the initial conditions. One can adopt smooth initial conditions by setting the initial temperature according to

[TABLE]

where $L$ is the length scale over which the initial discontinuity is smoothed, which we set to twice the initial particle spacing. The top row of Figure 4 demonstrates that the oscillations indeed vanish in this case (compare top to bottom row).

Price (2008) showed that, in general, one requires artificial dissipation terms whenever a discontinuity occurs in SPH. Hence an alternative solution is to add an artificial conductivity term at the discontinuity. For shock capturing the usual term is of the form

[TABLE]

where $\alpha_{u}\in[0,1]$ and $v_{\rm sig}^{u}$ is the signal speed. For hydrodynamics, Price (2008) proposed

[TABLE]

which is designed to smooth artificial pressure blips which can occur at contact discontinuities. For the purpose of this paper — considering the heat equation in isolation — we adopt a simpler signal speed of the form

[TABLE]

which could be easily generalised for diffusion of any quantity. Importantly, the addition of this term does not affect the convergence of the method, since the artificial diffusion with $v^{u}_{\rm sig}$ set as above is $\mathcal{O}(h^{2})$. Given this it is safe to simply adopt $\alpha_{u}=1$.

Figures 4 and 5 demonstrate that both using the artificial conduction term (34) or starting with the smoothed initial conditions effectively eliminate the ‘carbuncle mode’ caused by the discontinuity. For this figure we also placed the particles on a more realistic particle distribution, by initially placing particles randomly in the domain and relaxing them using the usual SPH equations with a damping term. We refer to this as a ‘glass-like’ particle arrangement.

3.2 3D diffusion with constant heat conduction tensor

For our next test, we consider a domain $x,y,z\in[-1,1]$ with Dirichlet boundary conditions. Assuming an initial heat distribution in the form of a delta function located at the origin and using Green functions, we can find that the solution for the isotropic equation in 3D is given by

[TABLE]

As it is impossible numerically to start with this solution, we set the temperature at $t=0$ to be zero everywhere except a sphere of radius $\epsilon=0.1$ around the origin where it is

[TABLE]

We also considered an anisotropic case with diffusion acting only along the $x$ axis with $\kappa_{xx}=1$. As the solution is then still the product of Green functions acting in each direction, the time evolution occurs only in the function corresponding to the nonzero component of the heat conduction tensor. The exact solution is then

[TABLE]

Figure 6 shows the results for isotropic diffusion using the Brookshaw method and for anisotropic diffusion using the Espanol-Revenga method (bottom row). The corresponding exact solutions are shown in the top row, respectively. The results highlight the instability inherent in the Espanol-Revenga method when the diffusion is highly anisotropic. The solutions to the same problem performed with two first derivatives are stable (middle row).

It is possible to obtain accurate results at this resolution (642 particles) using the direct second derivative method as well, provided one employs the quintic kernel. However, the method does not guarantee stability and may become unstable, as demonstrated in Figure 7.

3.3 3D diffusion with variable heat conduction tensor

If we consider the previous problem in cylindrical coordinates $(\rho,\phi,z)$, as shown by Hopkins (2017), we can define the initial heat source according to

[TABLE]

where $\rho_{0}=0.3$ is a radial position of the source, $\delta\rho=0.05$ is a radial size of the source, and $\delta\phi=0.5$ is a azimuthal size of the source, $T_{0}=0$, $T_{1}=1$.

Then, the solution is given by

[TABLE]

If we set the heat conduction tensor to have only one nonzero component $\kappa_{\phi\phi}=1$, it means that there should be the only diffusion in the angular direction. The corresponding tensors in both cylindrical and cartesian coordinates are given by

[TABLE]

Another possible approach here would be to transform the formulae for SPH derivatives into cylindrical coordinates. We found similar results regardless of the method employed.

Figure 8 shows that two first derivative method can handle the situation where the direction of the heat flow is not aligned with the particle distribution. For the same problem, the Espanol-Revenga method results in numerical instability. The direct derivative method is an acceptable, but risky, choice, as it may become unstable depending on the particle arrangement, the number of neighbours, or kernel involved (we used the $M_{6}$ quintic spline here). By contrast, the proof of increasing entropy we gave in 2.2.4 for the two first derivatives method does not depend on any of the above.

3.4 Convergence and kernels

For each the preceding tests, we performed a convergence study, shown in Figures 9 and 10. These show that the kernels and initial particle distribution have a significant influence on the overall accuracy. The kernel bias is the region in which the accuracy cannot be increased by simply using more particles.

Figure 9 (left panel) compares the convergence properties of the Brookshaw method and two first derivatives for isotropic heat conduction with particles on a regular lattice with the Gaussian pulse as the initial temperature distribution and quintic spline (M6) as a smoothing kernel. Convergence is quadratic for both methods, even though the Brookshaw method is one order of magnitude more accurate overall. When the error reaches $\approx 10^{-4}$, it stops improving with higher resolution. The error then depends purely on the smoothing kernel. The central panel of the Figure 9 shows that the convergence rate and kernel bias stay the same when we switch to the anisotropic problem (where the Brookshaw method is no longer applicable). The right panel of Figure 9 shows that the choice of kernel and the number of neighbours sets the kernel bias. Using higher kernels in the spline series, we can obtain progressively more accurate results at high spatial resolution.

Figure 10 shows the convergence results for one dimensional heat diffusion at $t=0.025$ (Section 3.1). We adopted a glass-like lattice for different kernels and applied different methods to deal with the discontinuous initial conditions. Importantly, we see that the carbuncle mode is not an instability, since the introduced noise decreases with resolution even if left without any treatment (left panel). The convergence is linear. When we use the artificial thermal conduction term (Eq. 34) similar to the one used in hydrodynamics, the order of convergence increases to $\sim 1.5$ (central panel). Finally, from the right panel, the convergence of the same problem with continuous initial conditions is second order once again.

Comparing Figures 9 and 10 demonstrates that the kernel biases remain similar independent of the initial conditions, particle placement, whether or not diffusion is anisotropic, and whether or not a discontinuity is present. For the $M_{6}$ spline kernel it is $\sim~{}5\times 10^{-4}$.

4 Discussion

Our main finding is that the second law of thermodynamics is the most important consideration when assessing the stability of diffusion operators in SPH. This consideration leads to the conclusion that the only stable operator for anisotropic heat conduction in SPH is the ‘two first derivatives’ method described in Section 2.1.5. In the same time, we emphasised that it is crucial to alternate between differential and symmetric operators when constructing a two first derivatives method, in order to preserve positive increase of entropy. This conclusion was already reached by Price et al. (2017) in the context of physical viscosity and dust-gas mixtures. Similarly, the need for conjugate derivative operators to satisfy conservation laws is common in SPH (e.g. Monaghan, 1992; Cummins & Rudman, 1999; Tricco & Price, 2012).

Likewise, the two first derivatives method has already been used widely, e.g. for dust-gas mixtures Price & Laibe (2015), physical viscosity (Flebbe et al., 1994; Sijacki & Springel, 2006), resistive and ambipolar diffusion (Wurster et al., 2014) and the Hall effect (Wurster et al., 2016). In the above papers, the authors found the method to be the most reliable.

For the case of anisotropic diffusion, our findings give an alternative to the fix suggested by Petkova & Springel (2009). In particular, using two first derivatives does not require limiting the anisotropy of the flow in order to achieve stability. The Español & Revenga (2003) method should not be used for anisotropic diffusion.

The main caveat to using two first derivatives is the appearance of carbuncle modes if the initial conditions are discontinuous. The problem disappears with smooth initial conditions. Another solution is to introduce additional artificial conduction term that acts only when the solution (e.g. in pressure) is discontinuous, as proposed by Price (2008). These terms are similar to the diffusion that would arise in a Finite Volume scheme when fluxes are reconstructed at discontinuous interfaces. The ‘integral Godunov’ methods proposed by Hopkins (2017) indeed require flux limiters to ensure that the entropy is positive. Hopkins (2017) pointed out that symmetric operators, in general, can result in low-order convergence if the particles are disordered. We found this can be mitigated by the use of smoother spline kernels. We found the Español & Revenga (2003) method to be significantly less accurate.

An interesting follow-up would be to search for a kernel function that satisfies condition (26). If such a kernel exists, it may be possible to take direct second derivatives that are both stable and give non-oscillating solutions. An interesting application beyond the scope of this paper would be to the Magneto-Thermal Instability (see Hopkins, 2017).

5 Conclusions

We analysed the stability of methods for anisotropic diffusion in SPH. Our conclusions are:

In case of isotropic diffusion the Brookshaw method is stable and also the most accurate method. 2. 2.

Two first derivatives are the only method for anisotropic diffusion where stability is guaranteed. The only caveat is that smoothing or an artificial diffusion term is required for accurate treatment of discontinuities. 3. 3.

We recommend against the use of the Espanol-Revenga method for anisotropic diffusion because it is not only unstable under certain circumstances but also inaccurate. 4. 4.

We find that use of the M6 quintic spline kernel reduces the kernel bias by approximately one order of magnitude compared to the cubic spline kernel for both isotropic and anisotropic diffusion. 5. 5.

For the two first derivatives method one can use a timestep 3–8 times larger (depending on the choice of kernel) than for the Brookshaw method, while remaining stable. This offers a potentially large cost saving.

Acknowledgements

We thank the organisers of SPHERIC2018 for a useful and stimulating conference in Galway, Ireland. DP acknowledges funding from the Australian Research Council via FT130100034 and DP180104235.