\subsection{Gravitational softening} \label{ssec:potential_softening} To avoid artificial two-body relaxation, the Dirac $\delta$-distribution corresponding to each particle is convolved with a softening kernel of a given fixed, but time-variable, scale-length $H$. Instead of the commonly used spline kernel of \cite{Monaghan1985} (e.g. in \textsc{Gadget}), we use a C2 kernel \citep{Wendland1995} which leads to an expression for the force that is cheaper to compute and has a very similar overall shape. The C2 kernel has the advantage of being branch-free leading to an expression which is faster to evaluate using vector units available on modern architectures; it also does not require any divisions to evaluate the softened forces. We set $\tilde\delta(\mathbf{r}) = \rho(|\mathbf{r}|) = W(|\mathbf{r}|, 3\epsilon_{\rm Plummer})$, with $W(r, H)$ given by \begin{align} W(r,H) &= \frac{21}{2\pi H^3} \times \nonumber \\ &\left\lbrace\begin{array}{rcl} 4u^5 - 15u^4 + 20u^3 - 10u^2 + 1 & \mbox{if} & u < 1,\\ 0 & \mbox{if} & u \geq 1, \end{array} \right. \end{align} and $u = r/H$. The potential $\varphi(r,H)$ corresponding to this density distribution reads \begin{align} \varphi(r,H) = \left\lbrace\begin{array}{rcl} f(\frac{r}{H}) \times H^{-1} & \mbox{if} & r < H,\\ r^{-1} & \mbox{if} & r \geq H, \end{array} \right. \label{eq:fmm:potential} \end{align} with $f(u) \equiv -3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3$. These choices lead to a potential at $|\mathbf{x}| = 0$ equal to the central potential of a Plummer sphere (i.e. $\varphi(0) = 1/\epsilon_{\rm Plummer}$)\footnote{Note the factor $3$ in the definition of $\rho(|\mathbf{x}|)$ which differs from the factor $2.8$ used for the cubic spline kernel as a consequence of the change of the functional form of $W$.}. From this expression the softened gravitational force can be easily obtained: \begin{align} \mathbf{\nabla}\varphi(r,H) = \mathbf{r} \cdot \left\lbrace\begin{array}{rcl} g(\frac{r}{H}) \times H^{-3} & \mbox{if} & r < H,\\ r^{-3} & \mbox{if} & r \geq H, \end{array} \right. \label{eq:fmm:force} \end{align} with $g(u) \equiv f'(u)/u = -21u^5+90u^4-140u^3+84u^2-14$. This last expression has the advantage of not containing any divisions or branching (besides the always necessary check for $r