Paul Ricker Home Research Teaching Papers Interests Links
Background Fluids Particles Gravity Codes References

Old research pages

I'm keeping these old research pages here for the benefit of anyone who might arrive here from an external page. They are no longer updated. Please visit my new web site at Eventually some content from these pages will migrate to the new site.

Astrophysical Hydrodynamics

The behavior of matter under astrophysical conditions

The behavior of matter in space is very heavily influenced by one parameter: the ratio of the mean free path, or the average distance a particle travels between collisions, to the length scale characterizing a system. In addition to determining the dynamical properties of a system, this parameter determines what numerical methods are most efficient for studying the system.

Collisional vs. Collisionless
What happens when two opposing streams of matter strike one another, if the matter is (a) collisional or (b) collisionless.

If the mean free path is much smaller than length scales of interest, particles collide frequently, randomizing their initial velocities. The velocity field is characterized by a randomly fluctuating `thermal' component superimposed on a `bulk' or average component which is a single-valued function of position. Matter streams cannot interpenetrate, but instead stop when they collide, producing sound or shock waves if the matter is compressible. This type of matter is called `collisional' or `fluid' matter; the term `hydrodynamics' refers to the dynamics of collisional matter. Examples of collisional fluids include the gas in Earth's atmosphere or in the Sun, the material in accretion disks around black holes, and the gas in clusters of galaxies. Generally collisional flows in astrophysics are inviscid, turbulent, and supersonic (so that interacting matter streams produce shock waves). Often astrophysical flows are not truly collisional, but because of interactions between the particles and magnetic fields they behave as though they were.

If the mean free path is much larger than length scales of interest, particles tend to stream from one place to another, turning gradually in their mutual gravitational or electromagnetic field but not undergoing sudden changes of direction which would randomize the velocity distribution. Matter streams can interpenetrate, causing the bulk velocity to be a multivalued function of position. This type of matter is called `collisionless' matter. Examples include the `gas' of stars in a galaxy, the neutrinos escaping from the solar interior, and the `gas' of galaxies in a cluster. The dark matter needed to bind galaxies and galaxy clusters, which by hypothesis interacts only gravitationally with baryonic matter, is collisionless if the simplest hypothesis is true. Collisionless matter can have (and does develop, if given enough time) an isotropic velocity distribution like collisional matter, but when disturbed, its response is very different.

As mentioned, the numerical method used to study a system depends on whether the system is collisional or collisionless. The starting point for all methods is the Boltzmann equation, which describes the evolution in space and time of the particle distribution function f(x, p), which gives the probability of finding a particle with momentum p at location x. Collisional systems are best studied by taking velocity moments of the Boltzmann equation. This results in a hierarchy of conservation equations for the fluid density, bulk momentum, and thermal energy, which we close by assuming an equation of state (relationship between pressure and thermal energy). We then solve these equations by using techniques developed for other partial differential equations (by sampling solution values at points on a grid or at suitably chosen test particle positions). Collisionless systems are best studied by following the motion of test particles drawn from the particle distribution function. Cases where the mean free path is roughly comparable to the size of a system are the most difficult: these require us to solve the Boltzmann equation directly, which means solving a partial differential equation in six independent variables and assuming explicitly what happens at the `microscopic' level when two particles collide.

Collisional methods: PPM vs. SPH

Even when gas in an astrophysical setting is collisional, very frequently it is highly compressible, and its temperature is often such that typical interactions (say between clouds) occur supersonically. Under such conditions shock waves form. A shock wave is a discontinuity in the density, pressure, and bulk velocity of a fluid where bulk kinetic energy is suddenly transformed into thermal energy. Hence the density and pressure go up as a shock wave passes a fixed point. For an observer riding with the shock, material flows into the shock at a supersonic speed and exits the shock at a subsonic speed. Thus the preshock gas has no information about the coming shock, and its deceleration takes place within a few mean free path lengths of the shock position. As a result, shocks typically have a very fine structure which is impossible to resolve spatially in a simulation (if the gas is truly collisional).

Two methods are predominantly used in astrophysical hydro to study the evolution of collisional flows with shocks. The Piecewise-Parabolic Method (PPM; Colella & Woodward 1984) solves for the average values of the density, pressure, and velocity in cells of a grid. The flow of material from one cell to another is determined using the known analytical solution to a simple shock problem called the Riemann problem. In smooth flow the technique reduces to a simple second-order finite difference, but when shocks are present, they are typically represented by discontinuities no more than one or two zones across. This property makes PPM ideal for problems in which shocks are important, although like all grid-based methods it has difficulty covering a wide range of length scales without some form of adaptive gridding.

The second technique, Smoothed Particle Hydrodynamics (SPH; Gingold & Monaghan 1977, Lucy 1977), evolved from the N-body methods used for collisionless particle dynamics (see below). In SPH each particle is assigned a finite size which depends on the local particle density; in denser regions, the size is smaller. Each particle is also assigned a pressure corresponding to the value of the actual pressure field integrated against a smoothing kernel which is centered at the particle's location and has a `width' corresponding to the size of the particle. The pressure gradient, estimated using the pressures of nearby particles, is then used to determine a pressure force on the particle. SPH does not handle shocks particularly well, but because it can concentrate particles where they are needed most, it is well-suited to problems involving a wide range of length scales. The tradeoff here comes in mass resolution: the particle mass determines the least massive objects which can be resolved, as well as the particle size in very underdense regions, where the temperature can be poorly determined (since it depends on a small number of particles).

