© ESO, 2017

1. Introduction

Many interesting problems in Astrophysics involve the evolution of fluids and plasmas coupled with complex physics. For example, in core collapse supernova, magnetohydrodynamics meets with general relativity, nuclear processes and radiation transport. Other scenarios, such as neutron star mergers, Type Ia supernova, and planet- or star formation, face similar challenges in terms of complexity. Besides that, these phenomena often have a strong dependence on the dimensionality and they must be studied in three dimensions. This requires accurate numerical tools which translate to rather sophisticated hydrodynamic codes. Because of its adaptability to complex geometries and good conservation properties, the Smoothed Particle Hydrodynamics (SPH) method is a popular alternative to grid-based codes in the astrophysics community. SPH is a fully Lagrangian method, born forty years ago (Lucy 1977; Gingold & Monaghan 1977), that since then has undergone sustained development (Monaghan 1992, 2005; Rosswog 2015a; Springel 2010b; Price 2012). Recent years have witnessed a large range of improvements specially aimed at reducing the numerical errors inherent to the technique. These errors are known as E₀ errors (Read et al. 2010) and mainly appear due to the conversion of the integrals, representing local-averaged magnitudes of the fluid, into finite summations. The most simple and naive way to get rid of them would be to work closer to the continuum limit, which implies working with a number of particles and neighbors as big as possible (ideally, N → ∞ and N_nb → ∞, see Zhu et al. 2015). Unfortunately, this is not feasible in common applications of the technique because the total number of particles is limited by the computing power (both by speed and storage). Moreover, the number of neighbors of a given particle cannot be arbitrarily increased without suffering pairing-instability (Schuessler & Schmitt 1981). This is a numerical instability that acts as an attractive force which appears at scales slightly shorter than the smoothing length h, provoking artificial particle clumping and effectively decreasing the quality of the discretization, which eventually leads to unrealistic results.

In order to reduce the E₀ errors, another more practical possibility has been studied during recent years: finding interpolating functions that are less prone to particle clustering than the widely used M₄ or cubic-spline kernel (Monaghan & Lattanzio 1985). Among the various candidates, the most used (pairing-resistant) kernels come either from an extension of the M_n family to higher-order polynomials (Schoenberg 1946) or those based on the Wendland functions (Wendland 1995). In particular, the Wendland family is specially suited to cope with pairing instability (Dehnen & Aly 2012). Another possibility is the sinc family of kernels (Cabezón et al. 2008), which are functions of type S(x) = C(n)(sinx/x)ⁿ, and add the capability of dynamically modifying their shape simply by changing the exponent n. That adaptability of the sinc kernels makes the SPH technique even more flexible and can be used, in particular, to prevent particle clustering, as shown in Sect. 2.1.

Historically, the growth of subsonic hydrodynamical instabilities has been problematic for SPH simulations, as they damp them significantly. The Rayleigh-Taylor (RT) instability is a ubiquitous phenomenon which serves as a paradigmatic example. It appears wherever there is cold and dense fluid on top of a hot and diluted one in the presence of gravity (or any inertial force, in virtue of the principle of equivalence). The entropy inversion leads to the rapid overturn of the fluid layers. In the real world, the overturn is triggered by small perturbations at the separation layer between the light and dense fluids. The RT instability is one of the most important agents driving the thermonuclear explosion of a white dwarf, which gives rise to the Type Ia supernova (SNIa) explosions. Their correct numerical description is also crucial to understanding the structure of the supernova remnants (SNR) and to model core collapse supernova (CCSN) explosions. Additionally, the RT instability is also the source of other interesting phenomena such as the KH instability or turbulence. The numerical simulation of the Rayleigh-Taylor instability using SPH has traditionally been a drawback for the technique, especially for low-amplitude initial perturbations in the presence of a weak gravitational force. At present, for a similar level of resolution, the best SPH codes cannot yet compete with the state-of-art grid-based methods. For example, the finite-volume/difference Godunov methods such as ATHENA and PLUTO, AMR codes such as FLASH, the Meshless Finite Mass (MFM) and Volume methods (MFV), and especially the moving mesh methods based on Voronoi tessellations, as in the code AREPO, provide a good approach to the RT instability. This problem is partially overcome using a large number of neighbors and by adding an artificial heat diffusion term to the energy equation, as in the PSPH proposal by Saitoh & Makino (2013) and Hopkins (2013). However, these problems still persist when either the size of the initial perturbation or the gravity value are reduced (Valdarnini 2012), meaning that they are a symptom of another source of numerical error in SPH named tensile instability. This is an artificial surface tension that appears at contact discontinuities because of an insufficient smoothness of pressure between both sides of the discontinuity (Monaghan 2000). As a consequence, the integration of the momentum equation gives incorrect results. An excess of that tension provokes the damping of fluid instabilities, especially those with short wavelengths. Several techniques have been proposed to treat this problem, like averaging the pressure by means of the interpolating kernel itself, scheme PSPH (Hopkins 2015), volume element estimation bounded to pressure, the density independent scheme (DISPH; Hopkins 2013; Saitoh & Makino 2013, 2016), or by adding an artificial diffusion of heat to the energy equation, which helps to balance the pressure across the discontinuity (Price & Monaghan 2007). They have paved the road that led the SPH technique to a new standard within the last few years, and have helped to overcome this long lasting inconvenience. In particular, it has been proved that it is fundamental to increase the accuracy of the gradient estimation across the contact discontinuities via reducing its numerical noise. To achieve that, García-Senz et al. (2012) used an integral scheme to calculate spatial derivatives, which proved to be especially efficient at handling fluid instabilities. The validity of that approach was assessed in subsequent works by Cabezón et al. (2012) and Rosswog (2015a). In this last case, the integral approach to the derivatives (IAD) was used to extend the SPH scheme to the special-relativistic regime (see also Rosswog 2015b). See also the recent work of Valdarnini (2016), where the efficiency of the IAD scheme to reduce the E₀ errors is studied in detail.

Finally, a recent breakthrough in SPH was the emergence of the concept of generalized volume elements (Ritchie & Thomas 2001; Saitoh & Makino 2013; Hopkins 2013). In these works, it was shown that a clever choice of the volume element (VE) can reduce the tensile-instability, leading to a better description of hydrodynamic instabilities. In this paper, we present a novel estimator to the VE that preserves the normalization of the kernel. We show in this work that having a good partition of the unity is also connected to the tensile-instability problem. The calculations of the growth of the Kelvin-Helmholtz and Rayleigh-Taylor instabilities using these new VE are encouraging.

In this work, we also introduce the hydrodynamics code SPHYNX, that gathers together the latest advances in the SPH technique, including those new ones presented here. SPHYNX has already been used in production runs simulating type Ia and core collapse supernova and is publicly accessible¹.

The organization of this paper is as follows. In Sect. 2 we review the main properties of the sinc kernels as well as the integral approach to the derivatives, which are at the heart of SPHYNX. Section 3 is devoted to the choice of the optimal volume element and to the update of the smoothing-length h and of the kernel index n. Section 4 describes the structure of the hydrodynamics code SPHYNX: moment and energy equations and included physics. Sections 5 and 6 are devoted to describing and analyzing a variety of tests carried out in two-dimensions and three-dimensions, respectively. Finally, we present our main conclusions and prospects for the future in Sect. 7.

2. Interpolating kernels and gradients evaluation

2.1. The Sinc kernels

Interpolating kernels in SPH must fulfill several basic requirements in order to be considered suitable interpolators. Some of these features are automatically satisfied if the functional form of the kernel is spherically symmetric and approaches the Dirac−δ in the continuum limit. A functional representation of the Dirac-δ is (Couch 1997) (1)Writing the function above becomes (2)where the magnitude is a widely known function used in signal analysis. Nevertheless, the expression given by Eq. (2) does not have compact support, limiting its practical applications as SPH interpolation kernel. Moreover, for x/h> 2, the sinc function oscillates and eventually produces negative values. An obvious solution to this problem is to restrict the domain to | x/h | ≤ 2, but in that case the first derivative of this function does not go to zero at the limits of the interval. To fix that, we define the kernel set as: (3)where S_n(·) = sincⁿ(·), q = | r | /h, and d is the spatial dimension. B_n is a normalization constant, whereas n is a real number with n ≥ 2 to guarantee the nullity of the first derivative at q = 2. We generically refer to the set as sinc kernels; see Cabezón et al. (2008; where they were called harmonic kernels) and García-Senz et al. (2014). These interpolators fulfill the regular desirable characteristics for SPH interpolators. Namely, tending to a δ-function in the continuum limit, having compact support, and being spherically symmetric. But they have a number of interesting additional features: a) they are connected by definition to the δ-function; b) they add flexibility to the SPH calculations as the kernel profile can be changed by varying the kernel index n; c) working with a high index (n ≥ 5) not only ensures that high-order derivatives are well behaved, but also overcomes particle clustering (see Sect. 2.1); and d) they tend to fulfill separability when the index n increases (see Appendix A).

Unfortunately, there is not a general analytic expression for the normalization constant B_n, but this small problem can be circumvented by fitting B_n, numerically calculated for a series of values of n. A fitting formula for B_n, valid in the range 3 ≤ n ≤ 12, was provided in García-Senz et al. (2014). We reproduce it here for the sake of completeness, (4)where the values of coefficients b₀,b₁,b₂,b₃ as a function of the dimensionality are provided in Table 1.

In the limit of large n, the sinc kernels also display an interesting feature: separability, which has not been sufficiently emphasized in the extant literature concerning SPH. Standard kernels used in SPH are spherically symmetric functions which naturally lead to a second-order accurate scheme. Interestingly, the majority of these kernels do not fulfill the identity: (5)This property guarantees the consistency of simulations involving planar symmetry, which should render identical results when calculated in 1D, 2D, or 3D, if the resolution and the interpolating kernel are the same in all three cases. Exploring in detail the applications of this property is beyond the scope of this paper. Nevertheless, we added a short discussion in Appendix A.

Additionally, some particular values of the kernel index n can mimic the profile of the most notable kernels used in SPH. In Table 2 we show the correspondence between the sinc and both, the M_n and C_m families of interpolators. The value of n shown in Table 2 was obtained by minimizing the metric distance , where W1 and W2 are the interpolators to be compared². It is remarkable that n ≃ 5 is able to approach both the quintic M₆ B-spline and the C₄ Wendland kernel. Thus, the choice n = 5 seems a suitable default value when the number of neighbors is moderate, n_b ≃ 60–120 in 3D (see f.e. Figs. 4 and 5 of Rosswog 2015a). For n_b ≥ 120 it may be advisable to raise the index of the sinc kernel to avoid the pairing instability. Interestingly, the choice n = 6.315 provides a very similar profile to that of the Wendland C₆ kernel. Nevertheless, a similarity in the real space is not the only characteristic to take into account. As proved by Dehnen & Aly (2012), having a positive Fourier transform for a large range of modes is of utmost importance when dealing with pairing instability (see Sect. 2.1).

