The fine-grained phase-space structure of Cold Dark Matter halos

The fine-grained phase-space structure of Cold Dark Matter halos

Mark Vogelsberger, Simon D.M. White, Amina Helmi, Volker Springel
Max-Planck Institut fuer Astrophysik, Karl-Schwarzschild Strasse 1, D-85748 Garching, Germany
Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, Netherlands
Accepted ???. Received ???; in original form ???
Accepted ???. Received ???; in original form ???

We present a new and completely general technique for calculating the fine-grained phase-space structure of dark matter throughout the Galactic halo. Our goal is to understand this structure on the scales relevant for direct and indirect detection experiments. Our method is based on evaluating the geodesic deviation equation along the trajectories of individual DM particles. It requires no assumptions about the symmetry or stationarity of the halo formation process. In this paper we study general static potentials which exhibit more complex behaviour than the separable potentials studied previously. For ellipsoidal logarithmic potentials with a core, phase mixing is sensitive to the resonance structure, as indicated by the number of independent orbital frequencies. Regions of chaotic mixing can be identified by the very rapid decrease in the real space density of the associated dark matter streams. We also study the evolution of stream density in ellipsoidal NFW halos with radially varying isopotential shape, showing that if such a model is applied to the Galactic halo, at least streams are expected near the Sun. The most novel aspect of our approach is that general non-static systems can be studied through implementation in a cosmological N-body code. Such an implementation allows a robust and accurate evaluation of the enhancements in annihilation radiation due to fine-scale structure such as caustics. We embed the scheme in the current state-of-the-art code GADGET-3 and present tests which demonstrate that N-body discreteness effects can be kept under control in realistic configurations.

cold dark matter, phase-space structure, dynamics, N-body
pagerange: The fine-grained phase-space structure of Cold Dark Matter halosReferencespubyear: 2007

1 Introduction

Dark Matter (DM) is still to be directly detected. The first indirect indications of its existence came in the 1930s, with measurements of the velocities of galaxies in clusters. The cluster mass required to gravitationally bind the galaxies was found to be more than an order of magnitude larger than the sum of the luminous masses of the individual galaxies (Zwicky, 1933; Smith, 1936). The early detection of large amounts of unseen matter associated with the Local Group (Kahn & Woltjer, 1959) was followed in the 1970s by observations of the rotation curves of spiral galaxies which showed that these are flat, or even rising, at distances far beyond their stellar components (Rubin & Ford, 1970; Faber & Gallagher, 1979; Rubin et al., 1980). Studies of satellite systems suggested that the mass distributions of most galaxies might be an order of magnitude larger and more massive than their visible parts (Ostriker et al., 1974; Einasto et al., 1974). All these discoveries led to the conclusion that a large fraction of mass in the Universe is dark. This has also been supported by recent gravitational lensing studies that demonstrate the existence of extended massive DM halos (e.g. Mandelbaum et al. 2006).

As the dominant mass component of galaxies and large-scale structures, DM has necessarily become a key ingredient in theories of cosmic structure formation. The most successful of these theories is the hierarchical paradigm. In the current version of this model, the DM is composed of nonbaryonic elementary particles known as Cold Dark Matter (CDM) (Peebles, 1982). The term “cold” derives from the fact that the DM particles had negligible thermal motions at the time of matter-radiation equality. Their abundance was set when the interaction rate became too small for them to remain in thermal equilibrium with other species in the expanding Universe. The first objects to form in a CDM Universe are small galaxies. They then merge and accrete to give rise to the larger structures we observe today. Thus structure formation occurs in a “bottom-up” fashion (Blumenthal et al., 1984; Davis et al., 1985; Springel et al., 2006).

The crucial test of this paradigm undoubtedly consists in the determination of the nature of DM through direct detection experiments. Among the most promising candidates from the particle physics perspective are axions and neutralinos. Axions were introduced to explain the absence in Nature of strong-CP (Charge conjugation and Parity) violations (Peccei & Quinn, 1977). The cosmological population formed out of equilibrium as a zero momentum Bose condensate. They can be detected through their conversion to photons in the presence of a strong magnetic field, (e.g. Ogawa et al. 1996; Hagmann et al. 1998). Neutralinos are the lightest supersymmetric particles, and may be considered as a particular form of weakly interacting massive particles (WIMPs) (Steigman & Turner, 1985). The most important direct detection process for neutralinos is through elastic scattering on nuclei.

Today many experiments are searching for these particles, (Akerib & et al, 2004; Sanglard & et al, 2005; Schnee, 2006; Aprile et al., 2007; Spooner, 2007). The main challenge lies in the large background they encounter. Having an idea what the detector signal might look like can help substantially in fine-tuning the experiments in order to increase the chance of a detection. In addition, many experiments are attempting to detect WIMPs indirectly by searching for -ray emission from their annihilation (de Boer, 2005; de Boer et al., 2005; Bergström & Hooper, 2006; Hooper & Serpico, 2007). Predictions for this radiation are currently uncertain because very substantial enhancements are possible, at least in principle, from fine-scale structure such as caustics in the dark matter distribution (Dalcanton & Hogan, 2001; Mohayaee et al., 2007).

The differential detection rate for WIMPs in a given material (per unit detector mass) is sensitive to the local velocity distribution:


where is the energy deposited in the detector and is the differential cross section for WIMP elastic scattering with the target nucleus. is the WIMP mass which lies in the GeV to TeV range for currently preferred models. is the local mass density and the velocity distribution of WIMPs that reach the detector.

From Eq. (1) it is evident that the count rate in a direction-insensitive experiment depends on the velocity distribution of the incident particles and will be modulated by the orbital motion of the Earth around the Sun (Drukier et al., 1986). In most studies, an isotropic Maxwellian distribution relative to the galactic halo has been assumed, (e.g. Freese et al. 1988), although there are other models in the literature, discussing, for example, multivariate Gaussians (Evans et al., 2000). Some attempts at understanding the effect of fine-scale structure in the velocity distribution of DM particles have also been made (Sikivie, 1998; Stiff et al., 2001; Hogan, 2001).

A significant signal could come from what are known as streams of DM (Sikivie et al., 1995; Helmi et al., 2002; Natarajan & Sikivie, 2005). Both axions and WIMPs are cold. In the absence of clustering their present day velocity dispersion would be negligible ( for WIMPs and for axions). They are effectively restricted to a 3D hypersurface, a sheet in 6D phase-space. The growth of structure results in continual stretching and folding in phase-space of this initially almost uniform sheet. This process is called mixing. The more strongly the system mixes, the more streams of DM particles are present at a given location in configuration-space. Mixing stretches each sheet and so its density decreases. The existence of distinct streams is a direct consequence of the collisionless character and the coldness of CDM. The situation is sketched in Fig. 1. At the points where the number of streams changes by two, the local configuration-space density of dark matter becomes extremely high. These are caustics of the kind studied by catastrophe theory (Gilmore, 1982; Tremaine, 1999). Note that Liouville’s Theorem prevents the CDM phase-space sheet from ever tearing, although it can be arbitrarily strongly stretched.

Figure 1: Sketch of an idealised CDM phase-space sheet in the plane. The thickness of the line represents the local velocity dispersion within each stream. The small wiggles correspond to initial density perturbations and the multi-valued region reflects the multiple streams created by winding in non-linear regions. Depending on position an observer sees one or three streams. At points where the number of streams changes by two, a caustic with a very high DM density is present. Such caustics may be significant for the total annihilation flux. The number of streams at each point is a measure of the local amount of mixing. The dashed line represents the Hubble Flow. The cross marks the phase-space coordinates of a particular CDM particle embedded in the flow.

If the dark matter density in the solar neighbourhood is dominated by a small number of streams, its velocity distribution will effectively consist of a few discrete values, one for each stream. If, on the other hand, the number of streams is very large, the velocity distribution will be smooth and individual streams will be undetectable. This issue has so far been addressed only under simplified conditions and divorced from its proper cosmological context (Helmi et al., 2002; Natarajan & Sikivie, 2005). This is largely because the only tool capable of studying cosmological structure formation in full generality, namely N-body simulations, cannot resolve the relevant scales. For example, Natarajan & Sikivie (2005) estimate that of the order of particles would be required in a simulation of the Milky Way’s halo to begin to resolve the streams in the solar neighbourhood. Even the largest simulation so far published, the Via Lactea model (Diemand et al., 2007), is 4 orders of magnitude short of this minimum requirement. It is thus impossible to figure out the number of streams near the Sun and the properties of “typical” streams with current N-body capabilities.

A related issue that has been much discussed is whether caustics can affect direct or indirect detection experiments. Some authors claim that caustics play an important role as stable phenomena connected to any collapsing CDM halo and taking the form of relatively massive rings or shells (Sikivie, 1999; Natarajan & Sikivie, 2006). Other authors have argued that caustics should be more abundant, weaker and dynamically negligible (Helmi et al., 2002). Even if caustics negligibly affect the gravitational potential, they may have very substantial effects on the annihilation rate of dark matter (Hogan, 2001; Bergström et al., 2001; Pieri & Branchini, 2005; Mohayaee & Shandarin, 2006; Mohayaee et al., 2007; Natarajan, 2007). All these papers were able to evaluate the enhancements due to caustics only under restrictive and unrealistic assumptions about symmetry, formation history, etc. While they demonstrate that large enhancement factors may be possible, they do not provide reliable estimates for the actual enhancement expected. The method we present below is capable of providing such estimates for halos growing as predicted by the CDM model.

In this paper we will present a novel approach that directly analyses structure in the fine-grained phase-space distribution. We circumvent the “number of particle” problem by solving the geodesic deviation equation (GDE) for every DM particle. This allows us to calculate the local properties of the DM stream each simulation particle is embedded in, in particular, its configuration-space density and its local velocity distribution. The mass-weighted number of streams near any point is then the total local density divided by the mean density of the individual local streams. Caustic passages can be detected robustly from the properties of each particle’s local stream, and the particle’s instantaneous annihilation probability within this stream can be evaluated simply and integrated accurately through caustics.

This paper is the first of a series. Here we introduce our method and apply it to study the evolution of streams in static potentials that are too complex to be analysed using previous techniques. We also implement our scheme in a state-of-the-art N-body simulation code, and use simple test problems to demonstrate that N-body discreteness effects can be kept under control in realistic applications. Later papers will address issues associated with mixing, caustics and annihilation radiation in the full cosmological context.

