diff --git a/theory/Multipoles/bibliography.bib b/theory/Multipoles/bibliography.bib
new file mode 100644
index 0000000000000000000000000000000000000000..abdb5855f52aefd9b4562851946d16bd402db68a
--- /dev/null
+++ b/theory/Multipoles/bibliography.bib
@@ -0,0 +1,27 @@
+@ARTICLE{Dehnen2014,
+   author = {{Dehnen}, W.},
+    title = "{A fast multipole method for stellar dynamics}",
+  journal = {Computational Astrophysics and Cosmology},
+archivePrefix = "arXiv",
+   eprint = {1405.2255},
+ primaryClass = "astro-ph.IM",
+     year = 2014,
+    month = sep,
+   volume = 1,
+      eid = {1},
+    pages = {1},
+      doi = {10.1186/s40668-014-0001-7},
+   adsurl = {http://adsabs.harvard.edu/abs/2014ComAC...1....1D},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@article{Hardy2006,
+  title={Combinatorics of partial derivatives},
+  author={Hardy, Michael},
+  journal={Electron. J. Combin},
+  volume={13},
+   number={1},
+   pages={13},
+   year={2006}
+ }
+	      
diff --git a/theory/Multipoles/fmm_standalone.tex b/theory/Multipoles/fmm_standalone.tex
new file mode 100644
index 0000000000000000000000000000000000000000..6948708622dbf0cfc05daa9d63917f7a3017bfc2
--- /dev/null
+++ b/theory/Multipoles/fmm_standalone.tex
@@ -0,0 +1,22 @@
+\documentclass[fleqn, usenatbib, useAMS,a4paper]{mnras}
+\usepackage{graphicx}
+\usepackage{amsmath,paralist,xcolor,xspace,amssymb}
+\usepackage{times}
+
+\newcommand{\swift}{{\sc Swift}\xspace}
+
+%opening
+\title{FMM in SWIFT}
+\author{Matthieu Schaller}
+
+\begin{document}
+
+\maketitle
+
+\input{gravity_derivatives}
+
+\bibliographystyle{mnras}
+\bibliography{./bibliography.bib}
+
+
+\end{document}
diff --git a/theory/Multipoles/gravity_derivatives.tex b/theory/Multipoles/gravity_derivatives.tex
new file mode 100644
index 0000000000000000000000000000000000000000..913bdf752c666aac5913c6e6452bdb06f3fcdc86
--- /dev/null
+++ b/theory/Multipoles/gravity_derivatives.tex
@@ -0,0 +1,54 @@
+We use the multi-index notation of \cite{Dehnen2014} to simplify expressions.
+
+\subsection{Derivatives of the gravitational potential}
+
+The calculation of all the
+$D_\mathbf{n}(x,y,z) \equiv \nabla^{\mathbf{n}}\phi(x,y,z)$ terms up
+to the relevent order can be quite tedious and it is beneficial to
+automatize the whole setup. Ideally, one would like to have an
+expression for each of this term that is only made of multiplications
+and additions of each of the coordinates and the inverse distance. We
+achieve this by writing $\phi$ as a composition of functions
+$\phi(u(x,y,z))$ and apply the \textit{Fa\`a di Bruno}
+formula \citep[i.e. the ``chain rule'' for higher order derivatives,
+e.g.][]{Hardy2006} to construct our terms:
+
+\begin{equation}
+\label{eq:faa_di_bruno}
+\frac{\partial^n}{\partial x_1 \cdots \partial x_n} \phi(u)
+= \sum_{A} \phi^{(|A|)}(u) \prod_{B \in
+A} \frac{\partial^{|B|}}{\prod_{c\in B}\partial x_c} u(x,y,z),
+\end{equation}
+where $A$ is the set of all partitions of $\lbrace1,\cdots, n\rbrace$,
+$B$ is a block of a partition $A$ and $|\cdot|$ denotes the
+cardinality of a set. For generic functions $\phi$ and $u$ this
+formula yields an untracktable number of terms; an 8th-order
+derivative will have $4140$ (!)  terms in the sum\footnote{The number
+of terms in the sum is given by the Bell number of the same order}. \\
+We choose to write
+\begin{align}
+   \phi(x,y,z) &= 1 / \sqrt{u(x,y,z)}, \\
+   u(x,y,z) &= x^2 + y^2 + z^2.
+\end{align}
+This choice allows to have derivatives of any order of $\phi(u)$ that
+only depend on powers of $u$:
+
+\begin{equation}
+f^{(n)}(u) = \frac{\Gamma(\frac{1}{2})}{\Gamma(\frac{1}{2} -
+n)}\frac{1}{u^{n+\frac{1}{2}}}.
+\end{equation}
+More importantly, this choice of decomposition allows us to have
+derivatives of $u$ only up to second order in $x$, $y$ or $z$. The
+number of non-zero terms in eq. \ref{eq:faa_di_bruno} is hence
+drastically reduced. For instance, when computing
+$D_{(4,1,3)} \equiv \frac{\partial^8}{\partial x^4 \partial y \partial
+z^3} \phi$, $4100$ of the $4140$ terms will involve at least one
+zero-valued derivative (e.g. $\partial^3/\partial x^3$ or
+$\partial^2/\partial x\partial y$) of $u$. Furthermore, among the 40
+remaining terms, many will involve the same derivatives and can be
+grouped together, leaving us with a sum of six products of $x$,$y$ and
+$z$. This is generally the case for most of the $D_\mathbf{n}$'s and
+figuring out which terms are identical in a given set of partitions of
+$\lbrace1,\cdots, n\rbrace$ is an interesting exercise in
+combinatorics left for the reader \citep[see also][]{Hardy2006}.
+
diff --git a/theory/Multipoles/run.sh b/theory/Multipoles/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..34304a37dc1cb1fad3a7442544a8df793ea94d35
--- /dev/null
+++ b/theory/Multipoles/run.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+python kernels.py
+pdflatex -jobname=fmm fmm_standalone.tex
+bibtex fmm.aux
+pdflatex -jobname=fmm fmm_standalone.tex
+pdflatex -jobname=fmm fmm_standalone.tex