From f55069e8e38bd77fe2155a44e794e45f774d3927 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 21 Apr 2015 18:15:30 +0100
Subject: [PATCH] Explicit version of the FMM equations. First part:
 multipoles.

---
 examples/theory/multipoles.tex     | 240 ++++------------------------
 examples/theory/multipoles_old.tex | 245 +++++++++++++++++++++++++++++
 2 files changed, 274 insertions(+), 211 deletions(-)
 create mode 100644 examples/theory/multipoles_old.tex

diff --git a/examples/theory/multipoles.tex b/examples/theory/multipoles.tex
index 0ae9db4..54cafb8 100644
--- a/examples/theory/multipoles.tex
+++ b/examples/theory/multipoles.tex
@@ -10,7 +10,8 @@
 \newcommand{\acc}{\mathbf{a}}
 \newcommand{\muu}{\boldsymbol{\mu}}
 
-\title{Multipole expansion}
+\title{FMM and B-H equations up to quadrupole terms}
+\author{Matthieu Schaller}
 
 \begin{document}
 
@@ -18,228 +19,45 @@
 
 Bold quantities are vectors. The indices $\alpha,\beta$ run over the directions $x,y,z$.
 
-For a set of particles at position $\p{i}=(x_i, y_i, z_i)$ with mass $m_i$, we can construct the total mass of the set 
-$M$ and its centre of mass $\muu=(\mu_x, \mu_y, \mu_z)$:
+\section{Construction of multi-poles}
 
-\begin{equation}
- M = \sum_{i=1}^N m_i, \qquad \muu = \frac{1}{M} \sum_{i=1}^N m_i\p{i}.
-\end{equation}
-
-For each particle, we can define $r_i= \sqrt{(x_i-\mu_x)^2 + (y_i-\mu_y)^2 + (z_i-\mu_z)^2}$, the distance to the 
-centre of mass.\\
-The quadrupole moment $\overline{\overline{Q}}$ of this mass distribution is then written as
-
-\begin{eqnarray}
- Q_{xx} &=& \sum_i m_i \left( 3(x_i-\mu_x)^2 - r_i^2\right) \\
- Q_{xy} &=& \sum_i m_i 3(x_i-\mu_x)(y_i-\mu_y) \\
- Q_{xz} &=& \sum_i m_i 3(x_i-\mu_x)(z_i-\mu_z) \\
- ~& & \nonumber\\
- Q_{yx} &=& \sum_i m_i 3(y_i-\mu_y)(x_i-\mu_x) \\
- Q_{yy} &=& \sum_i m_i \left( 3(y_i-\mu_y)^2 - r_i^2\right) \\
- Q_{yz} &=& \sum_i m_i 3(y_i-\mu_y)(z_i-\mu_z) \\
- ~& &\nonumber\\
- Q_{zx} &=& \sum_i m_i 3(z_i-\mu_z)(x_i-\mu_x)  \\
- Q_{zy} &=& \sum_i m_i 3(z_i-\mu_z)(y_i-\mu_y) \\
- Q_{zz} &=& \sum_i m_i \left( 3(z_i-\mu_z)^2 - r_i^2\right)
-\end{eqnarray}
-Note that this matrix is symmetric and traceless ($Q_{xx}+Q_{yy}+Q_{zz}=0$), so only $5$ components need to be 
-computed. \\
-
-\pagebreak
-
-The potential created by this distribution of mass on a point at position $\rr=(r_x,r_y,r_z)$ from the centre of mass 
-$\muu$ is then:
+For a set of particles $A$ at position $\p{i}=(x_i, y_i, z_i)$ with mass $m_i$, 
+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}
- \phi(\rr) = -\frac{G}{|\rr|}M + \frac{1}{2}\frac{G}{|\rr|^5}\sum_{\alpha,\beta} r_\alpha r_\beta Q_{\alpha\beta}
+ M_A = \sum_{i\in A} m_i, \qquad \muu_A = \frac{1}{M_A} \sum_{i \in A}
+m_i\p{i}.
 \end{equation}
 
-Taking the gradient of this expression to get the acceleration, we obtain:
+The multi-poles can then be computed around the centre of mass $\muu_A$. We use 
+Dehnen's notation.
 
+Monopole:
 \begin{equation}
- -\nabla_\gamma\phi(\rr)= -\frac{GM}{|\rr|^3} r_\gamma + 
-\frac{1}{2}G\sum_{\alpha,\beta}\left(\frac{\delta_{\alpha\gamma}r_\beta}{|\rr|^5} + 
-\frac{\delta_{\gamma\beta}r_\alpha}{|\rr|^5} -5\frac{r_\alpha r_\beta r_\gamma}{|\rr|^7}\right)Q_{\alpha\beta}
+M_{(0,0,0)} = M_A 
 \end{equation}
 
