\subsection{Coupling the FMM to a mesh for periodic long-range forces}
\label{ssec:mesh_summary}

We truncate the potential and forces computed via the FMM using a
smooth function that drops quickly to zero at some scale $r_s$ set by
the top-level mesh. Traditionally, implementations have used
expressions which are cheap to evaluate in Fourier space
\citep[e.g.][]{Bagla2003, Springel2005}. This, however, implies a
large cost for each interaction computed within the tree as the
real-space truncation function won't have a simple analytic form that
can be evaluated efficiently even by modern architectures (typically
an $\mathrm{erf}()$ function). Since the FMM scheme involves to not
only evaluate the forces but higher-order derivatives, a more
appropiate choice is necessary. We use the sigmoid
$\sigma(w) \equiv \frac{e^w}{1 + e^w}$ as the basis of our truncation
function and write for the potential:

\begin{align}
  \varphi_s(r) &= \frac{1}{r} \times \chi(r, r_s) = \frac{1}{r}\times\left[2 - 2\sigma\left(\frac{2r}{r_s}\right)\right].% \nonumber\\
  %&= \frac{1}{r}\left[2 - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}}.\right] 
\end{align}
This function alongside the trunctation function used in
\gadget\footnote{For completeness, the \gadget expression reads:\\
 $\varphi_s(r) = \frac{1}{r} \times \mathrm{erfc}(\frac{1}{2}\frac{r}{r_s})$.} is shown
on Fig.~\ref{fig:fmm:potential_short}. This choice of $\sigma(w)$ can
seem rather cumbersome at first but writing
$\alpha(w) \equiv (1+e^w)^{-1}$, one can express all derivatives of
$\sigma(w)$ as simple polynomials in $\alpha(w)$ (with an identical
$e^w$ pre-factor), which are easy and cheap to evaluate (see Appendix
\ref{sec:pot_derivatives}). For instance, in the case of the direct
force evaluation between two particles, we obtain
\begin{align}
  |\mathbf{f}_s(r)| &= \left|\frac{\partial}{\partial r}\varphi_s(r)\right|
                      = \left|\frac{\partial}{\partial r}\left(\frac{1}{r} \chi(r, r_s)\right) \right|\nonumber \\
  &=  \frac{1}{r^2}\times\left[-\frac{4r}{r_s}\sigma'\left(\frac{2r}{r_s}\right) -
    2\sigma\left(\frac{2r}{r_s}\right) + 2\right] \nonumber \\
  &=
  % \frac{1}{r^2}\times 2 \left[x\alpha(x) - x\alpha(x)^2 - e^x\alpha(x) + 1\right],
  % \frac{1}{r^2}\times 2 \left[1 - e^x\alpha(x) - xe^x\alpha^2(x)\right],
    \frac{1}{r^2}\times 2 \left[1 - e^x\left(\alpha(x) - x\alpha(x)^2\right) \right]
\end{align}
with $x\equiv2r/r_s$. The truncated force is compared to the Newtonian
force and to the \gadget truncated forces\footnote{For completeness,
  the \gadget expression for the norm of the truncated forces is:
  $|\mathbf{f}_s(r)| = \frac{1}{r^2} \times
  \left[\mathrm{erfc}\left(\frac{1}{2}\frac{r}{r_s}\right) +
    \frac{1}{\sqrt{\upi}}\frac{r}{r_s}\exp\left(-\frac{1}{4}\frac{r^2}{r_s^2}\right)\right]$.}
on Fig.~\ref{fig:fmm:force_short}. At distance $r<r_s/10$, the
truncation term is negligibly close to one and the truncated forces
can be replaced by their Newtonian equivalent. We use this
optimization in \swift and only compute truncated forces between pairs
of particles that are in tree-leaves larger than $1/10$ of the mesh
size or between
two tree-leaves distant by more than that amount.\\

\textcolor{red}{MORE WORDS HERE.}\\

The truncation function in Fourier space reads

\begin{equation}
  \tilde\varphi_l(k) =
  \frac{1}{k^2}\left[\frac{\upi}{2}kr_s\textrm{csch}\left(\frac{\upi}{2}kr_s\right)
    \right]
\end{equation}

\begin{figure}
\includegraphics[width=\columnwidth]{potential_short.pdf}
\caption{Truncated potential used in \swift (green line) and \gadget
  (yellow line) alongside the full Newtonian potential (blue dasheed
  line). The green dash-dotted line corresponds to the same
  trunctation function where the exponential in the sigmoid is
  replaced by a sixth order Taylor expansion. At $r>4r_s$, the
  truncated potential becomes negligible.}
\label{fig:fmm:potential_short}
\end{figure}



\begin{figure}
\includegraphics[width=\columnwidth]{force_short.pdf}
\caption{Norm of the truncated forces used in \swift (green line) and
  \gadget (yellow line) alongside the full Newtonian force term (blue
  dasheed line). The green dash-dotted line corresponds to the same
  trunctation function where the exponential in the sigmoid is
  replaced by a sixth order Taylor expansion. At $r<r_s/10$, the
  truncated forces becomes almost equal to the Newtonian ones and can
  safely be replaced by their simpler form. The deviation between the
  exact expression and the one obtained from Taylor expansion at
  $r>r_s$ has a small impact since no pairs of particles should
  interact directly over distances of order the mesh size. }
\label{fig:fmm:force_short}
\end{figure}


\begin{figure}
\includegraphics[width=\columnwidth]{potential_long.pdf}
\caption{cc}
\label{fig:fmm:potential_long}
\end{figure}


\subsubsection{Normalisation of the potential}

The gravitational potential computed by the combination of the mesh
and multipole method needs to be corrected to obtain the zero point of
the peculiar potential. Learning from the Ewald summation technique
(see Sec.~\ref{ssec:exact_forces}), we notice that of the three terms
entering the calculation of the potential with periodic boundary
conditions (see Eq. 2.11 of \cite{Hernquist1991}), the last two are
computed by the mesh and tree respectively. The only reamining term is
the constant
\begin{equation}
  \frac{4\pi r_s^2}{L^3}\sum_a m_a,
\end{equation}
where the sum runs over all particles and $L$ is the comoving
side-length of the simulation volume. This constant is added to the
potential of all the particles in the simulation.