© ESO, 2012

1. Introduction

Smoothed particle hydrodynamics (SPH) is a Lagrangian, mesh-free particle method that is used to model fluid hydrodynamics in numerical simulations. The technique was originally developed in an astrophysical context (Lucy 1977; Gingold & Monaghan 1977), but since then it has been widely applied in many other areas (Monaghan 2005) of computational fluid dynamics.

The method has several properties that make its use particularly advantageous in astrophysical problems (Hernquist & Katz 1989; Rosswog 2009; Springel 2010a). Because of its Lagrangian nature, the development of high matter concentrations in collapse problems is followed naturally. Moreover, the method is free of both advection errors and Galilean invariant. Finally, the method naturally incorporates self-gravity and possesses very good conservation properties, its fluid equations being derived from variational principles.

The other computational fluid dynamical method commonly employed in numerical astrophysics is the standard Eulerian grid-based approach in which the fluid is evolved on a discretized mesh (Stone & Norman 1992; Ryu et al. 1993; Norman & Bryan 1999; Fryxell et al. 2000; Teyssier 2002). The spatial resolution of an Eulerian scheme based on a fixed Cartesian grid is often insufficient, however, to adequately resolve the high dynamic range frequently encountered in many astrophysical problems, such as galaxy formation. This has motivated the development of adaptative mesh refinement (AMR) methods, in which the spatial resolution of the grid is locally refined according to some selection criterion (Berger & Colella 1989; Kravtsov, Klypin & Khokhlov 1997; Norman 2005). Additionally, the order of the numerical scheme has been improved by adopting the parabolic piecewise method (PPM) of Colella & Woodward (1984), such as in the AMR Eulerian codes ENZO (Norman & Bryan 1999) and FLASH (Fryxell et al. 2000).

Application of these different types of hydrodynamical codes to the same test problem with identical initial conditions should in principle lead to similar results. However, there has been growing evidence over recent years that in a variety of hydrodynamical test cases there are significant differences between the results produced from the two types of methods (O’Shea et al. 2005; Agertz et al. 2007; Wadsley et al. 2008; Tasker et al. 2008; Mitchell et al. 2009; Read et al. 2010; Valcke et al. 2010; Junk et al. 2010).

For instance, Agertz et al. (2007) show that in the standard SPH formulations the growth of Kelvin-Helmholtz (KH) instabilities in shear flows is artificially suppressed when steep density gradients are present at the contact discontinuities. Moreover, the level of central entropy produced in binary-merger non-radiative simulations of galaxy clusters is significantly higher by a factor  ~2 in Eulerian simulations than in those made using SPH codes (Mitchell et al. 2009). The origin of these discrepancies has been recognized (Agertz et al. 2007; Read et al. 2010) as partly due to the intrinsic difficulty for SPH in properly modeling density gradients across fluid interfaces, which in turn implies that there is a surface tension effect that inhibits the growth of KH instabilities. A second problem is the Lagrangian nature of SPH codes, which prevents the mixing of fluid elements at the particle level and leads to entropy generation (Wadsley et al. 2008; Mitchell et al. 2009). In particular, Mitchell et al. (2009) argue that the main explanation of the different levels of central entropy found in cluster simulations is the different degree of entropy mixing present in the two codes. Unlike SPH, in Eulerian codes fluid mixing occurs by definition at the cell level and a certain degree of overmixing is certainly present in simulations made using mesh-based codes (Springel 2010b).

Given the advantages of the SPH codes highlighted previously, it appears worth pursuing any improvement in the SPH method capable of overcoming its present limitations. Along this line of investigation, many efforts have been made by a number of authors (Abell 2011; Price 2008; Wadsley et al. 2008; Valcke et al. 2010; Read et al. 2010; Heß & Springel 2010; Cha et al. 2010; Murante et al. 2011).

Price (2008) introduces an artificial heat conduction term in the SPH equations with the purpose of smoothing thermal energy at fluid interfaces. This artificial conductivity (AC) term in turn gives a smooth entropy transition at contact discontinuities with the effect of enforcing pressure continuity and removing the artificial surface-tension effect that inhibits the growth of KH instabilities at fluid interfaces. Similarly, Wadsley et al. (2008) suggest that in SPH the lack of mixing at the particle level can be alleviated by adding a heat diffusion term to the equations so as to mimic the effects of subgrid turbulence, thereby improving the amount of mixing.

Read et al. (2010) present an SPH implementation in which a modified density estimate is adopted (Ritchie & Thomas 1991), together with the use of a peaked kernel and a much larger number of neighbors. The authors show that the new scheme is capable of following the development of fluid instabilities in a more improved way than in standard SPH. Abell (2011) presents an alternative derivation of the SPH force equation that avoids the problem encountered by standard SPH in handling fluid instabilities, although the approach is inherently neither energy nor momentum conserving and is prone to large integration errors when the simulation resolution is low.

The method proposed by Inutsuka (2002) reformulates the SPH equations by introducing a kernel convolution so as to consistently calculate the density and hydrodynamic forces. The latter are determined using a Riemann solver (Godunov SPH). The method was revisited by Cha et al. (2010) and Murante et al. (2011), who show that the code correctly follows the development of fluid instabilities in a variety of hydrodynamic tests.

A deeper modification than those presented here so far was introduced by Heß & Springel (2010), who replaced the traditional SPH kernel approach with a new density estimate based on Voronoi tesselation. The authors show that the method is free of surface tension effects and therefore the growth rate of fluid instabilities is not as adversely affected as in standard SPH.

Finally, a radically new numerical scheme was developed by Springel (2010b) to retain the advantages of both SPH and mesh-based codes. In the new code, the hydrodynamic equations are solved on a moving unstructured mesh using a Godunov method with an exact Riemann solver. The mesh is defined by the Voronoi tesselation of a set of discrete points and allowed to move freely with the fluid. The method is therefore adaptative in nature and thus Galilean invariant but, at the same time, the accuracy with which shocks and contact discontinuities are described is that of an Eulerian code. Bauer & Springel (2012) argue that the standard formulation of SPH fails to accurately resolve the development of turbulence in the subsonic regime (but see Price 2012b, for a different viewpoint). The authors draw their conclusions by analyzing results from simulations of driven subsonic turbulence made using the new moving-mesh code, named AREPO, and a standard SPH code.

Similar conclusions were reached in a set of companion papers (Sijacki et al. 2012; Vogelsberger et al. 2011), in which the new code was used in galaxy formation studies to demonstrate its superiority over standard SPH. However, the code is characterized by considerable complexity making the use of the SPH scheme still appealing, and more generally, it is desirable that simulation results produced with a specific code should be reproduced with a completely independent numerical scheme when complex non-linear phenomena are involved. It appears worthwhile, therefore, to investigate, along the line of previous authors, the possibility of constructing a numerical scheme based on the traditional SPH formulation that is capable of correctly describing the development of fluid instabilities and at the same time incorporates the effects of self-gravity when present.

This is the aim of the present study, in which the SPH scheme is modified by incorporating an AC diffusion term into the equations as described by Price (2008). However, in Price (2008) the strength of the AC is governed by a signal velocity that is based on pressure discontinuities. For simulations where gravity is present, this approach is not applicable because hydrostatic equilibrium requires pressure gradients. An appropriate signal velocity for conductivity when gravity is present is then used to construct an AC-SPH code with the purpose of treating the growth of fluid instabilities self-consistently. The viability of the approach is tested using a suite of test problems in which results obtained using the new code are contrasted with the corresponding ones produced in previous work using different schemes. The code is very similar in form to that presented by Wadsley et al. (2008), but here the energy diffusion equation is implemented in a different manner.

The paper is organized as follows. Section 2 presents the hydrodynamical method and introduces the AC approach. In Sect. 3, we investigate the effectiveness of the method by presenting results from a suite of purely hydrodynamical test problems, namely the two-dimensional Kelvin-Helmholtz instability, the Sod shock tube, the point explosion or Sedov blast-wave test, and the blob test. In Sect. 4, we then discuss the results of hydrodynamic tests that include self-gravity. Specifically, we consider the cold gas sphere or Evrard collapse test, the Rayleigh-Taylor instability, and the cosmological integration of galaxy clusters. Finally, the main results are summarized in Sect. 5.

2. The hydrodynamic method

Here we present the basic features of the method and for reviews refer to Rosswog (2009), Springel (2010a), and Price (2012a).

2.1. Basic SPH equations

