A full treatment requires the dynamical coupling of the multiple time scales and multiple length scales of the external disturbance and the galaxian components discussed above. Relevant characteristic length and time scales may differ by an order of magnitude between satellite and halo or disk orbits. In addition, we will see that the halo disturbance may be relatively weak and a significant perturbation of the outer disk at the same time. These multiple-scale weak regimes are a challenging task for an n-body computation. However, this class of problems is ideally suited to linear techniques and the work here will use the expansion technique known as the matrix method. Although the matrix method is computationally intensive, it is no more so than n-body methods and is practical on current workstations. In this section, I will give a brief overview of the general method with details on posing and implementing the coupled response solutions in the references cited below and in the Appendix.
In short, the matrix method represents the response of a galaxy to an external perturbation by a truncated series of orthogonal functions, similar to those one would use to solve an electrostatics problem. The perturbation is also represented by this series and the temporal dependence of each coefficient is Fourier transformed to a (complex) frequency distribution. The response of the galaxy to one of the orthogonal functions at a particular forcing frequency is then computed in the continuum limit using the collisionless Boltzmann (Vlasov) equation. The entire procedure is analogous to signal processing in Fourier space. Pursuing the analogy, we now do the inverse transform. The response to any perturbation, the weighted superposition of the response to each basis function, is then a matrix equation. Finally, to get the full time dependence of the response, one resums the solutions to the matrix equation at each frequency weighted by the Fourier coefficients.
This method assumes that the perturbation is small enough that the overall change to the structure of the galaxy is small. In this limit, the method has the advantage of accuracy and sensitivity to the large scale structures of interest. For contrast, the n-body simulation determines the response of a galaxy to a perturbation by solving the equations of motion for a representative set of orbits. The orbit is advanced in a fixed potential for a short time interval and the gravitational potential or force is then recomputed. The simulation works well for large perturbations but because the simulation uses a finite number of particles, fluctuation noise limits the sensitivity to small amplitude distortions. The matrix method nicely complements the n-body simulations, excelling in the regimes where the n-body simulation are suspect.
Historically, the approach is related to the treatment of general eigenvalue problems described in the mathematical physics literature (e.g. Courant & Hilbert 1953, Chap V). The matrix method in stellar dynamics had varied applications beginning with Kalnajs (1977) who investigated the unstable modes of stellar disks. Polyachenko & Shukhman (1981) adapted the method to study a spherical system (see also Fridman & Polyachenko 1984, Appendix) and it was later employed by both Palmer & Papaloizou (1987) in the study of the radial orbit instability and by Bertin & Pegoraro (1989) to study the instability of a family of models proposed by Bertin & Stiavelli (1984). In addition to Paper I, Weinberg (1989, Paper II) used the matrix formulation to study the response of a spherical galaxy to an encounter with a dwarf companion and Saha (1991) and Weinberg (1991) investigated the stability of anisotropic galaxian models.