Collisionless methods: Particle-Mesh vs. Tree

Collisionless matter is simulated by sampling the particle distribution function with a small number (compared with the actual number of molecules, stars, or galaxies, which can be quite large) of `superparticles' which obey the same equations of motion as the real particles. Therefore collisionless solvers are also referred to as N-body codes. The time integration of particle positions in an N-body code is typically done using some sort of leapfrog technique. The computation of particle forces, since it usually includes a long-range force such as electromagnetism or gravity, is the most time-intensive part of a collisionless solver. Hence it is also the force computation which primarily distinguishes among methods.

The simplest approach to computing forces, direct summation or particle-particle (PP), is also by far the least efficient, although advances have been made, for example, by using separate timesteps for each particle (Aarseth 1972) or by implementing the summation in dedicated co-processor hardware for computer workstations (Ebisuzaki et al 1993). A more efficient method is the particle-mesh (PM) method (Hockney & Eastwood 1988). In this method, one uses an interpolating kernel to assign densities to a grid given a set of particle positions. One then solves for the gravitational potential on this grid using a technique such as the Fast Fourier Transform (FFT) or multigrid. The same interpolating kernel is then used to determine the force on each particle. What this gains in speed it loses to the spatial resolution of the grid. The particle-particle-particle-mesh (PPPM) method (Efstathiou & Eastwood 1981) gets around this by correcting the PM force with a direct summation for particles which lie within the same grid zone. However, for highly clustered matter distributions, the inefficient PP calculation begins to dominate the required computer time. Finally, tree methods (Hernquist 1987) dispense with grids altogether by approximating the force due to distant concentrations of matter using their monopole moments, reserving direct summations for nearby particles.

Linking hydrodynamics and particle dynamics: Poisson solvers

Frequently we wish to study the gravitationally coupled evolution of collisional and collisionless matter -- for example, stars and gas in a galaxy, or galaxies, dark matter, and gas in a cluster of galaxies. In these cases we need to solve for the combined gravitational field of both matter components and use it as a source term for each method. If SPH is used for the hydrodynamics, either PM, PPPM, or the tree method can be used to solve for the gravitational force. If a grid-based method such as PPM is used for the hydrodynamics, the best choice for the collisionless solver is probably particle-mesh, since the hydro grid already limits the spatial resolution.

For particle-mesh codes, two general types of potential solver are used: FFT and multigrid. The FFT method (Hockney & Eastwood 1988) convolves the appropriate Green's function with the Fourier transform of the density field, then uses an inverse transform to compute the potential. Multigrid methods (Brandt 1977, Hackbusch 1985) accelerate standard relaxation techniques by solving the Poisson equation on grids of increasing coarseness, then using these solutions to correct the solution on the simulation grid. Both methods achieve roughly the same execution time scaling with problem size, and both can be used with periodic or isolated boundary conditions (the two most frequently encountered in astrophysical problems). However, multigrid is more straightforward to apply to nonuniform and adaptive meshes.

The COSMOS N-body/hydrodynamics code

The COSMOS code joins together PPMnD, a PPM/multigrid-based hydrodynamic solver that I wrote as part of my Ph.D. thesis work, and a particle-mesh code written by Scott Dodelson of the Fermilab Theoretical Astrophysics Group. COSMOS is designed to handle fairly general 3D astrophysical problems with Cartesian symmetry, as well as 1D spherical problems without a collisionless component. Visit the COSMOS home page for information about the structure and usage of the code and test results.

The FLASH adaptive-mesh hydrodynamics code

FLASH is a more modular simulation code designed to handle general astrophysical problems, with a tilt toward problems involving nuclear burning. I am helping to develop it as part of the FLASH Code Group in the ASCI FLASH Center at the University of Chicago. FLASH includes adaptive mesh refinement and nuclear burning as well as compressible hydrodynamics. It shares a history with COSMOS in that FLASH version 1.x uses a more advanced version of the code framework that Scott and I developed for COSMOS.


Aarseth, S. J. in Gravitational N-body Problem, IAU Colloquium No. 10, ed. M. Lecar (Reidel), 1972

Brandt, A. Math. Comp. 31 333 (1977)

Colella, P. and Woodward, P. J. Comp. Phys. 54 174 (1984)

Ebisuzaki, T., et al. PASJ 45 269 (1993)

Efstathiou, G and Eastwood, J. W. MNRAS 194 503 (1981)

Gingold, R. A. and Monaghan, J. J. MNRAS 181 375 (1977)

Hackbusch, W. Multi-Grid Methods and Applications (Springer), 1985

Hernquist, L. Ap. J. Supp. 64 715 (1987)

Hockney, R. W. and Eastwood, J. W. Computer Simulation using Particles (IOP), 1988

Lucy, L. B. Astr. J. 82 1013 (1977)

[Home] [Research] [Teaching] [Papers] [Interests] [Links] disclaimer 09/01/16