In Anderson's method, the multipole expansion of the potential due to a clump of particles is effectively expressed in terms of the values of potentials on a sphere surrounding the particles. The potential outside the sphere is given by the surface integral on that sphere, which is then approximated by the sum over sampling points. This method, though elegant, appears to be rather indirect.
An alternative approach would be to use multiple particles to represent the multipole expansion. The basic idea here is to place small number of pseudoparticles which reproduce the multipole expansion of the original physical particles. In the following, I first present the theory in two dimensions, and then that in three dimensions.
In two dimensions, the multipole expansion of the gravitational field due to one particle is given by
where and z are position of the particle and position at which to evaluate the potential in the complex plane, and m is the mass of the particle. Here we use the system of units where the gravitational constant G is unity. This formula converges if .
If we have N particles with mass at locations (), the potential field outside the circle of radius a is expressed as
where M is the total mass of particles and is defined as
Our goal is to find an efficient way to place K pseudoparticles to approximate the potential field . In theory, the total number of freedoms we can attain with K particles is 3K. Therefore, with arbitrary assignment of mass and position, K particles should be able to represent multipole expansions of up to , where denotes the maximum integer which does not exceed x. However, in order to determine such an arbitrary distribution of pseudoparticles, we have to solve the system of nonlinear equations with 2p+1 variables, and the calculation cost would be at least . In addition, it is not clear whether or not an acceptable solution exists. In the following, we describe a more systematic approach in which we do not have to solve nonlinear equations.
In Anderson's method, potential is calculated at equispaced points on a circle. In the same splits, here we place particles in a ring, and only allow the masses of particles to change. This implies that we use only K out of 3K degrees of freedoms, and we need 2p+1 particles to express the multipole expansion coefficients. However, with this choice we can determine the mass of these 2p+1 particles with or calculation cost.
Consider the mass distribution of on a ring of radius r. The multipole expansion coefficients is given by
where is the line density of mass at polar coordinate . Thus, from the expansion coefficients , can be calculated by evaluating the Fourier series
When we approximate this continuous m by 2p+1 discrete points at , is given by
Because of the nature of the Fourier series, these express exact values of multipole expansions up to p-th order. The potential outside this circle can be calculated as the sum of the potentials by these particles as
In the case of the tree algorithm, we can use Eq. (11) to calculate the gravitational interaction between a node and a particle.
In the M2M part, which is the same for both the tree algorithm and FMM, we still have to construct the multipole expansion or particle representation around the center of a node from those of the child nodes. Since the child nodes are already represented by particles, we can use Eq. (7) to obtain the expansion coefficients of the parent node.
Alternatively, we can eliminate the use of multipole expansion coefficient by calculating the mass of pseudoparticles directly from mass of physical particles (or the mass of the pseudoparticles in child nodes). We can derive the formula to calculate directly from by combining Eqs. (7) and (10)
For tree algorithm, this is the end of the story. In the case of FMM, we still have to specify the algorithms for M2L part and L2L part. We can use either Anderson's method or standard harmonic expansion. For local expansion, Anderson's method is simpler to implement than the spherical harmonics.
The formulation for three dimensions is essentially the same as that for two dimensions, except that we need to use spherical harmonics instead of . The expansion coefficients is expressed as
where is the mass of particle i and is its polar coordinate. The function is the spherical harmonics of degree l, which is expressed as
Here, is the associated Legendre function of degree l and order m. Using these , the potential at position is given by
Our goal here is to obtain a mass distribution on a sphere of radius a which satisfy
where S denotes the surface of the sphere. Because the spherical harmonics comprise an orthonormal system, this is expressed as
If we use K points on a sphere, their masses are calculated by
The series expansion must be truncated at a finite value as we've seen in the case of two dimensions. As discussed by Anderson , the cutoff order must be , if K points form a spherical t-design.
As in the case of two dimensions, we can directly translate the positions and masses of physical particles to that of pseudoparticles. The formal expression is a triple summation over i, l, and m. However, we can simplify it using the addition theorem of spherical harmonics. The addition theorem is
where is the angle between two vectors with directions and . Using this addition theorem, is expressed as
where is the angle between the direction of physical particle i and pseudoparticle j.
The potential due to physical particles is approximated by the sum of potentials due to these pseudoparticles.