# The geometry of phase mixing

## 1 Introduction

The evolution of the phase-space density of stars in galaxies is determined by the collisionless Boltzmann equation, which states that phase-space flow is incompressible. Thus a cloud of stars in phase space becomes more and more distorted as it evolves; the local or fine-grained density around any point in phase space remains the same, but the coarse-grained density evolves towards a stationary state. This process, known as phase mixing, is the principal mechanism by which stellar systems reach coarse-grained equilibrium. The standard analogy is stirring a glass containing 20% rum and 80% Coke; eventually every finite volume element in the glass contains 20% rum and 80% Coke, even though infinitesimal volume elements are either 100% rum or 100% Coke (e.g. Arnold and Avez (1968)).

Galaxies are not thoroughly phase-mixed, since they are at most a few dynamical times old in their outer parts, and since mergers and star formation continuously add new stars. There is a variety of direct observational evidence for incomplete phase mixing: (i) Sharp-edged features (“shells”) in the outer parts of at least 10% of elliptical galaxies are believed to arise from the recent tidal disruption of small galaxies (e.g. Hernquist and Quinn (1988)). (ii) Proper-motion and radial-velocity surveys of the local halo reveal clumping in phase space (Majewski et al. (1996)). (iii) The tidally disrupted Sagittarius dwarf galaxy (Ibata et al. (1997)) provides a concrete example of a cloud in the early stages of phase mixing. (iv) Moving groups in the solar neighborhood may result from dissolution and mixing of star clusters and associations (Eggen (1965); Dehnen (1998)).

Incomplete phase mixing is likely to play a growing role in the interpretation of observations, for several reasons: radial-velocity and proper-motion surveys are rapidly improving in size and quality; large telescopes with improved image quality will permit us to examine the surface-brightness structure of elliptical galaxies with high spatial resolution and signal-to-noise ratio; the Hipparcos mission has dramatically improved the precision of measurements of the phase-space distribution of stars in the solar neighborhood, and future space-based astrometric missions such as SIM and GAIA will provide proper-motion and parallax measurements over much of the Galaxy; and the Sloan Digital Sky Survey will soon provide a stellar database of unprecedented size and uniformity ($`10^8`$ stars with five-color photometry).

Virtually all theoretical descriptions of phase mixing have focused on mixing in a two-dimensional phase space, which is the simplest to visualize (Fig. 1). This paper discusses some of the geometrical features of partially phase-mixed systems (we shall call these “phase structures”) with higher dimensionality, and the projection of this geometry into the observables that are measured in a survey. We shall argue that this geometry can be organized by two integers: the dimension of the phase structure and the dimension of the survey.

## 2 The dimension of phase structures

A cloud of stars orbiting in an integrable potential (the “host galaxy”) can be viewed as a surface or manifold of dimension $`DN`$ embedded in $`N`$-dimensional phase space (usually $`N=6`$). For example, a cold disk has $`D=2`$ (the disk occupies a two-dimensional surface in configuration space, and the velocity at each position is unique), collisionless cold dark matter has $`D=3`$ (three spatial dimensions but zero random velocity), and a hot galaxy has $`D=6`$. Since real stellar systems always have non-zero thickness and velocity dispersion, we are referring to an “effective” or “coarse-grained” dimension; loosely speaking, a structure in a phase space of $`N`$ dimensions has effective dimension $`D`$ if its extent in $`D`$ independent directions is much larger than its extent in the other $`ND`$ directions.

Because the flow of the cloud through phase space is incompressible, the dimension of the manifold does not change in the initial stages of phase mixing, although the shape of the manifold becomes more and more complicated. However, the effective dimension can change over longer times, through several distinct processes. The most obvious is that the phase-mixing scale grows smaller and smaller with time (see Figure 1), and eventually becomes smaller than the coarse-graining scale. At this point the phase structure has the same (coarse-grained) dimension as the phase space, unless all of the stars in the cloud have the same value of one or more isolating integrals of the motion. Gravitational scattering by small-scale irregularities (e.g. massive objects) can also lead to diffusion of the phase structure and an increase in effective dimension.