It is also worth noting that the family forms a continuous set, making it possible to locally change the index n during the runtime of the simulation so that one can, for example, raise the value of n in presence of sharp density gradients, (see Sect. 3.2).

Finally, the evaluation of a sinc function is indeed computationally more costly than the regular spline kernels. Nevertheless, it is very easy to circumvent this problem by interpolating the value of the sinc function from a relatively small pre-calculated table. In practical applications, though, the choice of the kernel has a negligible impact on the computational burden, being completely masked by efficient cache usage and other sections of the code, as neighbor search or gravity calculation.

Although SPHYNX uses the sinc kernels by default, it also incorporates the Wendland C₂,C₄ and C₆ kernels which can easily be selected by the user whenever necessary. On another note, the impact of the C_m kernels has recently been studied by other authors (Dehnen & Aly 2012; Rosswog 2015a) so we focus here almost exclusively on the set. Nevertheless, a comparison between the performance of the C₆ and sinc interpolators is provided in Sect. 5.5, where the evolution of the Gresho-Chan vortex is simulated and discussed.

Convergence to the continuum limit is only achieved when h → 0. According to Zhu et al. (2015), it can only be reached using both a large number of particles N, and a large number of neighbors n_b. The value of N chiefly depends on the available computer resources and on algorithmic details. However, working with a large number of neighbors is not only expensive but it may be totally impractical due to the tendency of the mass-particles to gather in clusters when n_b is high, a phenomenon known as pairing instability. A leap forward to alleviate this problem was the proposal of using Wendland functions (Wendland 1995) as interpolators by Dehnen & Aly (2012), which are much less sensitive to this instability than the low-order members of the M_n B-spline family. Another possible option to avoid particle clustering is to choose a sinc kernel with a large enough exponent, García-Senz et al. (2014).

As mentioned above, having a Fourier transform with positive modes is necessary for a kernel to be stable. In Fig. 1 we show the Fourier transform of several B-spline, Wendland, and sinc kernels. As expected, Wendland kernels show positive Fourier transform at very large modes. Nevertheless, it is clear that the sinc family systematically becomes negative at higher modes and their negative regions decrease in size very quickly as the exponent n increases. Therefore, particle clustering can be totally suppressed, in practice, simply by raising n. It is also worth noting that, in order to improve the performance of Wendland kernels, it is advisable to considerably increase the number of neighbors up to n_b ≃ 400 in 3D (Valdarnini 2016) which has a substantial increase on the computational burden of actual simulations³.

We tested the impact of the kernel and the number of neighbors on the emergence of pairing instability in a dynamical simulation. Results can be found in Sect. 5.1.

2.2. The integral approach to derivatives

In SPH, gradients are usually calculated simply by applying the nabla operator to the kernel function (Gingold & Monaghan 1977). Nevertheless, the same procedure is not applied to calculating higher-order derivatives, such as the Laplace operator, needed, for example, to evaluate the heat transport contribution to the energy equation. For these cases, an integral approach is preferred (Brookshaw 1985) because it is less noisy and provides excellent results. Working in that direction, an integral approach to calculating first derivatives, called IAD, was recently introduced by García-Senz et al. (2012). The IAD approach provides a more accurate estimation of the first derivative, thus improving the overall quality of the simulations. This was demonstrated through the study of several test cases in 1D, 2D, and 3D by Cabezón et al. (2012), Rosswog (2015a), and Valdarnini (2016), hence we refer the reader to these papers for technical details of the method. In García-Senz et al. (2012) it was shown that a restriction of the full integral approach (called IAD₀) leads to a conservative formulation of the SPH equations in the same way as the standard scheme⁴ does. This, however is not without cost. Some accuracy is lost in the gradient calculation, which is exact for linear functions when IAD is used. In his work, Rosswog (2015a) evaluated this accuracy loss in several tests, and the outcome was that even if IAD₀ was used, the accuracy loss, in comparison with IAD, was much smaller than the accuracy gain when comparing with standard gradient estimators. Furthermore, Valdarnini (2016) presents a series of tests that show the effectiveness of the IAD₀ method at removing sampling errors in sub-sonic flows. Additionally, the SPH equations linked to the IAD₀ scheme are formally similar to those of the standard method provided that the kernel gradient is computed as follows (for the sake of clarity, from now on we drop the dependence of the kernel on the particle inter-distance, r_b−r_a), (6)where, (7)c_ij,a being the coefficients of the inverse of matrix for particle a, whose elements are, (8)In García-Senz et al. (2012), it is shown that for the Gaussian kernel the IAD₀ scheme is equivalent to the gradient estimation with the analytical derivative of the kernel function, this being a particular case of the integral approach. It was also proven that IAD₀, in fact, exactly reproduces the gradient of linear functions provided that, (9)Equation (9) is, in general, fulfilled only approximately, and this is the main difference between IAD and IAD₀ accuracies. The best gradient estimation is therefore attained by IAD₀ only when ⟨ Δr ⟩ ≃ 0. In this respect, the initial model should be built so that Eq. (9) is approached as much as possible. As we see in the following section, it turns out that getting a small ⟨ Δr ⟩ is closely related to the choice of the adequate volume element for the integration. An analytical proof showing that an improvement in the partition of the unit translates to a better estimation of gradients within the IAD₀ paradigm is provided in the Appendix B.

A collateral, albeit positive, feature of IAD₀ is that it hinders the emergence of the pairing stability. This is shown in Sect. 5.1.

3. Preconditioning the calculation

3.1. The choice of the volume element

A recent breakthrough in SPH has been the development of a general method to assign the volume element to a particle (Saitoh & Makino 2013; Hopkins 2013, 2015; Rosswog 2015a). In SPH, the volume element of particle a has traditionally been V_a = m_a/ρ_a but, as noted in the seminal proposal by Ritchie & Thomas (2001), other options could be more convenient for handling specific problems. SPHYNX makes use of the scheme developed by Hopkins (2013), where an estimator X_a is introduced so that the particle volume is (10)The density of the particle is then calculated as ρ_a = m_a/V_a. Current choices for the estimator are X_a = 1, m_a, and (k ≤ 1), where P_a is the pressure. Taking X_a = m_a leads to the common choice V_a = m_a/ρ_a. Actually, the same volume element is obtained with X_a = 1 when the mass of the particles is the same. On the other hand, the choice helps to suppress the surface tension that appears at contact discontinuities.

Here we want to discuss another option for X_a, not considered in the aforementioned papers, which could also be of value for SPH users. We choose (11)and put it into Eq. (10): (12)which for p = 1 becomes, (13)where the summation underneath is simply the kernel normalization condition for the standard volume choice m_a/ρ_a. Therefore, the volume element V_a given by Eq. (13) comes after re-normalizing m_a/ρ_a. We note that if , then Eq. (13) leads to the identity ρ_a = ρ_a, as expected. Furthermore, one can take Eq. (12) as a recursive equation, which, jointly with ρ_a = m_a/V_a, allows to find the optimal ρ_a leading to an almost perfect kernel normalization after several iterations. That is, starting from an initial guess of the density, for example ρ_a = ∑ _bm_bW_ab, the value of V_a is computed in an explicit way using Eq. (12). This gives a new density ρ_a = m_a/V_a to be used in the next integration step, and so on. The consistency and robustness of this procedure has been checked with the tests described in Sects. 5, 6, and in Appendix C.

The impact of changing the volume element is highlighted in the hydrostatic square test described by Saitoh & Makino (2013) and that we reproduce in Sect. 5.2.

3.2. Choosing smoothing length and kernel exponent: equalization

The use of the sinc kernels allows us to work with continuum adaptive index n(r,t). Interestingly, the smoothing length h(r,t) and kernel index n(r,t) can be calculated jointly with the density estimation. Because of the large number of involved variables, ρ_a, h_a, n_a, V_a, and Ω_a, this point deserves some discussion.

Firstly, the volume estimator X_a is updated only when the global iteration has finished. Therefore, X_a is handled explicitly so that the coupling with other variables during the current iteration is avoided. Secondly, self-consistent values of h_a, n_a, ρ_a, and Ω_a are simultaneously calculated with a Newton-Raphson iterative scheme, that we named equalization (because it equalizes the resolution at the post-shock tail of shock waves) and is described in García-Senz et al. (2014). We summarize it here for the sake of completeness:

• 1.

  Choose a trial value of h_a, as well as the baseline kernel index n = n₀, at thebeginning of the integration step.

• 2.

  Calculate the density of each particle and its logarithm average over neighbors: .

• 3.

  Evaluate λ_a, the ratio between ρ_a and , which is taken as a local indicator of linearity and is quantified with the following:

(14)

• 4.

  Use λ_a to assign a new kernel index n_a according to:

(15)where Δn is the maximum allowed jump from n₀, λ_c ≃ 1 is a scaling parameter, and f(ξ): (16)

• 5.

  Solve Eq. (15) jointly with mass conservation:

(17)where is a constant set at the beginning of the simulation jointly with n₀, Δn, and λ_c. After this process, we obtain the values of h_a, n_a, ρ_a, and Ω_a. Obviously, setting Δn = 0 turns the equalization off and the evolution is computed with n_a = n₀. Typical values for the simulations presented in the following sections are n₀ = 5, 0 ≤ Δn ≤ 5, and λ_c = 0.5.

4. The hydrodynamics code SPHYNX

SPHYNX gathers all the advances on the SPH technique that have been presented in Sects. 2 and 3 and represents our flagship code for astrophysical applications. In the following, we present the most relevant details on its implementation and itemize the mathematical expressions of the SPH equations implemented in the Cartesian version of the code. The precise form of the equations as well as the notation follows the discussion made in the preceding sections. They are similar to those presented in Rosswog (2015b).

4.1. Work-flow and formulation

First of all, we specify the initial number of neighbors n_b and kernel index variables: n₀ (baseline value), Δn (maximum allowed index jump), and λ_c (scaling parameter). This is followed by the choice of the volume estimator X_a. Then, the preconditioning stage starts with the self-consistent calculation of the volume element V_a, smoothing-length h_a, kernel index n_a, density ρ_a, and grad-h/n term Ω_a (see Sect. 3.2).

