\subsection{The multipole acceptance criterion}

The main remaining question is to decide when two cells are far enough from
each others that the truncated Taylor expansion used as approximation for
the potential (eq. \ref{eq:fmm:expansion}) is accurate enough. The
criterion used to make that decision is called the \emph{multipole
  acceptance criterion} (MAC). \\
We know that (\ref{eq:fmm:expansion}) is converging towards the correct
answer provided $1>|\mathbf{r}_a + \mathbf{r}_b| / |\mathbf{R}|$. This is
hence the most basic (and always necessary) MAC that can be designed. If
this ratio is lower, the accuracy (at a fixed expansion order) is improved
and it is hence common practice to define a critical \emph{opening angle}
$\theta_{\rm cr}$ and allow the use of the multipole approximation between
two cells if

\begin{equation}
  \theta_{\rm cr} > \frac{\rho_A + \rho_B} {|\mathbf{R}|}.
  \label{eq:fmm:angle}
\end{equation}
This lets users have a second handle on the accuracy on the gravity
calculation besides the much more involved change in the expansion order
$p$ of the FMM method. Typical values for the opening angle are in the
range $[0.3, 0.7]$, with the cost of the simulation growing as $\theta_{\rm
  cr}$ decreases. \\
This method has the drawback of using a uniform criterion across the entire
simulation volume and time evolution, which means that the chosen value of
$\theta_{\rm cr}$ could be too small in some regions (leading to too many
operations for the expected accuracy) and too large in some other other
ones (leading to a lower level of accuracy than expected). \swift instead
uses a more adaptive criterion to decide when the multipole approximation
can be used. This is based on the error analysis of FMM by
\cite{Dehnen2014} and is summarised below for completeness. The key idea is
to exploit the additional information about the distribution of particles
that is encoded in the higher-order multipole terms.\\
We start by defining the scalar quantity $P_{A,n}$, the
\emph{power} of the multipole of order $n$ of the particles in cell $A$,
via
\begin{equation}
  P_{A,n}^2 = \sum_{|\mathbf{m}|=n} \frac{\mathbf{m}!}{|\mathbf{m}|!}\mathsf{M}_{A,\mathbf{m}}^2,
\end{equation}
where the sum runs over all the multipole terms of order $n$ in the
cell\footnote{Note that $P_{0} \equiv \mathsf{M}_{(0,0,0)}$ is
  just the mass of the cell and since \swift uses the centre of mass as the
  centre of expansion of the multipoles, $P_{1} = 0$.}. This
quantity is a simple upper bound for the amplitude of the multipole
($\mathsf{M}_{A, \mathbf{m}} < P_{A,|\mathbf{m}|}/|\mathbf{m}|!$)
and can hence be used to estimate the importance of the terms of a given
order in the Taylor series of the potential. Following \cite{Dehnen2014} we
then consider a sink cell $A$ and a source cell $B$ (figure \ref{fig:fmm:cells}) for which we evaluate
at order $p$ the scalar
\begin{equation}
  E_{BA,p} = \frac{1}{M_B|\mathbf{R}|^p} \sum_{n=0}^p \binom{p}{n} P_{B,n}
  \rho_A^{p-n},
  \label{eq:fmm:e_ab}
\end{equation}
with $M_B \equiv \mathsf{M}_{B,(0,0,0)}$, the sum of the mass of the
particles in cell $B$. Note that since $P_{B,n} \leq M_B
\rho_B^n$, we have $E_{BA, p} \leq \left((\rho_A +
\rho_B)/|\mathbf{R}|\right)^p$, where the right-hand side is the
expression used in the basic opening angle condition
(\ref{eq:fmm:angle}). We finally scale the $E_{BA,p}$'s by the relative
size of the two cells to define the error estimator $\tilde{E}_{BA,p}$:
\begin{equation}
  \tilde{E}_{BA,p} = 8\frac{\max(\rho_A, \rho_B)}{\rho_A + \rho_B}E_{BA,p}.
  \label{eq:fmm:e_ab_tilde}