In the SPH method, the fluid is described by a set of N particles with mass m_i, velocity v_i, density ρ_i, and a thermodynamic variable such as the specific thermal energy u_i or the entropic function A_i. The particle pressure is then defined as ${\mathit{P}}_{\mathit{i}}\mathrm{=}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{\rho }}_{\mathit{i}}{\mathit{u}}_{\mathit{i}}\mathrm{=}{\mathit{A}}_{\mathit{i}}{\mathit{\rho }}_{\mathit{i}}^{\mathit{\gamma }}$, where γ = 5/3 for a monoatomic gas. The density estimate ρ(r) at the particle position r_i is given by ${\mathit{\rho }}_{\mathit{i}}\mathrm{=}\sum _{\mathit{j}}{\mathit{m}}_{\mathit{j}}\mathit{W}\mathrm{\left(}\left|{r}\mathit{ij}\right|\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathit{,}$(1)where r_ij ≡ r_i − r_j, W( | r_ij | ,h_i) is the interpolating kernel that has compact support and is zero for  | r_ij |  ≥ ζh_i (Price 2012a). The kernel is normalized to the condition ^∫Wd^Dr = 1. The sum in Eq. (1) is over a finite number of particles and the smoothing length h_i is a variable that is implicitly defined by the equation ${\mathit{f}}_{\mathit{D}}{\left(\mathit{\zeta }{\mathit{h}}_{\mathit{i}}\right)}^{\mathit{D}}{\mathit{\rho }}_{\mathit{i}}\mathrm{=}{\mathit{N}}_{\mathrm{sph}}{\mathit{m}}_{\mathit{i}}\mathit{,}$(2)where D is the number of spatial dimensions, f_D = π, 4π/3 for D = 2, 3, and N_sph is the number of neighboring particles within a radius ζh_i. This equation can be rewritten by defining h_i in units of the mean interparticle separation ${\mathit{h}}_{\mathit{i}}\mathrm{=}\mathit{\eta }{\left({\mathit{m}}_{\mathit{i}}\mathit{/}{\mathit{\rho }}_{\mathit{i}}\right)}^{\mathrm{1}\mathit{/}\mathit{D}}\mathit{,}$(3)so that ${\mathit{N}}_{\mathrm{sph}}^{\mathrm{2}\mathit{D}}\mathrm{=}\mathit{\pi }\mathrm{(}\mathit{\zeta \eta }{\mathrm{)}}^{\mathrm{2}}$ and ${\mathit{N}}_{\mathrm{sph}}^{\mathrm{3}\mathit{D}}\mathrm{=}\mathrm{4}\mathit{\pi }\mathrm{(}\mathit{\zeta \eta }{\mathrm{)}}^{\mathrm{3}}\mathit{/}\mathrm{3}$. The smoothing length h_i is determined by solving the non-linear Eq. (2), and it should be kept in mind that N_sph does not necessarily need to be an integer but can take arbitrary values if η is used as the fundamental parameter determining h_i (Price 2012a). A kernel commonly employed is the M₄ (cubic spline), which is zero for ζ ≥ 2. In 3D, typical choices of N_sph lie in the range N_sph ~ 33 − 64, which for the M₄ kernel corresponds to η ≃ 1. − 1.25. The equation of motion for the SPH particles can be derived from the Lagrangian of the system (Springel 2010a; Price 2012a) and is given by $\frac{\mathrm{d}{v}\mathit{i}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\sum _{\mathit{j}}{\mathit{m}}_{\mathit{j}}\left[\frac{{\mathit{P}}_{\mathit{i}}}{{\mathrm{\Omega }}_{\mathit{i}}{\mathit{\rho }}_{\mathit{i}}^{\mathrm{2}}}\mathbf{\nabla }\mathit{i}{\mathit{W}}_{\mathit{ij}}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}\mathrm{+}\frac{{\mathit{P}}_{\mathit{j}}}{{\mathrm{\Omega }}_{\mathit{j}}{\mathit{\rho }}_{\mathit{j}}^{\mathrm{2}}}\mathbf{\nabla }\mathit{i}{\mathit{W}}_{\mathit{ij}}\mathrm{\left(}{\mathit{h}}_{\mathit{j}}\mathrm{\right)}\right]\mathit{,}$(4)where the coefficients Ω_i are defined as ${\mathrm{\Omega }}_{\mathit{i}}\mathrm{=}\left[\mathrm{1}\mathrm{-}\frac{\mathit{\partial }{\mathit{h}}_{\mathit{i}}}{\mathit{\partial }{\mathit{\rho }}_{\mathit{i}}}\sum _{\mathit{k}}{\mathit{m}}_{\mathit{k}}\frac{\mathit{\partial }{\mathit{W}}_{\mathit{ik}}\mathrm{\left(}{\mathit{h}}_{\mathit{i}}\mathrm{\right)}}{\mathit{\partial }{\mathit{h}}_{\mathit{i}}}\right]\mathrm{·}$(5)These terms are present in the momentum Eq. (4) because the smoothing length h_i itself is implicitly a function of the particle coordinates through Eq. (2).

2.2. Artificial viscosity

In SPH, the momentum equation must be generalized by adding an appropriate viscous force, which is aimed at correctly treating the effects of shocks. An artificial viscosity (AV) term is then introduced with the purpose of dissipating local velocities and preventing particle interpenetration at the shock locations. The new term is given by ${\left(\frac{\mathrm{d}{v}\mathit{i}}{\mathrm{d}\mathit{t}}\right)}_{\mathrm{AV}}\mathrm{=}\mathrm{-}\sum _{\mathit{i}}{\mathit{m}}_{\mathit{j}}{\mathrm{\Pi }}_{\mathit{ij}}\mathbf{\nabla }\mathit{i}\mathit{W̅}\mathit{ij}\mathit{,}$(6)where the term $\mathit{W̅}{}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\left(}\mathit{W}\mathrm{\right(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{i}}\mathrm{)}\mathrm{+}\mathit{W}\mathrm{(}{\mathit{r}}_{\mathit{ij}}\mathit{,}{\mathit{h}}_{\mathit{j}}\mathrm{\left)}\mathrm{\right)}$ is the symmetrized kernel and Π_ij is the AV tensor. In the SPH entropy formulation (Springel & Hernquist 2002), it is the entropy function per particle A_i that is integrated and its time derivative is calculated as follows $\frac{\mathrm{d}{\mathit{A}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\gamma }\mathrm{-}\mathrm{1}}{{\mathit{\rho }}_{\mathit{i}}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\sum _{\mathit{j}}{\mathit{m}}_{\mathit{j}}{\mathrm{\Pi }}_{\mathit{ij}}{v}\mathit{ij}\mathrm{·}{\mathrm{\nabla }}_{\mathit{i}}\mathit{W̅}\mathit{ij}\mathrm{\equiv }\frac{\mathit{\gamma }\mathrm{-}\mathrm{1}}{{\mathit{\rho }}_{\mathit{i}}^{\mathit{\gamma }\mathrm{-}\mathrm{1}}}{\left(\frac{\mathrm{d}{\mathit{u}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\right)}_{\mathrm{AV}}\mathit{,}$(7)where v_ij = v_i − v_j. For the AV tensor, we adopt here the form proposed by Monaghan (1997) in analogy with the Riemann solvers $\begin{array}{ccc}{\mathrm{\Pi }}_{\mathit{ij}}\mathrm{=}\mathrm{-}\frac{{\mathit{\alpha }}_{\mathit{ij}}}{\mathrm{2}}\frac{{\mathit{v}}_{\mathit{ij}}^{\mathrm{AV}}{\mathit{\mu }}_{\mathit{ij}}}{{\mathit{\rho }}_{\mathit{ij}}}{\mathit{f}}_{\mathit{ij}}\mathit{,}& & \end{array}$(8)where ${\mathit{v}}_{\mathit{ij}}^{\mathrm{AV}}\mathrm{=}{\mathit{c}}_{\mathit{i}}\mathrm{+}{\mathit{c}}_{\mathit{j}}\mathrm{-}\mathrm{3}{\mathit{\mu }}_{\mathit{ij}}$(9)is the signal velocity and μ_ij = v_ij·r_ij/ | r_ij |  if v_ij·r_ij < 0 but zero otherwise, so that Π_ij is non-zero only for approaching particles. Here scalar quantities with the subscripts i and j denote arithmetic averages, c_i is the sound speed of particle i, the parameter α_i regulates the amount of AV, and f_i is a controling factor that reduces the strength of AV in the presence of shear flows. The latter is given by Balsara (1995), as ${\mathit{f}}_{\mathit{i}}\mathrm{=}\frac{{\left|\mathbf{\nabla }\mathrm{·}{v}\right|}_{\mathit{i}}}{{\left|\mathbf{\nabla }\mathrm{·}{v}\right|}_{\mathit{i}}\mathrm{+}{\left|\mathbf{\nabla }\mathrm{\times }{v}\right|}_{\mathit{i}}}\mathit{,}$(10)where (∇·v)_i and (∇ × v)_i are the standard SPH estimates for divergence and curl (Monaghan 2005). For pure shear flows,  | ∇ × v | _i ≫  | ∇·v | _i, so that the AV is strongly damped.

Using the signal velocity (9), the Courant condition on the timestep of particle i is given by $\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}^{\mathrm{C}}\mathrm{\simeq }\mathrm{0.3}\frac{{\mathit{h}}_{\mathit{i}}}{\mathit{ma}{\mathit{x}}_{\mathit{j}}\left|{\mathit{v}}_{\mathit{ij}}^{\mathrm{AV}}\right|}\mathrm{·}$(11)In the standard SPH formulation, the viscosity parameter α_i controling the strength of the AV is given by α_i = constant ≡ α₀, with α₀ = 1 being a common choice (Monaghan 2005). This scheme has the disadvantage that it generates viscous dissipation in regions of the flow that are not undergoing shocks. To reduce any spurious viscosity effects, Morris & Monaghan (1997) proposed making the viscous coefficient α_i time-dependent so that it can evolve in time under certain local conditions. The authors proposed an equation of the form $\frac{\mathrm{d}{\mathit{\alpha }}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{{\mathit{\alpha }}_{\mathit{i}}\mathrm{-}{\mathit{\alpha }}_{\min }}{{\mathit{\tau }}_{\mathit{i}}}\mathrm{+}{\mathit{S}}_{\mathit{i}}\mathit{,}$(12)where S_i is a source term and τ_i regulates the decay of α_i to the floor value α_min away from shocks. For the source term, the adopted expression is ${\mathit{S}}_{\mathit{i}}\mathrm{=}{\mathit{S}}_{\mathrm{0}}{\mathit{f}}_{\mathit{i}}\left(\mathrm{-}{\mathrm{(}\mathbf{\nabla }\mathrm{·}{v}\mathrm{)}}_{\mathit{i}}\mathit{,}\mathrm{0}\right)\mathrm{(}{\mathit{\alpha }}_{\max }\mathrm{-}{\mathit{\alpha }}_{\mathit{i}}\mathrm{)}\mathit{,}$(13)which is constructed in such a way that it increases in the presence of shocks. The prefactor S₀ is unity for γ = 5/3 and the damping factor f_i is inserted to account for the presence of vorticity. The original form proposed by Morris & Monaghan (1997) was refined by the factor (α_max − α_i) (Rosswog et al. 2000), which has the advantage of being more sensitive to shocks than the original formulation and of preventing α_i from becoming higher than a prescribed value α_max. Cullen & Dehnen (2011) presented an improved version of the Morris & Monaghan (1997) AV scheme, employing as a shock indicator a switch based on the time derivative of the velocity divergence.

The decay parameter τ_i is of the form ${\mathit{\tau }}_{\mathit{i}}\mathrm{=}\frac{{\mathit{h}}_{\mathit{i}}}{{\mathit{c}}_{\mathit{i}} {\mathit{l}}_{\mathrm{d}}}\mathit{,}$(14)where l_d is a dimensionless parameter that controls the decay timescale. In a number of test simulations, Rosswog et al. (2000) found that appropriate values for the parameters α_max,α_min, and l_d are 1.5, 0.05, and 0.2, respectively. In principle, the effects of numerical viscosity in regions away from shocks can be reduced by setting l_d to higher values than l_d = 0.2. In practice, the minimum time necessary to propagate through the resolution length h_i sets the upper limit l_d = 1. However, the time evolution of the viscosity parameter can be affected if very short damping timescales are imposed. Neglecting variations in the coefficients, the solution to Eq. (12) at times t > t_in can be written as ${\mathit{\alpha }}_{\mathit{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{q}}_{\mathit{i}}\mathrm{+}\left({\mathit{\alpha }}_{\mathit{i}}\mathrm{\left(}{\mathit{t}}_{\mathrm{in}}\mathrm{\right)}\mathrm{-}{\mathit{q}}_{\mathit{i}}\right){\exp }^{\mathrm{-}\left(\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{in}}\right)\mathit{/}{\mathit{\tau }}_{\mathit{i}}^{\mathrm{\prime }}}\mathit{,}$(15)where ${\mathit{\tau }}_{\mathit{i}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{\tau }}_{\mathit{i}}}{\mathrm{1}\mathrm{+}{\mathit{S}}_{\mathit{i}}{\mathit{\tau }}_{\mathit{i}}}\mathit{,}$(16)and q_i is a modified source term ${\mathit{q}}_{\mathit{i}}\mathrm{=}\frac{{\mathit{\alpha }}_{\min }\mathrm{+}{\mathit{S}}_{\mathit{i}}{\mathit{\tau }}_{\mathit{i}}{\mathit{\alpha }}_{\max }}{\mathrm{1}\mathrm{+}{\mathit{S}}_{\mathit{i}}{\mathit{\tau }}_{\mathit{i}}}\mathit{.}$(17)From Eq. (15), it can be seen that α_i(t) ≃ α_max in the strong shock regime S_iτ_i ≫ 1 but this condition is not satisfied if S_iτ_i ≲ 1. Therefore, for mild shocks this implies that, if very short decay timescales are imposed, the peak value of α_i(t) at the shock front might be below the AV strength necessary to properly treat shocks.

To overcome this difficulty, a modification to Eq. (13) was adopted (Valdarnini 2011, hereafter V11) that, when l_d = 1, compensates for the smaller values of S_iτ_i with respect to the small l_d regimes. This is equivalent to considering a higher value for α_max, so that in Eq. (13) α_max is substituted by α_max → ξα_max.

The correction factor ξ was calibrated using the shock-tube problem as reference and requiring a peak value of  ~0.6−0.7 for the viscosity parameter at the shock front, as in the l_d = 0.2 case. The results indicate that a normalization of the form ξ = (l_d/0.2)^0.8 for l_d ≥ 0.2 satisfies these constraints. This normalization has been validated in a number of other test problems showing that it is able to produce a peak value of the viscosity parameter at the shock location that is independent of the chosen value of the decay parameter l_d.

The time-dependent AV formulation of SPH has been shown to be effective in reducing the damping of turbulence due to the effects of numerical viscosity in simulations of galaxy clusters (Dolag et al. 2005, V11). Moreover, it was used in Price (2012b) to argue that the conclusion of Bauer & Springel (2012) about the difficulty for SPH codes in properly modeling the development of subsonic turbulence is based on using only a SPH code in its standard AV formulation. In the following, unless otherwise stated, we describe simulations performed using a time-dependent AV, this will be fully specified by the set of parameters  { α_min,α_max,l_d } .

2.3. Artificial conductivity

As already outlined in the Introduction, different formulations have been proposed for overcoming the problems encountered by standard SPH in the treatment of fluid discontinuities. Here we follow the approach suggested by Price (2008), who proposed adding a dissipative term to the thermal energy equation for smoothing the energy across contact discontinuities. The presence of these dissipative terms introduces a smoothing of entropy at fluid interfaces with the effect of removing pressure discontinuities.

The motivations behind this approach were discussed by Monaghan (1997), who showed how the momentum and energy equations in SPH must contain a dissipative term related to jumps in the variables across characteristics, in analogy with the corresponding Riemann solutions.

The artificial conductivity (AC) term for the dissipation of energy takes the form ${\left(\frac{\mathrm{d}{\mathit{u}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\right)}_{\mathrm{AC}}\mathrm{=}\sum _{\mathit{j}}\frac{{\mathit{m}}_{\mathit{j}}{\mathit{v}}_{\mathit{ij}}^{\mathrm{AC}}}{{\mathit{\rho }}_{\mathit{ij}}}\mathrm{[}{\mathit{\alpha }}_{\mathit{ij}}^{\mathrm{C}}{\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{-}{\mathit{u}}_{\mathit{j}}\mathrm{)}}^{\mathrm{]}}{{e}}_{\mathit{ij}}\mathrm{·}\mathbf{\nabla }\mathit{i}\mathit{W̅}\mathit{ij}\mathit{,}$(18)where ${\mathit{v}}_{\mathit{ij}}^{\mathrm{AC}}$ is the signal velocity, which does not need to be the same as that used in the momentum equation, e_ij ≡ r_ij/r_ij, and ${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$ is the AC parameter, which is of order unity¹. Equation (18) represents the SPH analogue of a diffusion equation of the form ${\left(\frac{\mathrm{d}{\mathit{u}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\right)}_{\mathrm{AC}}\mathrm{\simeq }{\mathit{D}}_{\mathit{i}}^{\mathrm{AC}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathit{i}}\mathit{,}$(19)where ${\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathit{i}}\mathrm{=}\mathrm{2}\sum _{\mathit{j}}{\mathit{m}}_{\mathit{j}}\frac{{\mathit{u}}_{\mathit{i}}\mathrm{-}{\mathit{u}}_{\mathit{j}}}{{\mathit{\rho }}_{\mathit{j}}}\frac{{e}\mathit{ij}\mathrm{·}\mathbf{\nabla }{\mathit{W}}_{\mathit{ij}}}{{\mathit{r}}_{\mathit{ij}}}$(20)is the SPH expression for the Laplacian (Brookshaw 1985) and, in analogy with the analysis of Lodato & Price (2010) for defining an equivalent physical viscosity coefficient for the SPH numerical viscosity, we can define ${\mathit{D}}_{\mathit{i}}^{\mathrm{AC}}$ as a numerical heat-diffusion coefficient given by ${\mathit{D}}_{\mathit{i}}^{\mathrm{AC}}\mathrm{\simeq }\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}{\mathit{v}}_{\mathit{ij}}^{\mathrm{AC}}{\mathit{r}}_{\mathit{ij}}\mathit{.}$(21)An important issue concerns the choice of the AC parameter ${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$, which must be constructed so that diffusion of thermal energy is kept under control and is damped away from discontinuities. This can be achieved by introducing a switch similar to that devised for AV. The dissipation parameter is then evolved according to $\frac{\mathrm{d}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}\mathrm{-}{\mathit{\alpha }}_{\min }^{\mathrm{C}}}{{\mathit{\tau }}_{\mathit{i}}^{\mathrm{C}}}\mathrm{+}{{\mathit{S}}^{\mathrm{C}}}_{\mathit{i}}\mathit{,}$(22)the meaning of the terms being similar to that of Eq. (12). In particular, for the decay timescale ${\mathit{\tau }}_{\mathit{i}}^{\mathrm{C}}\mathrm{=}{\mathit{h}}_{\mathit{i}}\mathit{/}\mathrm{𝒞}{\mathit{c}}_{\mathit{i}}$ we here set and the floor value ${\mathit{\alpha }}_{\min }^{\mathrm{C}}$ to zero. For the source term ${\mathit{S}}_{\mathit{i}}^{\mathrm{C}}$, we use the expression ${{\mathit{S}}^{\mathrm{C}}}_{\mathit{i}}\mathrm{=}{\mathit{f}}_{\mathrm{C}}{\mathit{h}}_{\mathit{i}}\frac{\begin{array}{}\mathrm{|}\\ \mathrm{|}\\ \mathrm{|}\end{array}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathit{i}}\begin{array}{}\mathrm{|}\\ \mathrm{|}\\ \mathrm{|}\end{array}}{\sqrt{{\mathit{u}}_{\mathit{i}}\mathrm{+}\mathit{\epsilon }}}\mathrm{(}{\mathit{\alpha }}_{\max }^{\mathrm{C}}\mathrm{-}{{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}}^{\mathrm{)}}\mathit{,}$(23)where f_C is a dimensionless parameter of order unity and the parameter ε avoids divergences when u_i → 0. This equation differs from the corresponding source term proposed by Price & Monaghan (2005) and Price (2008) by the factor ${\mathit{\alpha }}_{\max }^{\mathrm{C}}\mathrm{-}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$, which has been inserted here in analogy with Eq. (13) for introducing a stronger response by the switch in the presence of discontinuities. The choice of values for the parameters f_C and ${\mathit{\alpha }}_{\max }^{\mathrm{C}}$ depends on the problem under consideration, for example Price & Monaghan (2005) proposed f_C = 0.1. Here a number of numerical experiments has shown that the best results in terms of the response of the AC parameter ${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$ to the presence of discontinuities are obtained by setting f_C = 1 and ${\mathit{\alpha }}_{\max }^{\mathrm{C}}\mathrm{=}\mathrm{1.5}$, which we assume henceforth as reference values. Moreover, it is found that significant benefits in terms of the sharpness of the AC parameter profile across the discontinuity are obtained by inserting the ${\mathit{\alpha }}_{\max }^{\mathrm{C}}\mathrm{-}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$ term. In principle, the choice of the derivative term used to detect discontinuities is arbitrary, but in practice a second derivative (Price & Monaghan 2005) term ensures higher sensitivity to sharp discontinuities in thermal energy than a first derivative (Price 2005a).

An example of a signal velocity specifically designed to remove pressure gradients at contact discontinuities is that originally introduced by Price (2008) for pure hydrodynamical simulations $\begin{array}{ccc}{\mathit{v}}_{\mathit{ij}}^{\mathrm{AC}}\mathrm{\left(}\mathit{P}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\left|{\mathit{P}}_{\mathit{i}}\mathrm{-}{\mathit{P}}_{\mathit{j}}\right|}{{\mathit{\rho }}_{\mathit{ij}}}}\mathit{,}& & \end{array}$(24)and further refined by Valcke et al. (2010). The ability of the new AC formulation of SPH to follow the development of KH instabilities using this expression for the signal velocity was verified in a number of tests (Price 2008; Valcke et al. 2010; Merlin et al. 2010). However, the disadvantage of using the signal speed (24) is that it cannot be applied when self-gravity is considered because in that case a pressure gradient is present at hydrostatic equilibrium. Imagine, for example, a self-gravitating gas sphere in hydrostatic equilibrium in which there is a negative radial temperature gradient, with the gas temperature decreasing radially outward from the center of the sphere. An application of the SPH equations, with the AC term of Eq. (22) now using the signal velocity (24), leads to a heat flux from the inner to the outer regions and, in the long term, to the development of an unphysical isothermal temperature profile. A signal velocity that avoids this difficulty is given by $\begin{array}{ccc}{\mathit{v}}_{\mathit{ij}}^{\mathrm{AC}}\mathrm{\left(}\mathit{grav}\mathrm{\right)}\mathrm{=}\left|\mathrm{(}{v}\mathit{i}\mathrm{-}{v}{\mathit{j}}^{\mathrm{)}}\mathrm{·}{r}\mathit{ij}\right|/{\mathit{r}}_{\mathit{ij}}\mathit{,}& & \end{array}$(25)which corresponds to the formulation proposed by Wadsley et al. (2008) in which a dissipative term is added to the evolution of the thermal equation with the purpose of modeling heat diffusion due to turbulence. The new term is constructed in analogy with the subgrid-scale model of Smagorinsky (1963) and is given by $\begin{array}{ccc}{\left(\frac{\mathrm{d}{\mathit{u}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\right)}_{\mathrm{AC}}\mathrm{=}\mathit{C}\sum _{\mathit{j}}\frac{{\mathit{m}}_{\mathit{j}}\left|{v}\mathit{i}\mathrm{-}{v}\mathit{j}\right|\mathrm{(}{\mathit{h}}_{\mathit{i}}\mathrm{+}{\mathit{h}}_{\mathit{j}}\mathrm{)}}{{\mathit{\rho }}_{\mathit{ij}}}\mathrm{(}{\mathit{u}}_{\mathit{i}}\mathrm{-}{\mathit{u}}_{\mathit{j}}\mathrm{)}{e}\mathit{ij}\mathrm{·}\mathbf{\nabla }\mathit{i}\mathit{W̅}\mathit{ij}\mathit{,}& & \end{array}$(26)where C is a coefficient of order unity, whose precise value depends on the problem under consideration. For the rising hot bubble problem considered by Wadsley et al. (2008), the closest agreement with the corresponding PPM results taken as reference is obtained by setting C = 0.1, higher values being too diffusive. Throughout this paper, the SPH simulations of the hydrodynamic tests are performed by adopting the expression (18) for the thermal energy dissipative term. With respect to the formulation of Wadsley et al. (2008), this approach has the advantage of using a diffusion parameter that is not constant but evolves in time according to a source term, so that the amount of diffusion away from discontinuities is minimized. For the AC signal velocity, we then use Eq. (25), whose performances in the AC formulation (18) has not been fully tested before in SPH simulations of hydrodynamic test problems.

We note that, using the signal velocity (25), the Von Neumann stability criterion becomes unimportant with respect the Courant condition (11). The Von Neumann constraint requires Δt ≤ 0.5Δx²/D, which in SPH is given by $\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}^{\mathrm{AC}}\lesssim \frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{h}}_{\mathit{i}}^{\mathrm{2}}}{{\mathit{D}}_{\mathit{i}}^{\mathrm{AC}}}\mathrm{·}$(27)Since α^C ≃ O(1), this condition implies that $\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}^{\mathrm{AC}}\gtrsim \mathrm{\Delta }{\mathit{t}}_{\mathit{i}}^{\mathrm{C}}$ in both the supersonic and subsonic regimes.

3. Hydrodynamic tests

In the following, simulation results obtained by applying the new AC-SPH code to a number of hydrodynamic test problems are discussed with the objective of validating the code and assessing its performance. To this end, the simulation results of the tests will be compared with the corresponding ones obtained by previous authors using different codes and/or numerical methods. The problems considered are usually presented in the literature in dimensional or complexity order, but here we follow a different approach. Since the KH instability is the hydrodynamic test that initially (Agertz et al. 2007) inspired the discussion about the difficulties of standard SPH in properly handling the development of instabilities, and it is also the most widely considered in this context (Abell 2011; Price 2008; Cha et al. 2010; Heß & Springel 2010; Read et al. 2010; Valcke et al. 2010; Murante et al. 2011), we here discuss first the two-dimensional (2D) version of the test in detail. The setup of the other runs is then considered in the light of the results obtained from the 2D KH test.

3.1. 2D Kelvin-Helmholtz instability

The KH instability arises in the presence of shear flow between two fluid layers when a small velocity perturbation is imposed in the direction perpendicular to the interface between the two fluids. The development of the instability is characterized by an initial phase, where the fluids interpenetrate each other, and then the forming of vortices, which become progressively more pronounced and turn into KH rolls at the onset of non-linearity. In the case of incompressible fluids for a sinusoidal perturbation of wavelength λ, this phase is reached on a timescale (Chandrasekhar 1961) ${\mathit{\tau }}_{\mathrm{KH}}\mathrm{=}\frac{\mathit{\lambda }\left({\mathit{\rho }}_{\mathrm{1}}\mathrm{+}{\mathit{\rho }}_{\mathrm{2}}\right)}{{\left({\mathit{\rho }}_{\mathrm{1}}{\mathit{\rho }}_{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{v}}\mathit{,}$(28)where ρ₁ and ρ₂ are the fluid densities and v = v₁ − v₂ is the relative shear velocity. A proper modeling of instabilities in numerical hydrodynamics is essential since the KH instability plays a crucial role in the development of turbulence that occurs in many hydrodynamical phenomena. In particular, the KH instability is relevant to many astrophysical processes, such as gas stripping from cold gas clouds occurring in galactic halos and the production of entropy in the intracluster gas of galaxy clusters due to the injection of turbulence.

The growth of the KH instability in hydrodynamic simulations has been addressed by various authors (Abell 2011; Price 2008; Wadsley et al. 2008; Read et al. 2010; Valcke et al. 2010; Heß & Springel 2010; Cha et al. 2010; Murante et al. 2011). These studies found that the development of the instability is artificially suppressed because of two distinct effects: the first problem is the so-called local mixing instability (LMI) due to entropy conservation, which inhibits mixing on the kernel scale and thus introduces pressure discontinuities; the second problem is caused by errors in the momentum equation, which cannot be reduced by increasing the number of neighbor particles without causing particle clumping.

Given the wide variety of numerical parameters and initial conditions with which the KH instability has been addressed, we choose here to perform the tests using as reference the simulations presented by Valcke et al. (2010). In particular, the authors implement in their SPH equations an AC term similar to that of Eq. (18), but with the AC parameter set to unity, and with a signal velocity given by Eq. (24). The authors performed a systematic analysis of the capabilities of SPH to capture the KH instability using different SPH formulations, kernels, numerical resolutions, and KH timescales. This choice allows us to assess the code capabilities by constructing a large suite of simulations whose numerical results can be contrasted against the corresponding ones discussed by Valcke et al. (2010).

3.1.1. Initial conditions set-up

The problem domain consists of a periodic box of unit length with Cartesian coordinates x ∈  { 0,1 } , y ∈  { 0,1 } . Within this domain, there is a fluid with adiabatic index γ = 5/3, which satisfies the conditions $\mathit{\rho ,} \mathit{T,}{\mathit{v}}_{\mathit{x}}\mathrm{=}\{\begin{array}{c}\\ {\mathit{\rho }}_{\mathrm{1}}\mathit{,}{\mathit{T}}_{\mathrm{1}}\mathit{,}{\mathit{v}}_{\mathrm{1}}& \mathrm{|}\mathit{y}\mathrm{-}\mathrm{0.5}\mathrm{|}\mathrm{\le }\mathrm{0.25}\\ {\mathit{\rho }}_{\mathrm{2}}\mathit{,}{\mathit{T}}_{\mathrm{2}}\mathit{,}{\mathit{v}}_{\mathrm{2}}& \mathrm{|}\mathit{y}\mathrm{-}\mathrm{0.5}\mathrm{|}\mathit{>}\mathrm{0.25.}\end{array}$(29)As in Valcke et al. (2010), we choose here ρ₁ = 10 and ρ₂ = 1, with the index 1 referring to the high density layer. This choice is motivated by the finding that the difficulties of SPH in reproducing KH instabilities increase as the density contrast between the two fluid layers gets higher. The two layers are in pressure equilibrium with P₁ = P₂ = 10, so that the sound velocities in the two layers are ${\mathit{c}}_{\mathrm{1}}\mathrm{=}\sqrt{\mathit{\gamma }{\mathit{P}}_{\mathrm{1}}\mathit{/}{\mathit{\rho }}_{\mathrm{1}}}\mathrm{=}\mathrm{1.29}$ and ${\mathit{c}}_{\mathrm{2}}\mathrm{=}\sqrt{\mathit{\gamma }{\mathit{P}}_{\mathrm{2}}\mathit{/}{\mathit{\rho }}_{\mathrm{2}}}\mathrm{=}\mathrm{4.08}$, respectively. The layers slide against each other with opposite shearing velocities v₁ =  −v₂ and to seed the KH instability a small single-mode velocity perturbation is imposed along the y-direction ${\mathit{v}}_{\mathit{y}}\mathrm{=}{\mathit{w}}_{\mathrm{0}}\sin \mathrm{(}\mathrm{2}\mathit{\pi x}\mathit{/}{\mathit{\lambda }}^{\mathrm{)}}\mathit{,}$(30)where w₀ = 0.025 and λ = 1/6. To restrict the perturbation to spatial regions in the proximity of the interfaces, the perturbation (30) is applied only if  | y − σ |  < 0.025, where σ takes the values of 0.25 and 0.75, respectively. With this choice of parameters, the Mach number of the high-density layer is M ≃ v₁/c₁ ≃ v₁/1.29 and the KH timescale is τ_KH ≃ 0.29/v₁. As in Valcke et al. (2010), we consider KH simulations with a range of five different Mach numbers; Table 1 reports the values of M together with the respective values of v₁ and τ_KH.

To implement the density set-up (29), a two-dimensional lattice of equal-mass particles is placed inside the simulation box. We adopt here an isotropic hexagonal-close-packed (HCP) configuration for the particle coordinates instead of the more commonly employed Cartesian grid. The advantage of this configuration is that, for a given number of neighbors, it gives a more robust density estimate owing to its symmetry properties. To construct the initial density configuration (29), the lattice spacing of the particles is varied until the SPH density estimate (1) satisfies the required values within a certain tolerance criterion ( ≲ 1%) for the relative density error. The simulations were run using a total number of N = 512² particles and the initial particle specific energies were assigned after the density calculation so as to satisfy pressure equilibrium. This particle number is larger than that used in the runs of Valcke et al. (2010) ( ≃ 200 K), but guarantees a density setup with the specified tolerance criterion.

As noticed by Valcke et al. (2010) and other authors, SPH is a numerical method that can only represent smoothed quantities, hence applying it to hydrodynamic problems where strong density gradients are present can lead to inconsistencies This is, in fact, the situation encountered by standard SPH in the treatment of KH instabilities, where the lack of entropy mixing induces an artificial pressure discontinuity at fluid interfaces with a jump in density.

Motivated by these difficulties, we consider here a set of runs in which the density discontinuity at the interfaces is replaced by a smooth transition. To allow us to make a consistent comparison with the corresponding runs of Valcke et al. (2010), we adopt the same smoothing profile $\mathit{\rho }\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\mathit{D}\mathrm{\pm }\mathit{A}\mathrm{atan}\left[\mathit{B}\left({\mathit{y}}^{\mathrm{\prime }}\mathrm{+}\mathit{C}\right)\right]\mathit{,}$(31)where the coefficients are given by $\begin{array}{ccc}\mathit{A}& \mathrm{=}& \\ \mathit{B}& \mathrm{=}& \\ \mathit{C}& \mathrm{=}& \\ \mathit{D}& \mathrm{=}& \end{array}$and y′ = y − σ + δ/2. In Eq. (31), the sign in front of the coefficient A refers to σ = 0.25, 0.75, respectively. The parameters β and δ determine the thickness of the density transition and take the values β = 10, δ = 0.5/0.75.

The following procedure is adopted to construct a lattice configuration that satisfies the density profile given by Eq. (31). An HCP lattice of particle is first constructed in the domain 0 ≤ x < 1 and 0 ≤ y ≤ 0.5 with a spacing a_down such that the SPH density is ρ = ρ₂. Proceeding upward from y = 0, the lattice spacing is progressively reduced according to a = a_down [ρ₂/ρ(y)] ^1/2 until y = 0.25 and N_down particles are left. The same procedure is applied to a high-density lattice that in the same domain satisfies the condition ρ = ρ₁: starting from y = 0.5 and proceeding downward, the lattice spacing is increased so that a = a_up [ρ₁/ρ(y)] ^1/2 until y = 0.25 and N_up particles remain of the original high-density lattice. The whole procedure is numerically iterated until the numbers N_down and N_up satisfy the conditions N_down + N_up = N/2 and N_up/N_down = ρ₁/ρ₂. The lattice is then replicated in the top half of the box (y → 1 − y) so that the initial conditions are fully symmetric around y = 0.5. In the following, simulations with initial conditions for which the particle positions are arranged in a unsmoothed HCP lattice are denoted with the suffix SPH, whereas all of the others adopt the smoothed density profile constructed according to the procedure described here.

Another issue that is found to have a significant impact on the ability of standard SPH to properly capture KH instabilities is the choice of the kernel function. According to Read et al. (2010), the accuracy of the momentum equation for particle i is governed by the leading error ${{E}}_{{i}}\mathrm{0}\mathrm{=}\sum _{\mathit{j}}\frac{{\mathit{m}}_{\mathit{j}}}{{\mathit{\rho }}_{\mathit{j}}}\left[\frac{{\mathit{\rho }}_{\mathit{i}}}{{\mathit{\rho }}_{\mathit{j}}}\mathrm{+}\frac{{\mathit{\rho }}_{\mathit{j}}}{{\mathit{\rho }}_{\mathit{i}}}\right]{\mathit{h}}_{\mathit{i}}\mathbf{\nabla }\mathit{i}\mathit{W̅}\mathit{ij}\mathit{,}$(36)which vanishes in the continuum limit. However, this does not hold for a finite number of irregularly distributed particles or, more specifically, in the proximity of a contact discontinuity where a density step is present. A possible solution consists of increasing the number of neighbors so as to improve kernel sampling, although this approach presents some difficulties when the commonly employed M₄ or cubic spline (CS) SPH kernel is used. This occurs because for a large number of neighbors the CS kernel is subject to the so-called clumping instability, in which pairs of particles with interparticle distance q = r/h < 2/3 remain close together because the CS kernel gradient tends to zero below this threshold distance. A stability analysis (Morris 1996; Børve et al. 2004; Price 2005a; Read et al. 2010) shows that for the CS kernel, a Cartesian lattice of particles is unstable for η ≃ 1.5 or N_sph ~ 28,110 when D = 2 and D = 3, respectively. The clumping degrades the spatial resolution because it reduces the effective number of neighbors with which integrals are sampled, thus still has large E₀ errors even when the resolution is increased. To overcome this problem, one can modify the kernel shape in order to have a non-zero kernel derivative at the origin. However, for a fixed number of neighbors, kernels of this kind have the drawback of providing a less accurate density estimate than that given by the CS kernel since near the origin the kernels have a steeper profile.

As an alternative to the CS kernel, we consider here the linear quartic (LIQ) kernel, introduced by Valcke et al. (2010), which is a quartic polynomial for q_s ≤ q ≤ ζ = 2 and is linear for 0 ≤ q < q_s. The choice of parameter q_s determines the quality of density estimates. From a set of 2D Sod shock tube runs, Valcke et al. (2010) recommend q_s = 0.6, which is the value adopted here. For the functional form and normalization of the kernel, we refer to Valcke et al. (2010).

Another kernel that has been introduced for the purpose of avoiding particle clumping is the core triangle (CRT) kernel (Read et al. 2010), which has a constant first derivative for 0 ≤ q < α and is similar to the CS kernel for α ≤ q ≤ ζ = 2. This kernel is of the form $\mathit{w}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\frac{\mathit{\sigma }}{{\mathit{h}}^{\mathit{D}}}\{\begin{array}{c}\\ \mathrm{(}\mathrm{-}\mathrm{3}\mathit{\alpha }\mathrm{+}\frac{\mathrm{9}}{\mathrm{4}}{\mathit{\alpha }}^{\mathrm{2}}\mathrm{)}\mathit{q}\mathrm{+}\mathrm{1}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\alpha }}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\alpha }}^{\mathrm{3}}\mathit{,}& \mathrm{0}\mathrm{\le }\mathit{q}\mathit{<}\mathit{\alpha ,}\\ \mathrm{1}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{q}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{4}}{\mathit{q}}^{\mathrm{3}}\mathit{,}& \mathit{\alpha }\mathrm{\le }\mathit{q}\mathit{<}\mathrm{1}\mathit{,}\\ \frac{\mathrm{1}}{\mathrm{4}}{\left(\mathrm{2}\mathrm{-}\mathit{q}\right)}^{\mathrm{3}}\mathit{,}& \mathrm{1}\mathrm{\le }\mathit{q}\mathit{<}\mathrm{2}\mathit{,}\\ \mathrm{0.}& \mathit{q}\mathrm{\ge }\mathrm{2}\mathit{,}\end{array}$(37)where $\mathrm{1}\mathit{/}\mathit{\sigma }\mathrm{=}\mathrm{2}\mathit{\pi }{}^{\mathrm{(}}\frac{\mathrm{7}}{\mathrm{20}}\mathrm{+}\frac{{\mathit{\alpha }}^{\mathrm{4}}}{\mathrm{8}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{20}}{\mathit{\alpha }}^{\mathrm{5}}{}^{\mathrm{)}}\mathit{,} \mathrm{4}\mathit{\pi }\left(\frac{\mathrm{1}}{\mathrm{4}}\mathrm{+}\frac{{\mathit{\alpha }}^{\mathrm{5}}}{\mathrm{20}}\mathrm{-}\frac{{\mathit{\alpha }}^{\mathrm{6}}}{\mathrm{16}}\right)$ for D = 2 and D = 3, respectively. The value of α is fixed by the condition of continuity for the second derivative, giving α = 2/3. For the grid of initial conditions previously described, the sample of KH simulations is then constructed by performing SPH runs with the same initial conditions but using different kernels. We consider the CS kernel, together with the LIQ and CRT kernels. For all of these runs, the number of neighbors is N_sph = 32   (η ≃ 1.5).

However, another solution for reducing sampling errors consists of keeping a B-spline kernel function but increasing its order (Price 2012a, Sect. 5.4). This approach has the advantage of retaining the bell shape of the kernel, a feature that provides good density estimates (Fulk & Quinn 1996). After the M₄ (CS) kernel, at the next order is the M₅ or quartic kernel $\mathit{w}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\frac{\mathit{\sigma }}{{\mathit{h}}^{\mathit{D}}}\{\begin{array}{c}\\ {\mathrm{(}\frac{\mathrm{5}}{\mathrm{2}}\mathrm{-}{\mathit{q}}^{\mathrm{)}}}^{\mathrm{4}}\mathrm{-}\mathrm{5}{\mathrm{(}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}{\mathit{q}}^{\mathrm{)}}}^{\mathrm{4}}\mathrm{+}\mathrm{10}{\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\mathit{q}}^{\mathrm{)}}}^{\mathrm{4}}\mathit{,}& \mathrm{0}\mathrm{\le }\mathit{q}\mathit{<}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}\\ {\mathrm{(}\frac{\mathrm{5}}{\mathrm{2}}\mathrm{-}{\mathit{q}}^{\mathrm{)}}}^{\mathrm{4}}\mathrm{-}\mathrm{5}{\mathrm{(}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}{\mathit{q}}^{\mathrm{)}}}^{\mathrm{4}}\mathit{,}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathit{q}\mathit{<}\frac{\mathrm{3}}{\mathrm{2}}\mathit{,}\\ {\mathrm{(}\frac{\mathrm{5}}{\mathrm{2}}\mathrm{-}{\mathit{q}}^{\mathrm{)}}}^{\mathrm{4}}\mathit{,}& \frac{\mathrm{3}}{\mathrm{2}}\mathrm{\le }\mathit{q}\mathit{<}\frac{\mathrm{5}}{\mathrm{2}}\mathit{,}\\ \mathrm{0.}& \mathit{q}\mathrm{\ge }\frac{\mathrm{5}}{\mathrm{2}}\mathit{,}\end{array}$(38)which is truncated to zero for ζ ≥ 2.5 and has σ = 96/1199π,1/20π for D = 2 and D = 3, respectively.