• •

  Volume element (18)with, (19)For the volume elements we use X_a = (m_a/ρ_a)^p. From the experiments presented below, we saw that taking p ≃ 1 leads to the best results in most cases, but it has the disadvantage of overshooting density interpolations in contact discontinuities when p> 0.7. We overcome this problem by using a SPH-averaged value of the VE estimator, making it robust and allowing us to safely rise the exponent to p = 1 (see Sect. 5.2 for more details.)

• •

  Density equation (20)which we calculate jointly with Ω_a including both contributions: the grad-h and grad-n terms (García-Senz et al. 2014), (21)Once the preconditioning is completed, the hydrodynamic equations are evaluated.

• •

  Momentum equation (22)

• •

  Energy equation (23)

where is given by Eq. (7), but with the τ_ij,a coefficients computed using the generalized volume elements, Eq. (18), (24)The terms with superscript AV refer to the artificial viscosity acceleration and energy contributions. As detailed in the following subsection, we have slightly changed the implementation of these terms in order to make them fully compatible with the generalized volume elements.

4.2. The artificial viscosity

Regarding the inclusion of physics, SPHYNX incorporates, by default, an artificial viscosity algorithm to handle shocks, as well as routines for the calculation of the equation of state (EOS) and gravitational force. Heat transport by conductive and diffusive means is also included as a basic code unit. More specific routines dealing with nuclear or chemical reactions, ionization or more complex transport schemes can be modularly added to the code (see, e.g., the application of SPHYNX to simulate core collapse supernova, including a spectral neutrino treatment, in Perego et al. 2016). In the following, we explain with some detail the AV algorithm.

First of all, we take the viscosity Π_ab as in Monaghan (1997): (25)where is an estimate of the signal velocity between particles a,b and w_ab = v_ab·r_ab/ | r_ab |. It is common to take but in that case the volume element is implicitly assumed to be m_a/ρ_ab. Thus, the VE of particle a is not an unequivocally defined magnitude because it also depends on the density of the neighbor particle b. A more compatible option between AV and VE, which can be generalized to other VE, is obtained using , which leads to the following viscous acceleration, (26)where (27)and f_a,f_b are the Balsara limiters (Balsara 1995): (28)Expression (26) can be written as (29)with V_a = m_a/ρ_a, V_b = m_b/ρ_b. Unlike Π_ab, the magnitude is not divided by the density and the viscous force, , is symmetric with respect to any pair of particles, a,b. We note that the volume elements are now unequivocally defined when the mass of the particles is the same. The Balsara limiters work more efficiently if they are averaged as, (30)because it gives a lower AV than the arithmetic average in regions with strong shear flows. By default SPHYNX sets these limiters to: (31)which retains a residual viscosity to damp the smallest numerical fluctuations. Equation (31) works remarkably well in regions where strong shocks and instabilities cohabit, such as in the Triple Point test described below. In the case of strong shocks with no shear, Eq. (30) could be replaced by the standard arithmetic mean if necessary.

Following Springel (2010b), we used a constant α = 4/3 so that Π_ab remains close to the classical SPH artificial viscosity introduced by Monaghan & Gingold (1983). The contribution of the AV to the energy equation (Eq. (23)) is then, (32)Additional information on the physics included in SPHYNX, such as gravity, heat transport or alternative energy equations is deferred to Appendix D.

5. Two-dimensional tests

In this section we present the outcome of applying SPHYNX to several two-dimensional (2D) tests that have traditionally been problematic for the SPH technique. We show that the sinc kernels help against the appearance of pairing instability, that the new VE provide a better treatment of discontinuities than the often-used m/ρ value, and that, when applied jointly with IAD₀, there is an overall improvement in the results of the shock treatment in the Sedov and wall-shock tests, and also in the development of subsonic instabilities; in particular, the challenging cases of Kelvin-Helmholtz with a density jump of 8 and Rayleigh-Taylor with very weak gravitational field (g = −0.1). We also found a convergence on the Gresho-Chan test closer to those of Eulerian codes. Finally, we prove that our implementation suppresses the tensile instability, enabling mixing in tests like the wind-cloud interaction and the triple-point shock test.

Our results below suggest that improving volume conservation produces better results than using volume elements that do not ensure volume equipartition. At first glance, it also renders equalization useless because the advantages of using variable adaptive kernel indices are apparently obscured by the new VE. Nevertheless, equalization and volume elements work on different bases. As demonstrated in García-Senz et al. (2014), the use of a variable kernel index improves the interpolation of non-linear functions, even when they are estimated using integrals instead of finite summations. It is therefore expected that the equalization will be useful when both the number of particles, N, and neighbors, n_b, are high enough so that the error in interpolations is dominated by non-linear effects rather than by the evaluation of integrals as finite summations (see also Sect. 6.2). For a moderate amount of neighbors the use of the new VE given by Eq. (12) makes equalization unnecessary. Thus, unless explicitly stated, many hydrodynamic tests presented in this work have been carried out taking Δn = 0.

5.1. Pairing instability

Here we show a numerical experiment, similar to that described in Zhu et al. (2015), used to study the relationship between the kernel index and pairing instability. We set a sample of particles in a regular 2D lattice leading to a homogeneous density distribution. Then, we seed a random small fluctuation of the internal energy, taking the system out of mechanical equilibrium. Afterwards, we add a frictional dissipative force proportional to the velocity in the momentum equation so that the system returns to equilibrium after a few sound crossing times. The final equilibrium distribution is sensitive to the number of neighbors and to the exponent of the sinc function. This is summarized in Fig. 2, which depicts the distribution of the minimum inter-particle distance (q), normalized to its initial value (q₀), as a function of pairs (n,n_b) at several elapsed times. As we can see, a careful choice of the exponent n keeps the distribution of normalized particle distances closer to 1, making particle clustering difficult, thus avoiding the pairing-instability even at large n_b. For example, n = 3 keeps q/q₀ close to 1 when the number of neighbors is small (top-left panel of Fig. 2), even at large times, while it completely fails when n_b increases (top-right panel). Nevertheless, increasing the exponent n of the sinc kernel suppresses again the pairing instability. This numerical experiment points to n = 3, 5, 7 and n ≥ 10 as indicative values to handle with n_b ≃ 15, 25, 50, and 100 neighbors, respectively (in 2D). For equivalent numbers in 3D, for example n_b ≃ 60–120, currently used in SPH calculations, a sinc kernel with n ≃ 4–6 will be sufficient to suppress clustering in many applications (see also Fig. 3 in García-Senz et al. 2014). We provide a comparison between S_n and the Wendland C₃ in the bottom-right panel of Fig. 2. As we can see, the C₃ kernel works better than S₁₀ but the difference is small.

The previous test was done using the standard gradient estimator. Figure 3 depicts the results for the same numerical experiment, but carried out with IAD₀. The comparison between Figs. 2 (standard gradient estimator) and 3 (IAD₀) suggests that the integral scheme also helps with the pairing problem, as it keeps the particle distribution closer to q/q₀ = 1.

5.2. Static square of test

We consider a system of two fluids with different density but identical pressure: (33)The system is isobaric with P = 2.5 and the EOS is that of a perfect gas with γ = 5/3.

To carry out the simulations we use Eqs. (22) and (23) from Sect. 4, with constant kernel index n = 5, n_b = 40 neighbors, and a variety of volume elements. We tried two different initial settings: using particles with the same mass (and uneven spaced grid), as well as particles with unequal mass spread on a uniform lattice. The outcome of the simulations at times t = 0, 0.4, 1, and 2 is depicted in Fig. 4. For equal-mass particles and X_a = (m_a/ρ_a)^p (with p = 0.9), the code was able to keep the hydrostatic equilibrium during several crossing times (first row). The evolution calculated using X_a = 1, X_a = m_a, and (k = 0.5) (second row) led to the same wrong results in much shorter times. Rows 3 and 4 in Fig. 4 summarize the evolution using a homogeneous lattice of particles with different mass. Hydrostatic equilibrium is very well preserved using X_a = 1, , and X_a = (m_a/ρ_a)^0.9 (third row), while X_a = m_a (fourth row) does not give a satisfactory result.

In Fig. 5 we show the magnitudes ∑ _bV_bW_ab (left panel) and ⟨ Δr ⟩ = ∑ _bV_b (r_b−r_a)W_ab (right panel) along a 1D cut taken around y = 0.5, from x = 0.1 to x = 0.5 (center of the system) in the hydrostatic system described above after 20 integration steps. As can be seen in the left panel, choosing X_a = (m_a/ρ_a)^p(p = 1) gives the best results, leading to an almost perfect volume normalization, while the case with p = 0 (equivalent to X_a = 1) shows a large oscillation around the contact discontinuity. The case with p = 0.9 is also quite satisfactory although, a small oscillation is still present. The right panel in Fig. 5 suggests that an adequate choice of the volume element may improve even further the accuracy of the IAD₀ method to estimate gradients because, as we can see, a value of p ≃ 1 helps to keep ⟨ Δr ⟩ ≃ 0, which is the necessary condition for IAD₀ to exactly reproduce the gradient of linear functions.

Unfortunately, taking p = 1 has the unwanted side effect that density tends to overshoot close to the contact discontinuity. The overshooting in density drops the value of the estimator X_a = m_a/ρ_a, so that in the next integration step the density undergoes a further overshooting. Occasionally, the feedback between ρ_a and X_a might rise the density to unrealistic values after several time steps. This problem can be avoided by reducing the value of the exponent p but in that case the adequate optimal value becomes problem dependent. For example, p = 0.9 works well in the hydrostatic test depicted in Fig. 4 but it produces worse results for the Sedov test described later in Sect. 5. For that test p ≃ 0.6–0.7 was a better choice. A way to circumvent the overshooting problem is to consider (34)where ⟨ m/ρ ⟩ _a is the SPH average of the standard volume element. Although less efficient than X = (m/ρ)^p, this simple recipe enhances the robustness of the scheme so that the exponent of the estimator can be safely raised to p ≃ 1, independently of the problem to be simulated. This procedure is especially well suited to handle contact discontinuities hosting large density contrasts, and it has been checked in the blob and Evrard tests described in Sects. 5.6 and 6.3. A more quantitative discussion on the convergence rate of the estimator X_a as a function of p and the density contrast is given in Appendix C.

5.3. Shocks

5.3.1. Point explosion

The study of the propagation of a point-like explosion in more than one dimension is of great usefulness for checking hydrodynamics codes. We have carried out several simulations with SPHYNX to explore the impact of different combinations of volume element and level of equalization on handling shock waves. The explosion was triggered by depositing a considerable amount of energy inside a Gaussian surface of characteristic width σ = 0.05 located at the center of a uniform lattice of size [0,1] × [0,1]. The initial density profile is ρ(r,0) = 1 and the total injected energy E = 1. All simulations used N = 256 × 256 particles and the initial number of neighbors was set to n_b = 50. The different parameters used in the simulations are summarized in Table 3.