Less obviously, phase mixing can also decrease the number of effective dimensions. For example, consider the cloud resulting from the disruption of a small hot galaxy. The motion of the stars in the cloud can be described by a Hamiltonian $`H(𝐉)`$, where $`(𝐰,𝐉)`$ are action-angle variables in the host potential. The equations of motion are

$$𝐉=\text{const},\frac{d𝐰}{dt}=\frac{H(𝐉)}{𝐉}𝛀(𝐉).$$
(1)   
The trajectory of the centroid of the cloud is $`(𝐉_c(t),𝐰_c(t))=(𝐉_0,𝐰_0+𝛀_0t)`$, where $`𝛀_0𝛀(𝐉_0)`$ and $`(𝐉_0,𝐰_0)`$ is approximately the location of the center of the satellite at the disruption time $`t=0`$.

Similarly, the trajectory of a star in the cloud is $`(𝐉_0+\mathrm{\Delta }𝐉,𝐰_0+\mathrm{\Delta }𝐰+𝛀(𝐉_0+\mathrm{\Delta }𝐉)t)`$, where $`|\mathrm{\Delta }𝐉|𝐉_0`$ and $`|\mathrm{\Delta }𝐰|1`$. This can be simplified to

$$(𝐉(t),𝐰(t))=(𝐉_c,𝐰_c(t))+(\mathrm{\Delta }𝐉,\mathrm{\Delta }𝐰+𝐇\mathrm{\Delta }𝐉t)+\mathrm{O}(\mathrm{\Delta }𝐉,\mathrm{\Delta }𝐰)^2,$$
(2)   
where

$$H_{ik}=^2H(𝐉_0)/J_iJ_k$$
(3)   
is the Hessian of the Hamiltonian. At large times the extent of the cloud is dominated by the terms $`t`$ in equation (2), so we may write approximately

$$𝐉(t)𝐉_0+\mathrm{O}(\mathrm{\Delta }𝐉),𝐰(t)𝐰_c(t)+𝐇\mathrm{\Delta }𝐉t+\mathrm{O}(\mathrm{\Delta }𝐰).$$
(4)   
Since $`𝐇`$ is symmetric, it is diagonalizable, that is, there exists an orthogonal matrix $`𝐀`$ such that

$$\mathrm{𝐀𝐇𝐀}^1=𝐃(𝝀),$$
(5)   
where $`D_{ij}(𝝀)=\lambda _i\delta _{ij}`$, $`𝝀=\{\lambda _i\}`$ are the eigenvalues of $`𝐇`$, and $`𝐀^t=𝐀^1`$. We now make a canonical transformation to new action-angle variables $`(𝐉^{},𝐰^{})`$ using the mixed generating function $`S(𝐉^{},𝐰)=𝐉^{}\mathrm{𝐀𝐰}`$; thus

$$𝐰^{}=\frac{S}{𝐉^{}}=\mathrm{𝐀𝐰},𝐉=\frac{S}{𝐰}=𝐀^t𝐉^{},𝐰=𝐀^t𝐰^{},𝐉^{}=\mathrm{𝐀𝐉}.$$
(6)   
Equation (4) now simplifies to