An additional set of SPH simulations was then performed using the M₅ spline as the kernel. For consistency with the other runs, we kept the same ratio of smoothing length to particle spacing (η ~ 1.5), so that for these runs the chosen number of neighbors is N_sph = 50. A non-trivial issue concerns the role of pairing instability for this class of kernels. Because the gradient of the M₅ kernel still goes to zero as q → 0, one would expect the instability still to occur for η ~ 1.5. Nonetheless, it will be seen that KH simulations with the M₅ kernel do not exhibit particle clumping, in contrast to corresponding simulations performed with the CS kernel. This suggests that the stability properties of the M₅ kernel are better than those of the CS one; this topic is discussed in a subsequent dedicated section (Sect. 3.1.3).

Table 2 summarizes the main simulation parameters used in the KH tests. For comparison purposes, a set of mirror runs was carried out for the LIQ simulations in which Eq. (25), for the signal velocity, was replaced in Eq. (18) by Eq. (24), which is based on pressure discontinuities.

3.1.2. Results

Figure 1 shows the density maps at t = τ_KH for some of the KH simulations listed in Table 2. The panels can be compared directly with the corresponding ones in Fig. 10 of Valcke et al. (2010). A visual inspection reveals that the expected features of the KH runs are correctly reproduced and of a general consistency between the results produced by the two codes. In particular, the LP simulations, which use the pressure-based AC signal velocities of Eq. (24), appear to produce results almost identical to the corresponding LIQ ones. Thus, showing how, for the KH tests examined here, the use of the two signal velocities can be considered equivalent, within the spanned range of Mach numbers.