Figure 6 shows the results of the Sedov test for several combinations of volume elements V_a and equalization Δn at t = 0.1. Models S₁, S₂, S₃, and S₄ were calculated using and different values for the exponent p. All models are quite satisfactory but model S₂, calculated with p = 0.7 and Δn = 0, leads to the best results. The density peak is closer to the analytic value and the pressure profile close to the origin behaves well. Model S₄, calculated using with p = 1, also provides good results while models S₁, S₃, with the choice p = 0 (standard VE) lead to lower peak values and oscillating pressure tails that depart from the analytic result.

The lower panels of Fig. 6 depict the profile of the moduli | ⟨ Δr ⟩ | _a = | ∑ _bV_b (r_b−r_a)W(h_a,n_a) | (left panel) and that of the normalization condition ∑ _bV_bW(h_a,n_a) (right panel). These results are encouraging because they strongly suggest that the new volume elements not only enhance volume conservation (lower right panel of Fig. 6) but also make gradient estimation more accurate, because Eq. (9) is better fulfilled, ensuring that the IAD₀ approximation is, in fact, close to the full IAD in terms of accuracy (lower left panel of Fig. 6). Again models S₂ and S₄ lead to the best results whereas S₁ and S₃ show worse volume conservation and lower fulfillment of Eq. (9).

It has been suggested that other choices of VE could be of interest for handling specific problems. In particular, it has been claimed that the pressure-based estimator (k< 1) could be well suited to handling contact discontinuities (Hopkins 2013). We carried out one simulation using that estimator, with the recommended value k = 0.05 (Rosswog 2015a). We found that the results with the pressure-based estimator were not as good as those of models S₂ and S₄ of Table 3. In Fig. 7 we compare the results of using the estimator with those of S₁ and S₂. In the right panel it is clear that the volume normalization is not so well preserved as in model S₂, obtaining values similar to those of using standard VE (model S₁). Moreover, it leads to a lower peak density value and a spurious precursor shock is clearly seen ahead of the shock in the left panel. Although small, such precursor shock is a numerical artifact which was also reported by other authors (Rosswog 2015a). For that reason, hereafter we have limited the choice of volume elements to those given by Eq. (12)⁵.

5.3.2. The wall-shock test

Along with the point explosion (Sedov test), the so called wall-shock (or Noh) test is also used to check the performance of multidimensional hydrodynamics codes. Unlike the Sedov test, the Noh experiment is an implosion towards a geometrical center. The interaction of the converging inflow of gas with the stagnated material at the central volume provokes the formation of a shock-wave moving outwards. The Noh problem has an analytic solution (Noh 1987) to compare with. It is, however, not an easy test owing to the large jump in density across the shock front, 16 in 2D, 64 in 3D, for γ = 5/3. SPH codes traditionally had difficulties in handling this test because: a) many particles are needed to correctly reproduce the shock region (especially in 3D); and b) close to the center the internal energy shows a pronounced spike and, consequently, the density abruptly drops to keep the pressure constant. In fact, both problems are connected to the use of the AV and it is not clear to what extent a change in the SPH formalism may alleviate these shortcomings (Brookshaw 2003).

A sample of 256² equal-mass particles was distributed in an bcc lattice with circular perimeter. The radius of the initial configuration is R = 1 and the density is homogeneous with ρ = 1 (except at the outer edge of the system). A radial profile of v_r = −1 was imparted to all particles at t = 0, so that the system implodes. The evolution was tracked with SPHYNX until t = 0.5, when a clear steady shock moving outwards is formed. The results, calculated with n = 5, p = 0, and p = 0.7 for the VE and n_b = 50, have been compared to both, the analytic predictions and to the output of GADGET-2 for the same initial model.

Table 4 and Figs. 8 and 9 summarize initial parameters of each simulated model and the main results of the simulations. The density profiles at t = 0.5 for the three calculated models, WHS₁, WHS₂, WHS_G are shown on the leftmost panel of Fig. 8. A sharp shock, moving towards the right, is clearly seen at r ≃ 0.17. The density jump across the shock front is Δρ ≃ 16, hence close to the analytic value. Close to the center of the configuration there is a pronounced dip in density caused, as expected, by the use of the AV. The SPHYNX calculation with p = 0.7 (model WHS₂ with green dots) gives slightly better results at the shock position than WHS₁ and WHS_G because the density profile is steeper and the maximum density is closer to the analytic value of 16. We note that all three calculations do show density oscillations at the shock front and that the GADGET-2 simulation is more blurred. The color-maps of density for models WHS_G and WHS₂ are also shown in the central and right panels of Fig. 8. As we can see, the shock region is sharp and well defined in the SPHYNX calculation whereas the simulation with GADGET-2 shows stronger oscillations close to the shock location and an inhomogeneous particle distribution close to the origin.

Additionally, model WHS₂ has both a better volume normalization than WHS₁ (i.e., closer to 1), as shown in Fig. 9 (left panel) and a better behavior of | Δr | (i.e., closer to 0) at the shock location (right panel).

5.4. Fluid instabilities

Instabilities play a central role in hydrodynamics because they are usually connected to shear mixing and turbulence. In the cosmos, instabilities are especially important because the large Reynolds numbers involved in astrophysical processes make these systems prone to turbulence. Particle- and grid-based codes have been applied with different levels of success to simulate the evolution of hydrodynamical instabilities such as Kelvin-Helmholtz (KH) or Rayleigh-Taylor (RT) instabilities. Here we re-visit them with SPHYNX.

5.4.1. Kelvin-Helmholtz

Grid-based codes are, for the most part, almost free of numerical viscosity, leading to a good match between simulations and the analytic predictions during the linear phase of growth of the KH instability (Junk et al. 2010). On the other hand, SPH codes are Galilean invariant and do not suffer from numerical diffusion, but intrinsically have more numerical viscosity. It has been pointed out that SPH may not be appropriate for handling fluid instabilities across contact discontinuities with large density jumps (Agertz et al. 2007). Although there have been proposals to enhance the modeling of these systems (Price 2008), the SPH calculation of shear flows with large density contrast is still an open question.

We ran two sets of models simulating the growth of the KH instability around the boundary layer separating two flows with moderate and large density contrasts, respectively. Each set was in turn calculated using two values for the volume elements namely p = 0 and p = 0.7 in Eq. (12). The initial setting was the same as in McNally et al. (2012), where three stratified fluid layers inside a 2D box of size [0,1] × [0,1] were considered. The fluid layers span for y ≤ y₁, y₁ ≤ y ≤ y₂ and y ≥ y₂ with densities ρ₁, ρ₂, and ρ₃, respectively. The adopted values for y₁, y₂, ρ₁, ρ₂, and ρ₃, as well as the number of particles used in the simulation, are shown in Table 5. Periodic boundary conditions were implemented at all sides of the box.

A velocity v_x = 0.5 is given to the central strip, whereas the rest of the box moves in the opposite direction with v_x = −0.5. Prior to the calculations, the density and v_x distributions were smoothed following the method described in McNally et al. (2012; their Eqs. (1) and (3)). Thus, the growth-rate of the instability obtained with SPHYNX can be compared to the templates obtained by McNally and collaborators using the PENCIL Code⁶, a state-of-the-art hydrodynamics code of Eulerian type (Brandenburg & Dobler 2002). The pressure is set to P = 2.5 everywhere, with γ = 5/3. The fluid layer is in almost vertical equilibrium except for a small seeded perturbation in v_y given by, (35)with w₀ = 0.01. The velocity perturbation has therefore a wavelength λ = 0.5 so that the box hosts two complete waveforms.

Models KH₁ and KH₂ in Table 5 were calculated assuming a density contrast of two, whereas an initial jump of 8 was imposed to models KH₃ and KH₄. To build the initial models, a sample of 256 × 256 particles were first spread in a squared lattice and then stretched in the vertical direction so that the ensuing density profile adapts to that of Eq. (1) in McNally et al. (2012). This setting naturally leads to a density contrast of two for equal-mass particles and provides the necessary smoothness around the separation layers between the fluids. To build the high-density contrast models, the initial lattice was stretched until the density ratio was . The mass of the particles was then made proportional to the density of the stretched grid, so that a smooth profile with a density contrast of 8 was achieved.

The Balsara limiter f_a given by Eqs. (28), (30), and (31) was applied to all models to reduce the shear viscosity. For these initial conditions, previous simulations with the traditional formulation of SPH predict the growth of the KH instability only for the low-density contrast case (Agertz et al. 2007; Junk et al. 2010). When the density contrast rises well above 2 it is necessary either to add an artificial heat flux to keep isobaricity (Price 2008) or to redefine the volume elements (Saitoh & Makino 2013; Hopkins 2013). In particular, it was shown that the formulation which uses volume elements linked to pressure is able to reproduce the KH instability across high density jumps. Also, codes which directly use a smoothed pressure in the SPH equations, such as the PSPH formulation described in the Appendix of Hopkins (2015), can also handle with high density contrasts. Here we show that the IAD approach combined with new volume elements, which preserve kernel normalization, is also able to successfully simulate the KH instability with high density jumps.

Figure 10 shows a color-map of density of the evolution of models KH₁ and KH₂ at times t = 1.5 and 2.5. We see that there is a growth of the instability, qualitatively similar to the results obtained using other novel formulations of the technique such as the PSPH scheme with ≃200 neighbors (see for example Fig. 19 in Hopkins 2015). The results of SPHYNX are also very similar to those obtained with the recent CKRSPH formulation (see for example the upper-row of figure 17 in Frontiere et al. 2017). Model KH₂, calculated using the new volume elements, appears to evolve slightly faster than KH₁, computed with the standard choice, VE = m/ρ, showing more evolved structures in the non-linear stage. Unlike the PSPH scheme, our method does not estimate the pressure by kernel smoothing nor does it incorporate any artificial flux of heat to smooth the pressure.

