We invented the pseudo-particle multipole method (P2M2), a new method to handle multipole expansion in high accuracy Barnes-Hut treecode. This method expresses multipole expansion using a small number of pseudo-particles. With this method, the implementation of high-order terms of multipole expansion is greatly simplified. Furthermore, we can use GRAPE to evaluate high-order expansion terms, since the forces due to high-order terms are expressed as forces from particles. Previously we could not handle these terms on GRAPE and thus the calculation speed was rather slow.

The basic idea of P2M2 is to approximate potential field of particles using pseudo-particles. We distribute pseudo-particles so that they have the same multipole expansion as the original distribution up to a given order. Pseudo-particles have the degree of freedom both in their positions and masses. However, if we fully use their degree of freedom, we need to solve non-linear equations in order to determine the distribution of pseudo-particles. It is difficult to solve such equations, and even worse, it is not clear whether or not an acceptable solution exists for arbitrary distribution of physical-particles. In order to avoid this difficulty, we fix the positions of pseudo-particles and change only their masses. When positions are given, we can calculate the masses of pseudo-particles for an arbitrary distribution by solving linear equations, although the necessary number of pseudo-particles increases.

In practice, we put pseudo-particles on a ring in two dimensional calculation, and on the surface of a sphere in three dimensional calculation. We use these position sets to simplify the calculation of pseudo-particle masses. In the case of two dimensions, when we use a position set equally spaced on a ring, the multipole expansion coefficients are expressed by the Fourier transform of the masses of pseudo-particles. Therefore we can express the masses by the inverse Fourier transform of the expansion coefficients. In the case of three dimensions, the masses of pseudo-particles are determined in the same way as in two dimensions, except for two differences. One difference is that we choose the positions on the surface of a sphere so that they form a "spherical design", which is a set of positions "equally spaced" in some sense. Mathematically, it is defined as a set of points on which the numerical integration of an arbitrary polynomial of a given degree is exact. The other difference is that we use spherical harmonics as base functions instead of trigonometric functions. For both two and three dimensions, we derived formulae to calculate the masses of pseudo-particles from the distribution of physical-particles.

We applied P2M2 to treecode and combined it with GRAPE-4. Numerical tests have shown that P2M2 achieves the same order of accuracy with the original multipole expansion, and is faster than standard Barnes-Hut treecode with monopole when high accuracy is required. Timing results on GRAPE-4 shows that one processor board of GRAPE-4 accelerates the calculation by a factor of about 10 for the calculation of low accuracy, and the factor increases more than 200 as accuracy becomes higher.