\documentclass[a4paper,10pt]{article}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}

%opening
\title{SPH kernels in SWIFT}
\author{Matthieu Schaller}

\begin{document}

\maketitle

\section{General Definitions}

The smoothing kernels used in SPH are almost always isotropic and can hence be written in 3D as

\begin{equation}
 W(\vec{x},h) = \frac{1}{h^3}f\left(\frac{|\vec{x}|}{h}\right),
\end{equation}

where $f(q)$ is a dimensionless function, usually a low-order polynomial, normalized to unity. For computational 
reasons, this kernel 
usually has a finite support of radius $H$. In other words,

\begin{equation}
 W(\vec{x},h) = 0\quad \forall\quad |\vec{x}| > H.
\end{equation}
 One can then define the weighted number of neighbours within $H$ as

\begin{equation}
 N_{ngb} = \frac{4}{3}\pi H^3 \sum_j W(\vec{x}_i - \vec{x}_j,h).
\end{equation}

The value of $N_{ngb}$ is often used in the codes to find the smoothing length of each particle via Newton iterations 
or a bissection algorithm. $H$ is defined as \emph{the smoothing length} in the GADGET code. This definition is useful 
for implementation reasons but does not really correspond to a true physical quantity. \\
The main question is the definition of the smoothing length. The function $W(\vec{x},h)$ is invariant under the 
rescaling $h\rightarrow \alpha h,~f(q)\rightarrow\alpha^{-3}f(\alpha q)$, which makes the definition of $h$ difficult. 
This ambiguity is present in the litterature with authors using different definition of the \emph{physical} smoothing 
length, $h=\frac{1}{2}H$ or $h=H$ for instance. \\
A more physically motivated estimate is the standard deviation of the kernel:

\begin{equation}
 \sigma^2 = \frac{1}{3} \int \vec{x}^2~W(\vec{x},h)~d^3\vec{x}
\end{equation}
 
which then allows us to set $h=2\sigma$. This definition of the smoothing length is more physical as one can 
demonstrate that the reconstruction of any smooth field $A(\vec{x})$ using interpolation of particles at the point 
$\vec{x}_i$ can be expanded as

\begin{equation}
A_i  \approx A(\vec{x}_i) + \frac{1}{2}\sigma^2 \nabla^2A(\vec{x}_i) + \mathcal{O}\left(\sigma^4\right).
\end{equation}

The quantity $H/\sigma$ is independant of the choice of $h$ made and is purely a functional of $f(q)$. The number of 
neighbours (used in the code to construct the neighborhood of a given particle) can then be expressed as a function of 
this \emph{physical} $h$ (or $\sigma$). Or to relate it even more to the particle distribution, we can write 
$h=\eta\Delta x$, with $\Delta x$ the mean inter-particle separation:

\begin{equation}
 N_{ngb} = \frac{4}{3}\pi \left(\frac{1}{2}\eta\frac{H}{\sigma}\right)^3 = \frac{4}{3}\pi 
\left(\eta\zeta\right)^3
\end{equation}

This definition of the number of neighbours only depends on $f(q)$ (via $\zeta$) and on the mean inter-particle 
separation. The problem is then fully specified by specifying a form for $f(q)$ and $\eta$. \\
Experiments suggest that $\eta \approx 1.2 - 1.3$ is a reasonnable choice. The bigger $\eta$, the better the smoothing 
and hence the better the reconstruction of the field. This, however, comes at a higher computational cost as more 
interactions between neighbours will have to be computed. Also, spline kernels become instable when $\eta>1.5$. 

\section{Kernels available in SWIFT}

The different kernels available are listed below.
\paragraph{Cubic Spline Kernel}
\begin{equation*}
 f(q) = \frac{1}{\pi}\left\lbrace \begin{array}{rcl}
                      \frac{3}{4}q^3 - 15q^2 + 1 & \mbox{if} & 0 \leq q < 1 \\
		      -\frac{1}{4}q^3 + \frac{3}{2}q^2-3q+2 & \mbox{if} & 1 \leq q < 2 \\
		      0 & \mbox{if} & q \geq 2 \\
                     \end{array}\right.
\end{equation*}
with $\zeta = \frac{1}{2}\sqrt{\frac{40}{3}} \approx 1.825742$. Thus, for a resolution of $\eta = 1.235$, this kernel 
uses $N_{ngb} \approx 48$. The code uses $h = \frac{1}{2}H = \zeta \sigma$ internally.