Figure 11 depicts the mode-amplitude evolution of v_y calculated taking the Fourier transform (FT) of the v_y field. The FT was calculated with the expressions given in McNally et al. (2012) so that a comparison with the amplitude growth calculated with the PENCIL code (continuum-red line) and the NDSPMHD code, described in Price (2012), (dashed lines) is straightforward. The results are encouraging as SPHYNX reproduces the KH growth closer to the reference simulation of the PENCIL code than the NDSPMHD code, the latter even having a factor 2 higher resolution than our simulations with SPHYNX (see also Fig. 7 in McNally et al. 2012, and comments therein). Our results are close to those obtained with the PSPH scheme with the same resolution. In Fig. 11 we also show the L₁ errors of each model with respect to the reference values, calculated as: , where A^r and A^s stand for the reference and simulated amplitudes, respectively. These results are an indication of the importance of using an accurate gradient evaluation, this being the most relevant difference between SPHYNX and NDSPMHD.

Model KH₂, calculated with p = 0.7, exhibits better volume normalization than KH₁, computed with p = 0, as shown in Fig. 12. That figure depicts the value of the normalization condition ∑ _bV_bW_ab as a function of the density of the particles. For both runs, most of the particles cluster around the expected value of 1, but the summation falls below one in the low-density region and exceeds it in the high-density region. Nevertheless, the dispersion around the correct value is considerably lower in the case with p = 0.7, showing its greater capability to achieve equipartition in disordered particle distributions.

Figure 13 summarizes the simulations with an initial density ratio ρ₂/ρ₁ = 8. The upper and lower panels were calculated with p = 0 and p = 0.7 in Eq. (12), respectively. The snapshots correspond to times t = 1.5 and 2.5. First of all, we see that despite the large jump in density, the instability is able to develop in both cases. These results, therefore, suggest that the use of the IAD scheme to calculate gradients improves by itself the quality of the simulations of hydrodynamic instabilities. Although the instability evolves slightly faster when p = 0.7, and shows more structure, it is probably affected by the noise introduced by the VE estimator (m/ρ)^p when the mass of the particles is not constant. As in the precedent tests, the model calculated with p = 0.7 has the best kernel normalization properties (not shown in the figures).

5.4.2. The Rayleigh-Taylor instability

In García-Senz et al. (2012) it was shown that the IAD₀ scheme is able to simulate the gross features of the growth of the RT instability for initial perturbations as low as Δv_y = 0.01 in the velocity field. In this section we present clear proof that SPHYNX is able to cope with the growth of the RT instability also when the gravitational force is small.

Our initial model is similar to that described in Springel (2010a), where the numerical experiment takes place in a box sizing [0.5 × 1.5]. The lower and upper halves of the box are filled with equal-mass particles. The particle distribution is then arranged in the vertical direction so that ρ_u = 2 (upper region) and ρ_d = 1 (lower region).

The contact discontinuity around the interface was smoothed using: (36)with Δy = 0.0083 and y₀ = 0.75. The integration of the hydrostatic equilibrium equation , where g is the gravitational acceleration, with an ideal EOS with γ = 1.4 gives the pressure profile in the vertical direction, (37)with P₀ = 2.5. We used 256² particles and an initial number of neighbor particles n_b = 50. The boundary conditions are periodic at the left and right sides of the box but reflective on the bottom and the top. The main features of the models used in this test are summarized in Table 6.

The fluid was initially at rest everywhere, except for a small, single-mode perturbation applied to the vertical velocity field, (38)where w₀ = 0.0025. In Fig. 14 we show the density color-map of the evolution of the instability when g = −0.5 is adopted. As we can see, the instability grows at a good rate and it enters into the non-linear regime at t> 3.5. The contact discontinuity between both fluids is well marked, but not totally sharp, owing to both the smooth initial conditions and to the spread of the physical magnitudes over the characteristic smoothing length h. During the non-linear stage, the calculation with p = 0.7 (model RT₂ in Table 5) displays more structure than that of model RT₂ calculated with p = 0. At t = 7.5 the vertical extension of the unstable region is nearly the same but model RT₂ shows a richer structure inside the mushroom-like instability.

The impact caused by the choice of VE is much more marked when a smaller gravitational force, g = −0.1, is adopted. This was the case reported in Springel (2010b) concerning the simulations of the RT instability with the AREPO code. As far as we know, there are no satisfactory calculations on this same scenario using the SPH technique. The combination of a small initial perturbation plus a tiny gravitational force makes the problem difficult because the E₀ errors and the tensile instability conspire to totally suppress the growth of the instability (Valdarnini 2012). The upper panels in Fig. 15 depict the evolution of the instability as simulated by SPHYNX taking standard VE (i.e., p = 0 in Eq. (13)). As we can see, the instability is able to grow, although the growth-rate is slow and there is a clear lack of structure during the non-linear stage. The lower row of panels in the same figure depict the evolution when VE are calculated with p = 0.7. In this case there is a clear boost in the development of the instability as it grows faster and shows a richer structure in the non-linear regime.

Quantitative numbers on this test can be extracted from the evolution of the points located on the tip/bottom of the bubbles/spikes. Additionally, the numerically inferred terminal velocity of the vertex of the bubble can be compared with the analytical results by Goncharov (2002). The results are summarized in the Fig. 16, where the left panel indicates the bubble/spike evolution of the farthest point achieved by the lighter and heavier fluids (bubble and spike respectively), and the right panel indicates the evolution of the velocity of the bubble. As we can see, the evolution of the bubble/spike sample is faster when the VE are calculated with p = 0.7 and p = 1 (in this last case the averaged estimator, Eq. (34), was used). We note that at t ≃ 15, the tip/bottom of the bubble/spike are close to colliding with the box limits. Consequently the rising velocity of the bubble velocity was computed in the interval 0 ≤ t ≤ 15. As we can see, the terminal velocity of the bubble matches well the analytical estimation by Goncharov (2002), (39)where A_t = (ρ₂−ρ₁)/(ρ₂ + ρ₁) = 1/3 is the Atwood number, k = 4π is the wave-number of the applied single perturbation and C_g = 3 in two-dimensions. The terminal velocity calculated using Eq. (39) (solid pink horizontal line in Fig. 16) is close to those obtained using SPHYNX.

Although our results are not as good as those obtained with the code AREPO for the same initial conditions (see Fig. 35 in Springel 2010b, for a comparison), they are encouraging, as they suggest that the simulation of the RT instability in a tiny gravitational field is also affordable for SPH codes when an E0-suppressing technique like IAD₀ is used jointly with equipartition-preserving volume elements.

5.5. The Gresho-Chan vortex

The simulation of a stationary vortex that is in stable equilibrium is a very demanding test for any hydrodynamics code. This is a fundamental test for the accuracy of the numerical scheme, the preservation of symmetry and the conservation of angular momentum. The Gresho-Chan vortex has been especially problematic for SPH, which, in spite of having excellent angular momentum conservation, shows very poor convergence (Springel 2010b) due to numerical deviations from the initial conditions that degenerate the stability of the vortex within short timescales. This is specially relevant in the simulation of self-gravitating disks, where the centrifugal force and pressure gradient should be (up to a certain extent) balanced. Therefore, numerical noise can trigger deviations from this equilibrium configuration that lead, for example, to an artificial fragmentation of the disk.

Several recent works (Dehnen & Aly 2012; Zhu et al. 2015; Rosswog 2015a) have managed to keep the vortex in steady state, with low dispersion, during more than one orbit-cycle. In his work, Rosswog (2015a) proves that the use of high-order kernels is crucial to obtain accurate results, and in combination with IAD₀ and an improved artificial viscosity with a noise dissipation trigger, leads to a roughly linear convergence rate.

In this test, we use the common initial conditions for the azimuthal velocity and pressure profile: (40)(41)where ψ = r/R₁, R₁ = 0.2, v₀ = 1, P₀ = 5, for an ideal gas with γ = 5/3. This setting corresponds to a low Mach number test with . We use N = 256² particles evenly distributed in a square lattice of size [0,1] × [0,1] so that the density is uniform with value ρ = 1. The artificial viscosity parameter in Eq. (25) is set to the current value α = 4/3 and the Balsara coefficients estimated with Eqs. (28), (30), and (31). We note that because the vortex is in steady conditions, the divergence of the velocity should vanish everywhere and f_a ≃ 0. As in other simulations, however, we allow a low level of AV to control the numerical noise and to reduce the dispersion of the particles. For this test, the particular choice of the volume element is not relevant because the particle distribution is homogeneous due to the constant density, hence the kernel is properly normalized at all times, independently of the exponent. Therefore, we take the nominal value of p = 0.7 in Eq. (12) in all calculated models.

We carried out several simulations of the vortex evolution for different choices of the kernel index n and initial number of neighbors⁷, n_b(2D). We also calculated models with the Wendland C₆ kernel which, according to Dehnen & Aly (2012) and Rosswog (2015a), was especially suited to handle this test. The quality of the simulations and convergence rate were checked with the parameter where v^t and v^s stand for the theoretical and the simulated value of velocity, respectively, and the summation runs over all particles.

In Fig. 17, we show the profile of the tangential velocity at time t = 1 (roughly the elapsed time needed by the vortex to complete one turn) for three kernel choices and different number of neighbors (n_b). Using a low n_b (≃ 30, in 2D) the best results were obtained with the sinc (n = 5) whereas the choice n = 6.315 or the Wendland C₆ leads to qualitatively similar results. The situation is reverted when n_b = 50, where the C₆ interpolator gives a slender velocity profile and the lowest L₁ value. However, the sinc(n = 6.315) also shows low dispersion profiles and L₁ values. Increasing the number of neighbors to n_b = 70 (≃ 500 in 3D) reduces the dispersion, as expected. The C₆ provides the lowest L₁ value here but closely followed by the sinc (n = 7).

Figure 18 summarizes the functional dependence of the estimator L₁ with respect to the index n (left panel) and number of neighbors n_b (right panel). Looking at the left panel we see that for n_b ≃ 20 the sinc (n = 4) (similar to the M₅ and C₂ kernels) gives the best results. Increasing the number of neighbors moves the minimum value of L₁ to the right (higher n values). At n_b ≃ 30, the index n = 5 becomes the optimal choice. Setting the initial number of neighbors to 70 and 80 demands larger kernel indices (n ≃ 7–8) at the position of minimum L₁. The right panel of Fig. 18 explores the convergence rate of the Gresho-Chan vortex test with respect to the number of neighbors at a constant total particle number N = 256². When n_b ≤ 30 it seems that the low-order interpolators with n ≤ 5 give the lowest L₁ values. Nevertheless, above n_b ≥ 40 it is necessary to raise the index of the sinc kernel, or use the Wendland C₆, to achieve convergence. The C₆ gives the lowest L₁ values when n_b ≥ 40. The curve referred as nmix (dotted black line) is a fit which takes advantage of the optimal n index at each n_b. Both, the nmix profile and the C₆ curve can be used to estimate the convergence rate in this test. These profiles show a rapid initial decline: , followed by a more gentle reduction of the convergence rate .