We note, however, the absence of KH rolls for the LP run with M = 0.4, which contrasts with the same simulation in Valcke et al. (2010) where the rolls have already been developed. Given the general agreement between the two codes, it is hard to ascertain the origin of the discrepancy for this specific run. However, the panels of Fig. 1 show that the AC-SPH formulation, as discussed in Sect. 4.2, is still unable to properly capture the development of KH instabilities for very subsonic shear flows. This suggests how small differences in the time integration procedure of the two codes can manifest themselves in the long-term evolution of cold flows.

In Fig. 2, the results for the same set of KH simulations of Table 1 are shown, but with the use of kernels LIQ, CRT, and M₅. An important feature is now the appearance in the CRT and M₅ runs of KH rolls in the M = 0.4 case. This improvement in the description of KH instabilities suggests that the integration errors that are present with the LIQ kernel are now absent or smaller in the CRT and M₅ runs. However, in the M = 0.2 case, the kernels are still unable to capture the development of the KH instability. For this test case, a high-resolution run (N = 1024²) using the M₅ kernel (not shown here) still reveals the absence of any KH roll at t = τ_KH, thus showing that in SPH the problem of an accurate description of KH instabilities in the very subsonic regimes is neither a resolution issue nor due to the AC implementation. We refer to McNally et al. (2012) for a recent discussion of this topic.

To test the effectiveness of the switch (13) in ensuring a fast response of the ${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$ parameters to the presence of thermal energy discontinuities, Fig. 3 renders the plots of the AC parameters that are shown for the test runs of Fig. 2. The maps are color-coded according to the range of values of ${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$. The highest values lie in the range  ~0.7−0.8 and are confined to a narrow strip around the interface layers, with the floor value  ~0 as the background value away from the discontinuities. We note that in general, and in particular for the M = 1 test case, the maximum values of the α^C for the M₅ runs are below those of the other runs. This result is due to the faster growth rate of the instability, which is ensured by an improved kernel accuracy in interpolating the variables.

The long-term evolution of the KH tests is shown at t = 2 in Fig. 4. The overall features of the density plots for different runs are similar to their counterparts displayed in Fig. 11 of Valcke et al. (2010). We note the tendency for the M₅ runs to show small-scale features at the layer contacts.

The aim of this section is to test the consistency of the present AC implementation by comparing results extracted from a suite of AC-SPH simulations of the 2D KH instability problem, with those of similar runs (Valcke et al. 2010). The results of the KH tests indicate a code behavior that is in accord with that expected and the simulations of Valcke et al. (2010).

However, Valcke et al. (2010) argue that the absence of the AC term is not the main reason for the SPH failure to develop KH instabilities; however the presence of AC is necessary for the long-term evolution of the instabilities, but this difficulty of SPH has two distinct reasons. The first is a general problem of consistency of the initial condition set-up, as SPH by definition can only deal with smoothed quantities and therefore its application to problems where sharp discontinuities are present leads to inconsistencies. This is part of the more general problem (Robertson et al. 2010) of properly smoothing in numerical simulations of hydrodynamic test problems the discontinuities initially present at the interfaces, in order to achieve convergence in the solution. This aspect of the KH test problem can be cured by properly stretching the initial particle lattice so as to introduce a smooth density transition at the interfaces. Results suggest that there is a significant improvement in the capability of SPH to properly capture the correct growth rate of the KH instability.

The other problem that causes the poor performance of SPH when handling KH instabilities is the leading error in the momentum equation due to the incomplete kernel sampling. This error can be reduced by removing particle clumping, which depends on the kernel stability properties. In fact, the kernels LIQ and CRT were introduced (Read et al. 2010; Valcke et al. 2010) with the aim of removing the clumping instability. Since the results of the simulations suggest performances for the M₅ kernel that are quite similar to those achieved by these kernels, it is therefore interesting quantify the relative performances of these kernels in a better way.

Following Junk et al. (2010), we then measure the growth rate of the KH instability for some of the runs and compare it with the linear theory growth rate expectation  ∝ e^t/τ_KH. This is achieved, using Fourier transforms, by measuring the time-evolution of the λ = 1/6 growing mode of the v_y velocity perturbation component. The details of the procedure are given in Junk et al. (2010) and so are not reported here. For the sake of clarity, the results of the LIQ runs are not displayed since they are quite similar to those of the CRT simulations. Moreover, we only display the growth rates for three different KH test cases, those with Mach number M = 0.2, 0.4, and M = 1, the results of the others being intermediate between these.

A number of distinct features are apparent in the panels of Fig. 5. The first is that simulations with smoothed IC (RHO) perform systematically better than the unsmoothed ones (SPH). The second is that only for M = 1 are the simulations with kernels CRT and M₅ able to correctly recover the expected growth rate. Finally, this capability progressively declines as one considers lower Mach numbers. While this behavior is in accordance with the visual impression of the maps previously displayed, and confirms that the present SPH implementation still has problems in describing instabilities in subsonic flows, we note that the performances of the M₅ runs are similar to those obtained using the CRT kernel. This behavior is particularly interesting, since the simulations were performed by setting the ratio of the smoothing lengths to particle spacing to η ≃ 1.5 so that the pairing instability, which is present in runs that use the CS kernel, should be present in the M₅ simulations as well.

How the clumping instability affects sampling errors can be assessed by computing the particle errors ${E}\begin{array}{}\mathrm{0}\\ \mathit{i}\end{array}$, according to Eq. (36). The distribution, throughout the simulation domain of the ${E}\begin{array}{}\mathrm{0}\\ \mathit{i}\end{array}$ errors versus y, is shown in the top panels of Fig. 6 for M = 1 simulations performed using different kernels. The solid lines represent the mean of the binned distributions. Similar plots were produced by Valcke et al. (2010) and their Fig. 3 can be used for comparative purposes.

As expected, the largest E⁰ errors are present in the SPH simulation, whereas better results are obtained by using the LIQ and CRT kernels, as indicated by the error distribution in their respective panels. This is a result of the absence of particle clumping for these simulations, owing to the specific stability properties of these kernels (Read et al. 2010; Valcke et al. 2010).

A striking feature is given by the E⁰ error distribution of the simulation performed using the M₅ kernel, which shows how the magnitude of the E⁰ errors are even smaller than those of the two runs CRT and LIQ for this kernel. This result is in accordance with what has been found previously by analyzing the growth rates, suggesting that, since all of the simulations were performed using the same η, the stability properties of the M₅ kernel are better than those of the CS one.

To further investigate this point, the bottom panels of Fig. 6 show the “nearest neighbor” map of the corresponding top panels. This is defined by interpolating the quantity ${\mathit{q}}_{\mathit{i}}^{\mathrm{2}}$ at the grid points, according to the SPH prescription, where for the k − th neighbor ${\mathit{q}}_{\mathit{i}}^{\mathit{k}}\mathrm{=}\mathrm{|}{x}{}_{\mathit{i}}\mathrm{-}{x}{}_{\mathit{k}}\mathrm{|}\mathit{/}{\mathit{h}}_{\mathit{i}}$ and the neighbors have been sorted according to their distance from the particle i itself.

The map of the M₅ kernel indicates a distribution of the second nearest-neighbor distribution ${\mathit{q}}_{\mathit{i}}^{\mathrm{2}}$ that is quite similar to that displayed by the CRT and LIQ kernels. The only difference occurs at the layer interfaces, where the distribution of the quantities ${\mathit{q}}_{\mathit{i}}^{\mathrm{2}}$ for the M₅ kernel is slightly shifted towards smaller values than for the other kernels. Conversely, we note that because of the pairing instability, for the CS kernel the distribution of the neighboring distances is flipped with respect to that of the other kernels.

The results of Fig. 6, however, clearly show the absence of clumping instability for the M₅ kernel. In the next section, we investigate, in more detail, the stability properties of this kernel.

3.1.3. Stability issues

The stability properties of SPH have been investigated by a number of authors (Swegle & Hicks 1995; Morris 1996; Monaghan 2000; Børve et al. 2004; Read et al. 2010; Dehnen & Aly 2012). The instabilities are studied analytically by analyzing the dispersion relation for sound waves of small amplitude, propagating in a uniform medium. We now derive the dispersion relation for the SPH equation of motion in a manner similar to the analysis of previous authors. A uniform lattice of equal-mass particles of mass m, density ρ₀, pressure P₀, and sound speed ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathit{\gamma }{\mathit{P}}_{\mathrm{0}}\mathit{/}{\mathit{\rho }}_{\mathrm{0}}$ is perturbed by a small wave ${x}\mathit{i}\mathrm{=}{x}\begin{array}{}\mathrm{0}\\ \mathit{i}\end{array}\mathrm{+}{a}\exp \mathrm{[}{k}\mathrm{·}{x}\begin{array}{}\mathrm{0}\\ \mathit{i}\end{array}\mathrm{-}\mathit{\omega }{\mathit{t}}^{\mathrm{]}}\mathit{,}$(39)where a is the perturbation, k is the wavevector, and ${x}\begin{array}{}\mathrm{0}\\ \mathit{i}\end{array}$ are the unperturbed particle positions. By linearizing the equation of motion for the perturbation, one has the dispersion relation $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}{\mathit{a}}_{\mathit{\mu }}& \mathrm{=}& [\frac{\mathrm{2}\mathit{m}{\mathit{P}}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{2}}}\sum _{\mathit{j}}{\mathit{H}}_{\mathit{\mu \nu }}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathrm{(}{k}\mathrm{·}{x}{\begin{array}{}\mathrm{0}\\ \mathit{ij}\end{array}}^{\mathrm{)}}{\mathit{a}}_{\mathit{\nu }}\\ & & \mathrm{+}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{2}\mathrm{)}\frac{{\mathit{m}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{3}}}\mathrm{(}{a}\mathrm{·}{q}{\mathit{i}}^{\mathrm{)}}{\mathit{q}}_{\mathit{\mu }}^{\mathit{i}}]\mathit{,}\end{array}$where the summations are over the neighbors j of particle i, the vector q_i is defined as ${q}\mathit{i}\mathrm{=}\sum _{\mathit{j}}\sin {k}\mathrm{·}{x}\begin{array}{}\mathrm{0}\\ \mathit{ij}\end{array}\mathbf{\nabla }\mathit{i}{\mathit{W}}_{\mathit{ij}}\mathit{,}$(40)and H(W) is the Hessian of the kernel ${\mathit{H}}_{\mathit{\mu \nu }}\mathrm{=}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{W}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{\partial }{\mathit{x}}_{\mathit{\mu }}\mathit{\partial }{\mathit{x}}_{\mathit{\nu }}}\mathrm{·}$(41)For a given smoothing length h and wavevector k, the stable solutions of Eq. (40) are defined by the condition ω² ≥ 0. Solutions for which ω² < 0 represent exponentially growing or decaying perturbations. Moreover, it is useful to define a numerical sound speed ${\mathit{C}}_{\mathrm{num}}^{\mathrm{2}}\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}\mathit{/}{\mathit{k}}^{\mathrm{2}}$ and a scaled numerical sound speed ${\mathit{c}}_{\mathrm{num}}^{\mathrm{2}}\mathrm{=}{\mathit{C}}_{\mathrm{num}}^{\mathrm{2}}\mathit{/}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$. The condition ${\mathit{c}}_{\mathrm{num}}^{\mathrm{2}}\mathrm{=}\mathrm{1}$ should clearly be satisfied for the sound speed to be correctly modeled.