-Writing this explicitly for each coordinate (using the symmetry of $\overline{\overline{Q}}$), we get:
-
+Dipole:
 \begin{eqnarray}
- a_x &=& -\frac{Gr_x}{|\rr|^3}M + \frac{1}{2} \frac{G}{|\rr|^5}\left[2r_x Q_{xx} + 2r_yQ_{xy} + 2r_z Q_{xz}\right] - 
-\frac{5}{2} \frac{G r_x}{|\rr|^7} \Xi, \\
- a_y &=& -\frac{Gr_y}{|\rr|^3}M + \frac{1}{2} \frac{G}{|\rr|^5}\left[2r_x Q_{yx} + 2r_yQ_{yy} + 2r_z Q_{yz}\right] - 
-\frac{5}{2} \frac{G r_y}{|\rr|^7} \Xi, \\
- a_z &=& -\frac{Gr_z}{|\rr|^3}M + \frac{1}{2} \frac{G}{|\rr|^5}\left[2r_x Q_{zx} + 2r_yQ_{zy} + 2r_z Q_{zz}\right] - 
-\frac{5}{2} \frac{G r_z}{|\rr|^7} \Xi, \\
+M_{(1,0,0)} &=& 0\\
+M_{(0,1,0)} &=& 0\\
+M_{(0,0,1)} &=& 0
 \end{eqnarray}
-with the common coefficient $\Xi$ defined as:
-
-\begin{equation*}
- \Xi = \rr\cdot(\overline{\overline{Q}}\cdot \rr) = r_x^2Q_{xx} + r_y^2Q_{yy} + r_z^2Q_{zz} + 2r_xr_yQ_{xy} + 
-2r_xr_zQ_{xz} + 2r_yr_zQ_{yz}
-\end{equation*}
-
 
