In order to estimate the evolution on large scales over long time periods, we need to suppress small scale noise as much as possible and this need is satisfied by the biorthogonal expansion technique (e.g. Clutton-Brock 1972, 1973, Kalnajs 1976, Fridman & Polyachenko 1984, Hernquist & Ostriker 1992). The approach uses the eigenfunctions of the Laplacian to construct a complete set of orthogonal functions that satisfies the Poisson equation. The projection of the ensemble particle positions on each member of the orthogonal series yields a set of coefficients, similar to determining a potential from a charge distribution in electrostatics. These coefficients may then be used to describe the density, gravitational potential or force in the expansion code. The most efficient implementations use analytically derived recursion relations to generate the functions. The expansion converges quickly and is most efficient when the lowest order function is a good fit to the underlying equilibrium. However, most astronomical distributions do not match the available sets of special functions. The problem described in this paper motivated the algorithm developed in Weinberg (1999) which allows the adaptive construction of both spherical and three-dimensional cylindrical bases and this algorithm is used here.
Another advantage of this force solver is its ability to separately follow distinct kinematic components. Each component may be tied to a basis tailored to its geometry; this helps remove the bottleneck in simultaneously resolving multiple spatial scales. In particular, we can assign the halo particles to a spherical basis and the disk particles to a cylindrical basis. These are gravitationally coupled through the force evaluation. This procedure is easily extended. For example, a simulation which follows the response of the Milky Way halo simultaneously can be implemented by specifying an appropriate basis for the halo and following the halo response and its back-reaction on the LMC directly. This would increase the run time of simulations described here by only about 30%. These simulations will be performed in the next phase of this project.
Experimentation reveals that multiple expansion centers may introduce numerical feedback and excite oscillations. In principle, each component may be uncoupled as long as feedback is suppressed. For simulations here, the halo expansion center is tied to the disk center. The force solver can easily resolve small offsets correctly so this is not a limitation for this application.
A parallel code with these features is implemented on a Linux-based cluster using MPI. Each node in the cluster has two processors; the algorithm was multi-threaded on each node to reduce the memory overhead. This allows the computation to be in core for all harmonic orders used here. The code performs load balancing but this is usually not required for the dedicated cluster. In practice, the average CPU load efficiency is approximately 90% for N=105 on 32 processors and improves for larger numbers of particles.