We now examine the stability properties of the CS and M₅ kernels. For simplicity, we consider a plane wave propagating along the x-axis, k = k(1,0,0), and assume γ = 5/3. The bottom panels of Fig. 7 show, for the longitudinal and the two transverse waves of the perturbation, the instability regions of the CS kernel, which are denoted by the gray areas and represent the solutions to Eq. (40) in the domain (h,k), for which ω² < 0. Similarly, the top panels are for the M₅ kernel.

The longitudinal wave perturbation is responsible for the clumping instability, whereas the traverse waves give rise to the so-called banding instability (Read et al. 2010). Unlike the clumping instability, the latter is relatively unimportant in causing sampling errors (Read et al. 2010) and is not considered further here. A comparison of the two stability plots for the longitudinal wave solution clearly shows that the stability properties of the M₅ kernel are much better than those of the CS one.

This is part of a more general result that was previously recognized by Morris (1996): the stability properties of SPH improve when higher-order spline kernels are used in place of the CS kernel. This occurs, basically, because the higher the order of the spline, the better it approximates a Gaussian kernel. In Eq. (40), one can see that the numerical sound speed ${\mathit{c}}_{\mathrm{num}}^{\mathrm{2}}$ depends on the first and second derivatives of kernel W. Ideally, one should have ${\mathit{c}}_{\mathrm{num}}^{\mathrm{2}}\mathrm{\simeq }\mathrm{1}$ to accurately describe the sound wave propagation and this condition is fulfilled when smoother kernels are used, so as to keep the numerical dependence of ${\mathit{c}}_{\mathrm{num}}^{\mathrm{2}}$ as weak as possible.

We note however that here, in contrast to the stability properties of the CRT kernel, the clumping instability is not entirely suppressed but instead it is present whenever η ≳ 2.5.

Finally, it must be stressed that the paring instability that occurs in the hydrodynamic tests described here, is unlikely to have a significant impact on many astrophysical problems of interest, where very cold flows are generally absent. Rasio & Lombardi (1999) estimate for instance, from SPH simulations of a stationary fluid, that lattice effects become important when velocity dispersions are below  ~3%−4% of the sound speed.

The results of this section and the previous one, therefore, suggest that to avoid the pairing instability, the M₅ kernel can be considered a viable alternative to the use of the otherwise steeper CRT and LIQ kernels, provided that the parameter η is consistently rescaled. In the next section, we then discuss the related performances of these kernels in a test simulation.

3.2. Sod’s shock tube

A classic test used to investigate the hydrodynamic capabilities of SPH codes is the Sod shock-tube problem (Hernquist & Katz 1989; Wadsley et al. 2004; Springel 2005; Price 2008; Tasker et al. 2008; Rosswog 2009). This test consists in a fluid, initially at rest, in which a membrane located at x = 0 separates the fluid on the right, of high density and pressure, from the fluid on the left, of relatively lower density and pressure. The membrane is removed at t = 0 and a shock wave develops propagating toward the left, followed by a contact discontinuity and a rarefaction wave propagating to the right.

A well-known problem with standard SPH codes in reproducing the analytic solution of the shock-tube problem is the presence of a pressure discontinuity that arises across the propagating contact discontinuity. Simulations incorporating an artificial thermal conductivity term in the SPH equations (Price 2008; Rosswog 2009; Price 2012a) show shock profiles in which density and thermal energy are resolved across the discontinuity, hence producing a continuous pressure profile. However, in these runs the AC formulation adopts the AC signal velocity (24) where jumps in thermal energy are smoothed in accordance with the presence of pressure discontinuities. This is in contrast to the AC signal velocity (25) employed here, for which in the absence of shear flows, contact discontinuities are unaffected and therefore cannot be used to remove the blip seen at the contact discontinuity in the pressure profile of the shock tube SPH runs.

The application of the AC-SPH code to this test is nonetheless of interest because it can still be used to validate the code’s performances. In particular, we consider a 3D setup of the shock-tube test and with these initial conditions we construct a set of AC-SPH simulations performed with different kernels and neighbor numbers. Shock tube profiles extracted from these simulations are compared with the aim of assessing the goodness of different kernels in reproducing, for a given test problem, profiles of known analytic solutions.

The initial condition setup consists of an ideal fluid with γ = 5/3, initially at rest at t = 0. An interface set at the origin separates the fluid on the right of density and pressure (ρ₁,P₁) = (4,1) from the fluid on the left with (ρ₂,P₂) = (1,0.1795). To construct these initial conditions, a cubic box of side-length unity was filled with 10⁶ equal mass particles, so that 800   000 were placed in the right-half of the cube and 200   000 in the left-half. The particles were extracted from two independent, uniform glass-like distributions contained in a unit box consisting of 1.6    ×    10⁶ and 4 × 10⁵ particles, respectively. This initial condition setup is the same as that previously implemented in Sect. 5.1 of V11, to test the time-dependent AV scheme described in Sect. 2.2.

For the same initial setup, to investigate the performances of different kernels in reproducing the analytic profiles of the shock-tube problem, we perform runs with different kernels and neighbor numbers. The kernels considered are CS, LIQ, CRT, and M5. The number of neighbors varies as a power of two between 32 and 128. The SPH runs were realized by imposing periodic boundary conditions along the y and z axes and the results were examined at t = 0.12. We show the results obtained using the standard AV formulation, the results of the other runs being unimportant from the viewpoint of the kernel performances.

Some of the simulation profiles extracted from the 3D SPH runs of the shock-tube test are shown in Fig. 8. A striking feature is the large scatter between the pressure profiles of the runs performed using different kernels or neighbor numbers. The same behavior is present for the thermal energy profiles, whilst very similar density profiles are exhibited by the same runs.

There are several conclusions that can be drawn from the results of Fig. 8. The performances of the M₅-64 run are quite similar to those of CS-32, although for the former simulation a closer inspection reveals a slightly better treatment of the thermal energy spike at the contact discontinuity, the spike being due to the initial condition set-up.

The worst results are obtained by the LIQ simulations when using N_sph = 32 or N_sph = 64 neighbors. The profiles of the LIQ-128 run are not shown here to avoid overcrowding in the plots, and are quite similar to those of the simulation CRT-128. This clearly demonstrates the need for this class of kernels to use a large (say  ≳ 128) number of neighbors, so as to compensate for the density underestimate due to the steeper profiles introduced to avoid the pairing instability.

To better quantify this density bias, Table 3 reports, for different kernels and neighbor numbers, the mean SPH density estimated from a glass-like configuration of N = 10⁶ particles, in a unit periodic box of total mass one. The results illustrate how the density estimate of the M₅ kernel with 64 neighbors is comparable to the CS one obtained using 32 neighbors, the value of η being the same (η ≃ 1). Conversely, for the LIQ and CRT kernels, only when N_sph = 128 does the mean density approach the M₅-64 estimate.

The density values given in Table 3 also help us to explain the rather poor performances of the LIQ kernel when using 32 neighbors. From Fig. 8, it can be seen that for the corresponding run the relative pressure error is ε_P ≃ 10% in the unperturbed right zone of the cube. This error is already present in the pressure profile at t = 0 and originates from the density error due to the adopted kernel and neighbor number, together with the use of an entropy-conserving code to perform the simulations. Particle entropies are initially assigned according to the initial conditions so that, for particles satisfying x_i > 0, ${\mathit{A}}_{\mathit{i}}\mathrm{\equiv }{\mathit{A}}_{\mathrm{1}}\mathrm{=}{\mathit{P}}_{\mathrm{1}}\mathit{/}{\mathit{\rho }}_{\mathrm{1}}^{\mathit{\gamma }}\mathrm{\simeq }\mathrm{0.1}$. During the integration, particle pressures are calculated according to ${\mathit{P}}_{\mathit{i}}\mathrm{=}{\mathit{A}}_{\mathit{i}}{\mathit{\rho }}_{\mathit{i}}^{\mathit{\gamma }}$, and for x > 0 a relative pressure error ε_P ≃ γε_ρ ≃ 10% is present when ε_ρ ≃ 6%.

Taken at face value, the results of Table 3 demonstrate that in 3D simulations a conservative lower limit for the kernels with a modified shape should be to assume at least N_sph ≳ 200 neighbors. In fact, Read & Hayfield (2012) present a new formulation of SPH where they adopted, as reference, the so-called high-order core triangle (HOCT, Price 2005a). The profile of this kernel is a generalization of that of the CRT kernel and the authors assumed N_sph = 442 as the reference value for their scheme.

However, the results of the previous sections suggest that the stability properties of the M₅ kernel can be profitably used, with an appropriate choice of the η parameter, to avoid the pairing instability when dealing with tests of hydrodynamic instabilities in which cold flows are present.

Finally, the density estimates of Table 3 suggest that great care should be taken when deciding the goodness of a particular kernel based on its relative performance in terms of density estimates. The results of the 3D Sod shock tube clearly indicate how there could be a large difference between the simulation and the expected solution profile of some hydrodynamic variable, such as pressure or thermal energy, and, at the same time, a much smaller difference in the corresponding density profile.

3.3. Sedov blast wave