The outline of our paper is the following: In Section 2 we present a detailed derivation of the GDE and show how it can be used to quantify mixing, to locate caustics, and to calculate stream densities and annihilation rates. Section 3 describes our code, DaMaFlow, which is designed to solve the GDE for single orbits in a wide variety of potentials. In Section 4 we analyse static, separable potentials and compare results from our method to previous work. The fifth section applies our scheme to non-integrable, but still static potentials, revealing their complex phase-space structure. In Section 6 we turn to more realistic non-spherical CDM halo potentials and discuss the influence of halo shape on stream density behaviour. Section 7 discusses the implementation of our method in an N-body code and presents results of simple tests of when discreteness effects compromise studies of stream densities and caustics. The final section summarises our results and gives some conclusions.

2 The geodesic deviation equation

Our scheme for calculating the evolution of the fine-grained phase-space distribution in the neighbourhood of a DM particle is based on the evolution of the distance between two infinitesimally separated particle trajectories. This evolution is described by the Geodesic Deviation Equation (GDE). We use the following notation to clearly distinguish between three- and six-dimensional quantities: an underline denotes a vector and two of them denote a matrix. An overline denotes a vector and two of them denote a matrix. Thus a general phase-space vector is composed of two three-dimensional vectors: .

To derive the GDE we first write down the equations of motion for a DM particle. These are simply


with initial conditions .

As the equation of motion in phase-space can be written:


with initial conditions .

We want to calculate the local stream density around the DM particle whose trajectory in phase-space is given by . To do so, we first ask how the displacement vector to a neighbouring DM particle in phase-space evolves with time:


Note that is the displacement between the reference DM particle, which was at at time , and another particle which was at at . We are interested in properties in the immediate neighbourhood of the reference particle, so is an infinitesimal displacement. This allows us to work to linear order:


Introducing the phase-space distortion tensor (note that this is a tensor)


we can rewrite Eq. (5) as a simple linear transformation from the starting phase-space displacement to the displacement at time :


Because is an arbitrary displacement vector, the distortion tensor describes how the complete local phase-space neighbourhood around the reference DM particle gets distorted while it is orbiting in the given potential. The time evolution of follows from the time evolution of the two trajectories. Again we can work this out in linear order111As can be chosen arbitrarily small, this is always possible.:


To derive the equation of motion for the distortion tensor itself we evaluate Eq. (8) for six unit vector phase-space displacements with where , and is the Kronecker delta. Taking into account Einstein’s sum convention this yields for each component of Eq. (8)


where we have introduced the phase-space tidal tensor


and is the configuration-space tidal tensor given by the second derivatives of the gravitational potential . As we are are only interested in linear order we replace by in Eq. (9) and get an equation of motion for the 6D distortion tensor:


Note that this first-order tensor differential equation represents a system of 36 coupled ordinary first-order differential equations. To solve them we need to specify initial conditions. These follow from the constraint


DM behaves like a collisionless fluid and its fine-grained phase-space density is described by the well-known Vlasov equation:


Using the Lagrangian derivative this reads , meaning that the local fine-grained phase-space density has to be conserved along the orbit of every particle in the system. Imagine particles that fill a small phase-space volume around the DM reference particle at time . At a later time these particles fill a volume around the reference particle at . Conservation of phase-space density implies that the two volumes are the same . The evolution of the initial displacement vectors from the reference particle to one of the other particles in the volume is described by the distortion tensor associated with the reference particle:


The change in volume due to this linear transformation is given by the determinant . As this volume has to be conserved, the determinant of the phase-space distortion tensor has to be conserved. Note that not only must the volume be conserved, but also the orientation of the volume element. This means that the sign of the determinant is also fixed. From the initial conditions Eq. (12) one gets =1 at all times. This fact can be used to check numerical solutions of the equations.

The structure of Eq. (11) allows the equations of motion for the distortion to be broken down to a set of equations that is more convenient to work with. Let us first rewrite Eq. (11) using blocks of tensors:222We suppress the argument to avoid confusion.

Writing down the equation for each matrix block yields four equations:




These can be combined to give:


Thus we get two identical differential equations of second-order for two tensors whose dynamics is driven by the ordinary tidal tensor. From the initial condition for the 6D distortion one can see that the only difference between and lies in the appropriate initial conditions: , and , in the two cases. From the solutions of these two initial condition problems the 6D distortion solution can then be constructed:


Up to this point we have worked out all equations in phase-space. As we are interested in the stream density in configuration-space, we need to project down to this space. As already mentioned, CDM lies on a thin sheet in phase-space. This sheet has a certain orientation at the starting point of the reference DM particle. Take to be the local parametrisation of the sheet surrounding this particle at time 333Such a parametrisation is always possible locally, but due to mixing there is, in general, no simple global parametrisation of the stream. This is only possible for very early times, where the sheet is dominated by the Hubble Flow, the relation is one-to-one, and the stream density is almost uniform.. Now we ask how an infinitesimal displacement in configuration-space is distorted by evolution. First we note that any displacement in implies a displacement in velocity-space due to the restriction of particles to the sheet:


The phase-space distortion describes how the corresponding phase-space displacement evolves. We are here interested in the configuration-space part of the phase-space displacement at time :


The evolution of the displacement in configuration-space can also be described by a linear transformation:

where we have introduced the configuration-space distortion tensor (note that this is a tensor)


This tensor can also be derived with the help of two projection operators:


As in the case of phase-space distortion, the change in volume due to the linear transformation in Eq. (2) is given by the determinant, so that the stream density in configuration-space is proportional to the inverse of this determinant:


At time the configuration-space distortion tensor equals unity. Thus, if we norm the stream density to its initial value, we get the following relation for the normed stream density:


In the rest of this paper we will almost always discuss this normed stream density.

In Fig. 2 we sketch the distortion of the infinitesimal cloud around the reference DM particle. Note the difference between the stream density evolution in configuration-space and in phase-space. The volume of the small cloud grows in Fig. 2 and is not constant anymore! This is a result of the projection from phase-space to configuration-space. Nearby trajectories spatially diverge in time.

Figure 2: The configuration-space distortion tensor describes how an initial small configuration-space displacement evolves in time. This reflects the stretching of an infinitesimally small cloud of virtual particles around the reference particle that is placed at at time . The stretching of the cloud is driven by the tidal field that the DM particle encounters as it orbits.

This formalism can be used to identify caustics in a very efficient way. If the DM particle passes through a caustic passes through zero, and the stream density goes to infinity (for perfectly cold DM). A small cloud surrounding the reference particle turns inside out as it passes through the caustic. This corresponds to a change in sign of the determinant, and thus can easily identified numerically. This property allows the location of caustics to be mapped accurately even in complex configurations. Note that the possibility of sign changes explains why we took the modulus of the determinant in Eq. (23).

The complex fine-grained phase-space structure and especially the caustics expected in CDM halos are likely to substantially enhance annihilation radiation. These effects have so far been analysed only under quite simplified conditions (Hogan, 2001; Bergström et al., 2001; Pieri & Branchini, 2005; Mohayaee & Shandarin, 2006; Mohayaee et al., 2007; Natarajan, 2007). As a result, it is unclear how strong such enhancement effects will be in the proper cosmological context. Previous studies of annihilation radiation from N-body halos had realistic formation histories but were unable to resolve caustics. They estimated emissivities from the local mean CDM density, thus effectively excluding contributions from single streams (Stoehr et al., 2003; Diemand et al., 2007). Hogan (2001) noted that this results in an underestimation of the annihilation rate, and suggested that annihilation might, in fact, be dominated by contributions from the neglected caustics (which he baptised as “micropancakes”).

Our formalism enables a robust and accurate calculation of the contribution to the annihilation radiation from individual streams. The annihilation rate for each particle due to encounters with other particles in its own stream can be evaluated directly from the local stream velocity distribution and density. Integrating these rates along the trajectories of all particles produces a Monte Carlo estimate of the intrastream annihilation rate for the system as a whole. This automatically includes the contributions from all caustics and is exactly the fine-scale contribution which is missing from standard N-body-based estimates of annihilation luminosities.

We now discuss briefly how this is done. Given the very small primordial velocity dispersion of CDM, and approximating the initial configuration-space density as a constant we can write the phase-space density around a reference particle at the initial time as follows: 444 denotes the transpose of a matrix


and . Note that this represents a Gaussian distribution in velocity-space and a constant density in configuration-space.

Using we obtain the phase-space density around the particle at the later time :


The configuration-space density around the reference particle at is simply the integral of the phase-space density over all velocities evaluated at :


where the velocity dispersions are given by and are the eigenvalues of the velocity submatrix of . The velocity distribution in the principal axis frame of the velocity ellipsoid centred on the particle’s position and velocity is given by:

where . Note that this velocity distribution is normalised, i.e. its integral over velocity space is unity.

These quantities allow us to calculate the instantaneous annihilation rate at each point on the particle’s trajectory:


where is the particle mass and the annihilation cross-section. We note that in many WIMP models the annihilation cross-section is inversely proportional to encounter velocity, and in this case is independent of the actual local velocity distribution (Jungman et al., 1996). An image of the system in annihilation radiation can be constructed by integrating all particles forward over a short time interval and summing their annihilation contributions into a pixel array. We note that equation (26) exhibits near-singular behaviour as particles pass through caustics and as a result special care is needed to obtain the correct contribution to the annihilation luminosity in this situation. This will be discussed more fully in later papers.

3 The DaMaFlow code

We have developed the code DaMaFlow to test our GDE scheme by analysing the behaviour of streams in a broad range of static potentials. DaMaFlow integrates the equations of motion and in parallel the GDE for a single orbit within user-specified potentials. This requires solving second-order differential equations in parallel. The integration algorithm can be chosen to be a symplectic second-order leapfrog (Drift-Kick-Drift or Kick-Drift-Kick formulation) or alternatively the energy-conserving and adaptive Dop853 algorithm (Hairer et al., 1993) of order that allows dense output and is very fast. Studies focusing on complex phase-space structures, especially in the field of chaos analysis, often use the Dop853 algorithm (or even higher order schemes) because of its high precision. On the other hand N-body codes often implement the leapfrog method because it is the best compromise between performance and accuracy. We find that with a moderate fixed time-step, both formulations of leapfrog are able to give comparable results to Dop853. This is an important point because it is not possible to run N-body simulations with slow but accurate high-order ODE solvers like Dop853.

