diff --git a/examples/theory/multipoles.tex b/examples/theory/multipoles.tex
index bfb2f7f97a86a066fddebe49319d40fbdfd4a94f..e89f4eae802287ba2e51c7d788945b5e47100eab 100644
--- a/examples/theory/multipoles.tex
+++ b/examples/theory/multipoles.tex
@@ -9,7 +9,10 @@
\newcommand{\p}[1]{\mathbf{p}_{#1}}
\newcommand{\acc}[1]{\mathbf{a}_{#1}}
\newcommand{\muu}{\boldsymbol{\mu}}
-
+\newcommand{\qq}[1]{\mathbf{Q}_{#1}}
+\newcommand{\jj}[1]{\underline{J_{#1}}}
+\newcommand{\ii}[1]{\underline{I_{#1}}}
+\newcommand{\identity}{\rm{Id_3}}
\title{B-H and FMM equations up to quadrupole terms}
\author{Matthieu Schaller}
@@ -18,7 +21,8 @@
\maketitle
Bold quantities are vectors, underlined quantities are matrices. The indices
-$\alpha,\beta$ run over the directions $x,y,z$.
+$\alpha,\beta$ run over the directions $x,y,z$. $\identity$ is the identity
+matrix in 3D.
\section{Construction of multipoles}
@@ -37,33 +41,33 @@ centre of mass $\muu_A$. The first LHS term uses Dehnen's notation.\\
Monopole:
\begin{equation}
-M_{(0,0,0)} = M_{{\rm tot}, A}
+M_{(0,0,0)}(\muu_A) = M_{{\rm tot}, A}
\end{equation}
Dipole:
\begin{eqnarray}
-M_{(1,0,0)} &=& P_{A,x}~= 0\\
-M_{(0,1,0)} &=& P_{A,y}~=0\\
-M_{(0,0,1)} &=& P_{A,z}~=0
+M_{(1,0,0)}(\muu_A) &=& P_{A,x}~= 0\\
+M_{(0,1,0)}(\muu_A) &=& P_{A,y}~=0\\
+M_{(0,0,1)}(\muu_A) &=& P_{A,z}~=0
\end{eqnarray}
Quadrupole:
\begin{eqnarray}
-M_{(2,0,0)} &=& \frac{1}{2}I_{A,xx}~= \frac{1}{2}\sum_{i\in A}m_i (
+M_{(2,0,0)}(\muu_A) &=& \frac{1}{2}I_{A,xx}~= \frac{1}{2}\sum_{i\in A}m_i (
p_{i,x}-\mu_{A,x})^2\\
-M_{(0,2,0)} &=& \frac{1}{2}I_{A,yy}~= \frac{1}{2}\sum_{i\in A}m_i (
+M_{(0,2,0)}(\muu_A) &=& \frac{1}{2}I_{A,yy}~= \frac{1}{2}\sum_{i\in A}m_i (
p_{i,y}-\mu_{A,y})^2\\
-M_{(0,0,2)} &=& \frac{1}{2}I_{A,zz}~= \frac{1}{2}\sum_{i\in A}m_i (
+M_{(0,0,2)}(\muu_A) &=& \frac{1}{2}I_{A,zz}~= \frac{1}{2}\sum_{i\in A}m_i (
p_{i,z}-\mu_{A,z})^2\\
-M_{(1,1,0)} &=& \frac{1}{2}I_{A,xy}~= \frac{1}{2}\sum_{i\in A}m_i (
+M_{(1,1,0)}(\muu_A) &=& \frac{1}{2}I_{A,xy}~= \frac{1}{2}\sum_{i\in A}m_i (
p_{i,x}-\mu_{A,x})( p_{i,y}-\mu_{A,y})\\
-M_{(0,1,1)} &=& \frac{1}{2}I_{A,yz}~= \frac{1}{2}\sum_{i\in A}m_i (
+M_{(0,1,1)}(\muu_A) &=& \frac{1}{2}I_{A,yz}~= \frac{1}{2}\sum_{i\in A}m_i (
p_{i,y}-\mu_{A,y})( p_{i,z}-\mu_{A,z})\\
-M_{(1,0,1)} &=& \frac{1}{2}I_{A,xz}~= \frac{1}{2}\sum_{i\in A}m_i (
+M_{(1,0,1)}(\muu_A) &=& \frac{1}{2}I_{A,xz}~= \frac{1}{2}\sum_{i\in A}m_i (
p_{i,x}-\mu_{A,x})( p_{i,z}-\mu_{A,z})
\end{eqnarray}
-\section{Recursive construction of the quadrupoles}
+\section{Recursive construction of the multipoles}
Given a set of multipoles $B$ expressed around their centre of masses
$\muu_{B}$, we can construct the total multipoles around the centre of mass
@@ -127,7 +131,7 @@ In the B-H approximation, the potential at position $\p{i}$ due to a set of
particles $A$ is given by
\begin{equation}
- \phi(\p{i}) = -\sum_{\bf n} M_{A, \bf n} D_{\bf n}(\muu_A - \p{i})
+ \phi(\p{i}) = -\sum_{\bf n} M_{\bf n}(\muu_A) D_{\bf n}(\muu_A - \p{i})
\end{equation}
Keeping only the terms up to second order (i.e. letting the sum run over
@@ -146,20 +150,101 @@ I_{A,\alpha} \left(\frac{3Gr_\alpha^2}{|\rr_{i,A}|^5} -
& & - \frac{1}{2}\sum_{\alpha,\beta} I_{A,\alpha\beta} \frac{3Gr_\alpha
r_\beta}{|\rr_{i,A}|^5} \\
&=& -G \left[ \frac{M_{{\rm tot},A}}{|\rr_{i,A}|} - \frac{1}{2} \frac{{\rm
-tr}(\underline{I_A})}{|\rr_{i,A}|^3} + \frac{3}{2}\frac{\rr_{i,A}^T \cdot
-\underline{I_A} \cdot \rr_{i,A}}{|\rr_{i,A}|^5}\right]
+tr}(\ii{A})}{|\rr_{i,A}|^3} + \frac{3}{2}\frac{\rr_{i,A}^T \cdot
+\ii{A} \cdot \rr_{i,A}}{|\rr_{i,A}|^5}\right]
\end{eqnarray}
The accelerations $\acc{i} = -\nabla_{\rr_{i,A}}\phi(\p{i})$ are then given by:
\begin{eqnarray}
\acc{i} = G\left[\frac{M_{{\rm tot},A} \rr_{i,A}}{|\rr_{i,A}|}
--\frac{3}{2}\frac{{\rm tr}(\underline{I_A})\rr_{i,A}}{|\rr_{i,A}|^5}
--\frac{3\underline{I_A}\cdot\rr_{i,A}}{|\rr_{i,A}|^5}
-+\frac{15}{2}\frac{(\rr_{i,A}^T \cdot \underline{I_A} \cdot
+-\frac{3}{2}\frac{{\rm tr}(\ii{A})\rr_{i,A}}{|\rr_{i,A}|^5}
+-\frac{3\ii{A}\cdot\rr_{i,A}}{|\rr_{i,A}|^5}
++\frac{15}{2}\frac{(\rr_{i,A}^T \cdot \ii{A} \cdot
\rr_{i,A})\rr_{i,A}}{|\rr_{i,A}|^7}\right]
\end{eqnarray}
+The last two expressions are used in the Quickshed example code as well as in
+Bonsai. Gadget only uses the first term in each expression, avoiding the
+construction of the matrices $\ii{A}$. In practice, the 2-nd order
+accurate B-H method requires the storage of 10 variables per cell ($M_{{\rm
+tot},A}, \muu_A, \ii{A}$).
+
+\section{FMM Field tensors}
+
+Instead of computing the potential and accelerations of each particle, field
+tensors (generated by the set of particles A) at position $\muu_B$ can
+be derived. The accelerations of each particle can then be obtained from
+these field tensors. These tensors are given by:
+
+\begin{equation}
+ F_{\bf n}(\muu_B) = \sum_{|{\bf m}| + |{\bf n}|\leq2} M_{\bf m}(\muu_A)D_{{\bf
+n}+{\bf m}}(\muu_B - \muu_A)
+\end{equation}
+
+Writing these out explicitly with $\rr_{BA} = \muu_B - \muu_A$, we obtain the
+following set of expressions.\\
+
+Monopole:
+\begin{eqnarray}
+ F_{(0,0,0)}(\muu_B) ~=~ N_B &=& \sum_{|{\bf m}| \leq 2} M_{\bf
+m}(\muu_A)D_{{\bf m}}(\rr_{BA}) \\
+&=& M_{(0,0,0)}(\muu_A)D_{(0,0,0)}(\rr_{BA})\\
+& & +M_{(2,0,0)}(\muu_A)D_{(2,0,0)}(\rr_{BA}) \\
+& & +M_{(1,1,0)}(\muu_A)D_{(1,1,0)}(\rr_{BA})\\
+& & + \dots \\
+&=&M_{{\rm tot}, A}\phi(\rr_{BA}) +
+\frac{1}{2}\sum_{\alpha,\beta}I_{A,\alpha\beta}\partial_{\alpha\beta}\phi(\rr_{
+BA }) \\
+&=& \frac{GM_{{\rm tot}, A}}{|\rr_{BA}|} -\frac{1}{2}\frac{{\rm
+tr}(\ii{A})}{|\rr_{BA}|^3} +
+\frac{3}{2}\frac{\rr_{BA}^T \cdot
+\ii{A} \cdot \rr_{BA}}{|\rr_{BA}|^5}
+\end{eqnarray}
+
+Dipole:
+\begin{eqnarray}
+ F_{(1,0,0)}(\muu_B) ~=~ Q_{B,x} &=& \sum_{|{\bf m}| \leq 1} M_{\bf
+m}(\muu_A)D_{{\bf
+m}+(1,0,0)}(\rr_{BA}) \\
+&=& M_{(0,0,0)}(\muu_A)D_{(1,0,0)}(\rr_{BA})\\
+&=& M_{{\rm tot}, A} \partial_x \phi(\rr_{BA}) \\
+&=& M_{{\rm tot}, A} \frac{G}{|\rr_{BA}|^3}r_{BA,x}
+\end{eqnarray}
+
+Quadrupole (diagonal term):
+\begin{eqnarray}
+ F_{(2,0,0)}(\muu_B)~=~ J_{B,xx} &=& \sum_{|{\bf m}| \leq 0} M_{\bf
+m}(\muu_A)D_{{\bf
+m}+(2,0,0)}(\rr_{BA}) \\
+&=& M_{(0,0,0)}(\muu_A)D_{(2,0,0)}(\rr_{BA}) \\
+&=& M_{{\rm tot}, A} \partial_{xx} \phi(\rr_{BA}) \\
+&=& M_{{\rm tot}, A} \left(\frac{3Gr_{BA,x}^2}{|\rr_{BA}|^5}
+- \frac{G}{|\rr_{BA}|^3}\right)
+\end{eqnarray}
+
+
+Quadrupole (off-diagonal term):
+\begin{eqnarray}
+ F_{(1,1,0)}(\muu_B)~=~J_{B,xy} &=& \sum_{|{\bf m}| \leq 0} M_{\bf
+m}(\muu_A)D_{{\bf
+m}+(1,1,0)}(\rr_{BA}) \\
+&=& M_{(0,0,0)}(\muu_A)D_{(1,1,0)}(\rr_{BA}) \\
+&=& M_{{\rm tot}, A} \partial_{xy} \phi(\rr_{BA}) \\
+&=& M_{{\rm tot}, A}\frac{3Gr_{BA,x}r_{BA,y}}{|\rr|^5}
+\end{eqnarray}
+
+All these terms can be written using a more compact notation.
+\begin{eqnarray}
+N_B &=& \frac{GM_{{\rm tot}, A}}{|\rr_{BA}|} -\frac{1}{2}\frac{{\rm
+tr}(\ii{A})}{|\rr_{BA}|^3} +
+\frac{3}{2}\frac{\rr_{BA}^T \cdot
+\ii{A} \cdot \rr_{BA}}{|\rr_{BA}|^5} \\
+\qq{B} &=& \frac{GM_{{\rm tot}, A}}{|\rr_{BA}|^3} \rr_{BA}\\
+\jj{B} &=& \frac{3GM_{{\rm tot}, A}}{|\rr_{BA}|^5} \rr_{BA}
+\otimes\rr_{BA} - \frac{GM_{{\rm tot}, A}}{|\rr_{BA}|^3} \identity
+\end{eqnarray}
+