Numerical+Methods

This page summarizes the no-slides talk of R. Spurzem on Wed Feb. 17 14:00 on methods (N-Body, Monte Carlo, Fokker-Planck). The hardware part will follow under a separate heading (GRAPE, GPU, FPGA) . Comparing Monte Carlo, Gaseous and N-Body Model (M. Giersz, R. Spurzem, S. Aarseth at the IAU Symposium 246 in Capri, 2007...)

There is the duality between direct N-Body Methods and Statistical Methods based on the Fokker-Planck (FP) equation. The latter comprises direct numerical or Monte Carlo solutions of the orbit-averaged FP equation, as well as so-called gaseous or momentum models using moment equations of the local FP equation, similar to hydrodynamical equations. A comparison of simple physical runs (single mass core collapse and post-collapse) for all three methods can be found e.g. in [|Giersz & Spurzem (1994)]. Similar successful comparisons exist in the literature for multi-mass and rotating axisymmetric systems, systems with central massive black holes. The N-Body method, solving just N Newtonian equations of motion (second order ordinary diff. eqs.) is often deemed most accurate in principle, because it involves the least physical approximations. However due to the inherent chaotic nature of the N-Body problem it is not really guaranteed that an N-Body solution follows a 'true' solution in phase space, even if some global quantities (energy, angular momentum) are very well conserved. Only by performing many simulations and interpreting averaged results, and by comparison with statistical methods we can gain confidende in results.
 * ==__Introduction__==

For Fokker-Planck methods we start with the local FP equation for the phase space distribution function f = f(r,v,t) ; the dependence on space coordinate x is dropped here for clarity, and we assume the system to be stationary on orbital timescales, so it does not explicitly depend on t. To arrive at this equation from 6N-dim. Liouville equation we have assumed that particles are indistinguishable, i.e. there are no correlations between their positions or velocities (no binaries):

Here we have used the Rosenbluth potentials h,g, introduced by [|Rosenbluth, McDonald & Judd 1957] for plasma dynamics: Eq. (A)

Note that for any distribution function f, which is given as, where A0 is a Maxwellian and Pl are Legendre polynomials, the Rosenbluth potentials are given analytically to any order desired, as was already noted by Rosenbluth, McDonald & Judd (1956). We have used this for multi-mass anisotropic systems [|(Spurzem & Takahashi 1995)] (see also below on gaseous models); hint: the trick is simply to expand 1/|v-vf| in terms of spherical harmonics, very much like we have seen it for the multipole series of 1/|r-r'| by Rosemary Mardling for the 3-body problem.

Note that the real diffusion coefficients, which are of the dimension velocity change per time unit (velocity change squared per time unit) are related to the Rosenbluth potentials by

Eq. (B)