\end{equation}
As shown by \cite{Dehnen2014}, these quantities are excellent estimators of
the error made in computing the accelerations between two cells using the
M2L and M2P kernels at a given order. We can hence use this property to
design a new MAC by demanding that the estimated acceleration error is no
larger than a certain fraction of the smallest acceleration in the sink
cell $A$. This means we can use the FMM approximation between to
approximate the accelerations in cell $A$ due to the particles in cell $B$ if
\begin{equation}
  \tilde{E}_{BA,p} \frac{M_B}{|\mathbf{R}|^2} < \epsilon_{\rm FMM} \min_{a\in
    A}\left(|\mathbf{a}_a|\right) \quad \rm{and} \quad \frac{\rho_A +
    \rho_B} {|\mathbf{R}|} < 1,
  \label{eq:fmm:mac}  
\end{equation}
where the $\mathbf{a}_a$ are the accelerations of the particles in cell $A$
and $\epsilon_{\rm FMM}$ is a tolerance parameter. Since this is self-referencing
(i.e. we need the accelerations to decide how to compute the
accelerations), we need to use a an estimator of $|\mathbf{a}_a|$. In
\swift, we follow the strategy used by \gadget and use the acceleration of
the previous time-step\footnote{On the first time-step of a simulation this
  value has not been computed yet. We hence run a fake 0th time-step with
  the simpler MAC (eq. \ref{eq:fmm:angle}), which is good enough to obtain
  approximations of the accelerations.}. The minimal norm of the
acceleration in a given cell can be computed at the same time as the P2M
and M2M kernels are evaluated in the tree construction phase. The second
condition in (\ref{eq:fmm:mac}) is necessary to ensure the convergence of the
Taylor expansion.\\
One important difference between this criterion and the purely
geometric one (\ref{eq:fmm:angle}) is that it is not symmetric in $A
\leftrightarrow B$ (i.e. $E_{AB,p} \neq E_{BA,p}$). This implies that
there are cases where a multipole in cell $A$ can be used to compute
the field tensors in cell $B$ but the multipole in $B$ cannot be used
to compute the $\mathsf{F}$ values of cell $A$ and vice versa. This
affects the tree walk by breaking the symmetry and potentially leading
to cells of different sizes interacting. \\
For the M2P kernel, the sink is a single particle $a$ and hence
$\rho_A = 0$, which simplifies some of the expressions above. In this
case, at order $p$, we get:
\begin{equation}
  E_{BA,p} = \frac{P_{B,p}}{M_B |\mathbf{R}|^p}, \qquad
  \tilde{E}_{BA,p} = 8E_{BA,p} \nonumber
\end{equation}
Note that, in this case, only the power term of the order of the
scheme appears; not a sum over the lower-order ones. This leads to the
following MAC for the M2P kernel:
\begin{equation}
  8\frac{P_{B,p}}{|\mathbf{R}|^{p+2}} < \epsilon_{\rm FMM} |\mathbf{a}_a| \quad
  \rm{and} \quad \frac{\rho_B} {|\mathbf{R}|} < 1.
    \label{eq:fmm:mac_m2p}  
\end{equation}
The value of $\epsilon_{\rm FMM}$ could in principle be different than the one
used for the M2L MAC. One special case is of particular interest to
link our expression to other results. Using the expression for order
$2$ and the approximation $P_{B,p} \approx M_B \rho_B^p$, we
get
\begin{equation}
  8\frac{M_B}{|\mathbf{R}|^2}\left(\frac{\rho_B}{|\mathbf{R}|}\right)^2
  < \epsilon_{\rm FMM} |\mathbf{a}_a| \nonumber
\end{equation}
for our MAC.  This is the same expression as the adaptive opening
angle used by \gadget \cite[see eq.18 of][]{Springel2005} up to
numerical factors and definition of the size of a multipole ($\rho$
vs. the cell edge). Note, however, that, in practice, since formally
$P_{B,p} \leq M_B \rho_B^p$, the dependence is slightly
different.\\
We conclude this section by noting that whilst the derivation of the
FMM equations and of the simple geometric MAC (eq. \ref{eq:fmm:angle})
do not make any assumptions about the functional form of $\varphi(r)$,
the more advanced MAC is valid in the specific case of the
gravitational potential $\varphi(r) = m/r$ as can be inferred from the
$m/r^2$ term appearing on the LHS of the criteria (\ref{eq:fmm:mac})
and (\ref{eq:fmm:mac_m2p}).