+Quadrupole:
+\begin{eqnarray}
+M_{(2,0,0)} &=& I_{xx}~= \sum_{i\in A}m_i ( p_{i,x}-\mu_{A,x})^2\\
+M_{(0,2,0)} &=& I_{yy}~= \sum_{i\in A}m_i ( p_{i,y}-\mu_{A,y})^2\\
+M_{(0,0,2)} &=& I_{zz}~= \sum_{i\in A}m_i ( p_{i,z}-\mu_{A,z})^2\\
+M_{(1,1,0)} &=& I_{xy}~= \sum_{i\in A}m_i ( p_{i,x}-\mu_{A,x})( 
+p_{i,y}-\mu_{A,y})\\
+M_{(0,1,1)} &=& I_{yz}~= \sum_{i\in A}m_i ( p_{i,y}-\mu_{A,y})( 
+p_{i,z}-\mu_{A,z})\\
+M_{(1,0,1)} &=& I_{xz}~= \sum_{i\in A}m_i ( p_{i,x}-\mu_{A,x})( 
+p_{i,z}-\mu_{A,z}
+\end{eqnarray}
 
-% 
-% Let's consider the gravitational acceleration $\acc=(a_x,a_y,a_z)$ that a set of point masses at position 
-% $\p{i}=(p_{i,x}, p_{i,y}, p_{i,z})$ with masses $m_i$ generate at position $\rr=(r_x, r_y, r_z)$:
-% 
-% \begin{equation}
-%  \acc(\rr) = \sum_i \frac{Gm_i}{|\rr - \p{i}|^3}(\rr - \p{i})
-% \end{equation}
-% 
-% This expression can split in one expression for each of the three spatial coordinates $u=x,y,z$:
-% 
-% \begin{equation}
-%  a_u (\rr) = \sum_i \frac{m_i G}{|\rr-\p{i}|^3} ( r_u - p_{i,u}) = \sum_i m_i f_u(\rr - \p{i}),
-% \end{equation}
-% 
-% with 
-% 
-% \begin{equation}
-%  f_u (\vv) = \frac{G}{|\vv|^3} v_u.
-% \end{equation}
-% 
-% We define two quantities to simplify the notations: the total mass $M$ of the set of point masses and the centre of 
-% mass $\muu=(\mu_x, \mu_y, \mu_z)$ of this set:
-% 
-% \begin{equation}
-%  M = \sum_{i=1}^N m_i, \qquad \muu = \frac{1}{M} \sum_{i=1}^N m_i\p{i}.
-% \end{equation}
-% 
-% 
-% We can now expand the functions $f_u$ around the vector linking the particle to the centre of mass $\rr-\muu$:
-% 
-% \begin{eqnarray}
-%  a_u(\rr) &\approx& \sum_i m_i f_u(\rr - \muu) \\
-%   & & + \sum_i m_i (\p{i}-\muu) \cdot \nabla f_u(\rr - \muu)\\
-%   & & + \frac{1}{2}\sum_i m_i (\p{i}-\muu)\cdot \nabla^2 f_u(\rr-\muu)\cdot (\p{i} - \muu).
-% \end{eqnarray}
-% 
-% The first order term is identically zero and can hence be dropped. Re-arranging some of the terms, introducing the 
-% vector $\dd=\rr-\muu=(d_x,d_y,d_z)$ and using the fact that the Hessian matrix of $f_u$ is symmetric, we get:
-% \begin{eqnarray}
-%  a_u(\rr) &=& Mf_u(\dd) \\
-%            & & + \frac{1}{2} \sum_i m_i \p{i}\cdot \nabla^2f_u(\dd)\cdot \p{i}  \\
-%            & & + \frac{1}{2} M \muu\cdot \nabla^2f_u(\dd)\cdot \muu \\
-%            & & - \sum_i m_i \p{i}\cdot \nabla^2f_u(\dd)\cdot \muu \\
-%            &=& Mf_u(\dd) \\
-%            & & + \frac{1}{2} \sum_i m_i \p{i}\cdot \nabla^2f_u(\dd)\cdot \p{i} \\
-%            & &- \frac{1}{2} M \muu\cdot \nabla^2f_u(\dd)\cdot \muu
-% \end{eqnarray}
-% 
-% The gradient of $f_u$ reads
-% \begin{equation}
-%  \nabla f_u(\dd) = \frac{-3Gd_u}{|\dd|^5}\dd + \frac{G}{|\dd|^3}\hat{\mathbf{e}}_u,
-% \end{equation}
-% 
-% with $\hat{\mathbf{e}}_u$ a unit vector along the $u$-axis. The different components of the Hessian matrix then read:
-% 
-% \begin{eqnarray}
-%  \nabla^2f_u(\dd)_{uu} &=& \frac{15Gd_u^3}{|\dd|^7} - \frac{9Gd_u}{|\dd|^5} \\
-%  \nabla^2f_u(\dd)_{uv} &=& \frac{15Gd_u^2d_v}{|\dd|^7} - \frac{3Gd_v}{|\dd|^5}\\
-%  \nabla^2f_u(\dd)_{vv} &=& \frac{15Gd_u^2d_v}{|\dd|^7} - \frac{3Gd_u}{|\dd|^5} \\
-%  \nabla^2f_u(\dd)_{vw} &=& \frac{15Gd_ud_vd_w}{|\dd|^7} 
-% \end{eqnarray}
-% 
-% Keeping only the $xx$ term of the Hessian matrix in the Taylor expansion and introducing $\sigma_{xx}^2 = 
-% \sum_im_ip_{i,x}^2$, we get for the accelerations:
-% 
-% \begin{equation}
-%  a_u(\rr) = Mf_u(\rr-\muu) + \frac{1}{2}(\sigma_{xx}^2 - M\mu_x^2)\nabla^2f_u(\dd)_{vv},
-% \end{equation}
-% 
-% with both $v=u$ or $v\neq u$. Expanding this coordinate by coordinate, we get:
-% 
-% \begin{eqnarray}
-%  a_x(\rr) &=& M\frac{G}{|\dd|^3} d_x + \frac{1}{2}\left(\sigma_{xx}^2 - M\mu_x^2\right)\left(\frac{15Gd_x^3}{|\dd|^7} - 
-% \frac{9Gd_x}{|\dd|^5}\right)\\
-%  a_y(\rr) &=& M\frac{G}{|\dd|^3} d_y + \frac{1}{2}\left(\sigma_{xx}^2 - 
-% M\mu_x^2\right)\left(\frac{15Gd_xd_y^2}{|\dd|^7}- 
-% \frac{3Gd_x}{|\dd|^5}\right) \\
-%  a_z(\rr) &=& M\frac{G}{|\dd|^3} d_z + \frac{1}{2}\left(\sigma_{xx}^2 - 
-% M\mu_x^2\right)\left(\frac{15Gd_xd_z^2}{|\dd|^7}- 
-% \frac{3Gd_x}{|\dd|^5}\right)
-% \end{eqnarray}
-% 
-% The quantities $M$, $\muu$ and $\sigma_{xx}^2$ can be constructed on-the-fly by adding particles to the previous total.
-
-%-------------------------------------------------------------------------------------------------
+\section{Recursive construction of the quadrupoles}
 
-% \begin{equation}
-%  \phi(\rr) = - \sum_{i=1}^N \frac{Gm_i}{|\rr - \p{i}|} = - \sum_{i=1}^N m_i f(\rr - \p{i}),
-% \end{equation}
-% 
-% with the function $f(\rr)=G/|\rr|$. The gradient and Hessian matrix of $f$ read:
-% 
-% \begin{equation}
-%  \nabla f(\rr) = \frac{G}{|\rr|^3}\rr, \quad \nabla^2 f(\rr) = G\left(
-%  \begin{array}{ccc}
-%   \frac{3r_x^2}{|\rr|^5} - \frac{1}{|\rr|^3} & \frac{3r_xr_y}{|\rr|^5} & \frac{3r_xr_z}{|\rr|^5} \\
-%   \frac{3r_yr_x}{|\rr|^5} & \frac{3r_y^2}{|\rr|^5} - \frac{1}{|\rr|^3} & \frac{3r_yr_z}{|\rr|^5} \\
-%   \frac{3r_zr_x}{|\rr|^5} &\frac{3r_zr_y}{|\rr|^5}&\frac{3r_z^2}{|\rr|^5} - \frac{1}{|\rr|^3} \\
-%  \end{array}
-%  \right)
-% \end{equation}
-% 
-% 
-% 
-% We define two quantities to simplify the notations: the total mass $M$ of the set of point masses and the centre of 
-% mass $\muu=(\mu_x, \mu_y, \mu_z)$ of this set:
-% 
-% \begin{equation}
-%  M = \sum_{i=1}^N m_i, \qquad \muu = \frac{1}{M} \sum_{i=1}^N m_i\p{i}
-% \end{equation}
-% 
-% Expanding the potential around $\muu$, we find
-% 
-% \begin{equation}
-%   \phi(\rr) \approx -\sum_{i=1}^N m_i f(\rr - \muu) -  \sum_{i=1}^N \frac{m_i}{2} (\p{i} - \muu) \cdot \nabla^2 
-% f(\rr - \muu) \cdot (\p{i} - \muu)
-% \end{equation}
-% 
-% Note that the first order term, the ``dipole'', is identically zero and not shown here. Re-arranging the terms, we 
-% get:
-% 
-% \begin{equation}
-%   \phi(\rr) \approx -Mf(\rr - \muu) - \frac{1}{2}\sum_{i=1}^N m_i \p{i}\cdot \nabla^2 
-% f(\rr - \muu) \cdot \p{i} + \frac{M}{2} \muu\cdot \nabla^2 
-% f(\rr - \muu) \cdot \muu \nonumber
-% \end{equation}
-% 
-% 
-% Let's now assume that on average $|x_{i,x} - \mu_x| \gg |x_{i,y} - \mu_y| \approx |x_{i,z} - \mu_z|$, then the only 
-% term in the matrix that needs to be computed is $\nabla^2 f(\rr-\muu)_{xx}$. The expression then reduces to
-% 
-% \begin{equation}
-%  \phi(\rr) \approx -Mf(\rr - \muu) - \frac{1}{2}\sum_{i=1}^N m_i x_{i,x}^2 \nabla^2f(\rr - \muu)_{xx} + 
-% \frac{M}{2} \mu_x^2 \nabla^2f(\rr - \muu)_{xx}
-% \end{equation}
-% 
-% We can introduce the quantity $d=\sum_{i=1}^N m_i x_{i,x}^2$ to simplify the expression even more:
-% 
-% \begin{equation}
-%  \phi(\rr) \approx -Mf(\rr - \muu) - \left(\frac{d}{2} - \frac{M\mu_x^2}{2}\right) \nabla^2f(\rr - \muu)_{xx}
-% \end{equation}
-% 
-% Using the definition of $f$, we get:
-% 
-% \begin{equation}
-%  \phi(\rr) \approx -\frac{GM}{|\rr - \muu|} - G\left(\frac{d}{2} - \frac{M\mu_x^2}{2}\right) \left(
-% \frac{3(r_x-\mu_x)^2}{|\rr-\muu|^5} - \frac{1}{|\rr-\muu|^3}\right).
-% \end{equation}
-% 
-% The acceleration created by the set of particles on a test particle at position $\rr$ is then
-% 
-% \begin{eqnarray}
-%  \mathbf{a}(\rr)&=& -\nabla\phi(\rr) \nonumber \\
-%  &\approx& - \frac{GM}{|\rr - \muu|^3}(\rr-\muu) \nonumber \\
-%  & &- G\left(\frac{d}{2} - \frac{M\mu_x^2}{2}\right)\left(\frac{15(r_x-\mu_x)^2}{|\rr-\muu|^7} - 
-% \frac{1}{|\rr-\muu|^5}\right) (\rr -\muu) \nonumber \\
-% & & + G\left(\frac{d}{2} - \frac{M\mu_x^2}{2}\right) \frac{6(r_x-\mu_x)^2}{|\rr-\muu|^7} \mathbf{e}_x,
-% \end{eqnarray}
-% 
-% where $\mathbf{e}_x$ is a unit vector along the x-axis. The quantities $M$, $\muu$ and $d$ can be constructed 
-% on-the-fly by adding particles to the previous total.
 \end{document}
diff --git a/examples/theory/multipoles_old.tex b/examples/theory/multipoles_old.tex
new file mode 100644
index 0000000..0ae9db4
--- /dev/null
+++ b/examples/theory/multipoles_old.tex
@@ -0,0 +1,245 @@
+\documentclass[a4paper,10pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage{amsmath}
+\usepackage{amssymb}
+
+\newcommand{\rr}{\mathbf{r}}
+\newcommand{\dd}{\mathbf{d}}
+\newcommand{\vv}{\mathbf{v}}
+\newcommand{\p}[1]{\mathbf{p}_#1}
+\newcommand{\acc}{\mathbf{a}}
+\newcommand{\muu}{\boldsymbol{\mu}}
+
+\title{Multipole expansion}
+
+\begin{document}
+
+\maketitle
+
+Bold quantities are vectors. The indices $\alpha,\beta$ run over the directions $x,y,z$.
+
+For a set of particles at position $\p{i}=(x_i, y_i, z_i)$ with mass $m_i$, we can construct the total mass of the set 
+$M$ and its centre of mass $\muu=(\mu_x, \mu_y, \mu_z)$:
+
+\begin{equation}
+ M = \sum_{i=1}^N m_i, \qquad \muu = \frac{1}{M} \sum_{i=1}^N m_i\p{i}.
+\end{equation}
+
+For each particle, we can define $r_i= \sqrt{(x_i-\mu_x)^2 + (y_i-\mu_y)^2 + (z_i-\mu_z)^2}$, the distance to the 
+centre of mass.\\
+The quadrupole moment $\overline{\overline{Q}}$ of this mass distribution is then written as
+
+\begin{eqnarray}
+ Q_{xx} &=& \sum_i m_i \left( 3(x_i-\mu_x)^2 - r_i^2\right) \\
+ Q_{xy} &=& \sum_i m_i 3(x_i-\mu_x)(y_i-\mu_y) \\
+ Q_{xz} &=& \sum_i m_i 3(x_i-\mu_x)(z_i-\mu_z) \\
+ ~& & \nonumber\\
+ Q_{yx} &=& \sum_i m_i 3(y_i-\mu_y)(x_i-\mu_x) \\
+ Q_{yy} &=& \sum_i m_i \left( 3(y_i-\mu_y)^2 - r_i^2\right) \\
+ Q_{yz} &=& \sum_i m_i 3(y_i-\mu_y)(z_i-\mu_z) \\
+ ~& &\nonumber\\
+ Q_{zx} &=& \sum_i m_i 3(z_i-\mu_z)(x_i-\mu_x)  \\
+ Q_{zy} &=& \sum_i m_i 3(z_i-\mu_z)(y_i-\mu_y) \\
+ Q_{zz} &=& \sum_i m_i \left( 3(z_i-\mu_z)^2 - r_i^2\right)
+\end{eqnarray}
+Note that this matrix is symmetric and traceless ($Q_{xx}+Q_{yy}+Q_{zz}=0$), so only $5$ components need to be 
+computed. \\
+
+\pagebreak
+
+The potential created by this distribution of mass on a point at position $\rr=(r_x,r_y,r_z)$ from the centre of mass 
+$\muu$ is then:
+
+\begin{equation}
+ \phi(\rr) = -\frac{G}{|\rr|}M + \frac{1}{2}\frac{G}{|\rr|^5}\sum_{\alpha,\beta} r_\alpha r_\beta Q_{\alpha\beta}
+\end{equation}
+
+Taking the gradient of this expression to get the acceleration, we obtain:
+
+\begin{equation}
+ -\nabla_\gamma\phi(\rr)= -\frac{GM}{|\rr|^3} r_\gamma + 
+\frac{1}{2}G\sum_{\alpha,\beta}\left(\frac{\delta_{\alpha\gamma}r_\beta}{|\rr|^5} + 
+\frac{\delta_{\gamma\beta}r_\alpha}{|\rr|^5} -5\frac{r_\alpha r_\beta r_\gamma}{|\rr|^7}\right)Q_{\alpha\beta}
+\end{equation}
+
+Writing this explicitly for each coordinate (using the symmetry of $\overline{\overline{Q}}$), we get:
+
+\begin{eqnarray}
+ a_x &=& -\frac{Gr_x}{|\rr|^3}M + \frac{1}{2} \frac{G}{|\rr|^5}\left[2r_x Q_{xx} + 2r_yQ_{xy} + 2r_z Q_{xz}\right] - 
+\frac{5}{2} \frac{G r_x}{|\rr|^7} \Xi, \\
+ a_y &=& -\frac{Gr_y}{|\rr|^3}M + \frac{1}{2} \frac{G}{|\rr|^5}\left[2r_x Q_{yx} + 2r_yQ_{yy} + 2r_z Q_{yz}\right] - 
+\frac{5}{2} \frac{G r_y}{|\rr|^7} \Xi, \\
+ a_z &=& -\frac{Gr_z}{|\rr|^3}M + \frac{1}{2} \frac{G}{|\rr|^5}\left[2r_x Q_{zx} + 2r_yQ_{zy} + 2r_z Q_{zz}\right] - 
+\frac{5}{2} \frac{G r_z}{|\rr|^7} \Xi, \\
+\end{eqnarray}
+with the common coefficient $\Xi$ defined as:
+
+\begin{equation*}
+ \Xi = \rr\cdot(\overline{\overline{Q}}\cdot \rr) = r_x^2Q_{xx} + r_y^2Q_{yy} + r_z^2Q_{zz} + 2r_xr_yQ_{xy} + 
+2r_xr_zQ_{xz} + 2r_yr_zQ_{yz}
+\end{equation*}
+
+
+
+% 
+% Let's consider the gravitational acceleration $\acc=(a_x,a_y,a_z)$ that a set of point masses at position 
+% $\p{i}=(p_{i,x}, p_{i,y}, p_{i,z})$ with masses $m_i$ generate at position $\rr=(r_x, r_y, r_z)$:
+% 
+% \begin{equation}
+%  \acc(\rr) = \sum_i \frac{Gm_i}{|\rr - \p{i}|^3}(\rr - \p{i})
+% \end{equation}
+% 
+% This expression can split in one expression for each of the three spatial coordinates $u=x,y,z$:
+% 
+% \begin{equation}
+%  a_u (\rr) = \sum_i \frac{m_i G}{|\rr-\p{i}|^3} ( r_u - p_{i,u}) = \sum_i m_i f_u(\rr - \p{i}),
+% \end{equation}
+% 
+% with 
+% 
+% \begin{equation}
+%  f_u (\vv) = \frac{G}{|\vv|^3} v_u.
+% \end{equation}
+% 
+% We define two quantities to simplify the notations: the total mass $M$ of the set of point masses and the centre of 
+% mass $\muu=(\mu_x, \mu_y, \mu_z)$ of this set:
+% 
+% \begin{equation}
+%  M = \sum_{i=1}^N m_i, \qquad \muu = \frac{1}{M} \sum_{i=1}^N m_i\p{i}.
+% \end{equation}
+% 
+% 
+% We can now expand the functions $f_u$ around the vector linking the particle to the centre of mass $\rr-\muu$:
+% 
+% \begin{eqnarray}
+%  a_u(\rr) &\approx& \sum_i m_i f_u(\rr - \muu) \\
+%   & & + \sum_i m_i (\p{i}-\muu) \cdot \nabla f_u(\rr - \muu)\\
+%   & & + \frac{1}{2}\sum_i m_i (\p{i}-\muu)\cdot \nabla^2 f_u(\rr-\muu)\cdot (\p{i} - \muu).
+% \end{eqnarray}
+% 
+% The first order term is identically zero and can hence be dropped. Re-arranging some of the terms, introducing the 
+% vector $\dd=\rr-\muu=(d_x,d_y,d_z)$ and using the fact that the Hessian matrix of $f_u$ is symmetric, we get:
+% \begin{eqnarray}
+%  a_u(\rr) &=& Mf_u(\dd) \\
+%            & & + \frac{1}{2} \sum_i m_i \p{i}\cdot \nabla^2f_u(\dd)\cdot \p{i}  \\
+%            & & + \frac{1}{2} M \muu\cdot \nabla^2f_u(\dd)\cdot \muu \\
+%            & & - \sum_i m_i \p{i}\cdot \nabla^2f_u(\dd)\cdot \muu \\
+%            &=& Mf_u(\dd) \\
+%            & & + \frac{1}{2} \sum_i m_i \p{i}\cdot \nabla^2f_u(\dd)\cdot \p{i} \\
+%            & &- \frac{1}{2} M \muu\cdot \nabla^2f_u(\dd)\cdot \muu
+% \end{eqnarray}
+% 
+% The gradient of $f_u$ reads
+% \begin{equation}
+%  \nabla f_u(\dd) = \frac{-3Gd_u}{|\dd|^5}\dd + \frac{G}{|\dd|^3}\hat{\mathbf{e}}_u,
+% \end{equation}
+% 
+% with $\hat{\mathbf{e}}_u$ a unit vector along the $u$-axis. The different components of the Hessian matrix then read:
+% 
+% \begin{eqnarray}
+%  \nabla^2f_u(\dd)_{uu} &=& \frac{15Gd_u^3}{|\dd|^7} - \frac{9Gd_u}{|\dd|^5} \\
+%  \nabla^2f_u(\dd)_{uv} &=& \frac{15Gd_u^2d_v}{|\dd|^7} - \frac{3Gd_v}{|\dd|^5}\\
+%  \nabla^2f_u(\dd)_{vv} &=& \frac{15Gd_u^2d_v}{|\dd|^7} - \frac{3Gd_u}{|\dd|^5} \\
+%  \nabla^2f_u(\dd)_{vw} &=& \frac{15Gd_ud_vd_w}{|\dd|^7} 
+% \end{eqnarray}
+% 
+% Keeping only the $xx$ term of the Hessian matrix in the Taylor expansion and introducing $\sigma_{xx}^2 = 
+% \sum_im_ip_{i,x}^2$, we get for the accelerations:
+% 
+% \begin{equation}
+%  a_u(\rr) = Mf_u(\rr-\muu) + \frac{1}{2}(\sigma_{xx}^2 - M\mu_x^2)\nabla^2f_u(\dd)_{vv},
+% \end{equation}
+% 
+% with both $v=u$ or $v\neq u$. Expanding this coordinate by coordinate, we get:
+% 
+% \begin{eqnarray}
+%  a_x(\rr) &=& M\frac{G}{|\dd|^3} d_x + \frac{1}{2}\left(\sigma_{xx}^2 - M\mu_x^2\right)\left(\frac{15Gd_x^3}{|\dd|^7} - 
+% \frac{9Gd_x}{|\dd|^5}\right)\\
+%  a_y(\rr) &=& M\frac{G}{|\dd|^3} d_y + \frac{1}{2}\left(\sigma_{xx}^2 - 
+% M\mu_x^2\right)\left(\frac{15Gd_xd_y^2}{|\dd|^7}- 
+% \frac{3Gd_x}{|\dd|^5}\right) \\
+%  a_z(\rr) &=& M\frac{G}{|\dd|^3} d_z + \frac{1}{2}\left(\sigma_{xx}^2 - 
+% M\mu_x^2\right)\left(\frac{15Gd_xd_z^2}{|\dd|^7}- 
+% \frac{3Gd_x}{|\dd|^5}\right)
+% \end{eqnarray}
+% 
+% The quantities $M$, $\muu$ and $\sigma_{xx}^2$ can be constructed on-the-fly by adding particles to the previous total.
+
+%-------------------------------------------------------------------------------------------------
+
+% \begin{equation}
+%  \phi(\rr) = - \sum_{i=1}^N \frac{Gm_i}{|\rr - \p{i}|} = - \sum_{i=1}^N m_i f(\rr - \p{i}),
+% \end{equation}
+% 
+% with the function $f(\rr)=G/|\rr|$. The gradient and Hessian matrix of $f$ read:
+% 
+% \begin{equation}
+%  \nabla f(\rr) = \frac{G}{|\rr|^3}\rr, \quad \nabla^2 f(\rr) = G\left(
+%  \begin{array}{ccc}
+%   \frac{3r_x^2}{|\rr|^5} - \frac{1}{|\rr|^3} & \frac{3r_xr_y}{|\rr|^5} & \frac{3r_xr_z}{|\rr|^5} \\
+%   \frac{3r_yr_x}{|\rr|^5} & \frac{3r_y^2}{|\rr|^5} - \frac{1}{|\rr|^3} & \frac{3r_yr_z}{|\rr|^5} \\
+%   \frac{3r_zr_x}{|\rr|^5} &\frac{3r_zr_y}{|\rr|^5}&\frac{3r_z^2}{|\rr|^5} - \frac{1}{|\rr|^3} \\
+%  \end{array}
+%  \right)
+% \end{equation}
+% 
+% 
+% 
+% We define two quantities to simplify the notations: the total mass $M$ of the set of point masses and the centre of 
+% mass $\muu=(\mu_x, \mu_y, \mu_z)$ of this set:
+% 
+% \begin{equation}
+%  M = \sum_{i=1}^N m_i, \qquad \muu = \frac{1}{M} \sum_{i=1}^N m_i\p{i}
+% \end{equation}
+% 
+% Expanding the potential around $\muu$, we find
+% 
+% \begin{equation}
+%   \phi(\rr) \approx -\sum_{i=1}^N m_i f(\rr - \muu) -  \sum_{i=1}^N \frac{m_i}{2} (\p{i} - \muu) \cdot \nabla^2 
+% f(\rr - \muu) \cdot (\p{i} - \muu)
+% \end{equation}
+% 
+% Note that the first order term, the ``dipole'', is identically zero and not shown here. Re-arranging the terms, we 
+% get:
+% 
+% \begin{equation}
+%   \phi(\rr) \approx -Mf(\rr - \muu) - \frac{1}{2}\sum_{i=1}^N m_i \p{i}\cdot \nabla^2 
+% f(\rr - \muu) \cdot \p{i} + \frac{M}{2} \muu\cdot \nabla^2 
+% f(\rr - \muu) \cdot \muu \nonumber
+% \end{equation}
+% 
+% 
+% Let's now assume that on average $|x_{i,x} - \mu_x| \gg |x_{i,y} - \mu_y| \approx |x_{i,z} - \mu_z|$, then the only 
+% term in the matrix that needs to be computed is $\nabla^2 f(\rr-\muu)_{xx}$. The expression then reduces to
+% 
+% \begin{equation}
+%  \phi(\rr) \approx -Mf(\rr - \muu) - \frac{1}{2}\sum_{i=1}^N m_i x_{i,x}^2 \nabla^2f(\rr - \muu)_{xx} + 
+% \frac{M}{2} \mu_x^2 \nabla^2f(\rr - \muu)_{xx}
+% \end{equation}
+% 
+% We can introduce the quantity $d=\sum_{i=1}^N m_i x_{i,x}^2$ to simplify the expression even more:
+% 
+% \begin{equation}
+%  \phi(\rr) \approx -Mf(\rr - \muu) - \left(\frac{d}{2} - \frac{M\mu_x^2}{2}\right) \nabla^2f(\rr - \muu)_{xx}
+% \end{equation}
+% 
+% Using the definition of $f$, we get:
+% 
+% \begin{equation}
+%  \phi(\rr) \approx -\frac{GM}{|\rr - \muu|} - G\left(\frac{d}{2} - \frac{M\mu_x^2}{2}\right) \left(
+% \frac{3(r_x-\mu_x)^2}{|\rr-\muu|^5} - \frac{1}{|\rr-\muu|^3}\right).
+% \end{equation}
+% 
+% The acceleration created by the set of particles on a test particle at position $\rr$ is then
+% 
+% \begin{eqnarray}
+%  \mathbf{a}(\rr)&=& -\nabla\phi(\rr) \nonumber \\
+%  &\approx& - \frac{GM}{|\rr - \muu|^3}(\rr-\muu) \nonumber \\
+%  & &- G\left(\frac{d}{2} - \frac{M\mu_x^2}{2}\right)\left(\frac{15(r_x-\mu_x)^2}{|\rr-\muu|^7} - 
+% \frac{1}{|\rr-\muu|^5}\right) (\rr -\muu) \nonumber \\
+% & & + G\left(\frac{d}{2} - \frac{M\mu_x^2}{2}\right) \frac{6(r_x-\mu_x)^2}{|\rr-\muu|^7} \mathbf{e}_x,
+% \end{eqnarray}
+% 
+% where $\mathbf{e}_x$ is a unit vector along the x-axis. The quantities $M$, $\muu$ and $d$ can be constructed 
+% on-the-fly by adding particles to the previous total.
+\end{document}
-- 
GitLab