$$J_i^{}(t)J_{0i}^{}+\mathrm{O}(\mathrm{\Delta }J^{}),w_i^{}(t)w_{ci}^{}(t)+\lambda _i\mathrm{\Delta }J_i^{}t+\mathrm{O}(\mathrm{\Delta }w^{}).$$
(7)   
This result shows that the shape of the expanding cloud is largely determined by the eigenvalues $`𝝀`$ of the Hessian $`𝐇`$, which are invariant under the orthogonal transformation $`𝐀`$. If two of the three eigenvalues are zero (say $`\lambda _2=\lambda _3=0`$) then the cloud expands only along the $`w_1^{}`$ direction in phase space, yielding a one-dimensional tidal streamer that grows linearly with time. If one of the eigenvalues is zero (say $`\lambda _3=0`$) then the cloud expands in two dimensions in phase space. However, if in addition $`|\lambda _1\mathrm{\Delta }J_1^{}||\lambda _2\mathrm{\Delta }J_2^{}|`$ or vice versa, then the expansion will effectively be one-dimensional. Thus the dimension $`D`$ has been reduced from 6 to 1 or 2, because the cloud expands much faster in some angles than in others, and not at all in action space.

We illustrate these remarks with some examples. The Hamiltonian for the triaxial harmonic oscillator may be written $`H(𝐉)=_{k=1}^3\omega _kJ_k`$, where $`\omega _k`$ is the frequency along axis $`k`$. The Hessian for this Hamiltonian is zero, so the disrupted system does not expand at all: the frequencies are independent of the actions and phase mixing does not occur. The Hamiltonian for the Kepler potential is $`H(𝐉)=\frac{1}{2}(GM)^2/J_1^2`$, where $`J_1=(GMa)^{1/2}`$ and $`a`$ is the semi-major axis. The eigenvalues of the Hessian are $`(3/a^2,0,0)`$, so the expanded cloud is one-dimensional, even if its original state was 6-dimensional. Simulations of disruptions of small galaxies in a spherical logarithmic potential yield one-dimensional clouds as well (Johnston et al. (1996)), but the precise conditions determining the effective dimension of the cloud resulting from tidal disruption of a small stellar system have not yet been studied.

## 3 Catastrophes

We now survey the phase space by measuring $`K`$ of the $`N`$ phase-space coordinates of each star (we call this $`K`$-dimensional space the “observable space”<sup>1</sup><sup>1</sup>1Other names include control space or external variables.), and ask how to detect the $`D`$-dimensional phase structure. For example, an image or set of sky positions has $`K=2`$; a data cube (positions and radial velocities) has $`K=3`$; a proper-motion survey has $`K=4`$, etc. A bolometric laboratory detector of cold dark matter measures only their energy and moreover is restricted to a single point in configuration space.

If $`D<K`$ the phase structure will appear as a distinct entity in the observable space, which can be detected by standard techniques: plotting the data on paper ($`K=2`$), visualization software ($`K=3`$), cluster-finding algorithms, etc. More specialized and powerful techniques can be used if the phase structure has known properties: for example, the disruption of distant satellite galaxies in a spherical potential leads to phase structures which are great circles on the sky (Lynden-Bell and Lynden-Bell (1995); Johnston et al. (1996)).

In this discussion we concentrate the more challenging case when $`D=K`$, so that the phase structure covers most or all of the observable space. In this case, the most prominent features in the survey will be the singularities of the phase structure, that is, the locations where the projection of the phase structure into the observable space leads to a singularity. For example, if in Figure 1 the observations yield only the horizontal coordinate, then near-singularities arise in the plot at $`t=50`$ when the curve representing the phase structure is vertical.

We can generalize this simple concept to singularities with higher dimension using catastrophe theory (Poston and Stewart (1978); Berry and Upstill (1980); Gilmore (1981)). We begin by setting up a coordinate system $`(u_1,\mathrm{},u_D)`$ on the phase structure that contains the cloud. We complete the coordinate system of the phase space by $`(u_{D+1},\mathrm{},u_N)`$, which are chosen so that $`u_{D+1}=\mathrm{}=u_N=0`$ on the phase structure. Thus the phase structure is specified by

$$\frac{\varphi (𝐮)}{u_j}=0,j=1,\mathrm{},N\text{where}\varphi (𝐮)=\underset{j=D+1}{\overset{N}{}}u_j^2.$$
(8)   
The advantage of specifying the phase structure in terms of the generating function $`\varphi `$ is that $`\varphi `$ is a single-valued function of the phase-space coordinates.

Next let $`(x_1,\mathrm{},x_D)`$ be the coordinates of the observable space, and complete the coordinate system of the phase space by $`(y_1,\mathrm{},y_{ND})`$ (we call this $`(ND)`$-dimensional space the “hidden space”<sup>2</sup><sup>2</sup>2Other names include state space or internal variables.). For given values of the observable coordinates, the location or locations of the phase structure are specified implicitly by the equations

$$\frac{\varphi (𝐲;𝐱)}{y_j}=0,j=1,\mathrm{},ND.$$
(9)   
Singularities in the observable space arise when there are displacements in the hidden space that leave (9) unaltered, that is, for which

$$\underset{k=1}{\overset{ND}{}}\frac{^2\varphi (𝐲;𝐱)}{y_jy_k}dx_k=0,j=1,\mathrm{},ND.$$
(10)   
These equations have a non-trivial solution if the determinant of the $`(ND)\times (ND)`$ Hessian matrix vanishes, that is, if

$$\text{det}\frac{^2\varphi (𝐲;𝐱)}{y_jy_k}=0.$$
(11)   
We can restrict our attention to structurally stable singularities (roughly speaking, ones for which small changes in $`\varphi `$ lead to small changes in the locus of the singularity); the reason is that unstable singularities represent a subset of measure zero and hence are atypical, at least for $`N7`$ (Zeeman (1977)).

Catastrophes are structurally stable singularities of gradient maps. They are classified by their codimension, which is the dimension of the observable space minus the dimension of the singularity. Thus codimension 1 corresponds to a singular point in a survey with 1 observable, a singular line in a survey with 2 observables, etc. The corank of the catastrophe is the minimum value of $`ND`$, the dimensionality of the hidden space, which is required for the singularity to occur. A third classification is by the degree $`k`$ of the singularity: the mean density in a neighborhood of the observable space of radius $`ϵ`$ around the singularity is $`ϵ^k`$.

### 3.1 Fold

The fold is the only catastrophe with codimension 1. Since the corank of the fold is also 1, the simplest example occurs in a phase space with two dimensions, one of which is observable. For example, the generating function

$$\varphi (y;x)=\frac{1}{3}y^3+xy$$
(12)   
implies that the phase structure is the parabolic curve $`x+y^2=0`$, and this manifold has a fold catastrophe at $`x=0`$. If the linear density on the curve is uniform, the observable density is

$$\rho (x)\left(1\frac{1}{4x}\right)^{1/2},x<0,$$
(13)   
and zero otherwise. The square-root divergence in density on one side of the singularity is characteristic of a fold, and implies that the degree $`k=\frac{1}{2}`$.

Examples of fold catastrophes occur in the projection of a spherical shell onto a plane, rainbows, sunlight sparkling on the sea, twinkling of starlight, gravitational lensing, etc.

The simplest examples of fold catastrophes in phase mixing are shell structures in elliptical galaxies, which are believed to arise from disrupted companion galaxies. Shells are one-dimensional structures that appear in two-dimensional observable space (the two coordinates on the sky plane) and hence if they are singularities they must have codimension 1 and hence must be folds. Hernquist and Quinn (1988) carried out simulations of the evolution of disrupted companions in spherical host galaxy potentials, and several of their conclusions can be interpreted in terms of the results we have derived so far. They find that shell formation requires the accretion of either a dynamically cold or a spatially compact companion; these are precisely the conditions required so that the cloud dimension $`D2`$2), which in turn is required so that singularities are present in a two-dimensional observable space. They distinguish “spatial folding” (projection from three-dimensional configuration space to two-dimensional space) from “phase-wrapping” (projection from phase space into three-dimensional space) but this distinction is not fundamental: both are projections from phase space into observable space. Finally, they point out that some shells arise from step functions in surface density rather than square-root singularities; the former are not singularities and hence are not described by catastrophe theory.

### 3.2 Cusp

The cusp is the only catastrophe with codimension $`2`$, and has corank 1. The cusp is a singular point in a 2-dimensional observable space, a line in 3-dimensional observable space, etc. For example, consider a phase space with $`N=3`$ dimensions, two of which are observable. The generating function

$$\varphi (y;𝐱)=\frac{1}{4}y^4+\frac{1}{2}x_2y^2+x_1y$$
(14)   
implies that the phase structure is given by the surface $`y^3+x_2y+x_1=0`$, which is singular (in the sense of eq. 11) along the lines

$$x_1=\pm \frac{2}{3^{3/2}}\left(x_2\right)^{3/2}.$$
(15)   
Each of these lines is a fold catastrophe, and their junction is a cusp catastrophe.

If the surface density of stars is uniform on the phase structure, the surface density in the observable space $`(x_1,x_2)`$ is given parametrically by

$$\mathrm{\Sigma }(x_1,x_2)\left[1+\frac{1+y^2}{(x_2+3y^2)^2}\right]^{1/2},\text{where}x_1=y^3x_2y.$$
(16)   
At a distance $`\mathrm{\Delta }x_2`$ from the fold ($`|\mathrm{\Delta }x_2||x_2||x_1|^{2/3}1`$), the density $`\mathrm{\Sigma }|x_1|^{2/3}`$ for $`\mathrm{\Delta }x_2>0`$ and $`\mathrm{\Sigma }|x_1|^{1/3}(\mathrm{\Delta }x_2)^{1/2}`$ for $`\mathrm{\Delta }x_2<0`$ (the characteristic fold behavior). The degree of the cusp singularity is $`k=\frac{2}{3}`$.

To illustrate the appearance of a cusp catastrophe, the first three panels of Figure 2 show a survey that samples the surface density (16) with 300, 1000, and 3000 data points in the interval $`x_1,x_2[1,1]`$. The fourth panel shows a higher-resolution survey with 3000 points in $`x_1,x_2[0.1,0.1]`$. The general cusp in a 2-dimensional observable space is formed from this special case by smooth distortion (i.e. diffeomorphism) of the observable space, and still consists of two fold catastrophes meeting with a common tangent, as in equation (15).

Examples of cusp catastrophes include the caustic curve seen on the surface of coffee in a cup, and the critical point of the van der Waals equation of state in the pressure-temperature plane.

One might speculate that some of the dwarf spheroidal satellite galaxies are cusp catastrophes rather than bound equilibrium stellar systems, but their isopleths are approximately elliptical (Irwin and Hatzidimitriou (1995)) and do not resemble the density contours that are expected near a cusp.

### 3.3 Catastrophes with codimension 3

There are three catastrophes with codimension $`3`$, the swallowtail, elliptic umbilic, and hyperbolic umbilic. The swallowtail has corank 1, and thus first appears in a phase space with $`N4`$ dimensions; the umbilics have corank 2 and thus require $`N5`$. All of these singularities are points in a 3-dimensional observable space, lines in 4-dimensional observable space, etc. Just as a cusp occurs at the junction of two fold curves, these occur at the junction of cusp curves, which are connected by fold surfaces.

#### 3.3.1 Swallowtail

Consider a phase space with $`N=4`$ dimensions, three of which are observable. The generating function

$$\varphi (y;𝐱)=\frac{1}{5}y^5+\frac{1}{3}x_3y^3+\frac{1}{2}x_2y^2+x_1y,$$
(17)   
implies that the phase structure is given by $`y^4+x_3y^2+x_2y+x_1=0`$, which is singular (in the sense of eq. 11) along the surfaces defined in terms of the parameters $`(y,x_3)`$ by

$$x_1=3y^4+x_3y^2,x_2=2(2y^3+x_3y).$$
(18)   
Each of these surfaces is a fold catastrophe; the folds meet at cusp lines defined by

