diff --git a/theory/Multipoles/bibliography.bib b/theory/Multipoles/bibliography.bib
index 193db42ea4947e49930b79cbd663562d971ec2d4..c3d1289584cab55cd8e0d4d0765d70e22f0fcf2e 100644
--- a/theory/Multipoles/bibliography.bib
+++ b/theory/Multipoles/bibliography.bib
@@ -161,4 +161,34 @@ keywords = "adaptive algorithms"
 }
 
 
+@ARTICLE{Springel2005,
+   author = {{Springel}, V.},
+    title = "{The cosmological simulation code GADGET-2}",
+  journal = {\mnras},
+   eprint = {astro-ph/0505010},
+ keywords = {methods: numerical, galaxies: interactions, dark matter},
+     year = 2005,
+    month = dec,
+   volume = 364,
+    pages = {1105-1134},
+      doi = {10.1111/j.1365-2966.2005.09655.x},
+   adsurl = {http://adsabs.harvard.edu/abs/2005MNRAS.364.1105S},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+@ARTICLE{Bagla2003,
+   author = {{Bagla}, J.~S. and {Ray}, S.},
+    title = "{Performance characteristics of TreePM codes}",
+  journal = {\na},
+   eprint = {astro-ph/0212129},
+     year = 2003,
+    month = sep,
+   volume = 8,
+    pages = {665-677},
+      doi = {10.1016/S1384-1076(03)00056-3},
+   adsurl = {http://adsabs.harvard.edu/abs/2003NewA....8..665B},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
 
diff --git a/theory/Multipoles/fmm_standalone.tex b/theory/Multipoles/fmm_standalone.tex
index d3030dc52c53eca421521023649d09522b39b7bf..65d6b522f3a6a1f9ede41091a39f4f5145cf041c 100644
--- a/theory/Multipoles/fmm_standalone.tex
+++ b/theory/Multipoles/fmm_standalone.tex
@@ -4,6 +4,7 @@
 \usepackage{times}
 \usepackage{comment}
 
+\newcommand{\gadget}{{\sc Gadget}\xspace}
 \newcommand{\swift}{{\sc Swift}\xspace}
 \newcommand{\nbody}{$N$-body\xspace}
 
@@ -32,7 +33,7 @@ Making gravity great again.
 
 \input{potential_softening}
 \input{fmm_summary}
-\input{gravity_derivatives}
+%\input{gravity_derivatives}
 \input{mesh_summary}
 
 \bibliographystyle{mnras}
diff --git a/theory/Multipoles/fmm_summary.tex b/theory/Multipoles/fmm_summary.tex
index 1ff9ab88ada6836d6118c7cfd74e39f4d1c504b3..051f01db866994b51e86d91059b8e3e06d317835 100644
--- a/theory/Multipoles/fmm_summary.tex
+++ b/theory/Multipoles/fmm_summary.tex
@@ -162,21 +162,21 @@ read:
 
 All the kernels (Eqs.~\ref{eq:fmm:P2M}-\ref{eq:fmm:L2L}) are rather
 straightforward to evaluate as they are only made of additions and
-multiplications (provided $\mathsf{D}$ can be evaluated quickly, see
-Sec.~\ref{ssec:grav_derivatives}), which are extremely efficient
-instructions on modern architectures. However, the fully expanded sums
-can lead to rather large and prone to typo expressions. To avoid any
-mishaps, we use a \texttt{python} script to generate C code in which
-all the sums are unrolled and correct by construction. In \swift, we
-implemented the kernels up to order $p=5$, as it proved to be accurate
-enough for our purpose, but this could be extended to higher order
-easily. This implies storing $56$ numbers per cell for each
-$\textsf{M}$ and $\textsf{F}$ plus three numbers for the location of
-the centre of mass. For leaf-cells with large numbers of particles, as
-in \swift, this is a small memory overhead. One further small
-improvement consists in choosing $\mathbf{z}_A$ to be the centre of
-mass of cell $A$ rather than its geometrical centre. The first order
-multipoles ($\mathsf{M}_{100},\mathsf{M}_{010},\mathsf{M}_{001}$) then
-vanish by construction. This allows us to simplify some of the
-expressions and helps reduce, albeit by a small fraction, the memory
-footprint of the tree structure.
+multiplications (provided $\mathsf{D}$ can be evaluated quickly),
+which are extremely efficient instructions on modern
+architectures. However, the fully expanded sums can lead to rather
+large and prone to typo expressions. To avoid any mishaps, we use a
+\texttt{python} script to generate C code in which all the sums are
+unrolled and correct by construction. In \swift, we implemented the
+kernels up to order $p=5$, as it proved to be accurate enough for our
+purpose, but this could be extended to higher order easily. This
+implies storing $56$ numbers per cell for each $\textsf{M}$ and
+$\textsf{F}$ plus three numbers for the location of the centre of
+mass. For leaf-cells with large numbers of particles, as in \swift,
+this is a small memory overhead. One further small improvement
+consists in choosing $\mathbf{z}_A$ to be the centre of mass of cell
+$A$ rather than its geometrical centre. The first order multipoles
+($\mathsf{M}_{100},\mathsf{M}_{010},\mathsf{M}_{001}$) then vanish by
+construction. This allows us to simplify some of the expressions and
+helps reduce, albeit by a small fraction, the memory footprint of the
+tree structure.
diff --git a/theory/Multipoles/mesh_summary.tex b/theory/Multipoles/mesh_summary.tex
index 3069257c8845804d9a307cc54fffec5e36e4ae8c..5e9a0a2fed2d95474a975aed39c17c117118970f 100644
--- a/theory/Multipoles/mesh_summary.tex
+++ b/theory/Multipoles/mesh_summary.tex
@@ -1,33 +1,82 @@
 \subsection{Coupling the FMM to a mesh for periodic long-range forces}
 \label{ssec:mesh_summary}
 
-\begin{equation}
-  S(x) = \frac{e^x}{1 + e^x}
-\end{equation}
+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 by computers. 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 $S(x) \equiv \frac{e^x}{1 + e^x}$
+as the basis of our truncation function write for the potential:
 
 \begin{align}
-  \varphi_s(r) &= \frac{1}{r}\left[2 - 2S\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] 
+  \varphi_s(r) &= \frac{1}{r} \chi(r, r_s) = \frac{1}{r}\times\left[2 - 2S\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 is
+shown on Fig.~\ref{fig:fmm:potential_short}. This choice of $S(x)$ can
+seem rather cumbersome at first but writing $\alpha(x) \equiv (1+e^x)^{-1}$,
+one can expressed all derivatives of $S(x)$ as simple polynomials in
+$\alpha(x)$, which are easy to evaluate. For instance, in the case of
+the direct force evaluation between two particles, we obtain 
+
+
 \begin{align}
-  |\mathbf{f}_s(r)| &= \frac{1}{r^2}\left[\frac{4r}{r_s}S'\left(\frac{2r}{r_s}\right) - 2S\left(\frac{2r}{r_s}\right) + 2\right] \nonumber \\
-  &= \frac{1}{r^2}\left[\frac{4r}{r_s}\frac{e^{\frac{2r}{r_s}}}{(1+e^{\frac{2r}{r_s}})^2} - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}} + 2\right]
+  |\mathbf{f}_s(r)| &=
+  \frac{1}{r^2}\times\left[\frac{4r}{r_s}S'\left(\frac{2r}{r_s}\right) -
+    2S\left(\frac{2r}{r_s}\right) + 2\right] \nonumber \\
+  %&=
+  %\frac{1}{r^2}\left[\frac{4r}{r_s}\frac{e^{\frac{2r}{r_s}}}{(1+e^{\frac{2r}{r_s}})^2}
+  %- \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}} + 2\right]
+  &=
+  \frac{1}{r^2}\times 2 \left[x\alpha(x) - x\alpha(x)^2 - e^x\alpha(x) + 1\right],
 \end{align}
+with $x\equiv2r/r_s$. The truncated force is compared to the Newtonian
+force 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.
+
+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]
+  \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{aa}
+\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{bb}
+\caption{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}
 