In this way the potentials h,g just describe the outcome of individual gravitative two-body encounter between two particles with velocities v and vf (well an average over all impact parameters p has been skipped here, as you see from the presence of the Coulomb logarithm (ln Lambda) above.Note that in this form the FP equation is a 3D integro-differential equation for f(x,v), which one could try to solve for every x with respect to v, remember that at this point f and also the D's still carry a dependence on position x, which is just dropped here in our notation. E.g. if we use local Maxwellian for f = A0, the normalisation factor contains the density and velocity dispersions at x..

The Monte Carlo Method in principle solves Eq. (A) by a Monte Carlo integral. Practically one directly determines a diffusion coefficient as in Eq. (B). So if we determine a representative cumulative encounter for one relaxation time in the MC method (see talk by J. Fregeau in this program) we have actually estimated the integrals in Eq. (A) by one value. We may only sum up the effect of many small angle encounters into one representative encounter over a relaxation time trx (or a MC step which should not too much differ from trx to avoid spurious errors and over-relaxation), if individual encounters are uncorrelated. Then, due to the central limit theorem, the effect of all encounters is a random variable with Gaussian distribution around the accumulated mean.
 * ==__Monte Carlo Method (MC)__==

But, this was only the integration over velocity space. To obtain a proper diffusion coefficient for a stellar orbit we should have to apply orbital changes at any location of an orbit with say given E and J, and integrate these over accessible coordinate to obtain a proper orbit-averaged diffusion coefficient. D(Delta E) is obtained from Eq. (B) above, after a variable transformation from the vi to E,J, and an angle. The "miracle" of the MC method is that it works by applying one scattering per stellar orbit at a given time, and the orbit average is realized through the ergodic principle - relaxing a star (whose position is picked randomly on its orbit) during many time steps means replacing the spatial average by a time average. That this works well has not been strictly proven as far as I know, and it obviously depends on the fact that the system evolves slowly.

Presently active MC methods all use the Hénon method [|(Hénon 1971)], which use this principle to realize the orbit average in time, i.e. applies the change for a star on its orbit, picking its position afterwards anywhere on the orbit for the next step. In contrast to this the MC method of Spitzer has applied MC changes to a particle locally on its orbit [|(Spitzer & Hart 1971)]. Presently, as to the knowledge of the author, only Hénon's method is in active use by stellar dynamics people (see again talk slides by John Fregeau).


 * ==__Orbit Averaged Fokker-Planck equation__==

Orbit Averaging proceeds in two steps: first apply Jeans' theorem, that f(x,v) is constant along a particle's trajectory through the system, so we can transform f into f(I1,I2,....) where the I's are the constant's of motion of a single particle in the gravitational potential Phi. Phi must be obtained as in the other methods from Poisson's equation using the density rho (integral of f over velocity space).

Generally we have in spherical symmetry f(E,J), i.e. f depends only on energy and the modulus of angular momentum. Often the dependence on J is neglected and only the isotropized (J-averaged) distribution function f(E) considered. In axisymmetric systems we have f(E,Jz,I3), i.e. f depends on energy, z-comp. of ang. momentum and a generally unknown third integral (if it exists). To show an unusual example we present here the orbit-averaged FP equation for an axisymmetric system, where the dependence on the 3rd integral is neglected (averaged away):

from [|Einsel & Spurzem (1999)]. This original pre-collapse single mass model has been extended to multi-mass, post-collapse and is used till the present day for more realistic models of globular and other star clusters, including the important degree of freedom of rotation (also with central black hole). Diffusion coefficients here cannot easily be computed from Rosenbluth potentials, rather they must be numerically computed by integrals over the accessible coordinate space for any given E and Jz, under the present potential Phi. The above equation can be directly numerically solved on a two-dimensional mesh of E and Jz values, a 2D FP equation. Note, that since the five diffusion coefficients above are known in principle they could be realized by a random choice of Delta E and Delta Jz for a new 2D MC method. That has not yet been done.


 * ==__Gaseous or Momentum Models__==

Building the hierarchy of moment equations in spherical symmetry from the FP equation yields

One realises well known continuity and Euler (linear momentum) equations in the first two lines. Note that in this case the so-called left-hand-side of the FP equation (not shown above, also known as Boltzmann equation) has been used, and the complicated FP equation above only shows up in the two "enc" terms on the right hand sides here. Therefore this model is the only one (except N-Body) which can treat dynamical phases of evolution. Note that this particular form is the anisotropic gaseous model by [|Louis & Spurzem (1991)] - it is closed by a heat flux equation with conductivity according to [|Lynden-Bell & Eggleton (1981)]. Closure equations and self-consistent encounter terms with analytical expressions have been generalized to multi-mass systems by [|Spurzem & Takahashi (1995)] and to systems with central massive black holes by [|Amaro-Seoane, Freitag & Spurzem (2004).]

Due to the long relaxation time we have to distinguish in the above equations radial and tangential pressure (temperatures or velocity dispersions), where the total pressure is 3p = pr + 2pt, and the energy flux 2F = Fr + Ft. The sum of the two last equations yields a standard energy equation, its difference an evolution equation for the anisotropy pressure pa = pr - pt. "bin3" denotes binary heating terms, which are not discussed here.