Figure 19 shows the convergence rate of L₁ at t = 1 when the one-dimensional (1D) equivalent number of particles, N_1D, is varied in the range 50 ≤ N_1D ≤ 400. In this study we have only considered the sinc kernels, but with several exponents. The convergence rate is explored using the following (n,n_b) pairs: (5,30), (6.315,50), (8,80) and (10,100). As can be seen, the high-order sincs with n = 8, n = 10 achieve a linear convergence of L₁, provided that a sufficient number of neighbors is taken. Moreover, the convergence rate is larger than those reported in Springel (2010b) () or in Dehnen & Aly (2012; ) and closer to the value obtained with Eulerian codes, (Springel 2010a); see also Valdarnini (2016) where the ability of the IAD₀ scheme to reduce sampling errors is proved in a similar subsonic test, but with Mach number as low as ℳ = 0.02.

5.6. Interaction of a supersonic wind with a cold cloud of gas

Popularly known as the “blob” test (Agertz et al. 2007), this problem has challenged SPH codes for a long time. The basic setting of this test gathers many pieces of physics, such as strong shocks and mixing due to the KH and RT instabilities in a multi-phase medium with a large density contrast. The initial configuration consists in a dense spherical cloud of cold gas embedded in a hotter ambient medium. The cloud is initially at rest while the ambient background (the wind) moves supersonically. The wind-cloud interaction generates a bow shock and, in short time, the cloud is fragmented and mixed with the background owing to the combined effect of ablation and KH and RT instabilities.

Modern grid-based codes are able to handle this scenario and all of them agree in the gross features of the evolution. Nevertheless, simulations using the SPH technique have historically had difficulty disrupting the cloud because of the poor development of the hydrodynamic instabilities that drive the erosion of the cloud (Agertz et al. 2007). As a matter of fact, the challenge introduced by this problem has led to interesting advances of the SPH technique in recent years, such as the pressure, PSPH, formulation of the SPH equations by Saitoh & Makino (2013) and Hopkins (2013), and the conservative reproducing kernel method, CRKSPH by Frontiere et al. (2017). As shown in these works, it is feasible to suppress the tensile-instability acting at the contact discontinuity between the two fluids, leading to a better agreement between SPH and grid-based codes.

Here we show that the new volume elements given by Eq. (12) also suppress the tensile-instability at the contact discontinuity. The code SPHYNX is thus able to cope with the blob test without leaving the density-formulation of SPH. Our initial setting is similar to that described in Agertz et al. (2007), although restricted to two dimensions where the dense cloud is no longer a sphere but a circle (a cylinder in 3D). The wind is simulated with N_w = 256² particles spread following a glass-like stable (previously relaxed) distribution in a box sizing [1,1/4]. The cloud is reproduced with N_C = 68² particles, spread in a regular square lattice tailored in a circle of radius R_C = 1/40, and centered at coordinates (1/8,1/8). All particles have the same mass. The features of the wind and cloud at t = 0 are ρ_w = 1,P_w = 1,u_w = 3/2,v_w = 2.7 and ρ_c = 10,P_c = 1,u_c = 3/20,v_c = 0, respectively, where the wind velocity v_w is supersonic with Mach number ℳ = 2.7. As the particles in the cloud are fully ordered we add a small random initial radial velocity, with .

We carried out two calculations of the blob test using SPHYNX: with the standard VE (model W₀) and with the new VE (model W₁), as well as with GADGET-2, (model W_G). Because of the large density contrast between the wind and the cloud the model W₁ was calculated with the averaged estimator X = ( ⟨ m/ρ ⟩ )^p and p = 1.

In Fig. 20 we show several snapshots of the simulations at normalized times (t/t_KH = 0.5,2,3), where t_KH is the characteristic KH time, as defined in Agertz et al. (2007). As we can see, the fragmentation and mixing of the cloud is only attained in model W₁, calculated with p = 1. Models W₀ and W_G give worse results, but we note that even in that case SPHYNX with p = 0 yields a slightly larger fragmentation of the cloud than GADGET-2 (see Fig. 21). Our results are also in good agreement with those obtained by Frontiere et al. (2017) using the CRKSPH method (see their Fig. 21). The large suppression of the tensile instability in the W₁ calculation is the characteristic that allows the fragmentation of the cloud. This is clearly visible in the detailed density snapshots shown in Fig. 22. It is therefore evident that the long-standing tensile instability problem in SPH is bounded to kernel normalization. Having a good partition of unity considerably reduces the strength of that instability.

The evolution of the surviving fraction of the cloud as a function of time is depicted in Fig. 21. The expected outcome is that the surviving fraction of the cloud reduces to 0 in about t/t_KH ~ 2.5, as obtained by grid codes (see, e.g., Fig. 25 of Hopkins 2015). As we can see, the W₁ calculation leads to almost complete destruction of the cloud in t/t_KH ≃ 4, while in the other two calculations hardly 50% of the cloud is mixed. Nevertheless, unlike in the grid-based calculations, the cloud is not totally mixed with the ambient gas as ≃ 10% of it remains unmixed in the W₁ simulation. The difference may arise from the different dimensionality of the calculations, because in 2D the cloud is a cylinder whereas in 3D the shape of the blob is spherical. However, a calculation of the same model in 3D with 3 × 10⁶ particles led to similar results. Actually, the complete dissolution of the cloud is only attainable if some heat-transport is included in the scheme (Hopkins 2015).

5.7. Shock plus vorticity: the triple-point shock test

The similarities among the results obtained by SPHYNX and the CRKSPH method, as found in the KH and the Blob tests above, are even more evident when considering the so called triple-point shock test. This is a shock-tube-like setting where three materials with different densities and pressures are put into contact through discontinuity lines. The first region contains a high-density, high pressure material, located in the vertical band on the left of the first snapshot in Fig. 23 (see also the sketch by Frontiere et al. 2017, shown in their Fig. 24). That region is in contact with two low-pressure horizontal bands, the upper band being considerably less denser than the lower band. The particular point where the three discontinuity lines intersect is the triple point. The vertical high-pressure band on the left launches a shock through the lower-pressure horizontal bands. Because of the pronounced difference in the speed of sound, the shock moves faster in the upper low-density, low-pressure region. As a consequence, a shear is induced around the line separating both low-pressure materials, and soon the KH instability develops, rolling-up the interface between these two regions. At late times, the reflected shocks at the boundaries of the box also induces the Richtmyer-Meskhov instability when they cross the interfaces between the different materials.

Our main aim is to simulate the triple-point shock scenario and to compare the results obtained with SPHYNX to those given by the CRKSPH method. A comparison among the performances of CRKSPH, the standard SPH formulation and the ReALE (Reconnecting arbitrary Lagrangian Eulerian) method by Loubère et al. (2010) is provided in the work by Frontiere et al. (2017) and we refer the reader to that paper. Their results show that for the triple-point shock scenario the CRKSPH scheme gives the best results.

The setting of the initial model, particle number in each region and EOS is exactly the same as described in Frontiere et al. (2017). The evolution of the system with SPHYNX and two choices for the VE, p = 0 and p = 1 in Eq. (34), was tracked until t = 8. We show in Fig. 23 the density color-maps for both cases at times t = 1, 3, 5, 7. As we can see, the results with SPHYNX and p = 1 are similar to those obtained with the state-of-the-art CRKSPH scheme. The main features of the shock wave are well simulated for both, the p = 0 and p = 1 cases, but the amount of structure seen in the calculation with p = 1 at t = 5 and t = 7 is larger than that of the standard VE choice. It is worth noting that our recipe for implementing the artificial viscosity is relatively basic and there is still some room for improvement in this area, such as, for instance, the use of time-dependent AV coefficients (Cullen & Dehnen 2010; Rosswog 2015a).

6. Three-dimensional tests

In this section we present three test cases where SPHYNX is applied to three-dimensional (3D) scenarios. We show here that the combined use of the new VE, the sinc kernels, and IAD₀ leads to better results when trying to obtain an homogeneous density field, decreasing discretization errors. In the following, we re-do the Sedov test, now in 3D, where the new VE stand out by better handling shock-waves and the equalization enhances the outcome in simulations with high number of neighbors. Finally, the collapse of an isothermal cloud tests the implementation of SPHYNX coupled with gravity evaluation. The results show that the new VE enhance the kernel normalization, and consequently the overall interpolations, leading to a better shock treatment when compared with standard implementations.

6.1. Approaching a uniform scalar field

The initial models used in many SPH simulations are often built so that the density and/or the pressure are as homogeneous as possible. These systems should ideally remain in equilibrium during many sound-crossing times. Although very common, arranging particles in ordered lattices is not the best option because the particles move off the lattice and introduce noise. In this section we discuss the ability of SPHYNX to build a uniform density/pressure field, which is obtained after relaxing a pseudo-ordered initial distribution of particles. Several features of the resulting glass-like structure are discussed as a function of the VE choice. To obtain a homogeneous density profile we follow a similar procedure such as that explained in Rosswog (2015a).

We first spread N = 128 000 particles in an ordered 3D squared bcc lattice. Afterwards, their position is displaced at random in each direction, with maximum amplitude 40% of the size of the smoothing length. Then, we allow the system to evolve keeping internal energy constant, at u₀ = 1, to regain equilibrium. During the relaxation, the velocities are set to zero but the particles are displaced with the simple recipe: Δr = 10^-3h (a_c/ | a_c |) where h is the smoothing-length and a_c stands for the acceleration vector. The value of h was the same for all particles, so that the number of neighbors is n_b ≃ 100 with low dispersion. During the evolution, the L₁ error estimators were monitored for several magnitudes, such as density, partition of unity, and ⟨ Δr ⟩. These are defined as: We note that as P = (γ−1)u₀ρ, with γ = 5/3, the error estimator L₁(ρ) also provides the dispersion in pressure.

The L₁ values after 800 iterations are summarized in Table 7. All settings of the VE lead to stable glass-like structures with small departures from homogeneity (L₁(ρ) ≤ 10^-3). Model D, calculated with X_a = ( ⟨ m_a/ρ_a ⟩ )^p; (p = 1), displays the lowest L₁(ρ) values. As expected, the calculation with X_a = (m_a/ρ_a)^0.7 leads to the best partition of unity but, paradoxically, is the most inhomogeneous of the sample. Model C, calculated with a negative exponent p = −0.7 is very similar to model A computed with the standard VE choice (p = 0.0). Because the internal energy was kept constant during the relaxation, a negative p value in Eq. (11) is equivalent to taking a VE linked to pressure. We therefore conclude that SPHYNX is able to handle homogeneous systems in pressure equilibrium for a wide range of the parameter p. For these systems, the best choice of VE is that given by Eq. (34) with p ≃ 1.