diff --git a/theory/Multipoles/potential_derivatives.tex b/theory/Multipoles/potential_derivatives.tex
index b000964f6ccee26891b6a0840d7ddef68411f2fb..9dc9244ca56eaa56047f561aa221dd6a761e7a5a 100644
--- a/theory/Multipoles/potential_derivatives.tex
+++ b/theory/Multipoles/potential_derivatives.tex
@@ -3,136 +3,165 @@
 
 For completeness, we give here the full expression for the first few
 derivatives of the potential that are used in our FMM scheme. We use
-the notation $\mathbf{r}=(r_x, r_y, r_z)$, $r = |\mathbf{r}|$ and
-$u=r/H$. We can construct the higher order derivatives by successively
-applying the "chain rule". We show representative examples of the
-first few relevant ones here split by order. We start by constructing
-common quantities that appear in derivatives of multiple orders.
+the notation $\mathbf{r}=(r_x, r_y, r_z)$, $r = |\mathbf{r}|$, $u=r/H$
+and $x=2r/r_s$. We also assume $H \ll r_s$. We can construct the higher order derivatives by
+successively applying the "chain rule". We show representative
+examples of the first few relevant ones here split by order. We start
+by constructing derivatives of the truncated potentials:
 
 \begin{align}
-  \mathsf{\tilde{D}}_{1}(r, u, H) =
+  \alpha(x) &= \left(1 + e^x\right)^{-1}  \nonumber \\
+  \chi(r, r_s) &= 2\left(1 - e^{2r/r_s}\alpha(2r/r_s) \right) \nonumber \\
+  \chi'(r, r_s) &= \frac{2}{r_s}\left(2\alpha(x)^2 - 2\alpha(x)\right) \nonumber \\
+  \chi''(r, r_s) &= \frac{4}{r_s^2}\left(4\alpha(x)^3 - 6\alpha(x)^2 + 2\alpha(x)\right) \nonumber \\
+  \chi^{(3)}(r, r_s) &= \frac{8}{r_s^3} \left(12\alpha(x)^4 - 24\alpha(x)^3 + 14\alpha(x)^2 -2 \alpha(x)\right) \nonumber \\
+  \chi^{(4)}(r, r_s) &= \frac{16}{r_s^4} \left(48\alpha(x)^5 - 120\alpha(x)^4 + 100\alpha(x)^3 -30 \alpha(x)^2 + 2\alpha(x)\right) \nonumber \\ 
+  \chi^{(5)}(r, r_s) &= \frac{32}{r_s^5} \left(240\alpha(x)^6 - 720\alpha(x)^5 + 780\alpha(x)^4 - 360\alpha(x)^3 + 62\alpha(x)^2 - 2\alpha(x) \right) \nonumber
+\end{align}
+We can now construct common quantities that appear in derivatives of
+multiple orders:
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{align}
+  \mathsf{\tilde{D}}_{1}(r, r_s, H) =
   \left\lbrace\begin{array}{rcl}
   \left(-3u^7 + 15u^6 - 28u^5 + 21u^4 - 7u^2 + 3\right)\times  H^{-1} & \mbox{if} & u < 1,\\
