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
where 
.
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
[1], 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.