6.2. Point explosion

The description of an explosion in the 3D space is usually a challenge for any hydrodynamics code. The lower resolution (in comparison with 2D calculations) makes it difficult to capture the density, pressure, and velocity at the peak of the blast and also to correctly reproduce the post-tail values of these variables. We have investigated the ability of our code to cope with a point-like explosion in 3D and compared the results with those obtained using GADGET-2.

We have run two sets of models labeled as { S₅,S₆,S₇,S₈ } and { S₉,S₁₀,S₁₁,S₁₂ } in Table 3. The first set consists of low-resolution calculations whereas the second set of simulations has not only more resolution, but also uses more neighbors to make the interpolations.

Figure 24 summarizes the results for the set of low-resolution models. As in the 2D simulations, models S₆ and S₈ in Table 3, with the new volume elements, lead to the best results. Both models give a higher density peak and a better post-shock evolution of density and pressure. They also exhibit an enhanced volume normalization and | Δr | has a better behavior across the wave front. The model calculated with p = 0 and Δn = 5 (S₇) is the third-best case, followed by model S₅ (standard VE setting p = 0,Δn = 0). Figure 24 also includes the density and pressure profiles (light blue lines) calculated using GADGET-2 for the same initial setting. Its density profile shows a similar post-shock tail evolution as model S₆, but lower peak value and the post-tail pressure profile diverges considerably from the analytical curve.

Figure 25 is the same as Fig. 24 but for models { S₉,S₁₀,S₁₁,S₁₂ } calculated with a larger number of particles and taking two times more neighbors in the summations. Therefore, these models are characterized by both a higher resolution and a less noisy estimation of gradients. According to our calculations, model S₁₂ computed using p = 0.7, Δn = 5 is the one closer to the analytic expectations. The density peak of the blast was the highest and at the correct position, while the values of density and pressure close to the origin behave well. Further, model S₁₂ also has the best volume normalization and minimum | Δr | at the position of the shock front. The profile of the sinc kernel index n of model S₁₂ is shown in Fig. 26. We can see that the highest value of the exponent (n ≃ 9.8) is reached at x = 0.3 behind the density peak. These results are closely followed by model S₁₀ (p = 0.7 and no equalization), and, in decreasing quality, by S₁₁ (p = 0 with equalization), and S₉ (p = 0 and no equalization).

The results above suggest that: 1) the use of the new volume elements is beneficial when handling shock-waves; 2) when using these VE, the equalization algorithm (i.e., considering n(r,t), variable kernel exponents) is only able to enhance the results when the number of neighbors is high (as in the cases depicted in Fig. 25).

6.3. Collapse of an isothermal cloud

A relevant test in Astrophysics is the numerical study of the gravitational collapse of a gaseous configuration, also known as the Evrard test (Evrard 1988). This test includes gravity and has been used many times to check hydrodynamics codes (see for example Springel 2010a and references therein).

We used here the initial configuration described in Springel (2010a) consisting of a sphere of gas with unit radius and mass. The density profile at t = 0 is (45)where R = 1 and M = 1. The initial configuration is assumed isothermal with u₀ = 0.05 per unit mass, obeying an ideal EOS with γ = 5/3. All particles are initially at rest. For this choice the internal energy is much lower than the gravitational energy (assuming G = 1) and the system collapses.

To build the initial model we have conveniently stretched a uniform grid with N = 40³ particles so that the profile given by Eq. (45) is obtained. This procedure reproduces well the 1D density profile, except near the surface layers. We have tracked the collapse of the sphere for two different values of the parameter p in both Eq. (11) and in its smoothed counterpart, Eq. (34). The results at time t = 0.8 are compared in Fig. 27 not only with an accurate 1D-PPM calculation (Steinmetz & Mueller 1993), but also with the results obtained using GADGET-2 for the same initial model.

The evolution of models calculated using Eq. (11) with p = 0.7, and with p = 1 in Eq. (34) were rather similar, so we describe the collapse of the configuration only for the latter case. The upper panels of Fig. 27 depict the density and radial velocity profiles at t = 0.8. We see only minor differences among the three calculated models. The density profiles are almost indistinguishable. Nevertheless, a close inspection of the region around the shock front (see the zoomed region in Fig. 27) reveals that the model calculated with p = 1 in Eq. (34) has a slightly better match to the exact solution than p = 0 and that obtained with GADGET-2. The profile of the radial velocity is similar in all three calculations. The evolutions calculated with GADGET-2 and with SPHYNX with p = 0 were almost identical and some minor differences may be due to the details in the implementation of the AV and gravity.

The new volume elements provide a better kernel normalization at the center and front locations (bottom-right panel in Fig. 27). Nevertheless, unlike in the point explosion tests, we have not seen a clear enhancement of | Δr | (Eq. (9)) when p ≃ 1 is used to compute the volume elements, although, overall, the profile is smoother (bottom-left panel).

7. Discussion and conclusions

In this paper, we present a new density-based SPH code, named SPHYNX, and test it in a series of traditionally problematic simulations for SPH codes in 2D and 3D. In particular, we have been able to perform a Rayleigh-Taylor simulation in a weak gravitational field g = 0.1. Additionally, the shock-blob interaction test proved that SPHYNX can efficiently suppress the tensile instability that prevents the rise of hydrodynamical instabilities and mixing in many scenarios simulated with SPH. Additionally, the outcome of other tests, such as the hydrostatic square, Kelvin-Helmholtz instability, Gresho-Chan vortex, Sedov explosion, Noh wall-shock, Evrard collapse and Triple point-shock, prove that our implementation provides results competitive with other state-of-the-art calculations. For these problems, SPHYNX produces better results than many of the extant density-based SPH codes, being qualitatively similar to those obtained with the recently developed CRKSPH scheme (Frontiere et al. 2017). But, unlike the CRKSPH method, our approach ensures the angular momentum conservation from the onset.

To achieve this, SPHYNX benefits from recent advances in the field and gathers together the latest methodologies to perform numerical simulations of astrophysical scenarios via the smoothed particle hydrodynamics technique. These methodologies include, as a novelty, a new generalized volume element estimator and a consistent update of the smoothing length and the sharpness of the interpolating kernel along with the particle density. Additionally, it counts with an integral approach to calculate gradients and a pairing-resistant family of interpolators. These features are summarized and discussed in the following.

The choice of non-standard volume elements to approximate the Euler integral equations as finite summations has a significant impact on the simulations. Following the works by Saitoh & Makino (2013) and Hopkins (2013), who generalized the VE so that they are not necessarily the trivial m/ρ choice, we postulate a new volume element which enhances the normalization of the kernel. As discussed in Sect. 3.1, the VE assigned to a particle is V_a = X_a/ ∑ _bX_bW_ab, where X_a = (m_a/ρ_a)^p is the weighting estimator of the kernel and 0 ≤ p ≤ 1 is a parameter chosen by the user. The value p = 0 reduces the VE to 1/ ∑ _bW_ab, which is the standard VE when the mass of the particles is the same. For p = 1, we have V_a = (m_a/ρ_a)/ ∑ _b(m_b/ρ_b)W_ab which is simply the re-normalized traditional volume element. As expected, a better kernel normalization (between a factor 2 and a factor 5) is obtained when these VE are used. A negative feature of the proposed VE is their tendency to overshoot the density estimation in the presence of sharp gradients when p ≃ 1. Actually, that is the fundamental reason for not taking p = 1 in the estimator X_a = (m_a/ρ_a)^p. The optimal value of p depends on the particular problem at hand, but the range 0 ≤ p ≤ 0.7, explored in this work seems to be safe. Nevertheless, a most robust implementation that allows taking p = 1 is to consider X_a = ( ⟨ m_a/ρ_a ⟩ )^p, where ⟨ . ⟩ is the SPH average of the magnitude. Although this last procedure requires the computation of the averages ⟨ m_a/ρ_a ⟩, it is the recommended default choice because of its robustness and ability to keep track of strong shocks and instabilities in the presence of sharp density gradients.

Another important feature is the dynamical choice of the interpolating kernel function. A large body of calculations carried out with SPH in the past made use of the M₄ cubic-spline function to perform interpolations. The M₄ polynomial has, however, a serious drawback: it is prone to the pairing-instability when the number of neighbors increases (e.g., exceeding n_b ≃ 60 in 3D calculations which uses the M₄ kernel). This is clearly a limitation, because in practical applications it is advisable to take as many neighbors as possible to reduce the E₀ errors in the SPH equations. A growing number of kernel candidates has been proposed during the last decade to alleviate this problem. For example, one option is to consider the natural extension of the M_n family to higher polynomial degrees, such as the quartic (M₅) or the quintic (M₆) kernels. More recently, a different family of interpolators has been proposed based on the Wendland functions, as discussed in Dehnen & Aly (2012), which shows a strong endurance against the pairing-instability. A third family of interpolators, called the sinc (harmonic-like) kernels, was introduced by Cabezón et al. (2008), which are also implemented in SPHYNX. As mentioned in Sect. 2.1, the definition of the sinc kernels is directly linked to that of the Dirac-δ function. Unlike the M_n family, which is discrete in the index n ∈ Z( +), the sinc kernels do form a continuous family, which depends on a leading exponent n ∈ ℜ( +). Actually, the M_n family could be considered as a subset of the sinc family (García-Senz et al. 2014). Using the sinc family of kernels endows the SPH technique with a flexible engine, as the shape of the kernel can be dynamically changed, in a continuous way, during run-time. This feature can be used, for example, to suppress the pairing instability (see Sect. 5.1) or to equalize the resolution behind a shock-wave (as shown in Fig. 26).

Additionally, SPHYNX estimates gradients by an integral approximation (IAD₀) which is more accurate than the traditional procedure based on the analytic derivative of the kernel function, and reduces the E0 errors caused by the particle sampling of the fluid. We fully confirm in this work the importance of this new approach, especially for handling hydrodynamic instabilities, in agreement with previous publications (García-Senz et al. 2012; Cabezón et al. 2012; Rosswog 2015a,b; Valdarnini 2016).

SPHYNX has been validated with several standard tests in two and three dimensions, ranging from strong shocks and subsonic fluid instabilities in boxes, to larger systems where the gravitational force takes over. From the analysis of these test cases we summarize the following conclusions.