The Sedov blast-wave test is used to validate, in three dimensions, the code capability in the strong shock regime. The test consists of a certain amount of energy E being injected at t = 0 into a very small volume of an ambient medium of uniform density ρ. The spherically symmetric shock propagates outward from the initial volume and at time t the shock front is located at the radius (Sedov 1959) $\mathit{R}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\simeq }\mathit{\beta }{\mathrm{(}\mathit{E}{\mathit{t}}^{\mathrm{2}}\mathit{/}{\mathit{\rho }}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{5}}\mathit{,}$(42)where β ~ 1.15 for γ = 5/3.

Previous investigations (Rosswog & Price 2007; Merlin et al. 2010) showed that, owing to the large discontinuities initially present in the thermal energy, incorporating an artificial conduction term in the energy equation greatly improves the description of the shock front in the simulations. Without the presence of this term, the initially large discontinuity in the thermal energy soon gives rise to a disordered particle distribution that degrades the shock profile (Rosswog & Price 2007).

The initial setting of the test was realized as follows. A HCP lattice of N = 2 × 64³ equal mass particles was arranged in a cubic box of side length unity. The particle masses were chosen so as to give ρ = 1 and periodic boundary conditions were imposed. The nearest particle to the position  { 0.5,0.5,0.5 }  was chosen as the center particle. To consistently represent a point-like explosion with the given numerical resolution, the particles j comprised within the kernel radius ζh_i of the center particle i were given an initial thermal energy such that the total injected energy was E = 1. This blast wave energy was distributed among the neighboring particles not uniformly but with a weight proportional to W_ij.

Using this initial condition set-up, we ran three test simulations, which differed in the choice of the adopted kernel and neighbor number. For the CS kernel, we used N_sph = 64 neighbors. We also ran two other test cases, now using the CRT and M₅ kernel and N_sph = 128 neighbors. All of the simulations were performed using the implemented AC scheme and the time-dependent AV formulation with parameters given by  {α_min,α_max,l_d}  =  {0.1,1.5,0.2}.

Figure 9 shows the radial density profiles at t = 0.06 for the different test runs. The solid black line represents the expected analytic solution, with the shock front being located at r ≃ 0.37. Radial simulation profiles were obtained by averaging, for each radial bin r_k, SPH densities calculated from the particle distributions over a set of (θ,φ) = (20,20) grid points uniformly spaced in angular coordinates: these were located at the surface of a sphere with radius r_k. The radial spacing was not uniform but was chosen so as to guarantee an accurate sampling of density in the proximity of the shock front, with about  ~40 radial bins between r ~ 0.34 and r ~ 0.42.

All of the simulations are in fair agreement with the analytical solution and the differences between the simulation profiles are negligible. At the shock front, the profiles exhibit a density jump of about  ~2, whereas the analytical solution gives a compression factor of γ + 1/γ − 1 = 4. These results are in accordance with previous findings (Rosswog & Price 2007; Springel 2010b; Heß & Springel 2010) and indicate that for 3D SPH simulations of the Sedov-Taylor point explosion problem the simulation profiles converge to the analytical solution as the resolution is increased, with approximately  ~345³ particles (Rosswog & Price 2007) being required to fully resolve the shock front.

Figure 9 can also be used to verify the behavior of the individual timestep algorithm. For problems involving very strong shocks, as demonstrated by Saitoh & Makino (2009), individual particle timesteps must be properly restricted so as to avoid in the proximity of the shock front their failure to satisfy the local Courant condition, thus leading to inaccuracies in the integration. The SPH simulations using individual timesteps, but without an appropriate limiter, fail to predict the expected solution profile (see Fig. 3 of Saitoh & Makino 2009). The profiles of Fig. 9 demonstrate the correctness of the algorithm used to update the timesteps, although different from the one devised by Saitoh & Makino (2009).

The present scheme adopts individual particle timesteps Δt_i that can vary in power-of-two subdivisions of the largest allowed timestep Δt₀ ≥ Δt_i (Hernquist & Katz 1989). At each step, particles whose time bin is synchronized with the current time are defined as “active” and their hydrodynamic quantities, as well as their smoothing lengths and densities are consistently updated. However, non-active particles j that are neighbors of an active particle i are defined as “low-order active” particles, and their hydrodynamic variables, as well as their time-step constraints but not their forces, are recalculated. This criterion is applied regardless of the local shock conditions and no particular conditions are imposed on Δt_j that depend on Δt_i.

3.4. The blob test

The blob test is another hydrodynamic test where the results of standard SPH differ significantly from those produced by grid-based simulations (Agertz et al. 2007; Read et al. 2010; Cha et al. 2010; Heß & Springel 2010; Murante et al. 2011). The test consists of a gas cloud of radius R_cl placed in an external medium ten times hotter and less dense than the cloud, so as to satisfy the pressure equilibrium. A large enough wind velocity v_W is given to the hot low-density medium so that a strong shock wave strikes the cloud. The interaction of the cloud with the supersonic medium produces a number of effects that are of interest in an astrophysical context, such as gas stripping and fragmentation.

The blob is initially perturbed by the development of Richtmyer-Meshkov and Rayleigh-Taylor instabilities (Agertz et al. 2007). Afterwards, large-scale (~R_cl) KH instabilities are created at the cloud surface because of the shear flows caused by the supersonic wind. This non-linear phase is supposed to develop over a timescale ${\mathit{\tau }}_{\mathrm{cr}}\mathrm{~}\mathrm{2}{\mathit{R}}_{\mathrm{cl}}\sqrt{\mathit{\chi }}\mathit{/}{\mathit{v}}_{\mathrm{W}}$ (Agertz et al. 2007), where χ is the density contrast, after which cloud disruption will take place.

To investigate the capability of the AC-SPH code to properly follow the hydrodynamics of the blob test, we compare results extracted from our set of SPH simulations realized with the same initial conditions but with different numerical parameters. The numerical setup of this test is the same as in Read et al. (2010), to which we refer for more details. A spherical cloud of radius R_cl = 197 kpc is placed in a periodic, rectangular box of size  { L_x,L_y,L_z }  =  { 2,2,8 }  Mpc. The cloud has density ρ_cl = 4.74 × 10^-33   gr   cm^-3 = χρ_ext and temperature T_cl = 10⁶   K = T_ext/χ. The ambient medium has density and temperature so that χ = 10. The cloud is initially located at  { 1,1,1 }  Mpc and the ambient medium is given a wind velocity v_W = 1000   km   s^-1, so that for an adiabatic index γ = 5/3 its Mach number is M = 2.7. An HCP lattice of equal mass particles is constructed to satisfy the above density requirements so as to use for this version of the test a total number of N ~ 1.1 × 10⁶ particles, as in Heß & Springel (2010). Finally, a velocity perturbation is imposed at the cloud surface in order to trigger the development of an instability; the amplitude and modes are given in Appendix B of Read et al. (2010).

We compare the results of SPH simulations where three different spline kernels were used CS, CRT, and M₅. We use 128 neighbors for the simulations with the CS and CRT kernels and 220 for the run with the M₅ kernel, so that the ratio η is the same for all the runs (η ~ 1.5). For the CS kernel, we run a simulation where the AC term given by Eq. (18) is absent in the equation of thermal energy evolution (CS-128 NOAC), and three other simulations (CS-128, CRT-128, M5-220) in which the AC term (18) is incorporated in the energy evolution equation. In the following, the simulation AV parameters are set as follows:  { α_min,α_max,l_d }  =  { 0.1,1.5,0.2 } .

For several runs, Fig. 10 shows the gas density maps at three different times t = 1, 2, 3, the time being in units of τ_KH ~ 1.6τ_cr (Agertz et al. 2007). The maps have been evaluated according to the SPH prescription on a 2D grid of (200x800) points in the central yz plane located at x = L_x/2. The top panels of Fig. 10 are for the standard NOAC run CS-128. It can be seen that in this case, unlike the results of mesh-based simulations (Agertz et al. 2007), the cloud survives the impact of the supersonic wind striking its surface and there is no disruption. This is due to the absence of fluid instabilities developing at the cloud surface. If they were present, they would in turn lead to the stripping of material and the cloud break-up. This code behavior is similar to what was found in previous findings (Agertz et al. 2007; Read et al. 2010; Heß & Springel 2010).

In contrast, incorporating the AC diffusion term in the SPH equations leads to a significant improvement in the code capability to properly model the cloud evolution. This can be seen from the middle and bottom row panels of Fig. 10, in which the density maps are displayed for the CS-128 and M5-220 runs. The simulation performed using the CRT kernel is not shown since its performances are quite similar to those obtained using the CS kernel.

This behavior is consistent with the expectation that introducing AC removes the numerical effects that in SPH prevent the treatment of contact discontinuities when large density jumps are present, thus the inconsistencies that suppress the growth of the instabilities. The AC-SPH formulation presented here can therefore, at least qualitatively, correctly follow the time evolution of the cloud as in the other SPH schemes that have been proposed (Read et al. 2010; Heß & Springel 2010; Murante et al. 2011; Saitoh & Makino 2012).

To quantify the code performances in a more quantitative way, Fig. 11 (left panel) shows the mass loss of the cloud as a function of time for the different runs. We follow previous definitions (Agertz et al. 2007) and a gas particle is defined to still be a member of the cloud if at the considered epoch its gas and temperature fulfill the conditions ρ > 0.64ρ_cl and T < 0.9T_ext. The standard SPH results show a cloud that is not disrupted and whose mass, after an initial transient, stays constant at about half of the original value. A striking result is instead given by the runs employing the AC term. For these simulations, there is now a high mass-loss rate occurring at early times, followed by complete cloud disruption at t ≳ 3τ_KH.

We note that the mass loss rate of the AC-SPH runs does not depend in a significant way on the choice of the kernel. This suggests that the cloud disruption is driven by large-scale instabilities and is relatively insensitive to small-scale perturbations. Given the similarities displayed in Fig. 11, left panel, by the mass loss rates of the AC-SPH simulations employing different kernels, when referring to the left panel we adopt the term mass loss rate to generically indicate the behavior of these curves.

In the test considered here, the new AC scheme is clearly now capable of properly removing the surface effects, present across the contact discontinuity in the standard SPH version, which artificially suppress the growth of hydrodynamic instabilities. However, a comparison of the mass loss rate with the corresponding one produced using the new moving-mesh code AREPO in simulations of equivalent resolution (Springel 2011, Fig. 10) shows that for the AC runs presented here the mass depletion of the cloud occurs much faster. This discrepancy suggests that the processes of heat diffusion, which in the adopted numerical scheme are mediated via the AC parameter α^C, should be somehow constrained by a physically motivated mechanism, which was not considered in the discussion of Sect. 2.3. This mechanism should be introduced with the purpose of avoiding a heat transfer mechanism, as governed in the code by Eq. (18), which is overly diffusive.

To evaluate the relative effectiveness of the heat and momentum transport, in the theory of heat transfer the Prandtl number Pr is defined as the ratio of the kinematic viscosity ν to the thermal diffusion coefficient D, Pr = ν/D (Blundell & Blundell 2006). For gases, the transport coefficients for the transport of heat and momentum are nearly equal and the Prandtl number is of order unity, with $\mathit{Pr}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}$ when γ = 5/3 (Blundell & Blundell 2006). This suggests that a constraint on the particle AC parameter ${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$ can be obtained by setting $\mathit{Pr}\mathrm{\simeq }\frac{{\mathit{\nu }}_{\mathrm{AV}}}{{\mathit{D}}_{\mathrm{AC}}}\mathrm{~}\frac{\mathrm{1}}{\mathrm{5}}\frac{{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{AV}}{\mathit{v}}_{\mathit{ij}}^{\mathrm{AV}}}{{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{AC}}{\mathit{v}}_{\mathit{ij}}^{\mathrm{AC}}}\gtrsim \frac{\mathrm{2}}{\mathrm{3}}\mathit{,}$(43)where the equivalent kinematic viscosity coefficient ν_AV due to AV is given approximately by ${\mathit{\nu }}_{\mathit{i}}^{\mathrm{AV}}\mathrm{~}\frac{\mathrm{1}}{\mathrm{10}}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{AV}}{\mathit{v}}_{\mathit{ij}}^{\mathrm{AV}}{\mathit{r}}_{\mathit{ij}}$ (Lodato & Price 2010), and the numerical heat diffusion coefficient ${\mathit{D}}_{\mathit{i}}^{\mathrm{AC}}$ has been estimated according to Eq. (21). We note that the reverse of Eq. (43) is relatively unimportant, since the source term S^C_i is driven by a second order derivative, while for AV it is a first order derivative that determines S_i.

It can be seen from Eq. (43) that, for the AC signal velocity adopted here, the constraint on the numerical heat diffusion becomes important in the presence of supersonic flows. Moreover, the condition Pr ~ O(1) is valid only for gases, where the transport of momentum and energy occurs simultaneously. For liquids, according to the dominant mechanism of heat conduction, the value of the Prandtl number can vary by several orders of magnitude among different substances (Dimotakis 2005). The condition given by Eq. (43) should, therefore, be considered problem-dependent.

A rigorous procedure for deriving the upper limits on the set of parameters ${}^{\mathrm{\right\{}}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}^{\mathrm{\left\}}}$ at any given timestep would require us to first define for the particle i the numerical kinematic viscosity and heat diffusion coefficients $\{\begin{array}{c}\\ \mathrm{⟨}{{\mathit{\nu }}_{\mathit{i}}^{\mathrm{AV}}}^{\mathrm{⟩}}& \mathrm{=}& \left|{\mathrm{(}{\frac{\mathrm{d}{v}\mathit{i}}{\mathrm{d}\mathit{t}}}^{\mathrm{)}}}_{\mathrm{AV}}\right|\hspace{0.17em}/\hspace{0.17em}\left|{\mathrm{\nabla }}^{\mathrm{2}}{v}\mathit{i}\mathrm{+}\mathrm{2}\mathbf{\nabla }\mathrm{(}\mathbf{\nabla }\mathrm{·}{v}{\mathit{i}}^{\mathrm{)}}\right|\\ \mathrm{⟨}{{\mathit{D}}_{\mathit{i}}^{\mathrm{AC}}}^{\mathrm{⟩}}& \mathrm{=}& {\mathrm{(}{\frac{\mathrm{d}{\mathit{u}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}}^{\mathrm{)}}}_{\mathrm{AC}}/{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathit{i}}\mathit{,}\end{array}$(44)where the bulk viscosity ζ is 5/3 of the shear viscosity η ≡ νρ (Lodato & Price 2010) and it is understood that the SPH expressions should be used for the operators at the denominators.

The upper limits on the parameters ${}^{\mathrm{\right\{}}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}^{\mathrm{\left\}}}$ are then given by the conditions $\mathrm{⟨}{{\mathit{D}}_{\mathit{i}}^{\mathrm{AC}}}^{\mathrm{⟩}}\mathrm{\le }\frac{\mathrm{3}}{\mathrm{2}}\mathrm{⟨}{{\mathit{\nu }}_{\mathit{i}}^{\mathrm{AV}}}^{\mathrm{⟩}}\mathit{,}$(45)which must be simultaneously satisfied by all the parameters. In place of the procedure just described, we adopt here a simpler approach and using Eq. (43) we estimate the upper limits on ${}^{\mathrm{\right\{}}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}^{\mathrm{\left\}}}$ as ${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{AC}}\lesssim \frac{\mathrm{3}}{\mathrm{10}}{\mathit{\alpha }}_{\mathit{i}}^{\mathrm{AV}}\mathit{MA}{\mathit{X}}_{\mathit{j}\mathit{<}}\left\{{\left|{\mathit{v}}_{\mathit{ij}}^{\mathrm{AV}}\right|}^{\mathrm{􏼮}}\left|{\mathit{v}}_{\mathit{ij}}^{\mathrm{AC}}\right|\right\}\mathit{,}$(46)where the notation j <  means that the maximum is taken for all the neighbors j that satisfy the condition v_ij·r_ij < 0. This is to be consistent with the definition of AV, for which Π_ij is non-zero only for approaching particles. Clearly, this definition does not guarantee an accurate constraint of the ${\mathit{\alpha }}_{\mathit{i}}^{\mathrm{C}}$ parameter, but the use of the maximum among all the approaching pairs should provide a floor value for the effective constraint.

For the CS kernel, we then performed a simulation identical to the one previously discussed, but now imposing the additional constraint given by Eq. (46) on the α^C parameters. This simulation is labeled as MAX N6 and the corresponding mass-loss curve is displayed in the right-hand panel of Fig. 11. To assess resolution effects, we ran a mirror simulation now using N = 10⁵ particles (MIN N5). Moreover, the uncertainties associated with the use of Eq. (46) can be estimated by looking at the mass loss rate of a simulation (MIN N6) where instead of Eq. (46) the same equation was used but the constraint was derived using the MIN operator instead of the MAX one.

For comparative purposes, in the right-hand panel of Fig. 11 the mass loss of the cloud (black squares), as found in a similar blob test performed by Springel (2011) using the new moving-mesh code AREPO, has been also inserted. The number of resolution elements is 64 × 64 × 128, so that the resolution is approximately equivalent to that of the tests shown here.

A striking result seen from the behavior of the mass loss rates is that, introducing the constraint (46) on the numerical heat transfer in the simulations, the mass depletion of the cloud is now in better agreement with previous results and in particular with the cloud mass evolution produced in the blob test by the Voronoi-based code AREPO.

Uncertainties associated with the procedure adopted to estimate the constraint on the α^C parameters are of limited impact, as can be assessed by the differences between the mass loss rates exhibited by the runs MAX N6 and MIN N6. Moreover, for the run MAX N5 the mass loss rate is much stronger than in the case MAX N6. This is indicative of the numerical diffusion effects that dominate the blob evolution. A numerical simulation with N = 10⁷ particles is clearly required to clarify this issue; however, the agreement with previous results of the mass loss for the run MAX N6 suggests that convergence is achieved using as few as N = 10⁶ particles.

