The SPH Density: Kernels, Sums, and Special Tricks

A good description of the particular SPH algorithm implemented in our code can be found in Rasio & Shapiro (1992). We summarize here the salient features of the method. In what follows, the letters $i$ and $j$ refer to hydrodynamical quantities defined for each particle, and

The key advance of SPH is to associate with each particle a ``smoothing length'' $h$, representing the finite spatial extent of the particle, which can differ in value for separate particles, as well as vary in time. The smoothing length is used in calculating all hydrodynamic terms, since rather than viewing the density of the particle as a spatial delta-function, we represent the spatial density distribution of particle ``$i$'' by a kernel function, such that

$$\rho_i(\vec{r})=m_i W(\vert\vec{r}-\vec{r}_i\vert,h),$$ (1)

where $\vec{r}_i$ is the actual position of particle $i$, and $W$ is a kernel function.

For the method to work, it is necessary that this kernel function have compact support, and the kernel function be at least singly differentiable. In practice, we define the kernel function $W$ such that it drops to zero at a radius equal to 2 smoothing lengths, and we choose a function which is actually $C^2$. The choice we use is defined as follows

$$W(r,h)={1\over\pi h^3}\cases{1-{3\over2}\left({r\over h}\rig... ...\right]^3,& $1\le{r\over h}<2$,\cr 0,& ${r\over h}\ge2$.\cr}$$ (2)

To save time, we calculate this function at the beginning of a run and store it in tabulated form.

As the final results of the SPH evolution calculation should not depend on the prescription for choosing the smoothing length $h$, we are free to choose a convenient expression. In practice, we let $h$ evolve over time such that the particle overlaps a fixed number of neighbors, a parameter set prior to starting the calculation.

By using this formulation, our realization of a fluid configuration into particles results in a well-defined, continuous and differentiable density field which can be evaluated at every point in space. The density at an arbitrary point is given in terms of a sum over all neighboring particles which overlap the point, denoted here by ``$j$'', as:

$$\rho(\vec{r})=\sum_j m_j W(\vert\vec{r}-\vec{r_j}\vert,h_j).$$ (3)

The SPH method allows us to calculate the density at each particle position following the formula above, but in practice we alter the formula slightly, by taking the average smoothing kernel function computed for the particle and its neighbor, finding

$$\rho_i = \sum_j m_j \frac{W(\vert\vec{r}_i-\vec{r}_j\vert,h_i)+W(\vert\vec{r}_i-\vec{r}_j\vert,h_j)}{2}$$

$$\equiv \sum_j m_j W_{ij},$$ (4)

where we define the quantity $W_{ij}$ to be the averaged smoothing kernel function. Furthermore, we modify even this equation slightly during computations, in order to ensure momentum conservation. It turns out that the property of ``neighborhood" is not exactly reflexive, i.e., particle $j$ can be a neighbor to particle $i$, but particle $i$ is not a neighbor to particle $j$, since the two particles have different smoothing lengths. Because of this, we implement a nice trick when computing the density according to the the formula above. If particle $j$ is a neighbor to particle $i$, we calculate the first half of $W_{ij}$, i.e. the term involving $h_i$, and add it to the sum for both particle $i$ and $j$. After all, this value is just the second half of $W_{ij}$ for particle $j$. If particle $i$ is in turn a neighbor to particle $j$, we pick up the second half of the density expression involving $h_j$ later. Although this may sound complicated, we have tested it numerous times, and it does in the end produce the right answer. If you will, you can think of it as a way to enforce the balances demanded by Newton's laws on a numerical level. In the comments to the code, this method is referred to as a ``gather-scatter'' technique.