The use of the Integral Approach to calculate gradients along with the traditional volume elements, V_a = m_a/ρ_a (p = 0 in Eq. (12)) and a sinc kernel with n = 5 improves the simulation of hydrodynamics instabilities subjected to small initial perturbations with respect the standard SPH. The quantitative amplitude growth-rate of the KH instability is closer to the correct growth-rate (as computed with state-of-the-art Eulerian codes) than current density-based SPH codes (with smaller L₁ errors by a factor 1.5–4), being similar to the results of the modern PSPH formulation. It is also able to reproduce the KH instability in stratified fluids with high density contrasts (ρ₂/ρ₁ ≃ 8). In the case of the RT instability, the scheme is also able to cope with small perturbations (w₀ = 0.0025) and tiny gravity values (g = −0.1), although in the latter case the non-linear evolution scarcely shows structure. In shocks, the results are similar to those provided by the standard method in identical conditions.

When the new VE are switched-on there is, in general, an increase of the quality of the simulations. We have monitored the volume normalization condition ∑ _bV_bW_ab = 1 in all calculated models and, without exception, it is better fulfilled (usually in the range of 10–20% closer to unity) with the new volume elements. This change has an impact on the overall evolution of the simulation, considerably improving the results of the simulations. A paradigmatic case is the RT instability, where the use of the VE leads to an increase of the growth-rate of the instability and to a richer evolution in the linear stage, even for the low-gravity simulation. In shock-waves, the front of the blast becomes steeper and the density peak is 10–25% higher, even in 3D. Regarding the Sedov test, the post-shock evolution of density and pressure is 5–10% closer to the analytic expectations. It is also worth noting that the VE also improve the condition | Δr | = 0 which, according to Eq. (9), is a necessary condition to exactly compute the gradient of linear functions when the IAD₀ scheme is used.

During the course of the simulations we did not see any sign of pairing instability, even when working with ≃ 50 neighbors in the 2D tests. In any case, to avoid the instability it is enough to raise the exponent of the sinc kernel above the adopted default value n = 5. We stress that, unlike other recent SPH schemes, the simulations of the KH and RT instabilities were carried out without including any artificial flux of heat or any other procedure to smooth the pressure.

Among the several improvements left for future work, we plan to improve the calculation of gravity by including a better treatment of the gravitational softening on short distances. The best way to do that is to include the gravity into the discretized SPH Lagrangian as described in Price & Monaghan (2007) and Springel (2010b). Also, the implementation and validation of switches to ensure that the AV is only added in regions where there are shocks (Cullen & Dehnen 2010), as well as noise triggers to control the velocity in subsonic flows (Rosswog 2015a) could be done with moderate effort. A more ambitious goal would be to directly calculate the volume elements solving implicitly the equation ∑ _bV_bW_ab = 1 on each particle of the system. Even though the strong coupling between particles renders any implicit calculation computationally expensive, it will probably solve the density overshooting problem seen in our explicit approach.

---

Acknowledgments

The authors want to thank S. Rosswog, M. Liebendörfer, and R. Käppeli for useful discussions and comments. We also thank M. Steinmetz for providing the PPM calculations of the Evrard test and C. McNally for giving us the reference model describing the growth-rate of the KH instability. This work has been supported by the Swiss Platform for Advanced Scientific Computing (PASC) project DIAPHANE and by the European Research Council (FP7) under ERC Advanced Grant Agreement No. 321263 – FISH. (R. Cabezón). It has also been supported by the MINECO-FEDER Spanish projects AYA2013-42762-P, AYA2014-59084-P and by the AGAUR (D. García-Senz and J. Figueira). The authors acknowledge the support of sciCORE (http://scicore.unibas.ch/) scientific computing core facility at University of Basel, where some of the calculations were performed.

References

Appendix A: Kernel separability

Spherically symmetric kernels sacrifice tensorial features in order to preserve the second-order accuracy. The notable exception is the Gaussian kernel, which is both symmetric and separable in each axis direction. But, unlike the Gaussian function, the compact-supported spherically symmetric kernels are not separable, hence neither is the sinc kernels (S_n) used in this work. Nevertheless, the S_n family approaches separability as n → ∞. This property is shown in Fig. A.1 (left panel) which depicts the profiles of (continuum lines) and that of the direct product (dotted lines) for n = 3, n = 6, and n = 12, in green, pink, and black, respectively. As we can see, the dispersion in the kernels product becomes narrower as the kernel index n increases. As a comparison, we show the same profiles for some Wendland kernels in Fig. A.1 (right panel). In this case, it is clear that Wendland kernels are not separable as the product of the unidimensional components departs from the direct multidimensional calculation. Furthermore, the result does not improve as we increase the order of the Wendland kernel. The study of the consequences that this property may have is beyond the scope of this paper, but it is undoubtedly worthy to explore in future works.

Appendix B: A better partition of unit improves the estimation of the gradient of a linear function

We consider the following integral in 1D: (B.1)Because the integrand is an odd function, I₀ is trivially equal to zero. Nevertheless, approaching I₀ with finite summations does not necessarily ensure I₀ = 0. We know that a necessary and sufficient condition to exactly reproduce the gradient of a linear function in SPH is that I₀ = ∑ _bV_bx_bW_ab = 0, (Eq. (9)), where V_b is the volume element of particle b.

Integrating I₀ by parts, (B.2)where G(x,h) is the primitive integral of W(x,h): G(x,h) = ^∫W(x,h) dx.

We now consider a very simple kernel W(x,h) whose primitive G(x,h) can be obtained analytically: (B.3)This kernel (Fulk & Quinn 1996) is spherically symmetric but currently it is not the preferred one because it is not a flat-top kernel. The normalization constant in 1D is C = 0.5. Nevertheless, it is an adequate kernel for this proof, as it admits an analytical primitive, (B.4)The first term of the RHS of Eq. (B.2) is zero, because G(x,h) vanishes very quickly when x/h increases. For spherically symmetric kernels I₀ = 0, thus the integral giving the second term should also be zero. The integral I₀ can be estimated as I₀ = I₁ + I₂ where, (B.5)(B.6)with W(x,h) given by Eq. (B.3). The analytical calculation of these integrals gives I₁ = −Ch, I₂ = + Ch, and I₀ = I₁ + I₂ = 0, which is obvious because Eq. (B.1) is an odd function. Things are different when we approach these integrals by finite summations: (B.7)On the other hand, if the V_b are well chosen, they should fulfill the equipartition condition: (B.8)Now we make the reasonable ansatz that improving V_b is actually enhancing each term inside the summation in Eq. (B.8). Thus, (B.9) and, according to Eq. (B.7), I₀ vanishes (or, in any case, it becomes very small) even when it is approached using summations. We note that this proof is strictly valid only for exponential kernels. Nevertheless, the main conclusion also holds (at least qualitatively) for other centrally-peaked, spherically-symmetric interpolating functions. The positive feedback between IAD₀ and the VE is also supported by many of the test cases presented in this work.

Appendix C: Partition of the unity: convergence rate of the estimator X = (m/ρ)^p

We discuss here a simple 2D static test which explores the effect of changing the exponent p of X_a = (m_a/ρ_a)^p into the normalization of the kernel. We considered a shock-tube filled with a two-density fluid separated by a contact discontinuity. The jump of the density across the frontier is estimated for several values of 0 ≤ p ≤ 1. The SPH particles were arranged in an ordered grid according to the density value at both sides of the discontinuity. The contact between both fluids was not smoothed. We have considered two density ratios: ρ₁/ρ₂ = 2 and ρ₁/ρ₂ = 8, and calculated the magnitude ∑ _bV_bW_ab for different values of p. The results are depicted in Fig. C.1 which shows the maximum relative error in the kernel normalization | ∑ _bV_bW_ab−1 | as a function of p, the density contrast and the particular procedure to estimate X_a, either X_a = (m_a/ρ_a)^p (continuum red lines) or X_a = ( ⟨ m_a/ρ_a ⟩ )^p (green dashed lines), where ⟨ . ⟩ is the SPH average. The smoothing length was adjusted so that there were n_b = 50 neighbors contributing to the summations.

For not too large density ratios, raising the exponent p leads to a linear improvement in kernel normalization. The convergence rate is slower when X_a = ( ⟨ m_a/ρ_a ⟩ )^p, as expected. Things are different for large density ratios where the convergence rate follows a parabolic line with a local minimum around p ≃ 0.7 (uppermost red line in Fig. C.1). Above that value it seems that any increase of p makes the convergence worse. The reason is that, in static configurations with sharp contact discontinuities, the density under/over-shoots at the two sides of the discontinuity. Taking both a large density ratio and a high p value induces a catastrophic feedback between the overshooting and X_a. Such behavior is due to the explicit nature of the implementation of the volume elements in SPH and would probably be overcome using an implicit approach. A simple recipe to circumvent the problem is to take X_a = ( ⟨ m_a/ρ_a ⟩ )^p, which leads to a slower, albeit safer, convergence rate (uppermost green line in Fig. C.1).

Appendix D: Included physics and alternative formulations of the energy equation

Most astrophysical applications need to calculate gravity, and to that extent SPHYNX incorporates an octal-tree structure (Hernquist & Katz 1989) with several levels of particle clustering. When the particle inter-distance r_ab is shorter than h_a + h_b we apply a simple smoothing to the gravitational force, (D.1)This softening usually gives satisfactory results but it is not fully compatible with the Euler-Lagrange derivation of the SPH equations (Price & Monaghan 2007).

SPHYNX also incorporates a thermal conductive transport equation compatible with the IAD formulation. That equation was described and checked in Cabezón et al. (2012) and it is reproduced here for completeness, (D.2)where d is the dimension of the space and the tilde symbol means the arithmetic average of the magnitude.

In high density plasmas with finite temperature, which characterize compact objects such as white dwarfs and neutron stars, it is often preferable to directly compute the temperature instead of the internal energy. SPHYNX can switch the energy equation, Eq. (23), to the temperature equation, (D.3)where C_v is the specific heat. This formulation of the temperature equation also leads to an almost perfect energy conservation. The internal energy can be obtained through the EOS when necessary. Equation (D.3) was recently used to calculate explosion models of Type Ia supernova with the SPHYNX code (García-Senz et al. 2016).

Also, the energy equation may be substituted without much effort by an entropy equation which ensures perfect entropy conservation in pure adiabatic fluxes (Springel & Hernquist 2002)⁸. An isentropic flow is characterized by (D.4)where A_a is a constant determined by the initial conditions of the flux. The entropy equation was incorporated as default in the GADGET-2 code (Springel 2005), bringing excellent results for ideal EOS. However, it has to be adapted to more complex and realistic EOS where the adiabatic index γ can be time-dependent and not so straightforward to know.

All Tables

All Figures