Finally, both the runs MAX N6 and MIN N6 exhibit at t ≳ 3τ_KH a remnant cloud mass that is  ~20−30% of the initial mass. Although this result can still be due to the use of a constraint on the numerical heat diffusion that still needs to be refined, we note that these remaining masses are equally present in SPH blob tests that have adopted, with the purpose of avoiding the problems of standard SPH, completely different approaches (Heß & Springel 2010; Murante et al. 2011).

We therefore suggest that this underestimate in the stripping rate of the blob mass is not due to the AC scheme, but depends rather on the SPH formulation and is associated with the intrinsic errors in gradient estimates. These errors, in turn, lead to the suppression of the small-scale instabilities. Such a issue will be discussed further in detail in Sect. 4.2, which is dedicated to the Rayleigh-Taylor instability.

To summarize, SPH simulations that incorporate the AC scheme can successfully be used to accurately model the blob evolution. However, the results of the tests indicate that it is necessary to properly constrain the AC parameter α^C in order to avoid unphysical heat diffusion when strong supersonic flows are present.

4. Gravity tests

To validate the performances of the AC signal velocity (25), we investigate several hydrodynamical test problems in which the self-gravity must be necessarily taken into account to properly model the system evolution. We first considered the 3D collapse of a cold gas sphere initially at rest. This test is customarily used in SPH to judge the code capability when strong shocks are present, such as those occurring during the formation of self-gravitating structures. We then examined a 2D version of the classic Rayleigh-Taylor instability test and finally the new code was used to assess the behavior of cluster entropies in cosmological structure formation.

4.1. Cold gas sphere

A standard hydrodynamical test for SPH codes in which gasdynamics is modeled including self-gravity is the 3D collapse of a cold gas sphere, also commonly referred to as the “Evrard” collapse test (Evrard 1988). The test follows in time the adiabatic collapse of a initially cold gas sphere and has been widely used by many authors (Hernquist & Katz 1989; Steinmetz & Mueller 1993; Wadsley et al. 2004; Springel 2005; Wetzstein et al. 2009; Springel 2010b; Heß & Springel 2010, V11) as a standard test for SPH codes.

The gas cloud is spherically symmetric and initially at rest with mass M, radius R, and density profile $\mathit{\rho }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{M}}{\mathrm{2}\mathit{\pi }{\mathit{R}}^{\mathrm{2}}}\frac{\mathrm{1}}{\mathit{r}}\mathrm{·}$(47)The gas obeys the ideal gas equation of state with γ = 5/3 and the thermal energy per unit mass is initially set to u = 0.05GM/R. The SPH simulations are performed using units for which G = M = R = 1 and the chosen time unit is the cloud free-fall timescale ${\mathit{t}}_{\mathrm{ff}}\mathrm{=}\mathrm{(}{\mathit{\pi }}^{\mathrm{2}}\mathit{/}\mathrm{8}\mathrm{)}\sqrt{{\mathit{R}}^{\mathrm{3}}\mathit{/}\mathit{GM}}\mathrm{=}{\mathit{\pi }}^{\mathrm{2}}\mathit{/}\mathrm{8}$.

With these initial conditions, the pressure support of the gas sphere is negligible and the cloud begins to collapse until a bounce occurs in the core with a subsequent shock wave propagating outward. Most of the kinetic energy is converted into heat at the epoch of maximum compression of the gas, which occurs at t ~ 1.1. The initial conditions setup is realized by stretching the radial coordinates of a glass-like uniform distribution of N = 88   000 equal mass particles located within a sphere of unit radius, so as to generate the density profile of Eq. (47).

To construct the set of SPH simulations, four tests were run using the kernels CS and M₅. The number of neighbors is set to 100 for the M₅ runs and 50 for the CS ones. To assess the convergence properties of the radial profiles of the considered hydrodynamic variables, the M₅ run was replicated using 200   000 particles (M₅ − 200 k). All of these runs were performed by incorporating in the SPH equations the AC term (18). For the CS kernel, a reference run with the AC term disabled was considered (CS-NOAC). The gravitational softening length was assumed to be ε_g = 0.02.

For each test case, we performed simulations with different settings of the AV parameter l_d. However, introducing a time-dependent AV scheme in SPH reduces the severity of the numerical viscosity effects and in particular, for the test investigated here, produces a radial entropy profile at the shock front in better agreement with the 1D PPM reference solution. For this reason, the results for different test cases are presented for runs where a standard AV formulation was used (runs AV₀ of V11). This is done to highlight differences in the simulation results due to the use of different kernels and the inclusion of the AC term in the SPH thermal energy equation. This choice of simulation parameters allows one to compare the test calculations performed here with previous results presented in Sect. 5.2 of V11.

The average radial profiles at t = 0.8 of density, entropy, time-dependent AC parameter, and radial velocity are displayed in Fig. 12. The black solid lines indicate the profiles of the 1D PPM calculation of Steinmetz & Mueller (1993). Broadly speaking, one expects to see some differences between the profiles of the NOAC simulation and those of the other runs. This should be valid in particular for the entropy profiles. However, at the considered epoch, Fig. 12 shows that the entropy profile of the NOAC run is still very similar to the others. This occurs because at t = 0.8 the shock front propagation is still in the early phase and the smoothing in the thermal energy due to the AC term is not yet significant.

The behavior of the α^C profiles is in accordance with that expected. In particular, the profiles exhibit a peak that corresponds to the shock front location, which occurs approximately at r ~ 0.18, and a decay as one moves away from it. The profile shapes depend on the kernel shape and weakly on resolution. The former is interpreted as a consequence of the improved mixing capability provided by the higher order of the kernel. We note that there is an increase in the profiles as one moves toward the sphere center, which occurs because of the presence of a temperature gradient. However, in Eq. (18) this radial dependence is of no consequence because behind the shock front the sphere quickly achieves a hydrostatic equilibrium and ${\mathit{v}}_{\mathit{ij}}^{\mathrm{AC}}\mathrm{\left(}\mathit{grav}\mathrm{\right)}\mathrm{\simeq }\mathrm{0}$ for r ≲ 0.1.

The epoch examined in Fig. 12 was chosen so as to compare the simulation profiles with the corresponding ones presented in Fig. 40 of Springel (2010b) and realized using the new moving-mesh code AREPO. In particular, a visual inspection shows that ahead of the shock front the accuracy of the entropy profiles shown in Fig. 12 can be considered to be in-between that of the profile produced using the AREPO code and the one realized with the SPH code Gadget-2 (right panels in Fig. 40 of Springel 2010b, test cases “moving-mesh” and “SPH”, respectively). We note, however, that a strict comparison between the different simulation profiles is difficult because the number of resolution elements employed by Springel (2010b) in the cold-gas sphere tests are lower by about  ~1/3 than the one used here. In particular, the smaller shock broadening found here can be ascribed to a resolution effect.

Finally, Fig. 13 shows a closer view of the radial density and entropy profiles in the proximity of the shock front at t = 0.9. The NOAC simulation exhibits the worst accord with the PPM reference solution, while for the high-resolution run M₅ − 200k the improvement in accuracy is minimal with respect to the other AC runs. This suggests that to obtain a significant improvement in the adherence of the simulation profiles with the PPM solution one must use a much higher number of simulation elements (say  ≳ 10⁶, Springel 2010b).

A noteworthy feature of the profiles of Fig. 13 is that now the AC simulation profiles are in much better agreement with the PPM profiles. This is interpreted as a numerical effect in which the pre-shock entropy, which is generated by the AV implementation ahead of the shock front owing to the converging flow, is now strongly reduced because the AC term in the SPH energy equation now removes this excess of internal energy. This in turn implies that, behind the shock front, a thermal energy profile is closer to the PPM solution profile.

Comparisons of the profiles of Fig. 13 with the corresponding ones displayed in Fig. 4 of V11 show that the dependence of the simulation profiles on the adopted AC scheme is stronger than that due to the time-dependent AV formulation. In fact, these results are almost unchanged when one performs SPH simulations using a time-dependent AV scheme in place of the standard one.

Finally, two runs were performed using the CRT and LIQ kernels and a number of neighbors set to N_sph = 50. The profiles of the two simulations have not been shown here to avoid overcrowding in the panels. The performances of the CRT run are quite similar to those of the CS one, whilst the profiles of the LIQ simulation are largely inconsistent with the reference PPM solution profiles.

4.2. 2D Rayleigh-Taylor instability

The Rayleigh-Taylor (RT) instability arises when a heavier fluid is placed on top of a lighter fluid (Chandrasekhar 1961) in the presence of an external gravitational field. The fluids are in pressure equilibrium with the external field, and in this configuration the system is unstable in the presence of small perturbations at the interface. The lighter fluid will then begin to rise and the denser one to fall. This process leads to the development of characteristic finger-like structures before the fluids enter the non-linear phase where they completely mix together. To evaluate the ability of the AC-SPH code to properly describe the evolution of RT instabilities, we consider a 2D version of the test. The initial conditions were chosen to be similar to the numerical test implemented by Abell (2011) to validate its new rpSPH formulation of SPH equations, so as to compare Abell’s results consistently. The computational domain consists in a 2D box with coordinates x ∈  {0,1/2}, y ∈  { 0,1 } . The boundaries are periodic in x and reflecting in y. The density is ρ₁ = 2 at the top and ρ₂ = 1 at the bottom, with a density profile $\mathit{\rho }\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{2}}\mathrm{+}\frac{\mathrm{(}{\mathit{\rho }}_{\mathrm{1}}\mathrm{-}{\mathit{\rho }}_{\mathrm{2}}\mathrm{)}}{\mathrm{[}\mathrm{1.}\mathrm{+}\exp {\mathrm{\{}\mathrm{-}\mathrm{2}\mathrm{(}\mathit{y}\mathrm{-}\mathrm{0.5}\mathrm{)}\mathit{/}{\mathrm{\Delta }}_{\mathit{y}}\mathrm{\}}}^{\mathrm{]}}}\mathit{,}$(48)where Δ_y = 0.05. This ensures a smoothing in density at the interface y = 0.5 that allows a consistent numerical behavior (Robertson et al. 2010). This density profile is realized by constructing a HCP lattice of N = 620² equal mass particles in which the spacing is varied according to the procedures described in Sect. 3.1.1 until the profile of Eq. (48) is satisfied.

The pressure at the interface is set to P₀ = ρ₁/γ = 10/7, where γ = 1.4 and varies with y according to P(y) = P₀ − gρ(y)(y − 1/2), g = 1/2 being the external acceleration, so that the system is initially in hydrostatic equilibrium. For the particles that satisfy the condition  | y − 0.5 |  < 0.2, a velocity perturbation is applied in the y direction given by $\begin{array}{ccc}{\mathit{v}}_{\mathit{y}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}& \mathrm{=}& \mathit{\delta }{\mathit{v}}_{\mathit{y}}\left\{\mathrm{1}\mathrm{+}\cos \left[\mathrm{8}\mathit{\pi }\mathrm{(}\mathit{x}\mathrm{+}\mathrm{1}\mathit{/}{\mathrm{4}}^{\mathrm{)}}\right]\right\}\\ & & \end{array}$(49)where δv_y = 0.1.

We present results for the CS, CRT, M₅, and, for reasons that we discuss later, we consider also the M₆ kernel (Price 2012a, Sect. 2.3).

The number of neighbors for these runs is set to N_sph = 32, 32, 50, 110, respectively. As in the other test cases considered here, for comparative purposes, we ran a simulation (CS-NOAC) in which the AC term in the thermal energy equation was disabled. The following settings were used for the AV parameters  { α_min,α_max,l_d }  =  { 0.1,1.5,0.2 } .

In Fig. 14, at the time t = 0.3, we present the 2D density maps of the different RT tests. As expected, standard SPH is unable to correctly capture the development of the RT instability, as indicated by the leftmost panel of Fig. 14 (CS kernel, NOAC run). However, from this viewpoint the improvement is minimal even for the corresponding AC version of the considered run (CS-32). On the contrary, a significant improvement is obtained if one uses the M₅ kernel or the next in order M₆, which is a quintic polynomial (Price 2012a). In fact, in the rightmost panel of Fig. 14 (M₆ run), the center of mass of the RT spikes appears at a lower position than in the M₅ run and therefore the convergence might still not have been reached yet, even for the M₆ run.

This strong dependence of the code performances on the kernel order demonstrates that, for the RT test, accuracy in pressure forces is a fundamental issue. These results are consistent with those of Sect. 3.1 and indicate that the poor performances of SPH when handling KH or RT instabilities, in particular for very subsonic flows, are mainly due to the leading errors in the momentum equation. These errors decrease as the order of the kernel increases, which implies that the accuracy in pressure force estimates is higher and thus the velocity noise is lower. This, in turn, indicates a better capability to capture the growth of the instability.

Similar results were obtained by McNally et al. (2012), who ran a suite of 2D KH simulations using carefully crafted initial conditions with the goal of assessing different hydrodynamic code capabilities. The authors conclude that for SPH the code performances are strongly related to the order of the kernel employed in the simulation and thus to the accuracy of the gradient estimates.

Moreover, the simulation behavior of Fig. 14 is consistent with the findings of García-Senz et al. (2012), who present a new formulation of SPH in which gradients are estimated using a tensorial approach that conserves both linear and angular momentum. The authors demonstrate that when using this higher-order gradient estimator the interpolation of physical quantities is significantly more accurate than in the standard method. In particular, several tests showed that the code is now able to successfully capture the development of KH and RT instabilities.

To summarize, the results of this section demonstrate that the performance of SPH in modeling hydrodynamic subsonic instabilities depend critically on the accuracy of the gradient estimator. The formulation proposed by García-Senz et al. (2012) looks particularly promising in this aspect in overcoming the present SPH difficulties.

Finally, it must be stressed that the simulation results of Fig. 14 were obtained by setting δv_y = 0.1. If the same test had been performed using δv_y = 0.01, as in the RT test of the new moving-mesh code AREPO (Springel 2010b, Sect. 8.8), the velocity noise would have suppressed the growth of the instability even for the M₆ run. We note that in their RT test, García-Senz et al. (2012) were able to recover the growth of the instability using the same amplitude of the velocity perturbation, δv_y = 0.01.

4.3. Cluster comparison

We now discuss a set of cosmological cluster simulations that we constructed using the standard SPH formulation as well as the new AC-SPH version. We only considered non-radiative, or “adiabatic” simulations, in which the hydrodynamics were modeled according to the formulation presented in Sect. 2. Previous investigations (Frenk et al. 1999; Wadsley et al. 2004, 2008; Mitchell et al. 2009; Springel 2010b) showed that for the same initial conditions there are systematic differences between results extracted from cluster simulations produced using standard SPH and Eulerian AMR codes. Specifically, the level of central entropy is found to be lower in SPH simulations than in the corresponding AMR runs. The latter simulations are characterized by the presence of a flat entropy core, whereas the radial entropy profile produced by SPH runs increases steadily with radius. As outlined in the Introduction, the origin of this discrepancy has been the subject of intense debate (Agertz et al. 2007; Mitchell et al. 2009; Read et al. 2010) that has given as the main explanation the different degrees of numerical mixing present in the two codes.

To assess the capability of the AC-SPH code in solving this issue, we then ran several cluster simulations. A comprehensive study of the differences between the hydrodynamic variables of cluster simulated using the standard and the AC version of the SPH code will be presented in a forthcoming paper. Here, we only show the final radial entropy profile of the simulated clusters as the main variable that can be used to test the effectiveness of the AC-SPH formulation in solving the entropy problem in cluster cores.