DaMaFlow is also set up to do a Numerical Analysis of Fundamental Frequencies (NAFF) (Laskar, 1988, 1990; Papaphilippou & Laskar, 1996; Laskar, 2003) and of the resonances associated with the chosen orbit. The fundamental frequencies are revealed by an integer programming routine. This is needed so that we can study the relation between the existence of well-defined fundamental frequencies or resonances and stream density behaviour. The NAFF method determines a quasi-periodic approximation to the orbital motion. For ordinary Fast-Fourier-Transforms (FFT) the accuracy of the determination of the frequencies is of the order of , where is the sampling interval. The NAFF method has an accuracy of . Thus it makes spectral analysis a lot faster compared to classical Fourier techniques, for example, those used in Binney & Spergel (1982).

To scan large parts of phase-space, DaMaFlow can be run in parallel on batch systems in order to integrate a large number of different orbits simultaneously. A fast automated stream density fitting routine was built in to facilitate efficient analysis of the underlying phase-space without user interaction. Before this fitting can be done, the stream density has to be smoothed to remove the large density spikes produced by caustics. DaMaFlow does this by extracting and fitting the lower envelope of the stream density time series. This envelope is constructed while the orbit integration is running by an iterative on-the-fly minimum finder.

From a numerical point of view, solving the GDE is quite difficult in chaotic regions of phase-space. In these regions the infinitesimal phase-space volume around the reference particle gets distorted very strongly. This produces large numerical values in the phase-space distortion tensor. And this can lead to overflows and round-off errors in numerical computations. In chaos analysis it is an established method to do some kind of re-norming to suppress these problems. We can do something similar to follow the evolution of phase-space density evolution. We use the transitivity of the phase-space distortion tensor:


where with is the solution of the GDE with evaluated at time . So the phase-space volume can be written as:


Thus dividing the time integration interval and re-initialising the distortion after each interval avoids large numerical values. This is very similar to the re-norming techniques used for calculating the largest Lyapunov exponents of chaotic systems, where the re-norming frequency is chosen to be of the order of the dynamical time-scale (Lichtenberg & Lieberman, 1983; El-Zant, 2002). Although this approach works nicely for such applications, it does not help us when calculating the stream density evolution, because we need the entire phase-space distortion information from initial to final time. It is not possible to separate the configuration-space part of the phase-space distortion and to do a similar re-norming. Thus one cannot avoid large numbers during the calculation. DaMaFlow therefore calculates all quantities in double-precision ( bits= bytes). Even in chaotic regions this is enough to follow the system for a substantial amount of time. We note that the phase-space density calculation involves the determinant of a matrix, whereas the stream density only involves the determinant of a matrix. As a result stream density calculations are less strongly affected by large numbers and overflows. We note that special software libraries can provide even higher precision (e.g. the GMP library 555

Since we wish to implement the GDE formalism also into an N-body code, execution speed and memory consumption are important considerations. For each particle we need to store the full 6D phase-space distortion tensor and the particle’s tidal tensor. This results in extra numbers per particle. Thus a particle simulation in double-precision needs already about  Gbytes of random access memory (RAM) just for the GDE calculation, assuming we store all information for every particle.

4 Integrable potentials

Figure 3: Stream density evolution for different initial sheet orientations in a static Hernquist potential. The general shape of the four curves is very similar. The long-term behaviour does not in general depend on the initial sheet orientation.

Before presenting some results for the evolution of stream densities in integrable potentials we want briefly to discuss the choice of the initial sheet orientation that goes into in equation (19). This depends, of course, on the starting time and the problem that is to be studied. For example, in a cosmological context is given by the linear initial conditions, so by the coupled initial density and velocity fields. We note that, to zeroth order, the sheet orientation is simply given by the Hubble Flow: .

What is the impact of the initial orientation on later evolution? In Fig. 3 we show the evolution of the normed stream density, as calculated from eq. (24), for a single orbit with four different choices of initial stream orientation. For this test we have used a spherical Hernquist potential:


with and kpc. The reference particle begins its orbit with kpc and . The initial sheet orientations were chosen to be (in units of ): , ,

It is striking that all four curves have very similar shape and very similar caustic spacings, although the caustic locations vary. The long-term behaviour does not depend on initial sheet orientation, at least in this case. Orientation produces lower densities than the others because the scale of is larger, but the shape of the lower envelope is very similar.

Natarajan & Sikivie (2006) show that the caustic shape in configuration-space is, in general, affected by the relative size of the matrix elements. A detailed analysis of caustic shape thus requires choosing . For example, in their model for the Milky Way halo Natarajan & Sikivie (2006) initialised trajectories at the turnaround sphere with a loosely motivated by tidal torque theory. This restricted the form and scale of the matrix, but still left a lot of freedom. Here our main motivation is not to analyse caustic shapes, but rather the long-term behaviour of the fine-grained phase-space distribution, in particular of stream densities. The initial sheet orientation is thus not an important issue for us. In the following we will consider quite general orbits, but will arbitrarily set unless otherwise stated 666Orbits starting on axes of symmetry in phase-space can show non-generic stream density behaviour for as we will discuss in the section on non-integrable potentials (section 5).. Note that the choice of does not influence the dynamical evolution of the distortion tensor as it is not part of the initial conditions for the GDE. Only the final projection to configuration-space is affected by initial sheet orientation.

From a dynamical point of view, static, integrable potentials are very simple. The motion within them can be described in terms of action/angle variables and their Hamiltonian can be expressed solely as a function of the actions . All motion in these potentials is regular, so there are no chaotic regions in their phase-space. In action-angle space the orbits lie on tori and are characterised by a fixed number of fundamental frequencies. DM particles in integrable potentials will experience only phase mixing.

Because of these simple properties Helmi & White (1999), hereafter HW, were able to develop an analytic linearised treatment based on action-angle variables to derive results for the stream density behaviour. In their paper they did not specifically focus on CDM, but rather analysed how Gaussian clouds in action-angle space evolve with time.

As an example of an integrable potential, we apply our method to several Eddington potentials (Lynden-Bell, 1962). These are separable in spherical coordinates. The third integral for this type of potential is . We chose the following example of an Eddington potential,


with , kpc, 777These values do not have any specific meaning. We have chosen them just in a convenient way., and studied an orbit which starts at kpc with a velocity of .

In Fig. 4 we show the evolution of the stream density for this orbit, i.e. the projection from phase-space to configuration-space for an initial condition with . Here and elsewhere (unless otherwise stated) we define ‘the ‘orbital period” as the radial oscillation period for the purpose of making such plots. The late-time behaviour of the stream density can be fitted by an analytic formula derived by HW:


with only one fitting parameter . Comparing this to the result for the long-term behaviour in HW (Eq. 37) it is clear that just reflects the initial phase-space distribution and the orbital parameters (the derivatives of the fundamental frequencies with respect to the actions). The results in HW are calculated for a Gaussian cloud in phase-space, not a cold sheet as in our case. Since the initial distribution only affects , the long-term behaviour of the two configurations is the same, as shown in Fig. 4.

Figure 4: The stream density evolution calculated using the GDE integrator DaMaFlow is compared to the analytic result obtained by a linearised treatment in action/angle variables (Helmi & White, 1999). The results agree essentially perfectly. Notice how well the numerical calculation matches the caustics. The upper panel clearly shows that the numerical result also has the correct behaviour at late times. The initial quasi-exponential decay is also visible.

One can clearly see that our method produces caustics at the correct positions and is also able to recover the secular evolution of the stream density, the density decrease. We note also the initial quasi-exponential stream density decrease that is often referred to as Miller’s instability (Miller, 1964). Recently Helmi & Gomez (2007), hereafter HG, showed this is a generic feature of Hamiltonian dynamics and not, as long believed, an artifact specific to N-body integrations (Hemsendorf & Merritt, 2002; Kandrup & Sideris, 2003). HG discuss the effect in detail for spherical potentials. All our tests with spherical, axisymmetric and triaxial potentials show a similar quasi-exponential initial decay. Since DaMaFlow integrates the equations of motion for a single particle in a perfectly smooth potential, it is evident that this behaviour has nothing to do with N-body effects.

A NAFF frequency analysis of particle orbits in the Eddington potential reveals, as expected, three linearly independent frequencies. It is this number that dictates the speed with which stream densities decrease in static, separable potentials. One can see this very clearly from Fig. 5. Here we show the density decrease in a simple Kepler-like toy-model:


For the orbit was started from and with an energy . This orbit has only one fundamental frequency. Changing to increases the number of frequencies to two; the loops of the orbit no longer close. The starting point for this second case has been set to and , corresponding also to . The increased number of frequencies results in a more rapid decrease in stream density: for one fundamental frequency, and for two. In this sense the long-term behaviour of streams in static, integrable systems is very simple and is determined only by the number of fundamental frequencies.

Figure 5: Stream density evolution for the normal Kepler potential () and for a modified potential (). Orbits in the Kepler potential have a single fundamental frequency. Changing the potential exponent from to increases the number of fundamental frequencies to two. As a result, the long-term stream density behaviour changes from to . For integrable potentials the long-term stream density decrease on each orbit is dictated by the number of fundamental frequencies.

5 Non-integrable potentials

Analytic methods like the HW formalism are only able to deal with integrable potentials. This is a serious limitation since realistic potentials are rarely integrable. To demonstrate that the GDE method can also deal with more complex phase-space structure we now discuss the well-known ellipsoidal logarithmic potential with a core,


analysing its stream density behaviour and its phase-space structure. There are two reasons why we have chosen this potential: first there has been substantial previous work on its phase-space structure (Binney & Spergel, 1982; Papaphilippou & Laskar, 1998), so we can compare directly with these earlier results. Second this potential is often considered as a good model for galactic halos because it reproduces a flat rotation curve and its shape can easily be tuned by two parameters that correspond to the axial ratios of the ellipsoidal isopotential surfaces. It has been used, for example, for dynamical studies of the debris streams of the Sagittarius dwarf galaxy (Helmi, 2004). The ellipsoidal logarithmic potential without a core () has also been used to study the influence of halo shape on the annual modulation signal in dark matter detectors (Evans et al., 2000).

It is known that, depending on the degree of triaxiality, the phase-space of the logarithmic potential can be occupied to a large extent by chaotic orbits (Papaphilippou & Laskar, 1998). In Fig. 6 we show how the stream density evolves along one of these chaotic orbits. This orbit was integrated in a potential with and started at ,. Here we apply the same system of units as Papaphilippou & Laskar (1998) and write all quantities as dimensionless numbers. It is evident from this plot that the system mixes very rapidly along this orbit. This is chaotic mixing, and contrasts with the phase mixing that we found before for regular motion in separable potentials. We note that chaotic orbits are difficult to handle from a numerical point of view because of the rapid spread in phase-space that characterises them. To check whether we can rely on our stream density values, we also calculated the 6-D phase-space density along this orbit. Over the full integration range shown in Fig. 6 ( orbits) it remained constant to an accuracy of .

Figure 6: Stream density evolution along a chaotic orbit in the ellipsoidal logarithmic potential with . The density decreases very rapidly, reflecting the chaotic mixing along this orbit. Note that the decrease is no longer a power-law, as in the case of regular motion.

Using DaMaFlow we have scanned the phase-space of box orbits within the logarithmic potential, identifying chaotic and regular regions by calculating the stream density evolution for about different orbits. Each orbit was integrated for a fixed time interval of using time-steps. This corresponds to about orbital periods. We chose this very long integration time in order to distinguish between chaotic and regular regions888For chaotic orbits elements of the distortion matrices can become very large. Here we are only interested in separating chaotic and regular orbits reliably. Thus we do not care about round-off errors in the calculation of these matrices..

Our results can be compared directly to previous work where the same potential was analysed using frequency shifts (Papaphilippou & Laskar, 1998). This method is based on the fact that chaotic motion, contrary to regular motion, has no stable fundamental frequencies, so that frequency estimates shift with time. By looking for such shifts one can distinguish between chaotic and regular motion. Fig. 7 shows maps of orbit type for two different sets of axial ratios . This figure can be compared directly to Fig. 6(c) and Fig. 4(b) in Papaphilippou & Laskar (1998). For the first of these calculations we adopted the following values: and . For the second case we changed the potential shape by instead taking .

The maps of Fig. 7 are constructed as follows. We start each orbit at the centre of the potential to be sure to get a box-orbit with zero angular momentum. Each individual orbit can be labelled by its initial and velocity components. The value of is then determined by the chosen value of the energy, . This is, of course, a very special point within the potential, and it turned out that our standard choice of initial stream orientation, , produces highly non-generic stream density behaviour in this case; the stream density remains constant! In order to get properly representative behaviour we therefore took when producing Fig. 7. With this set-up, we scanned the whole plane and saved the stream density of each orbit after a fixed amount of time. The greyscale in the plots corresponds to the stream density decrease after that fixed time. Black points denote a very rapid stream density decrease, thus regions of chaotic mixing. Grey regions show a slower density decrease, reflecting phase mixing and regular motion. Much of the box-orbit phase-space is chaotic for . Reducing the asphericity to results in a much larger fraction of the boxes being regular.

A comparison of Fig. 7 to Fig. 6(c) and Fig. 4(b) in Papaphilippou & Laskar (1998) shows excellent and detailed agreement. Regions of high frequency shift correspond, as expected, to those of rapid stream density decrease and chaotic mixing. Thus the GDE and the frequency shift technique work equally well for delineating regions of chaotic and normal phase mixing. We note that the Lyapunov exponent technique for identifying chaotic behaviour is closely related to the GDE (Lichtenberg & Lieberman, 1983; El-Zant, 2002), since these exponents are obtained from the eigenvalues of the 6D distortion tensor. Identifying and characterising chaotic behaviour is not the main goal of our work here, so we will not pursue this connection further in this paper.

Figure 7: Chaos maps for box orbits in the logarithmic ellipsoidal potential for two different sets of axial ratios. The upper plot is for a highly aspherical potential with , whereas the lower plot is for a rounder potential with . The greyscale indicates the stream density decrease after a fixed time interval (2000 time units or about 100 orbital periods). Black regions mark the very low final stream densities found for chaotic orbits, while grey regions mark the higher stream densities found for regular orbits. Densities decay quasi-exponentially in the former case, but only as a power law of time in the latter. These plots can be directly compared to Fig. 6(c) and Fig. 4(b) in Papaphilippou & Laskar (1998), where a frequency shift analysis of the same system reveals exactly the same structures.

So far we have classified orbits as either regular or chaotic, but the regular part of phase-space has substructure in the form of resonances (Carpintero & Aguilar, 1998; Wachlin & Ferraz-Mello, 1998; Merritt & Valluri, 1999). These are regions where the frequencies of motion are commensurate where the are integers and the are the three frequencies of the regular motion. As shown above, resonance influence the stream density behaviour over long timescales (Helmi & White, 1999; Siegal-Gaskins & Valluri, 2007). This is because they restrict the motion to a lower dimensional region in phase-space, implying that the orbit does not fill its KAM (Kolmogorov-Arnold-Moser) torus densely. In simple terms, the system cannot spread as fast as it would do in the non-resonant case because it occupies a space of lower dimension.

Figure 8: Stream density evolution on resonant orbits. Stream density drops at a rate which depends on the number of independent orbital frequencies. Non-resonant orbits have three independent frequencies, and their stream density decreases like at late times. Resonances reduce the number of independent frequencies and slow the decrease of stream densities. With one resonance there are two independent frequencies; the stream density then falls as the inverse square of time. With two resonances, only one independent frequency remains and density falls as the inverse first power of time. The number of resonances also strongly affects the orbit shape in configuration-space. As is visible for the examples in the lower plot, non-resonant orbits (red) fill a 3D volume, whereas orbits with two resonances (green) are restricted to a line.

Fig. 8 shows the stream density evolution along three different box orbits for . The initial conditions for these orbits were chosen so that they have different numbers of orbital resonances (non-resonant, one resonance, two resonances): , , , and , . It is clear that the number of resonances has a major effect on the final stream density decrease over timescales similar to those shown in this plot. The difference in stream density between the non-resonant case and the case with two resonances is about 3 orders of magnitude after 350 orbits! Resonances also have a strong influence on the shape of the orbit in configuration-space, as shown in the lower panel of Fig. 8. Two resonances restrict the orbit to a line. From the shape of the orbits it is evident why the stream density changes so much with the resonances. The particles cannot spread over a large region if they are confined to a space of small dimension.

Figure 9: Resonance structure of a logarithmic potential as revealed by the GDE. Green indicates regions with no resonance (i.e. three independent orbital frequencies). Stream densities for these orbits decrease as at late times. Red indicates orbits with one resonance, for which stream densities decreases as . Orbits with two resonances are coloured blue. This map was constructed by integrating orbits, each for about orbital periods, using DaMaFlow in parallel. Each orbit required time-steps using a KDK leapfrog algorithm.

The regular region of phase-space for the logarithmic potential is occupied by resonant and non-resonant regions. We used DaMaFlow to scan the regular region (as for the chaos maps above) and fitted the stream density decrease by a power law in time. A non-resonant motion then gives an exponent of , whereas regions with two resonances should give . We binned these power law exponents (bin size 0.1) and coloured them. The result of this procedure is shown in Fig. 9. For this map we integrated a total orbits for about orbital periods. Each integration required KDK leapfrog time-steps. The chaotic regions are shown in grey, while the regular regions are shown in three different colours depending on their stream density behaviour. We note that the chaos detection here was carried out by imposing a threshold on the stream density decrease. Every orbit that crosses this threshold during its evolution is considered as chaotic and marked as grey in the map. This is the reason why the chaotic pattern here is not identical to that in Fig. 7.

Most of the regular phase-space in this figure is occupied by non-resonant orbits shown in green. Superposed on these is a fine network of resonance lines shown in red. In blue regions the stream density decreases linearly with time because there are two resonances. Note how well our method locates the resonance lines in the initial condition space spanned by and . Papaphilippou & Laskar (1998) analysed resonances with frequency maps by plotting the rotation numbers defined as and , where are the fundamental frequencies along the long (L), short (S) and middle (M) axes. It is straightforward to identify the resonance lines in Fig. 9 with those in the frequency map of Papaphilippou & Laskar (1998) by applying a NAFF frequency analysis to the orbit corresponding to any specific initial condition, for example, one on a given resonance line. It turns out that we can identify all resonance lines in Fig. 9 with similar lines in the frequency map. At the intersection of these lines we have periodic orbits (satisfying two resonance conditions) with the same rotation numbers as those found in Papaphilippou & Laskar (1998). For example, the line going from the upper left corner to the lower right corresponds to the resonance, meaning that .

We conclude that our method can resolve the structure of phase-space equally as well as the standard frequency mapping technique.

6 Triaxial Dark Matter halos

Figure 10: Isopotentials for the outer and inner parts of one of our triaxial NFW halos. It is obvious that the halo becomes rounder as one moves outwards. In this case the transition scale was chosen to be equal to the scale radius of the NFW profile.

Figure 11: Comparison of the radial variation of isopotential axial ratios for N-body halos (Hayashi et al., 2007) to that for our simple triaxial NFW model. The N-body values are the average of those found in Hayashi et al. (2007). The transition scale of our model has been set equal to the scale radius of the underlying NFW profile.

In CDM cosmologies dark matter halos are not spherical. Furthermore, simulations suggest that their shape should vary with radius, both equidensity and equipotential surfaces being rounder (on average) at larger radii. Several studies have tried to constrain the shape of the Milky Way’s halo by analysing the properties of observed tidal streams like that of the Sagittarius dwarf galaxy (Ibata et al., 2001; Helmi, 2004). Recently, Hayashi et al. (2007) analysed the radial variation in potential shape of simulated halos that might correspond to that of the Milky Way. Although there is substantial object-to-object scatter, on average they found a relatively rapid transition from aspherical to almost spherical which occurs near the scale radius of the best fitting NFW profile. They provide a simple fitting formula for this mean behaviour,


for the principal axial ratios and . (Note that they actually provide two different sets of fitting parameters for Eq. (34) depending on the principal axial ratios.) They also propose a modified NFW potential that takes into account the variation in shape, but this potential is not very convenient because it is not straightforward to derive the corresponding equations of motion. The examples given above show that potential shape can have a substantial effect on stream density evolution, so it is interesting to see how strong such effects can be in a realistic model.

To analyse this we have built a simple extension of the NFW model that qualitatively reproduces the shape variation found by Hayashi et al. (2007) but which has simple equations of motion that can easily be implemented in DaMaFlow. (For another similar model, see Adams et al. (2007).)

We model the variable shape of the NFW halo by replacing the euclidean radius in the formula for the potential of a spherical NFW halo by a more general “radius” given by: