
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
http://sipapu.astro.illinois.edu/~ricker/.
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.

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 singlevalued
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.
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 PiecewiseParabolic 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 secondorder 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 gridbased 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 Nbody 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 wellsuited 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 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 Nbody codes.
The time
integration of particle positions in an Nbody code
is typically done using some sort of
leapfrog technique. The computation of particle forces, since it usually
includes a longrange force such as electromagnetism or gravity, is the
most timeintensive 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
particleparticle (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
coprocessor hardware for computer workstations
(Ebisuzaki et al 1993). A more efficient
method is the particlemesh (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 particleparticleparticlemesh (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.
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 gridbased method
such as PPM is used for the hydrodynamics, the best choice for the
collisionless solver is probably particlemesh, since the hydro grid
already limits the spatial resolution.
For particlemesh 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 code joins together PPMnD, a PPM/multigridbased hydrodynamic
solver that I wrote as part of my Ph.D. thesis work, and a
particlemesh 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.
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 Nbody 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. MultiGrid 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)

