multipoles.tex

\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}[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}

\begin{document}

\maketitle

Bold quantities are vectors, underlined quantities are matrices. The indices 
$\alpha,\beta$ run over the directions $x,y,z$. $\identity$ is the identity 
matrix in 3D.

\section{Construction of multipoles}

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}
 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}

For a set of particles $A$, the multipoles can be computed around the 
centre of mass $\muu_A$. The first LHS term uses Dehnen's notation.\\

Monopole:
\begin{equation}
M_{(0,0,0)}(\muu_A) = M_{{\rm tot}, A} 
\end{equation}

Dipole:
\begin{eqnarray}
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)}(\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)}(\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)}(\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)}(\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)}(\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)}(\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 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
$\muu_{A}$ of the system consisting of all the particles contained in each 
individual (disjoint) sub-sets $B$. This allows to construct the multipoles in 
the tree recursively.\\

Monopole:
\begin{equation}
 M_{{\rm tot},A} = \sum_{B\in A} M_{{\rm tot},B}
\end{equation}

Dipole:
\begin{equation}
 P_{A,\alpha} = 0
\end{equation}

Quadrupole:
\begin{equation}
 I_{A,\alpha\beta} = \sum_{B\in A}\left( I_{B,\alpha\beta} + 
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}


\section{Derivatives of the potential}

The gravitational potential is given by $\phi(\rr) = G/|\rr|$. Its derivatives 
are given by the following expressions. The first LHS term uses Dehnen's 
notation.\\

0-th order:
\begin{equation}
 D_{(0,0,0)}(\rr) = \phi(\rr) = \frac{G}{|\rr|}
\end{equation}

1-st order:
\begin{eqnarray}
D_{(1,0,0)}(\rr) &=& \partial_x \phi(\rr)~=-\frac{G}{|\rr|^3}r_x\\
D_{(0,1,0)}(\rr) &=& \partial_y \phi(\rr)~=-\frac{G}{|\rr|^3}r_y\\
D_{(0,0,1)}(\rr) &=& \partial_z \phi(\rr)~=-\frac{G}{|\rr|^3}r_z
\end{eqnarray}

2-nd order:
\begin{eqnarray}
D_{(2,0,0)} &=& \partial_{xx} \phi(\rr)~= \frac{3Gr_x^2}{|\rr|^5} - 
\frac{G}{|\rr|^3}\\
D_{(0,2,0)} &=& \partial_{yy} \phi(\rr)~= \frac{3Gr_y^2}{|\rr|^5} - 
\frac{G}{|\rr|^3}\\
D_{(0,0,2)} &=& \partial_{zz} \phi(\rr)~= \frac{3Gr_z^2}{|\rr|^5} - 
\frac{G}{|\rr|^3}\\
D_{(1,1,0)} &=& \partial_{xy} \phi(\rr)~= \frac{3Gr_xr_y}{|\rr|^5} \\
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}

3-rd order:
\begin{eqnarray}
 D_{(3,0,0)} &=& \partial_{xxx} \phi(\rr)~= -\frac{15Gr_x^3}{|\rr|^7} +
\frac{9Gr_x}{|\rr|^5}\\
D_{(1,2,0)} &=& \partial_{xyy} \phi(\rr)~= -\frac{15Gr_xr_y^2}{|\rr|^7} +
\frac{3Gr_x}{|\rr|^5}\\
D_{(1,0,2)} &=& \partial_{xzz} \phi(\rr)~= -\frac{15Gr_xr_z^2}{|\rr|^7} +
\frac{3Gr_x}{|\rr|^5}\\
D_{(2,1,0)} &=& \partial_{xxy} \phi(\rr)~= -\frac{15Gr_x^2r_y}{|\rr|^7} +
\frac{3Gr_y}{|\rr|^5}\\
D_{(2,0,1)} &=& \partial_{xxz} \phi(\rr)~= -\frac{15Gr_x^2r_z}{|\rr|^7} +
\frac{3Gr_z}{|\rr|^5}\\
D_{(1,1,1)} &=& \partial_{xyz} \phi(\rr)~= -\frac{15Gr_xr_yr_z}{|\rr|^7}
\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_{\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 
all vectors with $|{\bf n}|\leq2$) and writing $\rr_{i,A} = \muu_A - \p{i}$, we 
get:

\begin{eqnarray}
 \phi(\p{i}) &=& -M_{{\rm tot},A} \frac{G}{|\rr_{i,A}|} -\frac{1}{2}\sum_\alpha 