$$x_1=3y^4,x_2=8y^3,x_3=6y^2.$$
(19)   
The fold surfaces intersect one another along the parabola $`x_3^2=4x_1`$. The degree of the swallowtail is $`k=\frac{3}{4}`$.

Figure 3 shows the fold surfaces and cusp lines for the swallowtail catastrophe.

#### 3.3.2 Elliptic umbilic

Consider a phase space with $`N=5`$ dimensions, three of which are observable. The generating function

$$\varphi (y_1,y_2;𝐱)=y_1^33y_1y_2^2x_3(y_1^2+y_2^2)x_2y_2x_1y_1$$
(20)   
implies that the phase structure is given by the surface

$$3y_1^23y_2^22x_3y_1x_1=0,6y_1y_2+2x_3y_2+x_2=0.$$
(21)   
which is singular along the surfaces defined parametrically by

$$x_1=3y_1^23y_2^2\pm 6y_1(y_1^2+y_2^2)^{1/2},x_2=6y_1y_2\pm 6y_2(y_1^2+y_2^2)^{1/2},x_3=(y_1^2+y_2^2)^{1/2}.$$
(22)   
Each of these surfaces is a fold catastrophe; the folds meet at three parabolic cusp lines defined by $`x_1^2+x_2^2=x_3^4`$ and confined to planes at angles of $`120\mathrm{deg}`$ that intersect along the $`x_3`$-axis (see Figure 4). The degree is $`k=1`$.

#### 3.3.3 Hyperbolic umbilic

The generating function

$$\varphi (y_1,y_2;𝐱)=y_1^3+y_2^3x_3y_1y_2x_2y_2x_1y_1$$
(23)   
implies that the phase structure is given parametrically by

$$x_1=3y_1^2+ϵ_1(y_1y_2^3)^{1/2},x_2=3y_2^2+ϵ_1(y_1^3y_2)^{1/2},x_3=ϵ_2(y_1y_2)^{1/2},$$
(24)   
where $`ϵ_i=\pm 1`$ and $`y_1,y_20`$. Each of these surfaces is a fold catastrophe; the folds intersect along the coordinate axes $`x_1>0,x_2=x_3=0`$ and $`x_1=0,x_2>0,x_3=0`$ (Fig. 4). The cusp line is defined parametrically by

$$x_1=x_2=9y^2,x_3=6ϵ_2y,y>0.$$
(25)   
The degree is $`k=1`$.

These exhaust the catastrophes with codimension $`3`$. Catastrophes with higher codimension are only relevant when the observable space has dimension $`D4`$. Catastrophes of higher codimension become increasingly complicated, and are less important since we have argued in §2 that typical phase structures will often have dimension 3 or less.

## 4 Summary

It is important to understand the geometry of phase structures before searching for them. We have focused on two specific aspects of this geometry: the expected dimension of the phase structure in (usually 6-dimensional) phase space, and the classification of singularities when a $`D`$-dimensional phase structure is projected into a $`D`$-dimensional survey.

Phase mixing of disrupted small, hot galaxies generally leads to 3-dimensional phase structures but the dimensionality can be smaller if the potential of the host galaxy has special symmetries (e.g. spherical or Keplerian), or if the eigenvalues of the Hessian (eq. 3) are very different in magnitude. Phase mixing of large, cold galaxies leads to 2-dimensional phase structures.

The structurally stable singularities in $`D`$-dimensional surveys are folds ($`D1`$), cusps ($`D2`$), swallowtails, elliptic umbilics and hyperbolic umbilics ($`D3`$).

Even though phase mixing is a simple process, there remain many unresolved theoretical issues. How fast do small-scale gravitational irregularities or large-scale orbital chaos disrupt phase structures? Do phase structures provide a significant source of relaxation in galaxies (Tremaine and Ostriker (1998))? What are the properties of phase structures in galaxies that we expect from standard models of structure formation? What statistical measures can we use to characterize phase structures and their projections?