A Very Fast And Angular Momentum Conserving Tree Code

TL;DR

This paper presents a modification to the Cartesian fast multipole method that ensures conservation of both angular and linear momentum in astrophysical gravitational simulations, improving accuracy over traditional methods.

Contribution

The authors introduce a novel modification to the fast multipole method that guarantees conservation of angular momentum, a feature not present in previous implementations.

Findings

The modified tree code conserves angular momentum to machine precision.

Linear momentum conservation is maintained alongside angular momentum.

The method is efficient and applicable to astrophysical simulations.

Abstract

There are many methods used to compute the classical gravitational field in astrophysical simulation codes. With the exception of the typically impractical method of direct computation, none ensure conservation of angular momentum to machine precision. Under uniform time-stepping, the Cartesian fast multipole method of Dehnen (also known as the very fast tree code) conserves linear momentum to machine precision. We show it is possible to modify this method in a way that conserves both angular and linear momenta.

Equations

$$|{\mathbf{Z}}_{\mathbf{A}}-{\mathbf{Z}}_{\mathbf{B}}|\ge \frac{1}{\theta }\left(R_{A,\max }+R_{B,\max }\right),$$

$$M:=\sum _{{\mathbf{X}}_{\mathbf{n}}\in \mathrm{cell}}{\mu }_n,$$

$$M_i:=\sum _{{\mathbf{X}}_{\mathbf{n}}\in \mathrm{cell}}{\mu }_n{x}_{n,i},$$

$$M_{ij}:=\sum _{{\mathbf{X}}_{\mathbf{n}}\in \mathrm{cell}}{\mu }_n{x}_{n,i}{x}_{n,j}.$$

$$D:=-\frac{1}{R},$$

$$D_i:=\frac{R_i}{R^3},$$

$$D_{ij}:=-\frac{3R_i{R}_j-{\delta }_{ij}{R}^2}{R^5},$$

$$D_{ijk}:=\frac{15R_i{R}_j{R}_k-3\left({\delta }_{ij}{R}_k+{\delta }_{jk}{R}_i+{\delta }_{ki}{R}_j\right)R^2}{R^7}.$$

$$\begin{gathered} {\Phi }_{B\to A}\left(\mathbf{X}\right)\approx \left(M_{B}D-M_{Bi}{D}_i+\frac{1}{2}{M}_{Bij}{D}_{ij}\right)+ \\ {x}_i\left(M_B{D}_i-M_{Bj}{D}_{ij}+\frac{1}{2}{M}_{Bjk}{D}_{ijk}\right)+\frac{1}{2}{x}_i{x}_j\left(M_B{D}_{ij}-M_{Bk}{D}_{ijk}\right)+\frac{1}{6}{x}_i{x}_j{x}_k{M}_B{D}_{ijk}. \end{gathered}$$

$$\begin{gathered} {\mathbf{g}}_{B\to A}\left(\mathbf{X}\right)=-[\left(M_B{D}_i-M_{Bj}{D}_{ij}+\frac{1}{2}{M}_{Bjk}{D}_{ijk}\right)+ \\ {x}_j\left(M_B{D}_{ij}-M_{Bk}{D}_{ijk}\right)+\frac{1}{2}{x}_j{x}_k{M}_B{D}_{ijk}]{\hat{\mathbf{e}}}_{\mathbf{i}}. \end{gathered}$$

$$\begin{gathered} {\mathbf{g}}_{A\to B}\left(\mathbf{Y}\right)=-[\left(-M_A{D}_i-M_{Aj}{D}_{ij}-\frac{1}{2}{M}_{Ajk}{D}_{ijk}\right)+ \\ {y}_j\left(M_A{D}_{ij}+M_{Ak}{D}_{ijk}\right)-\frac{1}{2}{y}_j{y}_k{M}_A{D}_{ijk}]{\hat{\mathbf{e}}}_{\mathbf{i}}. \end{gathered}$$

$${\mathbf{g}}_{\mp \to \pm }=\mp \left(D_i+\left(x_j-y_j\right)D_{ij}+\left(x_j-y_j\right)\left(x_k-y_k\right)D_{ijk}\right){\hat{\mathbf{e}}}_{\mathbf{i}}.$$

$${\tau }_{\mathbf{AB}}=\frac{1}{2}{ϵ}_{piq}\left(M_{Apjk}{M}_B-M_A{M}_{Bpjk}\right)D_{ijk}{\hat{\mathbf{e}}}_{\mathbf{q}},$$

$$M_{pjk}:=\sum _{{\mathbf{X}}_{\mathbf{n}}\in \mathrm{cell}}{\mu }_n{x}_{n,p}{x}_{n,j}{x}_{n,k}.$$

$${\mathbf{F}}_{\mathbf{c}}=-\frac{1}{2}\left(M_{Ajkl}{M}_B-M_A{M}_{Bjkl}\right)D_{ijkl}{\hat{\mathbf{e}}}_{\mathbf{i}},$$

$${\mathbf{g}}_{C,B\to A}\left(\mathbf{X}\right)=+\frac{1}{M_A}{\mathbf{F}}_{\mathbf{c}}$$

$${\mathbf{g}}_{C,A\to B}\left(\mathbf{Y}\right)=-\frac{1}{M_B}{\mathbf{F}}_{\mathbf{c}}.$$

$$\begin{gathered} {D^{\prime }}_{ijkl}:= \\ +\frac{15\left({\delta }_{ij}{R}_k{R}_l+{\delta }_{ik}{R}_j{R}_l+{\delta }_{il}{R}_j{R}_k\right)}{R^7} \\ -\frac{3\left({\delta }_{ij}{\delta }_{kl}+{\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)}{R^5}, \end{gathered}$$

$${\mathbf{F}}_{\mathbf{c}}^{\prime }=-\frac{1}{2}\left(M_{Ajkl}{M}_B-M_A{M}_{Bjkl}\right){D^{\prime }}_{ijkl}{\hat{\mathbf{e}}}_{\mathbf{i}}.$$

$$\sqrt{\frac{{\sum }_d{\left({\sum }_n{\mu }_n{A}_{n,d}\right)}^2}{{\sum }_d{\left({\sum }_n{\mu }_n\left|A_{n,d}\right|\right)}^2}},$$

$${ϵ}_n=\frac{\left|{\mathbf{g}}_{\mathbf{n}}-{\mathbf{g}}_{\mathbf{PP},\mathbf{n}}\right|}{\left|{\mathbf{g}}_{\mathbf{PP},\mathbf{n}}\right|},$$

$${ϵ}_n=\frac{\left|\mathbf{X}\times {\mathbf{g}}_{\mathbf{n}}-\mathbf{X}\times {\mathbf{g}}_{\mathbf{PP},\mathbf{n}}\right|}{\left|\mathbf{X}\times {\mathbf{g}}_{\mathbf{PP},\mathbf{n}}\right|},$$

$$D_{l_1{l}_2...l_n}:=\frac{\partial }{\partial r_{l_1}}\frac{\partial }{\partial r_{l_2}}...\frac{\partial }{\partial r_{l_n}}G\left(R\right).$$

$${\Phi }_{i,B\to A}=\sum _{{\mathbf{Y}}_{\mathbf{j}}\in \mathrm{cellB}}{\mu }_j\sum _{m=0}^P\sum _{n=0}^{P-m}\frac{{\left(-1\right)}^n}{n!m!}{x}_{i,l_1}{x}_{i,l_2}...x_{i,l_m}{D}_{l_1{l}_2...l_m{q}_1{q}_2...q_n}{y}_{j,q_1}{y}_{j,q_2}...y_{j,q_n}.$$

$${\Phi }_{j,A\to B}=\sum _{{\mathbf{X}}_{\mathbf{i}}\in \mathrm{cellA}}{\mu }_i\sum _{m=0}^P\sum _{n=0}^{P-m}\frac{{\left(-1\right)}^n}{n!m!}{x}_{i,l_1}{x}_{i,l_2}...x_{i,l_m}{D}_{l_1{l}_2...l_m{q}_1{q}_2...q_n}{y}_{j,q_1}{y}_{j,q_2}...y_{j,q_n}.$$

$${\mathbf{g}}_{\mathbf{i},\mathbf{B}\to \mathbf{A}}=+{\mathbf{e}}_k\sum _{{\mathbf{Y}}_{\mathbf{j}}\in \mathrm{cellB}}{\mu }_j\sum _{m=0}^{P-1}\sum _{n=0}^{P-1-m}\frac{{\left(-1\right)}^n}{n!m!}{x}_{i,l_1}{x}_{i,l_2}...x_{i,l_m}{D}_{l_1{l}_2...l_{m}k{q}_1{q}_2...q_n}{y}_{j,q_1}{y}_{j,q_2}...y_{j,q_n}.$$

$${\mathbf{g}}_{\mathbf{j},\mathbf{A}\to \mathbf{B}}=-{\mathbf{e}}_k\sum _{{\mathbf{X}}_{\mathbf{i}}\in \mathrm{cellA}}{\mu }_i\sum _{m=0}^{P-1}\sum _{n=0}^{P-1-m}\frac{{\left(-1\right)}^n}{n!m!}{x}_{i,l_1}{x}_{i,l_2}...x_{i,l_m}{D}_{l_1{l}_2...l_{m}k{q}_1{q}_2...q_n}{y}_{j,q_1}{y}_{j,q_2}...y_{j,q_n}.$$

$$\sum _{{\mathbf{X}}_{\mathbf{i}}\in \mathrm{cellA}}{\mu }_i{\mathbf{g}}_{\mathbf{i},\mathbf{B}\to \mathbf{A}}+\sum _{{\mathbf{Y}}_{\mathbf{j}}\in \mathrm{cellB}}{\mu }_j{\mathbf{g}}_{\mathbf{j},\mathbf{A}\to \mathbf{B}}=0.$$

$$\begin{array}{rl}{\tau }_{\mathbf{AB}}:=& \sum _{{\mathbf{Y}}_{\mathbf{j}}\in \mathrm{cellB}}{\mu }_j{\mathbf{R}}_{\mathbf{j}}\times {\mathbf{g}}_{\mathbf{j},\mathbf{A}\to \mathbf{B}}+\sum _{{\mathbf{X}}_{\mathbf{i}}\in \mathrm{cellA}}{\mu }_i{\mathbf{x}}_{\mathbf{i}}\times {\mathbf{g}}_{\mathbf{i},\mathbf{B}\to \mathbf{A}}+\sum _{{\mathbf{Y}}_{\mathbf{j}}\in \mathrm{cellB}}{\mu }_j{\mathbf{y}}_{\mathbf{j}}\times {\mathbf{g}}_{\mathbf{j},\mathbf{A}\to \mathbf{B}}.\end{array}$$

$$\begin{array}{rl}\sum _{{\mathbf{Y}}_{\mathbf{j}}\in \mathrm{cellB}}{\mu }_j{\mathbf{R}}_{\mathbf{j}}\times {\mathbf{g}}_{\mathbf{j},\mathbf{A}\to \mathbf{B}}=& \\ {\mathbf{e}}_r\sum _{{\mathbf{X}}_{\mathbf{i}}\in \mathrm{cellA}}\sum _{{\mathbf{Y}}_{\mathbf{j}}\in \mathrm{cellB}}& {\mu }_i{\mu }_j\sum _{m=0}^{P-1}\sum _{n=0}^{P-1-m}\frac{{\left(-1\right)}^n}{n!m!}{ϵ}_{pkr}{R}_p{x}_{i,l_1}{x}_{i,l_2}...x_{i,l_m}{D}_{l_1{l}_2...l_{m}k{q}_1{q}_2...q_n}{y}_{j,q_1}{y}_{j,q_2}...y_{j,q_n}\end{array}$$

Full text

1 Introduction

Angular momentum plays an important role in a plethora of astrophysical phenomena. These phenomena include, but are not limited to, the formation of galactic discs, the accretion of matter in interacting binary star systems, the formation of proto-stars and proto-planetary discs, the dynamics of planetary orbits, and rapidly rotating neutron stars and black holes. In the limit that Newtonian gravity provides an accurate description of the gravitational field of any isolated astrophysical system, we expect the total angular momentum change due to the gravitational interaction to be zero.

Classical self gravitating astrophysical simulation codes compute the gravitational field in a variety of ways. Grid based hydrodynamics codes (e.g. Fryxell et al., 2000; Almgren et al., 2010; Stone et al., 2008; Dupuy & Liu, 2012; D’Souza et al., 2006) solve the discretized Poisson’s equation for the gravitational potential using iterative techniques, Fourier transforms, or a combination of the two. N-body codes and smoothed particle hydrodynamics (SPH) codes (e.g. Hernquist & Katz, 1989; Springel, 2005; Vanaverbeke et al., 2009; Lorén-Aguilar et al., 2010; Yokota & Barba, 2012) may use the tree code of Barnes & Hut (1986), particle-mesh methods (Bagla, 2002), or fast multipole methods (FMM) (Greengard & Rokhlin, 1997; Warren & Salmon, 1995; Dehnen, 2000; Dehnen & Read, 2011). Some of the gravity solvers in the aforementioned codes conserve linear momentum, but none ensure conservation of angular momentum.

Due to the symmetry of the equations for the multi-pole interactions, the method of Dehnen (2000) (hereafter “D2000”) naturally conserves linear momentum between any pair of particles. We have developed a modification to D2000 that preserves this property, while simultaneously conserving angular momentum between any pair of interacting multi-poles. This kind of conservation is not of the same quality as the conservation of linear momentum. It introduces artificial torques, however, the added torques are within the error bound of the original scheme. In §2 we describe our modification to D2000. In §3 we provide a numerical test of the method. In §4 we make the case for using this technique to model double white dwarfs (DWDs), as well as discuss some of the method’s shortcomings. In §A we provide a more general derivation of the method that applies to higher orders.

2 Method

The algorithm presented by D2000 decomposes the set of particles into an oct-tree structure, with each cell in the oct-tree containing a predetermined maximum number of particles, $n_{\mathrm{crit}}$. The code presented in this paper was run with $n_{\mathrm{crit}}=25$. Two cells, cell “A” and cell “B”, are considered “well separated” if they satisfy the “opening criterion”,

[TABLE]

where ${\bf{Z_{A}}}$ and ${\bf{Z_{B}}}$ are the respective centers of mass of cells A and B, $R_{A,\mathrm{max}}$ and $R_{B,\mathrm{max}}$ are the maximum distances from a particle within the cells to the centers of mass of their respective cells, and $\theta$ is an adjustable parameter called the “opening angle”, where $0<\theta\leq 1$. Forces within a cell are computed using multipole interactions and Taylor expansions for all cells that are well separated from it, while the force contributions from any remaining nearby particles are computed directly. (Note that Dehnen (2014) has recently developed a more complex selection criteria that uses an error estimate to select interaction pairs in a manner that maximizes execution speed for a given error. The development we present here is also applicable to that method.) Here we will present only what is necessary to describe the modifications we have made. Refer to D2000 for the full description of the original method.

Let ${\bf{R}}={\bf{Z_{B}}}-{\bf{Z_{A}}}$. Within cell A there are $N_{A}$ particles located at positions ${\bf{X_{1}}},{\bf{X_{2}}},...,{\bf{X_{N_{A}}}}$ and with masses $\mu_{A1},\mu_{A2},...,\mu_{AN_{A}}$. Similarly, within cell B there are $N_{B}$ bodies located at positions ${\bf{Y_{1}}},{\bf{Y_{2}}},...,{\bf{Y_{N_{B}}}}$ and with masses $\mu_{B1},\mu_{B2},...,\mu_{BN_{B}}$. Let ${\bf{x_{n}}}={\bf{X_{n}}}-{\bf{Z_{A}}}$ and ${\bf{y_{n}}}={\bf{Y_{n}}}-{\bf{Z_{B}}}$. The monopole, dipole, and quadrupole moments are

[TABLE]

[TABLE]

and

[TABLE]

Gradients of the Green’s function for the gravitational potential, ${\bf{\nabla^{n}}}g\left(R\right)=-{\bf{\nabla^{n}}}\frac{1}{R}$, are

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

The approximated potential generated by the particles in cell B at a position ${\bf{x}}:=x_{i}{\bf{\hat{e}_{i}}}$ in cell A is then

[TABLE]

The Cartesian FMM described by equations (2) - (9) is the same as in D2000, except that : (1) we have opted to express the quadrupole moments in the extensive form , and (2) we have added terms that involve the dipole moment. These terms drop out in the case that cell coordinate centers coincide with cell centers of mass. As in D2000, we have dropped the octupole moment from equation (9). For a given interaction, this term is constant in space and hence does not contribute to the force calculation.

Using equation (9), the gravitational acceleration caused by the particles in cell B at a point ${\bf{X}}$ within cell A can be expressed as

[TABLE]

Similarly, the gravitational acceleration caused by the particles in cell A at a point ${\bf{Y}}$ within cell B can be expressed as

[TABLE]

Using equations (2) - (4), (10), and (11), we can express the force between two unit masses as

[TABLE]

Although the computed force is an approximation of the force on an individual particle, equation (12) shows that the sum of forces between any two particles is exactly zero. This implies the sum of linear momentum changes due to gravitation over all the masses in pairs of interacting cells is zero, and therefore the change over the entire computational domain is zero.

The same result does not generally hold for the sum of the torques generated between pairs of cells. Referring to the more general derivation of the method in §A, we write equation (A12) to expansion order $P=3$ and find the sum of all torques to be

[TABLE]

where we define the octupole moments,

[TABLE]

For many evolution methods employing the FMM, such as SPH or N-body, the net torque found in equation (13) (or the equivalent expression for a higher expansion order) is the sole source of angular momentum non-conservation. Eliminating this efficiency would, therefore, guarantee angular momentum conservation to machine precision.

We seek a correction to the Cartesian FMM of D2000 that (1) balances the net torque found in equation (13), (2) produces an equal and opposite force on each cell, and (3) is within the error bounds of the computed force. One possible solution satisfying these requirements uses the corrective force

[TABLE]

where $D_{ijkl}$ is the fourth derivative of the Green’s function. Proof that ${\bf{F_{c}}}$ cancels the torque imbalance is found in §A. The correction for cell A is

[TABLE]

and the correction for cell B is

[TABLE]

The corrective accelerations, ${\bf{g}}_{C,A\rightarrow B}$ and ${\bf{g}}_{C,B\rightarrow A}$, are added to ${\bf{g}}_{A\rightarrow B}$ and ${\bf{g}}_{B\rightarrow A}$, respectively, to obtain the total acceleration. This correction produces an extra acceleration that is constant over each cell. The torque produced is equal in magnitude but opposite in direction to the torque imbalance found in equation (13). The sum of the corrective forces on one cell is equal and opposite to that on the other cell, preserving the force balance of the original method. Because it uses a higher order Green’s function derivative, the corrective force is within the error bounds of the original method.

Another possible solution is to replace $D_{ijkl}$ in equation (15) with the non-symmetric tensor,

[TABLE]

yielding the alternative force correction,

[TABLE]

Note that ${D^{\prime}}_{ijkl}$ is simply ${D}_{ijkl}$ with any terms that do not contribute to equation (13) removed. As shown below in §3, equation (19) yields a faster implementation for a given opening angle while slightly increasing the solution error.

It is important to note that the quality of torque conservation in our modified FMM is not the same as the quality of force conservation. As shown by equation (12), the force between any two individual particles sums to zero. An analogous relation does not hold for the torques. The torques satisfy the less strict requirement that the sum of torques between all the masses in two interacting cells is zero. The correction also introduces unphysical torques between particles in the same cell, however, these corrections are within the error bounds of the original scheme.

3 Numerical Test

To test our new method, we have implemented a minimalistic version of the method of D2000, with options to use the corrections described by equations (15) or (19). This code is written in C++ for serial execution on a single processing core. The version of the code used in this paper is available through the Zenodo repository,

https://doi.org/10.5281/zenodo.571523. Note that the code is intended only to illustrate our method, and is not intended for production purposes.

For our test problem, we have chosen a binary star system for which the net torque imbalance can be relatively high. As can be seen in equation (13), the torque imbalance grows with the difference between octupole moments of the stellar components of a binary. Therefore, the larger and less centrally condensed one star is compared to its companion, the larger the net torque imbalance. One such system is a high mass ratio double white dwarf (DWD) with the larger, less massive star filling its Roche lobe. A system like this, if stable to mass transfer, is a potential progenitor of an AM Canum Venaticorum (AM CVn) type cataclysmic variable binary star (Marsh et al., 2004; Kilic et al., 2016). Our test problem is an approximation of such a system. The accretor has a mass of $1M_{\odot}$ and the donor a mass of $0.2M_{\odot}$, with the donor’s volume equal to the volume of its Roche lobe. Each component is taken to be a spherical Lane-Emden polytrope. In realistic systems, the donor will be tidally distorted, however, the spherical approximation is sufficient to demonstrate the usefulness of our method. The donor has a polytropic index of $\frac{3}{2}$ and the accretor has a polytropic index of $3$, approximating the cold white dwarf equation of state in the low and high mass limits, respectively. The test problem consists of $10^{6}$ equal mass particles, chosen by sampling the density distribution computed from integrating the Lane-Emden equation for each component.

Our test was executed on a single core of a 2.8 GHz E5-2680v2 Intel Xeon Processor on the QB2 cluster of the Louisiana Optical Network Initiative (LONI). The code was compiled using the GNU C++ compiler version 4.9.0. The gravitational solution was generated, using opening angle $\theta=0.2,0.3,...,1.0$, for the original uncorrected D2000 method, the torque corrected method (using equation (15)), and the torque corrected and optimized method (using equation (19)). We refer to these three methods, respectively, as the “UC”, “TC”, and “TCO” variants.

The net force and torque balances were computed using the formula

[TABLE]

where $A_{n,d}$ is the computed linear or angular acceleration, respectively, for the $n^{\mathrm{th}}$ particle and $d^{\mathrm{th}}$ dimension, and $\mu_{n}$ is the mass of the $n^{\mathrm{th}}$ particle. As seen in 1(a), the TC and TCO variants preserve torque balance to many orders of magnitude greater precision than the UC variant, consistent with machine precision for double precision floating point arithmetic. In all cases, net force balance is preserved to machine precision. In 1(c) we plot the compute time for each variant. The TC variant takes approximately $2$ times the compute time of the UC variant, while the TCO variant takes $1.7$ times as long. The significant increase in computational time is due to the added expense of computing $D_{ijkl}$ (or ${D^{\prime}}_{ijkl}$ )and ${\bf{F_{c}}}$. The computation of ${\bf{F_{c}}}$ adds $60$ terms ($20$ for each direction) to the $44$ terms present in equation (9) when dipole moments are removed. The computation of $D_{ijkl}$ adds $83$ terms to the $27$ terms present in all of the lower derivatives, while the optimized ${D^{\prime}}_{ijkl}$ adds only $21$ more terms. The relative force error is defined as the average of

[TABLE]

over all particles, where ${\bf{g_{PP,n}}}$ is the exact, directly computed force on the the $n^{\mathrm{th}}$ particle. The relative torque error is similarly defined using

[TABLE]

where here ${\bf{X}}$ is the distance to the coordinate origin.We plot these errors in 1(b) and 1(d). Both the force and torque errors are virtually identical between the UC and TC variants. The torque correction in the TC variant does not result in a force or torque error higher than in the original scheme of the UC variant. Both errors are higher for the TCO variant, therefore, for a given error, the TCO variant requires a smaller $\theta$ than the TC variant, resulting in more interactions to compute.

4 Discussion

The method of Dehnen (2000) conserves linear momentum to machine precision, but not angular momentum. We have presented two modifications to this method that each enable it to also conserve angular momentum to machine precision. This extra feature comes at computational expense, requiring approximately twice the compute time.

Whether or not the extra computational effort is worth the benefit of conserving angular momentum to machine precision will depend on the particular astrophysical system under investigation. One example of such a system would be a DWD at the onset of stable mass transfer. Past simulations of interacting DWDs have found angular momentum is artificially either added or removed from the system as the simulation progresses. Motl et al. (2002) found a normalized gain rate of $\sim 10^{-4}/\mathrm{orbit}$ for polytropic binaries of mass ratios $1.0$ and $0.8436$. Dan et al. (2011) found a $\sim 10^{-3}/\mathrm{orbit}$ normalized violation rate using SPH to simulate $84$ orbits of an interacting 0.8 $M_{\odot}$ accretor and 0.2 $M_{\odot}$ donor. These loss rates are sufficient that over many hundreds or more orbits, the violation of angular momentum conservation may cause systems that should be stable to become unstable (or vice versa) . One possible way to avoid this problem is to increase the resolution to the point the artificial angular momentum gain or loss rate is small compared to changes in the orbital and spin angular momenta of the system, however, it is difficult to determine what resolution is needed a priori. Increased resolution also comes at significant computational cost. The method described in this paper provides a remedy without increasing resolution.

We also note that conservation of neither linear nor angular momentum holds when the time-stepping is not uniform throughout the entire domain. In practice many SPH and N-body codes use individual time-steps for particles or groups of particles, resulting in a faster computation speed (Ahmad & Cohen, 1973). In order to fully realize the benefits of the method presented here, one has to abandon individual time-stepping and the speed-up that comes with it. Another benefit of individual time-stepping is that the non-conservation of momentum is often used as a proxy for the measure of the force error. With exact conservation of momentum and angular momentum, this is no longer possible, necessitating the choice of a different proxy. One possibility is to sum the magnitudes of the highest order expansion terms over the entire domain.

A higher order extension to this method is presented in the §A. The method is also applicable for any Green’s function that is solely a function of the scalar distance between points, such as softened gravitational potentials.

Acknowledgements

We wish to acknowledge the support from the National Science Foundation through CREATIV grant AST- 1240655.

Portions of this research were conducted with high performance computing resources provided by Louisiana State University (http://www.hpc.lsu.edu).

Portions of this research were conducted with high performance computational resources provided by the Louisiana Optical Network Initiative (http://www.loni.org).

The author would like to thank Geoffrey C. Clayton, Patrick M. Motl, and Joel E. Tohline for their help with this article. The author also thanks the referee for help with the revision process. The referee’s insights and concerns enabled publication of an article of much greater quality than the original submission.

Appendix A Appendix

The $n^{\mathrm{th}}$ derivative of a Green’s function dependent only on $R$, $D_{l_{1}l_{2}...l_{n}}$, can be written

[TABLE]

where the position vector, ${\bf{r}}:={\bf{e_{l}}}r_{l}$, is with respect to the coordinate center of the entire system. The coordinate origin of cell $A$ is located at ${\bf{X}}$ and the coordinate origin of cell $B$ is located at ${\bf{Y}}$. The distance between the cells is ${\bf{R}}:={\bf{Y}}-{\bf{X}}$. For the $i^{\mathrm{th}}$ particle in cell $A$, we define ${\bf{x_{i}}}\left({\bf{r}}\right):={\bf{r}}-{\bf{X}}$, and for the $j^{\mathrm{th}}$ particle in cell $B$, we define ${\bf{y_{j}}}\left({\bf{r}}\right):={\bf{r}}-{\bf{Y}}$. The potential on the $i^{\mathrm{th}}$ particle of cell $A$ caused by all particles in cell $B$, expanded to order $P$, is

[TABLE]

Equation (A2) is Equation 3 from Dehnen (2002) expressed in tensor notation. The potential on the $j^{\mathrm{th}}$ particle of cell $B$ caused by particles in cell $A$, expanded to order $P$, is

[TABLE]

Taking the negative of the derivative of $\Phi_{i,B\rightarrow A}$ with respect to $x_{k}$, we obtain the acceleration on the $i^{\mathrm{th}}$ particle of cell $A$,

[TABLE]

Similarly, for the $j^{\mathrm{th}}$ particle of $B$,

[TABLE]

As expected, the sum of the forces over all particles is zero,

[TABLE]

The total torque, ${\bf{\tau_{AB}}}$, about the coordinate origin of cell $A$ is,

[TABLE]

The total torque can be thought of as the sum of a bulk torque,

[TABLE]

and the spin torques of each cell,

[TABLE]

Using fact that

[TABLE]

we can express the spin torques as

[TABLE]

We see that the RHSs of equation (A8) and equation (A11) differ only in sign and the range of summation indices. The spin torques of expansion order $m$ are canceled by the bulk torque of expansion order $m+1$. The spin torques that result from the highest expansion order do not have a bulk torque to cancel them, resulting in a net torque. Using equation (A7), equation (A8), and equation (A11), we can express the net torque

[TABLE]

If we apply a constant corrective force, ${\bf{F_{c}}}$, to the particles in cell $A$, and $-{\bf{F_{c}}}$ to the particles in cell $B$, the balance of force remains unaltered. The contribution to the torque is

[TABLE]

When the coordinate centers for cells $A$ and $B$ are coincident with the centers of mass for the respective cells, dipole moments vanish and the sum of corrective torques for the last two terms on the RHS of equation (A13) vanish. Comparing equation (A12) with the first term on the RHS of equation (A13), we find that if we set

[TABLE]

the sum of the original FMM torque and the corrective torque vanishes,

[TABLE]

Summing over all masses in each cell, the total corrective force, ${\bf{F_{C}}}$, is

[TABLE]

Here we have defined the generalized moments for each cell,

[TABLE]

and

[TABLE]

Making the definition,

[TABLE]

we can define an alternative corrective force,

[TABLE]

This corrective force also results in a balanced torque. Depending on the choice of Green’s function, equation (A20) may result in fewer terms to compute.