\subsubsection{Modifications for softened and truncated gravity}

\begin{figure}
\includegraphics[width=\columnwidth]{mac_potential.pdf}
\caption{The gravitational forces $f_{\rm SWIFT}$ computed by SWIFT
  (green line) including the force softening on the smallest scales
  and the long-range periodic mesh truncation on the largest scales
  for a simulation box of size $L$, a mesh scale-length $r_s$ and
  Plummer-equivalent softening $\epsilon_{\rm Plummer}$. The
  approximate fast estimator of the forces used in the MAC $f_{\rm
    MAC}$ is shown using yellow dash-dotted lines. Note that, by
  construction, $f_{\rm SWIFT} \leq f_{\rm MAC} \leq 1/r^2$ for all
  distances $r$.}
\label{fig:fmm:mac_potential}
\end{figure}

One drawback of using expression (\ref{eq:fmm:mac}) in the case of a
softened potential (or a potential truncated to apply long-range
forces from a mesh (Sec. \ref{ssec:mesh_summary}) is that the $M/R^2$
term will overestimate the expected contribution from the multipole to
the filed tensors, sometimes by large factors. This difference is
shown on fig. \ref{fig:fmm:mac_potential}, with for instance a ratio
of $3$ between the true forces and the Newtonian values reached a the
scale of the Plummer softening. Using the simple expression
(\ref{eq:fmm:mac}) will make the MAC too aggressive by preventing it
from using a given multipole as it will be difficult to make the large
term $M/R^2$ be below the fixed fraction $\epsilon_{\rm FMM}$ of the
total acceleration of the receiving cell. This implies more
computation as it will force the tree-walk algorithm to use more
interactions by going to the daughter cells. The estimation of the
contribution of the multipole in the MAC should hence be replaced by a
more realistic term, closer to the one actually used in the
interactions (eq. \ref{eq:fmm:force_norm}). In simulations with
periodic boundary conditions, the same reasoning applies to the
truncated force at the radii overlapping with the scale $r_s$ of the
mesh forces.

However, both the short- and long-range truncation functions are
expensive to evaluate in the context of the MAC which is called a
large number of times during a tree walk. We hence, construct a
cheaper to evaluate estimator $f_{\rm MAC}$ that is closer to the true
forces than the purely Newtonian term:
\begin{align}
f_{\rm MAC}(r) =
\left\lbrace\begin{array}{rcl}
  \left(\frac{9}{5}\right)^2 H^{-2} & \mbox{if} & r <
  \frac{5}{9}H,\\
  r^{-2} & \mbox{if} & \frac{5}{9}H \leq r < \frac{5}{3}r_s, \\
  \left(\frac{5}{3}\right)^2 r_s^2 r^{-4} & \mbox{if} & \frac{5}{3}r_s \leq r. \\
\end{array}
\right.
\label{eq:fmm:f_mac}
\end{align}
Since it is made of constants and even powers of the distance,
computin this term is much cheaper than the true forces.  This
esimator is shown as a dot-dashed line on
Fig. \ref{fig:fmm:mac_potential} and obeys the relation $f_{\rm
  SWIFT}(r) \leq f_{\rm MAC}(r) \leq 1/r^2$, with $f_{\rm SWIFT}(r)$
being the true truncated and softened norm of the gravity forces the
code solves for (green line). We use this expression in the multipole
acceptance criterion instead of the $1/|\mathbf{R}|$ term:
\begin{equation}
  \tilde{E}_{BA,p} M_Bf_{\rm MAC}(|\mathbf{R}|) < \epsilon_{\rm FMM} \min_{a\in
    A}\left(|\mathbf{a}_a|\right).
  \label{eq:fmm:mac_f_mac}  
\end{equation}
The same change is applied to the MAC used for the M2P kernel
(eq. \ref{eq:fmm:mac_m2p}). In the non-truncated un-softened case,
their expressions reduce to the \citep{Dehnen2014} one. Using this
$f_{\rm MAC}$ instead of the simpler purely-Newtonian one only makes a
difference in simulations where a lot of particles cluster below the
scale of the softening, which is often the case for hydrodynamical
simulations including radiative cooling processes. The use of this
term over the simpler $1/r^2$ estimator is a runtime parameter.