A detailed description of the procedures used to construct the simulated cluster sample is given in Sect. 2 of V11, and we provide a brief summary here. The simulations were carried out assuming a spatially flat ΛCDM model, with matter density parameter Ω_m = 0.3, vacuum energy density Ω_Λ = 0.7, baryonic density Ω_b = 0.0486, and h = H/100   km   s^-1   Mpc^-1 = 0.7. The scale-invariant power spectrum is normalized to σ₈ = 0.9 on an 8   h^-1 Mpc scale at the present age t₀. An N-body cosmological simulation involving only gravity was first run with a comoving box of size L₂ = 200   h^-1 Mpc, starting from an initial time t_in and ending at the final epoch t₀. At this epoch, clusters of galaxies were identified in the simulation box as groups of particles that are associated with overdensities approximately in excess of $\mathrm{~}\mathrm{200} {\mathrm{\Omega }}_{\mathrm{m}}^{-0.6}$. Several of these clusters were then resimulated individually using the AC-SPH code described here, coupled with a treecode gravity solver.

Initial conditions for a specified cluster were generated as follows: a spherical region with origin located at the cluster center was populated with a high-resolution (HR) grid, the radius of the HR sphere being such that it contains all of the original particles identified at t = t₀ as cluster members. Both a gas particle and a dark matter particle were associated with each grid node, whose positions were perturbed according to the random realization of the original cosmological simulations and only particles located within the HR sphere were kept for the hydro run. The HR particles were surrounded by a low-resolution shell of dark matter particles, extracted from a grid with spacing twice that of the HR grid, the shell being introduced for the purpose of mimicking the effects of tidal forces. The grid spacing was chosen such that at the end of the procedure a cluster was simulated with N_gas ~ 220   000 gas particles and N_dm ~ N_gas dark matter particles in the inner HR sphere, whereas ${\mathit{N}}_{\mathrm{dm}}^{\mathrm{ext}}\mathrm{~}{\mathit{N}}_{\mathrm{dm}}$ were used in the low-resolution shell. Particle masses were assigned according to the values of Ω_m and Ω_b.

The clusters selected to be re-simulated hydrodynamically were chosen with the following criterion. Originally (V11), the procedure previously described was repeated two more times to construct two new cluster samples extracted from cosmological simulations with box sizes L₄ = 2L₂ and L₈ = 2L₄. The three cluster samples were then combined to construct a final sample S_all of  ~ 160 clusters covering nearly a decade in cluster masses. Cluster members of sample S_all were then ranked according to their dynamical state, as measured by an appropriate statistical indicator. Two sub-samples of four test clusters each, denoted respectively by Q and P, were then constructed by extracting from sample S_all clusters with membership criterion for sample Q (P) of being in a fully relaxed (perturbed) state. These two sub-samples were then used to investigate the effects of AV in hydrodynamic SPH simulations of galaxy clusters (V11). Here the four test clusters of sub-sample Q were chosen to investigate the performances of the AC-SPH code. This choice was motivated by the need to proper disentangle the amount of entropy mixing produced during merger processes from that specific to the numerical method in the final cluster entropy. The clusters of sub-sample Q were carefully selected on the basis of their dynamical state and are among the most relaxed of sample S_all. Therefore, for these clusters, differences in the final radial entropy profile between standard and AC runs can be safely ascribed to the new AC term implemented in the code. The AV simulation parameters of the simulations were  { α_min,α_max,l_d }  =  { 0.1,1.5,0.2 } .

The final entropy profiles for the four test clusters are displayed in Fig. 15, where the cluster entropies have been normalized to ${\mathit{S}}_{\mathrm{200}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\left[\frac{\mathrm{2}\mathit{\pi }}{\mathrm{15}}\frac{{\mathit{G}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{200}}}{{\mathit{f}}_{\mathrm{b}}\mathit{H}}\right]}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathit{,}$(50)where f_b = Ω_b/Ω_m is the global baryon fraction and ${\mathit{M}}_{\mathrm{\Delta }}\mathrm{=}\mathrm{(}\mathrm{4}\mathit{\pi }\mathit{/}\mathrm{3}\mathrm{)}\hspace{0.17em}\mathrm{\Delta }\hspace{0.17em}{\mathit{\rho }}_{\mathrm{c}}\hspace{0.17em}{\mathit{r}}_{\mathrm{\Delta }}^{\mathrm{3}}$ denotes the mass contained in a sphere of radius r_Δ with mean density Δ times the critical density ρ_c = 3H²/8πG. As is customary, the cluster mass is defined by setting Δ = 200 and r₂₀₀ denotes the cluster radius.

The results of Fig. 15 indicate that, although there is a certain degree of scatter between the entropy profiles of individual clusters, all of the AC-SPH runs consistently exhibit much shallower entropy profiles and higher core entropies than in the standard SPH runs at r/r₂₀₀ ≲ 0.1. For a given cluster, the difference in the levels of central entropies produced by the two codes is about a factor  ~4, and now the values of central entropies are comparable with those produced using AMR codes in cosmological cluster simulations (Voit et al. 2005) and idealized cluster binary mergers (Mitchell et al. 2009). It should, however, be stressed that in Eulerian hydrodynamics it is the numerical scheme that forces the fluid to be mixed below the minimum cell size. This suggests that AMR simulations of galaxy clusters might overestimate the correct level of core entropy because of fluid mixing (Springel 2010b). To summarize, the development of entropy cores in the AC-SPH runs is clearly a consequence of the heat diffusion term described by Eq. (18) now present in the energy equation. This term acts to redistribute the internal energy produced in the shocks during dissipative processes, so that the subsequent entropy mixing leads, in the end, to higher levels of core entropies.

Similar results were obtained by Wadsley et al. (2008), who similarly added a heat diffusion term to the SPH energy equation, but using the prescription of Eq. (26) instead of the Eq. (18) adopted here. According to Wadsley et al. (2008), the value of the coefficient C tends to be problem dependent and for the cluster simulations the best agreement with the entropy profiles extracted from the mesh code runs is obtained by setting C ≃ 0.1−1. Here, the time-dependent formulation (22) presents the advantage of minimizing thermal diffusion away from shocks, as demonstrated in the cold gas sphere test by Fig. 12, where the radial profile of the AC parameter α^AC is peaked at the shock location.

The simulations presented in this section were realized using adiabatic gas physics. However, a realistic modeling of the intracluster medium physics must incorporate the possibility that the gas cools radiatively. Previous simulations (V11) showed that the presence of cooling leads to the development in the inner cluster regions of dense compact cool-gas cores and subsequently of high levels of turbulence, the latter being produced by the hydrodynamical instabilities generated by the interaction of the compact cool cores with the ambient medium.

How this scenario is affected by incorporating into the SPH equations a numerical heat diffusion term is an issue that can be properly addressed only by resorting to numerical simulations. Here, we outline that two competing effects are expected to influence the final level of turbulence in a cluster core: the improved capability of the code to describe the development of hydrodynamical instabilities, and thus an increase in the amount of turbulence, and the reduced availability of cool gas because of the presence of material of higher entropy, which has been brought in the inner core by the enhanced fluid-mixing properties of the code.

Finally, a resolution study previously performed in V11 found for the simulation resolution employed here that there is a substantial stability in the entropy profiles of the simulations. Although the simulations of V11 do not incorporate the AC term, we do not expect strong variations in the resolution dependence of the profiles, given the similarity in the level of core entropies with those found using mesh-based codes.

5. Summary and conclusions

We have presented an SPH numerical scheme that incorporates an artificial conductivity term and uses an appropriate signal velocity for simulations including gravity. The AC formulation was introduced by Price (2008) as a solution to the problems encountered by standard SPH to correctly follow the development of KH instabilities, owing to the inconsistencies of the standard formulation in the description of density at a contact discontinuity. Here, a suite of hydrodynamic test problems has been investigated with the purpose of validating the new AC-SPH code and assessing its performances when using the specifically adopted signal velocity.

The results of the KH instability test were presented in Sect. 3.1, in which the code capabilities were tested by considering SPH simulations of the KH test performed using a large variety of initial condition set-up and SPH kernels. The set of initial conditions was chosen as in Valcke et al. (2010), so as to consistently compare the results. These are in accordance with the corresponding ones of Valcke et al. (2010) and indicate that, for the version of the KH test analyzed here, the AC implementation is important for the long-term behavior of the simulation.

Moreover, the results of the KH test indicate (Valcke et al. 2010; Read et al. 2010; McNally et al. 2012) that the poor performances of standard SPH in properly treating KH instabilities can be explained in terms of two distinct effects. The first is a general problem of consistency, which for the problem under consideration requires that smooth interfaces should be present at contact discontinuities, in order to obtain numerical convergence. The second effect, which in standard SPH suppresses the growth of KH instabilities, is the leading error in the momentum equation, which is caused by incomplete kernel sampling (Read et al. 2010) and quantified by the norm E⁰ defined by Eq. (36).

To circumvent this problem, several proposals have been made (Valcke et al. 2010; Read et al. 2010) in which the standard SPH cubic spline M₄ kernel is replaced by the new LIQ or CRT kernels with steeper central profiles. These kernels present the advantage of being stable against particle clumping, so that the number of neighbors can be safely increased to reduce the error E⁰.

In this paper, we have considered the possibility of reducing sampling errors by considering higher-order B-spline kernels, specifically the M₅ or quartic spline kernel. A striking result of the KH runs presented in Sect. 3.1 is that, for a given value of the ratio η of the smoothing length to the mean interparticle spacing, simulations of the KH test performed using the M₅ kernel have amplitudes of the E⁰ error substantially smaller and in line with the behavior of the same quantity for simulations employing either the LIQ or CRT kernels. A linear stability analysis reveals that this result follows owing to the very good stability properties of the M₅ kernel. The analysis indeed suggests that the clumping instability is absent for values of η up to η ≲ 2, which in 3D corresponds to N_sph ~ 520 neighbors.

These findings are consistent with the results of Dehnen & Aly (2012), who proposed to adopt, with the specific purpose of avoiding the clumping instability, the Wedland (1995) functions as a new class of kernels. The authors showed that in terms of stability and accuracy properties, the quartic spline performs extremely well when compared with the proposed Wedland functions. These results are strictly connected to the properties of the kernel Fourier transform, which according to the authors must be non-negative to avoid particle clumping, and are consistent with the findings of Sect. 3.1.3 since when increasing the kernel order both the B-spline and the Wedland functions approach the Gaussian.

We therefore propose, as a compromise between the need to reduce sampling errors while keeping the computational cost to a minimum, the use of the M₅ kernel with a neighbor number in the range N_sph ~ 60−120 as the standard combination that guarantees sufficient accuracy in many SPH simulations of astrophysical problems. Moreover, the results of the Sod shock-tube test of Sect. 3.2 demonstrate that to obtain simulation profiles in accordance with the analytic solution, for simulations employing kernels with a modified shape the use of a much larger number of neighbors than in the case of the M₅ runs is necessary.

The results of the gravity tests show that the adoption of the AC-SPH scheme significantly reduces, at the level tested in this paper, the differences seen in the hydrodynamics between standard SPH and grid-based simulations of self-gravitating structures. For the cold gas sphere, the entropy profile is in closer agreement with the PPM reference solution and for the cosmological cluster simulations, a key result is the final level of core entropies, which are consistent with those of the central entropies produced using AMR codes. Thus, it appears that in hydrodynamic simulations where self-gravity is important the AC term, accompanied by the proposed signal velocity, plays a key role as a mechanism for redistributing thermal energy and hence as a source of entropy mixing.

To summarize, results extracted from simulations of hydrodynamic tests where self-gravity is the dominant factor are in much better agreement with the corresponding ones obtained using mesh-based codes. These results then demonstrate the capability of the implemented AC-SPH scheme in properly following the formation of cosmic structures.

We note that the artificial heat conduction term was originally proposed by Price (2008) with the purpose of avoiding the inconsistencies encountered by standard SPH in the presence of density steps at contact discontinuities. A complementary view was proposed by Wadsley et al. (2008), who introduced the same term, albeit in a different numerical formulation, with the aim of modeling the level of diffusion due to turbulence. The two interpretations are not mutually inconsistent, although the results of the self-gravity tests presented in this paper support the view of a heat diffusion term that in SPH is capable of mimicking the diffusion due to turbulence. In a similar fashion, Violeau & Issa (2007) presented an SPH scheme to model a free-surface incompressible flow that, in analogy with 3D Large Eddy Simulations, assumes a Smagorinsky (1963) model for the filtered Navier-Stokes equations.

The results of the blob test simulations demonstrate that the instabilities leading to the expected cloud disruption can develop only when the SPH energy equation incorporates the AC term. A particularly interesting result is that an appropriate limiting condition must be implemented on the AC coefficients α^C in order to avoid an unphysical amount of heat diffusion, which in turn leads to a cloud disruption that occurs too early. This limiter was identified with the Prandtl number and, for the AC signal velocity adopted here, severely limits the amplitude of the AC coefficients in the regime of strong supersonic flows.

The AC-SPH simulations of the blob test incorporating now the new constraint support this view, since Fig. 11 shows for the new runs a cloud mass-loss rate that is in closer agreement with the rates obtained from simulations realized using a completely independent numerical scheme (Springel 2011). However, it must be stressed that the condition specified by Eq. (43) was calculated for a perfect monoatomic gas with γ = 5/3, so that a physically motivated constraint on the artificial heat diffusion of the simulated medium should be considered problem-dependent.

However, the code still has several problems that render its use problematic if the development of hydrodynamic instabilities need to be followed in the regime of subsonic flows. The results of the simulations indicate that these shortcomings are not due to the AC implementation, but rather are intrinsic to the standard formulation with which gradients are calculated in SPH and the related errors are subsequently introduced in the momentum equation. In particular, simulations of the RT test show that increasing the kernel order alleviates the problems but does not solve them. These results are in accordance with recent findings (McNally et al. 2012; García-Senz et al. 2012) and clearly demonstrate that for very subsonic flows, the poor performances of SPH in modeling hydrodynamic instabilities are strictly connected to the code accuracy in gradient estimates.

The formulation of García-Senz et al. (2012), which was proposed with the aim of calculating SPH gradients with accuracy as high as possible while keeping the benefits of a Lagrangian formulation, looks in this respect very promising and encourages further investigations along the numerical approach proposed in this paper. These would be particularly relevant in the light of the results of Bauer & Springel (2012), who claim that standard SPH fails to properly model the regime of subsonic turbulence. They reach their conclusions by comparing results extracted from simulations of driven subsonic turbulence realized using the moving-mesh code AREPO and GADGET SPH. Their results were criticized by Price (2012a), for whom the use of a time-dependent AV is critical in SPH simulations of subsonic turbulence.

A re-analysis of the Bauer & Springel (2012) simulations using the AC-SPH code presented here, augmented with improved gradient operators, would then be fundamental for achieving a deeper understanding of the capability of the different numerical methods to model subsonic turbulence, which is expected to have a significant impact on shaping the thermodynamic properties of baryons in cosmological haloes and, subsequently, the process of galaxy formation.

---

References

All Tables

All Figures