-  r^{-1} & \mbox{if} & u \geq 1, 
+  %r^{-1} & \mbox{if} & u \geq 1,
+  \chi(r, r_s) \times r^{-1} & \mbox{if} & u \geq 1,
   \end{array}
   \right.\nonumber
 \end{align}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{align}
-  \mathsf{\tilde{D}}_{3}(r, u, H) =
+  \mathsf{\tilde{D}}_{3}(r, r_s, H) =
   \left\lbrace\begin{array}{rcl}
   -\left(21u^5 - 90u^4 + 140u^3 -84u^2 +14\right)\times  H^{-3}& \mbox{if} & u < 1,\\
-  -1 \times r^{-3} & \mbox{if} & u \geq 1, 
+  %-1 \times r^{-3} & \mbox{if} & u \geq 1,
+  \left(r\chi'(r, r_s) - \chi(r, r_s)\right) \times r^{-3} & \mbox{if} & u \geq 1, 
   \end{array}
   \right.\nonumber
 \end{align}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{align}
-  \mathsf{\tilde{D}}_{5}(r, u, H) =
+  \mathsf{\tilde{D}}_{5}(r, r_s, H) =
   \left\lbrace\begin{array}{rcl}
   \left(-105u^3 + 360u^2 - 420u + 168\right)\times  H^{-5}& \mbox{if} & u < 1,\\
-  3\times r^{-5} & \mbox{if} & u \geq 1, 
+  %3\times r^{-5} & \mbox{if} & u \geq 1,
+  \left(r^2\chi''(r, r_s) - 3r\chi'(r, r_s) + 3\chi(r, r_s) \right)\times r^{-5} & \mbox{if} & u \geq 1, 
   \end{array}
   \right.\nonumber
 \end{align}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{align}
-  \mathsf{\tilde{D}}_{7}(r, u, H) =
+  \mathsf{\tilde{D}}_{7}(r, r_s, H) =
   \left\lbrace\begin{array}{rcl}
   -\left(315u - 720 + 420u^{-1}\right)\times  H^{-7} & \mbox{if} & u < 1,\\
-  -15\times r^{-7} & \mbox{if} & u \geq 1, 
+  %-15\times r^{-7} & \mbox{if} & u \geq 1,
+  \left(r^3\chi^{(3)}(r, r_s) - 6r^2\chi''(r, r_s)+15r\chi'(r, r_s)-15\chi(r, r_s)\right) \times r^{-7} & \mbox{if} & u \geq 1, 
   \end{array}
   \right.\nonumber
 \end{align}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \begin{align}
-  \mathsf{\tilde{D}}_{9}(r, u, H) =
+  \mathsf{\tilde{D}}_{9}(r, r_s, H) =
   \left\lbrace\begin{array}{rcl}
   \left(-315u^{-1} + 420u^{-3}\right)\times  H^{-9}& \mbox{if} & u < 1,\\
-  105\times r^{-9} & \mbox{if} & u \geq 1.
+  %105\times r^{-9} & \mbox{if} & u \geq 1.
+  \left(r^4\chi^{(4)}(r, r_s) - 10r^3\chi^{(3)} + 45r^2\chi''(r, r_s) - 105r\chi'(r, r_s) + 105\chi(r, r_s) \right) \times r^{-9} & \mbox{if} & u \geq 1
   \end{array}
   \right.\nonumber
 \end{align}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{align}
+  \mathsf{\tilde{D}}_{11}(r, r_s, H) =
+  \left\lbrace\begin{array}{rcl}
+  -\left(315u^{-3} - 1260u^{-5}\right)\times  H^{-11}& \mbox{if} & u < 1,\\
+  %-945\times r^{-11} & \mbox{if} & u \geq 1.
+  \left(r^5\chi^{(5)}(r, r_s) - 15r^4\chi^{(4)}(r, r_s) + 105r^3\chi^{(3)}(r, r_s) - 420r^2\chi''(r, r_s) + 945r \chi'(r, r_s) - 945\chi(r, r_s)\right) \times r^{-11} & \mbox{if} & u \geq 1.
+  \end{array}
+  \right.\nonumber
+\end{align}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Starting from the potential (Eq. \ref{eq:fmm:potential},
 reproduced here for completeness), we can now build all the relevent derivatives
 \begin{align}
-  \mathsf{D}_{000}(\mathbf{r}) = \varphi (\mathbf{r},H) =
-    \mathsf{\tilde{D}}_{1}(r, u, H) \nonumber
+  \mathsf{D}_{000}(\mathbf{r}) = \varphi (\mathbf{r}, r_s, H) =
+    \mathsf{\tilde{D}}_{1}(r, r_s, H) \nonumber
 \end{align}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \noindent\rule{6cm}{0.4pt}
 \begin{align}
-  \mathsf{D}_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \varphi (\mathbf{r},H) =
-    r_x \mathsf{\tilde{D}}_{3}(r, u, H) \nonumber
+  \mathsf{D}_{100}(\mathbf{r}) = \frac{\partial}{\partial r_x} \varphi (\mathbf{r}, r_s, H) =
+    r_x \mathsf{\tilde{D}}_{3}(r, r_s, H) \nonumber
 \end{align}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \noindent\rule{6cm}{0.4pt}
 \begin{align}
-\mathsf{D}_{200}(\mathbf{r}) = \frac{\partial^2}{\partial r_x^2} \varphi (\mathbf{r},H) = 
-r_x^2 \mathsf{\tilde{D}}_{5}(r, u, H) +
-\mathsf{\tilde{D}}_{3}(r, u, H)\nonumber
+\mathsf{D}_{200}(\mathbf{r}) = \frac{\partial^2}{\partial r_x^2} \varphi (\mathbf{r}, r_s, H) = 
+r_x^2 \mathsf{\tilde{D}}_{5}(r, r_s, H) +
+\mathsf{\tilde{D}}_{3}(r, r_s, H)\nonumber
 \end{align}
 
 \begin{align}
-\mathsf{D}_{110}(\mathbf{r}) = \frac{\partial^2}{\partial r_x\partial r_y} \varphi (\mathbf{r},H) = 
-   r_x r_y \mathsf{\tilde{D}}_{5}(r, u, H) \nonumber
+\mathsf{D}_{110}(\mathbf{r}) = \frac{\partial^2}{\partial r_x\partial r_y} \varphi (\mathbf{r}, r_s, H) = 
+   r_x r_y \mathsf{\tilde{D}}_{5}(r, r_s, H) \nonumber
 \end{align}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \noindent\rule{6cm}{0.4pt}
 \begin{align}
-\mathsf{D}_{300}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \varphi (\mathbf{r},H) = 
-  r_x^3 \mathsf{\tilde{D}}_{7}(r, u, H)
-  + 3 r_x \mathsf{\tilde{D}}_{5}(r, u, H) \nonumber
+\mathsf{D}_{300}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^3} \varphi (\mathbf{r}, r_s, H) = 
+  r_x^3 \mathsf{\tilde{D}}_{7}(r, r_s, H)
+  + 3 r_x \mathsf{\tilde{D}}_{5}(r, r_s, H) \nonumber
 \end{align}
 
 \begin{align}
-\mathsf{D}_{210}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^2 r_y} \varphi (\mathbf{r},H) = 
-r_x^2 r_y \mathsf{\tilde{D}}_{7}(r, u, H) +
-r_y \mathsf{\tilde{D}}_{5}(r, u, H) \nonumber
+\mathsf{D}_{210}(\mathbf{r}) = \frac{\partial^3}{\partial r_x^2 r_y} \varphi (\mathbf{r}, r_s, H) = 
+r_x^2 r_y \mathsf{\tilde{D}}_{7}(r, r_s, H) +
+r_y \mathsf{\tilde{D}}_{5}(r, r_s, H) \nonumber
 \end{align}
 
 \begin{align}
-\mathsf{D}_{111}(\mathbf{r}) = \frac{\partial^3}{\partial r_x\partial r_y\partial r_z} \varphi (\mathbf{r},H) = 
-  r_x r_y r_z \mathsf{\tilde{D}}_{7}(r, u, H) \nonumber
+\mathsf{D}_{111}(\mathbf{r}) = \frac{\partial^3}{\partial r_x\partial r_y\partial r_z} \varphi (\mathbf{r}, r_s, H) = 
+  r_x r_y r_z \mathsf{\tilde{D}}_{7}(r, r_s, H) \nonumber
 \end{align}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \noindent\rule{6cm}{0.4pt}
 \begin{align}
   \mathsf{D}_{400}(\mathbf{r}) &= \frac{\partial^4}{\partial r_x^4}
-  \varphi (\mathbf{r},H) =
-  r_x^4 \mathsf{\tilde{D}}_{9}(r, u, H)+
-  6r_x^2 \mathsf{\tilde{D}}_{7}(r, u, H) +
-  3 \mathsf{\tilde{D}}_{5}(r, u, H)
+  \varphi (\mathbf{r}, r_s, H) =
+  r_x^4 \mathsf{\tilde{D}}_{9}(r, r_s, H)+
+  6r_x^2 \mathsf{\tilde{D}}_{7}(r, r_s, H) +
+  3 \mathsf{\tilde{D}}_{5}(r, r_s, H)
   \nonumber
 \end{align}
 
 \begin{align}
   \mathsf{D}_{310}(\mathbf{r}) &= \frac{\partial^4}{\partial r_x^3
-    \partial r_y} \varphi (\mathbf{r},H) =
-  r_x^3 r_y \mathsf{\tilde{D}}_{9}(r, u, H) +
-  3 r_x r_y \mathsf{\tilde{D}}_{7}(r, u, H)
+    \partial r_y} \varphi (\mathbf{r}, r_s, H) =
+  r_x^3 r_y \mathsf{\tilde{D}}_{9}(r, r_s, H) +
+  3 r_x r_y \mathsf{\tilde{D}}_{7}(r, r_s, H)
   \nonumber
 \end{align}
 
 \begin{align}
   \mathsf{D}_{220}(\mathbf{r}) &= \frac{\partial^4}{\partial r_x^2
-    \partial r_y^2} \varphi (\mathbf{r},H) =
-    r_x^2 r_y^2 \mathsf{\tilde{D}}_{9}(r, u, H) +
-    r_x^2 \mathsf{\tilde{D}}_{7}(r, u, H) +
-    r_y^2 \mathsf{\tilde{D}}_{7}(r, u, H) +
-    \mathsf{\tilde{D}}_{5}(r, u, H)
+    \partial r_y^2} \varphi (\mathbf{r}, r_s, H) =
+    r_x^2 r_y^2 \mathsf{\tilde{D}}_{9}(r, r_s, H) +
+    r_x^2 \mathsf{\tilde{D}}_{7}(r, r_s, H) +
+    r_y^2 \mathsf{\tilde{D}}_{7}(r, r_s, H) +
+    \mathsf{\tilde{D}}_{5}(r, r_s, H)
   \nonumber
 \end{align}
 
 \begin{align}
   \mathsf{D}_{211}(\mathbf{r}) &= \frac{\partial^4}{\partial r_x^2
-    \partial r_y   \partial r_z} \varphi (\mathbf{r},H) =
-    r_x^2 r_y r_z \mathsf{\tilde{D}}_{9}(r, u, H) +
-    r_y r_z \mathsf{\tilde{D}}_{7}(r, u, H)
+    \partial r_y   \partial r_z} \varphi (\mathbf{r}, r_s, H) =
+    r_x^2 r_y r_z \mathsf{\tilde{D}}_{9}(r, r_s, H) +
+    r_y r_z \mathsf{\tilde{D}}_{7}(r, r_s, H)
   \nonumber
 \end{align}
 