\paragraph{Quartic Spline Kernel}
\begin{equation*}
 f(q) = \frac{1}{20\pi}\left\lbrace \begin{array}{rcl}
                      6q^4 - 15q^2 + \frac{115}{8} & \mbox{if} & 0 \leq q < \frac{1}{2} \\
		      -4q^4 + 20q^3-30q^2 + 5q + \frac{55}{4} & \mbox{if} & \frac{1}{2} \leq q < \frac{3}{2} \\
		      q^4-10q^3+\frac{75}{2}q^2-\frac{125}{2}q+\frac{625}{16} & \mbox{if} & \frac{3}{2} \leq q < 
\frac{5}{2} \\
		      0 & \mbox{if} & q \geq \frac{5}{2} \\
                     \end{array}\right.
\end{equation*}
with $\zeta = \frac{1}{2}\sqrt{\frac{375}{23}} \approx 2.018932$. Thus, for a resolution of $\eta = 1.235$, this kernel 
uses $N_{ngb} \approx 64.9$. The code uses $h = \frac{2}{5}H =\frac{4}{5}\zeta \sigma$ internally.

\paragraph{Quintic Spline Kernel}
\begin{equation*}
 f(q) = \frac{1}{120\pi}\left\lbrace \begin{array}{rcl}
                      -10q^5 + 30q^4 - 60q^2 + 66 & \mbox{if} & 0 \leq q < 1 \\
		      5q^5 - 45q^4 + 150q^3 - 210q^2 + 75q + 51 & \mbox{if} & 1 \leq q < 2 \\
		      -q^5 + 15q^4 - 90q^3 + 270q^2 - 405q + 243 & \mbox{if} & 2 \leq q < 3 \\
		      0 & \mbox{if} & q \geq 3 \\
                     \end{array}\right.
\end{equation*}
with $\zeta = \frac{1}{2}\sqrt{\frac{135}{7}} \approx 2.195775$. Thus, for a resolution of $\eta = 1.235$, this kernel 
uses $N_{ngb} \approx 83.5$. The code uses $h = \frac{1}{3}H = \frac{2}{3}\zeta \sigma$ internally.

\paragraph{Wendland $C2$ Kernel}
\begin{equation*}
 f(q) = \frac{21}{2\pi}\left\lbrace \begin{array}{rcl}
                      4 q^5-15 q^4+20 q^3-10 q^2+1 & \mbox{if} & 0 \leq q < 1 \\
		      0 & \mbox{if} & q \geq 1 \\
                     \end{array}\right.
\end{equation*}
 with $\zeta = \frac{1}{2}\sqrt{15} \approx 1.93649$. Thus, for a resolution of $\eta = 1.235$, this kernel 
uses $N_{ngb} \approx 57.3$. The code uses $h = H = 2\zeta \sigma$ internally.


\paragraph{Wendland $C4$ Kernel}
\begin{equation*}
 f(q) = \frac{495}{32\pi}\left\lbrace \begin{array}{rcl}
                      \frac{35}{3} q^8-64 q^7+ 140 q^6-\frac{448}{3} q^5+70 q^4-\frac{28}{3} q^2+1 & \mbox{if} & 0 
\leq q < 1 \\
		      0 & \mbox{if} & q \geq 1 \\
                     \end{array}\right.
\end{equation*}
 with $\zeta = \frac{1}{2}\sqrt{\frac{39}{2}} \approx 2.20794$. Thus, for a resolution of $\eta = 1.235$, this kernel 
uses $N_{ngb} \approx 84.9$. The code uses $h = H = 2\zeta \sigma$ internally.

\paragraph{Wendland $C6$ Kernel}
\begin{equation*}
 f(q) = \frac{1365}{64\pi}\left\lbrace \begin{array}{rcl}
                      32 q^{11}-231 q^{10}+704 q^9-1155 q^8+1056 q^7-462 q^6+66 q^4-11 q^2+1 & \mbox{if} & 0 
\leq q < 1 \\
		      0 & \mbox{if} & q \geq 1 \\
                     \end{array}\right.
\end{equation*}
 with $\zeta = \frac{1}{2}\sqrt{24} \approx 2.44949$. Thus, for a resolution of $\eta = 1.235$, this kernel 
uses $N_{ngb} \approx 116$. The code uses $h = H = 2\zeta \sigma$ internally.

\section{Kernel Derivatives}

The derivatives of the kernel function have relatively simple expressions:

\begin{eqnarray*}
 \vec\nabla_x W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|} \\
 \frac{\partial W(\vec{x},h)}{\partial h} &=&- \frac{1}{h^4}\left[3f\left(\frac{|\vec{x}|}{h}\right) + 
\frac{|\vec{x}|}{h}f'\left(\frac{|\vec{x}|}{h}\right)\right]
\end{eqnarray*}

Note that for all the kernels listed above, $f'(0) = 0$. 

\end{document}