I_{A,\alpha} \left(\frac{3Gr_\alpha^2}{|\rr_{i,A}|^5} - 
\frac{G}{|\rr_{i,A}|^3}\right) \\ 
& &+ \frac{1}{2}\sum_{\alpha,\beta} (1-\delta_{\alpha\beta})I_{A,\alpha\beta} 
\frac{3Gr_\alpha r_\beta}{|\rr_{i,A}|^5} \\
             &=& -M_{{\rm tot},A} \frac{G}{|\rr_{i,A}|} +
\frac{1}{2}\left( I_{A,xx} + I_{A,yy} +I_{A,zz}\right)\frac{G}{|\rr_{i,A}|^3}  
\\
 & & - \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}(\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}(\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}$).\\

The accelerations can also written using Dehnen's short notation:

\begin{eqnarray}
 a_{i,x} &=& \sum_{\bf n} M_{\bf n}(\muu_A) D_{{\bf n}+(1,0,0)}(\muu_A - 
\p{i}) \\
 a_{i,y} &=& \sum_{\bf n} M_{\bf n}(\muu_A) D_{{\bf n}+(0,1,0)}(\muu_A - 
\p{i}) \\
 a_{i,z} &=& \sum_{\bf n} M_{\bf n}(\muu_A) D_{{\bf n}+(0,0,1)}(\muu_A - 
\p{i})
\end{eqnarray}


\section{FMM Field tensors for potential}

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_{BA} &=& \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{G}{2}\frac{{\rm 
tr}(\ii{A})}{|\rr_{BA}|^3} + 
\frac{3G}{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_{BA,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_{BA,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_{BA,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_{BA} &=&  \frac{GM_{{\rm tot}, A}}{|\rr_{BA}|} -\frac{G}{2}\frac{{\rm 
tr}(\ii{A})}{|\rr_{BA}|^3} + 
\frac{3G}{2}\frac{\rr_{BA}^T \cdot 
\ii{A} \cdot \rr_{BA}}{|\rr_{BA}|^5} \\
\qq{BA} &=&  \frac{GM_{{\rm tot}, A}}{|\rr_{BA}|^3} \rr_{BA}\\
\jj{BA} &=& \frac{3GM_{{\rm tot}, A}}{|\rr_{BA}|^5} \rr_{BA} 
\otimes\rr_{BA} - \frac{GM_{{\rm tot}, A}}{|\rr_{BA}|^3} \identity
\end{eqnarray}

Note that $\jj{B}$ is symmetric, meaning that $10$ additional numbers per cell 
only have to be stored. That is a total of $10$ for the multipoles and $10$ for 
the field tensors.

\section{Gravitational potential in FMM}

The potential at the position $\p{i}$ in the vicinity of $\muu_B$ due to the 
mass clustered around $\muu_A$ can be computed using the field tensors created 
by the mutipoles at position $\muu_A$ on the position $\muu_B$.

\begin{equation}
 \phi_A(\p{i}) = \sum_{|{\bf n}| \leq2} \frac{1}{ {\bf n}!} (\p{i} - 
\muu_B)^{\bf n} F_{\bf n}(\muu_B)
\end{equation}

Setting $\dd_{i,B}= \p{i} - \muu_B$ and $\rr_{BA} = \muu_B - \muu_A$, this sum 
expands into

\begin{eqnarray}
 \phi_A(\p{i}) &=& N_{BA} + \dd_{i,B}\cdot\qq{BA} + \dd_{i,B}^T \cdot \jj{BA} 
\cdot \dd_{i,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} \\
& &+ \frac{GM_{{\rm tot}, A}}{|\rr_{BA}|^3} \rr_{BA}\cdot\dd_{i,B} \\
& &+\frac{3}{2}\frac{GM_{{\rm tot}, A}}{|\rr_{BA}|^5} \left(\dd_{i,B}^T \cdot 
\rr_{BA}\otimes\rr_{BA}\cdot \dd_{i,B}\right) \\
& & - \frac{1}{2}\frac{GM_{{\rm tot}, A}}{|\rr_{BA}|^3} \left(\dd_{i,B}^T 
\cdot\identity \cdot \dd_{i,B}\right)
\end{eqnarray}

Note that if $\p{i}=\muu_B$ (i.e. $d_{i,B}={\bf 0}$), the expression for the 
potential in the B-H formalism is recovered.

\end{document}