@@ -140,47 +169,47 @@ r_y \mathsf{\tilde{D}}_{5}(r, u, H) \nonumber
 \noindent\rule{6cm}{0.4pt}
 \begin{align}
   \mathsf{D}_{500}(\mathbf{r}) &= \frac{\partial^5}{\partial r_x^5}
-  \varphi (\mathbf{r},H) =
-  r_x^5 \mathsf{\tilde{D}}_{11}(r, u, H) +
-  10r_x^3\mathsf{\tilde{D}}_{9}(r, u, H) +
-  15r_x\mathsf{\tilde{D}}_{7}(r, u, H)
+  \varphi (\mathbf{r}, r_s, H) =
+  r_x^5 \mathsf{\tilde{D}}_{11}(r, r_s, H) +
+  10r_x^3\mathsf{\tilde{D}}_{9}(r, r_s, H) +
+  15r_x\mathsf{\tilde{D}}_{7}(r, r_s, H)
   \nonumber
 \end{align}
 
 \begin{align}
   \mathsf{D}_{410}(\mathbf{r}) &= \frac{\partial^5}{\partial r_x^4
-    \partial r_y} \varphi (\mathbf{r},H) =
-  r_x^4 r_y \mathsf{\tilde{D}}_{11}(r, u, H) +
-  6 r_x^2 r_y \mathsf{\tilde{D}}_{9}(r, u, H) + 
-  3 r_y \mathsf{\tilde{D}}_{7}(r, u, H)
+    \partial r_y} \varphi (\mathbf{r}, r_s, H) =
+  r_x^4 r_y \mathsf{\tilde{D}}_{11}(r, r_s, H) +
+  6 r_x^2 r_y \mathsf{\tilde{D}}_{9}(r, r_s, H) + 
+  3 r_y \mathsf{\tilde{D}}_{7}(r, r_s, H)
   \nonumber
 \end{align}
 
 \begin{align}
   \mathsf{D}_{320}(\mathbf{r}) &= \frac{\partial^5}{\partial r_x^3
-    \partial r_y^2} \varphi (\mathbf{r},H) =
-  r_x^3 r_y^2 \mathsf{\tilde{D}}_{11}(r, u, H) +
-  r_x^3 \mathsf{\tilde{D}}_{9}(r, u, H) +
-  3 r_x r_y^2 \mathsf{\tilde{D}}_{9}(r, u, H) + 
-  3 r_x \mathsf{\tilde{D}}_{7}(r, u, H)
+    \partial r_y^2} \varphi (\mathbf{r}, r_s, H) =
+  r_x^3 r_y^2 \mathsf{\tilde{D}}_{11}(r, r_s, H) +
+  r_x^3 \mathsf{\tilde{D}}_{9}(r, r_s, H) +
+  3 r_x r_y^2 \mathsf{\tilde{D}}_{9}(r, r_s, H) + 
+  3 r_x \mathsf{\tilde{D}}_{7}(r, r_s, H)
   \nonumber
 \end{align}
 
 \begin{align}
   \mathsf{D}_{311}(\mathbf{r}) &= \frac{\partial^5}{\partial r_x^3
-    \partial r_y \partial r_z} \varphi (\mathbf{r},H) =
-  r_x^3 r_y r_z \mathsf{\tilde{D}}_{11}(r, u, H) +
-  3 r_x r_y r_z \mathsf{\tilde{D}}_{9}(r, u, H)
+    \partial r_y \partial r_z} \varphi (\mathbf{r}, r_s, H) =
+  r_x^3 r_y r_z \mathsf{\tilde{D}}_{11}(r, r_s, H) +
+  3 r_x r_y r_z \mathsf{\tilde{D}}_{9}(r, r_s, H)
   \nonumber
 \end{align}
 
 \begin{align}
   \mathsf{D}_{221}(\mathbf{r}) &= \frac{\partial^5}{\partial r_x^2
-    \partial r_y^2 \partial r_z} \varphi (\mathbf{r},H) =
-  r_x^2 r_y^2 r_z \mathsf{\tilde{D}}_{11}(r, u, H) +
-  r_x^2 r_z \mathsf{\tilde{D}}_{9}(r, u, H) +
-  r_y^2 r_z \mathsf{\tilde{D}}_{9}(r, u, H) +
-  r_z \mathsf{\tilde{D}}_{y}(r, u, H)
+    \partial r_y^2 \partial r_z} \varphi (\mathbf{r}, r_s, H) =
+  r_x^2 r_y^2 r_z \mathsf{\tilde{D}}_{11}(r, r_s, H) +
+  r_x^2 r_z \mathsf{\tilde{D}}_{9}(r, r_s, H) +
+  r_y^2 r_z \mathsf{\tilde{D}}_{9}(r, r_s, H) +
+  r_z \mathsf{\tilde{D}}_{y}(r, r_s, H)
   \nonumber
 \end{align}