diff --git a/examples/theory/multipoles.tex b/examples/theory/multipoles.tex
index 8e938e12c6e7588877a2040681496e6249f1003b..edae2f85018c87526c254f4996d8e7f72bfc5406 100644
--- a/examples/theory/multipoles.tex
+++ b/examples/theory/multipoles.tex
@@ -6,8 +6,8 @@
 \newcommand{\rr}{\mathbf{r}}
 \newcommand{\dd}{\mathbf{d}}
 \newcommand{\vv}{\mathbf{v}}
-\newcommand{\p}[1]{\mathbf{p}_#1}
-\newcommand{\acc}{\mathbf{a}}
+\newcommand{\p}[1]{\mathbf{p}_{#1}}
+\newcommand{\acc}[1]{\mathbf{a}_{#1}}
 \newcommand{\muu}{\boldsymbol{\mu}}
 
 \title{B-H and FMM equations up to quadrupole terms}
@@ -17,7 +17,8 @@
 
 \maketitle
 
-Bold quantities are vectors. The indices $\alpha,\beta$ run over the directions $x,y,z$.
+Bold quantities are vectors, underlined quantities are matrices. The indices 
+$\alpha,\beta$ run over the directions $x,y,z$.
 
 \section{Construction of multipoles}
 
@@ -26,7 +27,8 @@ we can construct the total mass of the set
 $M_A$ and its centre of mass $\muu_A=(\mu_{A,x}, \mu_{A,y}, \mu_{A,z})$:
 
 \begin{equation}
- M_A = \sum_{i\in A} m_i, \qquad \muu_A = \frac{1}{M_A} \sum_{i \in A}
+ M_{{\rm tot},A} = \sum_{i\in A} m_i, \qquad \muu_A = \frac{1}{M_{{\rm tot},A}} 
+\sum_{i \in A}
 m_i\p{i}.
 \end{equation}
 
@@ -35,7 +37,7 @@ centre of mass $\muu_A$. The first LHS term uses Dehnen's notation.\\
 
 Monopole:
 \begin{equation}
-M_{(0,0,0)} = M_A 
+M_{(0,0,0)} = M_{{\rm tot}, A} 
 \end{equation}
 
 Dipole:
@@ -47,15 +49,18 @@ M_{(0,0,1)} &=& P_{A,z}~=0
 
 Quadrupole:
 \begin{eqnarray}
-M_{(2,0,0)} &=& I_{A,xx}~= \sum_{i\in A}m_i ( p_{i,x}-\mu_{A,x})^2\\
-M_{(0,2,0)} &=& I_{A,yy}~= \sum_{i\in A}m_i ( p_{i,y}-\mu_{A,y})^2\\
-M_{(0,0,2)} &=& I_{A,zz}~= \sum_{i\in A}m_i ( p_{i,z}-\mu_{A,z})^2\\
-M_{(1,1,0)} &=& I_{A,xy}~= \sum_{i\in A}m_i ( p_{i,x}-\mu_{A,x})( 
-p_{i,y}-\mu_{A,y})\\
-M_{(0,1,1)} &=& I_{A,yz}~= \sum_{i\in A}m_i ( p_{i,y}-\mu_{A,y})( 
-p_{i,z}-\mu_{A,z})\\
-M_{(1,0,1)} &=& I_{A,xz}~= \sum_{i\in A}m_i ( p_{i,x}-\mu_{A,x})( 
-p_{i,z}-\mu_{A,z})
+M_{(2,0,0)} &=& \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 ( 
+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 ( 
+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 ( 
+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 ( 
+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 ( 
+p_{i,x}-\mu_{A,x})( p_{i,z}-\mu_{A,z})
 \end{eqnarray}
 
 \section{Recursive construction of the quadrupoles}
@@ -68,7 +73,7 @@ the tree recursively.\\
 
 Monopole:
 \begin{equation}
- M_{A} = \sum_{B\in A} M_B
+ M_{{\rm tot},A} = \sum_{B\in A} M_{{\rm tot},B}
 \end{equation}
 
 Dipole:
@@ -79,7 +84,7 @@ Dipole:
 Quadrupole:
 \begin{equation}
  I_{A,\alpha\beta} = \sum_{B\in A}\left( I_{B,\alpha\beta} + 
-M_B\mu_{B,\alpha}\mu_{B,\beta} \right) - \frac{1}{M_A} 
+M_{{\rm tot},B}\mu_{B,\alpha}\mu_{B,\beta} \right) - \frac{1}{M_{{\rm tot},A}} 
 \mu_{A,\alpha}\mu_{A,\beta}
 \end{equation}
 
@@ -115,4 +120,52 @@ D_{(0,1,1)} &=& \partial_{yz} \phi(\rr)~= \frac{3Gr_yr_z}{|\rr|^5}\\
 D_{(1,0,1)} &=& \partial_{xz} \phi(\rr)~= \frac{3Gr_xr_z}{|\rr|^5}
 \end{eqnarray}
 
+
+\section{B-H potential and accelerations}
+
+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}(\p{i} - \muu_A)
+\end{equation}
+
+Keeping only the terms up to second order (i.e. letting the sum run over 
+all vectors with $|{\bf n}|\leq2$) and writing $\rr = \p{i} - \muu_A$, we get:
+
+\begin{eqnarray}
+ \phi(\p{i}) &=& -M_{{\rm tot},A} \frac{G}{|\rr|} -\frac{1}{2}\sum_\alpha 
+I_{A,\alpha} \left(\frac{3Gr_\alpha^2}{|\rr|^5} - 
+\frac{G}{|\rr|^3}\right) \\ 
+& &- \frac{1}{2}\sum_{\alpha,\beta} \delta_{\alpha\beta}I_{A,\alpha\beta} 
+\frac{3Gr_\alpha r_\beta}{|\rr|^5} \\
+             &=& -M_{{\rm tot},A} \frac{G}{|\rr|} +
+\frac{1}{2}\left( I_{A,xx} + I_{A,yy} +I_{A,zz}\right)\frac{G}{|\rr|^3}  \\
+ & & - \frac{1}{2}\sum_{\alpha,\beta} I_{A,\alpha\beta}  \frac{3Gr_\alpha 
+r_\beta}{|\rr|^5} \\
+&=& -M_{{\rm tot},A} \frac{G}{|\rr|} + \frac{G}{2} \frac{{\rm 
+tr}(\underline{I_A})}{|\rr|^3} - \frac{3G}{2}\frac{\rr^T \cdot 
+\underline{I_A} \cdot \rr}{|\rr|^5}
+\end{eqnarray}
+
+The accelerations $\acc{i} = -\nabla_{\rr}\phi(\p{i})$ are then given by:
+
+\begin{eqnarray}
+ a_{i,x} &=& \left[M_{{\rm tot},A} \frac{G}{|\rr|^3} - \frac{3G}{2} 
+\frac{{\rm tr}(\underline{I_A})}{|\rr|^5} - \frac{3G\underline{I_A}}{|\rr|^5} + 
+\frac{15G}{2}\frac{(\rr^T \cdot 
+\underline{I_A} \cdot \rr)}{|\rr|^7}\right] r_x \\
+ a_{i,y} &=& \left[M_{{\rm tot},A} \frac{G}{|\rr|^3} - \frac{3G}{2} 
+\frac{{\rm tr}(\underline{I_A})}{|\rr|^5} - \frac{3G\underline{I_A}}{|\rr|^5} + 
+\frac{15G}{2}\frac{(\rr^T \cdot 
+\underline{I_A} \cdot \rr)}{|\rr|^7}\right] r_y\\
+ a_{i,z} &=& \left[M_{{\rm tot},A} \frac{G}{|\rr|^3} - \frac{3G}{2} 
+\frac{{\rm tr}(\underline{I_A})}{|\rr|^5} - \frac{3G\underline{I_A}}{|\rr|^5} + 
+\frac{15G}{2}\frac{(\rr^T \cdot 
+\underline{I_A} \cdot \rr)}{|\rr|^7}\right] r_z
+\end{eqnarray}
+
+
+
+
 \end{document}