N-body simulations of collisionless collapse have offered important clues to
the construction of realistic stellar dynamical models of elliptical galaxies.
Such simulations confirm and quantify the qualitative expectation that rapid
collapse of a self-gravitating collisionless system, initially cool and
significantly far from equilibrium, leads to incomplete relaxation. In this
paper we revisit the problem, by comparing the detailed properties of a family
of distribution functions derived from statistical mechanics arguments to those
of the products of collisionless collapse found in N-body simulations.Comment: 4 pages, 2 figures, to appear in the proceedings of the 39th
  Rencontres de Moriond, "Exploring the Universe," La Thuile, Italy, March 28 -
  April 4, 200

Bertin, G.

Trenti, M.

English

arXiv

N-body simulations of collisionless collapse have offered important clues to the construction of realistic stellar dynamical models of elliptical galaxies. Such simulations confirm and quantify the qualitative expectation that rapid collapse of a self-gravitating collisionless system, initially cool and significantly far from equilibrium, leads to incomplete relaxation. In this paper we revisit the problem, by comparing the detailed properties of a family of distribution functions derived from statistical mechanics arguments to those of the products of collisionless collapse found in N-body simulations

M. Trenti

G. Bertin

AIR Universita degli studi di Milano

ar
X
iv
:a
str
o-
ph
/0
40
62
36
v1
  9
 Ju
n 
20
04
DYNAMICAL MODELS AND NUMERICAL SIMULATIONS OF
INCOMPLETE VIOLENT RELAXATION
M. TRENTI1 & G. BERTIN2
(1) Scuola Normale Superiore, piazza dei Cavalieri 7, I-56126 Pisa, Italy;
(2) Dipartimento di Fisica, Universita` di Milano, via Celoria 16, I-20133 Milano, Italy
N-body simulations of collisionless collapse have offered important clues to the construction of
realistic stellar dynamical models of elliptical galaxies. Such simulations confirm and quantify
the qualitative expectation that rapid collapse of a self-gravitating collisionless system, initially
cool and significantly far from equilibrium, leads to incomplete relaxation. In this paper we
revisit the problem, by comparing the detailed properties of a family of distribution functions
derived from statistical mechanics arguments to those of the products of collisionless collapse
found in N-body simulations.
1 Introduction
The collapse of a dynamically cold cloud of stars can lead to the formation of realistic stellar
systems, with projected density profiles well represented by the R1/4 law 13. The theoretical
framework for the mechanism of incomplete violent relaxation that governs this scenario of
structure formation was proposed by Lynden-Bell 9, who argued that fast fluctuations of the
potential during collapse would lead to the formation of a well-relaxed isotropic core, embedded
in a radially anisotropic, partially relaxed halo. This general picture served as a physical justi-
fication for the construction of the so-called f∞ models
3, which indeed recovered the R1/4 law
and, suitably extended to the case of two-component systems (to account for the coexistence of
luminous and dark matter), led to a number of interesting applications to the observations 4.
The application of statistical mechanics to this formation scenario 10 also led to the deriva-
tion of a separate family of distribution functions, the f (ν) models, which was recently shown
to possess interesting thermodynamic properties in the context of the gravothermal catastro-
phe 5. The key ingredient for the construction of the f (ν) models is the conjecture that a third
quantity Q (defined as Q =
∫
Jν |E|−3ν/4fd3qd3p), in addition to the total mass M and the
total energy Etot, is approximately conserved during the process of collisionless collapse. This
quantity is introduced to model the process of incomplete violent relaxation, ensuring a radially
biased pressure tensor and a 1/r4 density profile in the outer parts of the system. A preliminary
inspection of the general characteristics of the f (ν) models then convinced us that, with signif-
icant advantage over the f∞ models, they might also serve as a good framework to interpret
the results of simulations of collisionless collapse not only qualitatively, but also in quantitative
detail, which is the main point of this paper.
2 Models
The extremization of the Boltzmann entropy under the three constraints described in the pre-
vious section leads to the distribution function f (ν) = A exp [−aE − d(J2/|E|3/2)ν/2], where
a, A, d, and ν are positive real constants; here E (E < 0) and J denote single-star specific
energy and angular momentum. At fixed value of ν, one may think of these constants as pro-
viding two dimensional scales (for example, M and Q) and one dimensionless parameter, such
as γ = ad2/ν/(4piGA). In the following we will focus on values of ν ranging from 3/8 to 1. The
corresponding models are constructed by solving the Poisson equation for the self-consistent
mean potential Φ(r) generated by the density distribution associated with f (ν). At fixed value
of ν, the models thus generated make a one-parameter family of equilibria, described by the con-
centration parameter Ψ = −aΦ(0), the dimensionless depth of the central potential well. The
projected density profile of concentrated models is well fitted by the R1/n law, with n close to
4, and the models provide reasonable fits to the surface brightness and to the kinematic profiles
of bright elliptical galaxies 12.
3 Simulations of collisionless collapse
In principle, to study the process of collisionless collapse we have two options: either a tree
code1 or a particle-mesh13 algorithm. We are interested in the large scale structure of the end-
products of collisionless collapse, for systems that do not exhibit large deviations from spherical
symmetry. The natural choice thus appears to be that of a particle-mesh code, based on a
spherical grid and an expansion in spherical harmonics. The code used in the present study is
thus a new version 11 of the van Albada code. For completeness, we have also run a number
of comparison simulations with the fast code developed by Dehnen 7, which confirmed that our
results do not depend on the numerical scheme employed.
During the process of collisionless collapse, violent relaxation is expected to wipe out much of
the details that characterize the initial conditions, so that the end-products are expected to have
properties that depend only or mostly on the initial value of the virial ratio u = (2K/|W |)t=0,
which sets the relevant collapse factor. In reality, violent relaxation is incomplete. Therefore,
the final state is that of an approximate dynamical equilibrium characterized by an anisotropic
distribution function, different from a Maxwellian (which would correspond to thermodynamic
equilibrium). Because of such incomplete relaxation, the end-products of the simulations do
conserve some memory of the initial state.
Earlier investigations13,8 compared “clumpy” to “homogeneous” initial conditions, showing
that clumpy initial conditions lead to end-states with projected density distributions well fitted
by the R1/4 law. We argue that, if the degree of symmetry in the initial conditions is excessive,
little room is left for the mechanism of violent relaxation; this is confirmed by the fact that little
or no mixing is observed in the single-particle angular momentum for homogeneous simulations,
as reported in Fig. 1. We have thus chosen to focus this paper on clumpy initial conditions, for
which an efficient mixing in phase space is guaranteed.
Figure 1: Single-particle angular momentum scatter (correlation between the initial and final values of J) for
a typical clumpy simulation (left panel) and for its symmetrized (homogeneous) version (right panel). The
symmetrization process in the initial state is performed by accepting the radius and the magnitude of the velocity
of each particle, following the procedure for initial clumpy conditions, but by redistributing uniformly the angular
variables.
We performed several runs varying the number of clumps and the initial virial ratio u in
the range from 0.05 to 0.25. In most cases the clumps are cold, i.e. their kinetic energy is
all associated with the motion of their centers of mass. The simulations are run for several
dynamical times beyond the violent collapse phase, ensuring that the system has settled into a
quasi-equilibrium. The final configurations are quasi-spherical, with shapes that resemble those
of E2 − E3 galaxies. The final equilibrium half-mass radius rM is basically independent of the
value of u. The central concentration achieved correlates with u, as can be inferred from the
conservation of maximum density in phase space 8.
4 Fits and phase-space properties
In order to study the output of the simulations we compare the density and the anisotropy
profiles, ρ(r) and α(r) (here the local anisotropy is defined as α(r) = 2− (〈p2θ〉+ 〈p
2
φ〉)/〈p
2
r〉), of
the end-products with the theoretical profiles of the f (ν) family of models. Smooth simulation
profiles are obtained by averaging over time in the last few dynamical times of the simulation.
For the fitting models, the parameter space explored is that of an equally spaced grid in (ν,Ψ),
with a subdivision of 1/8 in ν, from 3/8 to 1, and of 0.2 in Ψ, from 0.2 to 14.0; the mass and
the half-mass radius of the models are fixed by the scales set by the simulations. A minimum
χ2 analysis is performed as described elsewhere 6.
As illustrated in Fig. 2, the density of the simulations is well represented by the best-fit f (ν)
profile over the entire radial range. The fit is satisfactory not only in the outer parts, where
the density falls by nine orders of magnitude with respect to the center, but also in the inner
regions. Depending on the initial virial parameter u, the end-products possess a density profile
which, projected along the line-of-sight, may exhibit an R1/n behavior with different values of
n and yet it is well fitted by the f (ν) models. In turn, this may be interpreted in the framework
of the proposed weak homology of elliptical galaxies 2.
To some extent, the final anisotropy profiles for clumpy initial conditions are found to be
sensitive to the detailed choice of initialization, in particular to the number and the size of the
clumps. The agreement between end-products of the collapse and models (see Fig. 2) seems to
be best for simulations with 10 to 20 clumps, for a clump size such that the sum of their volumes
fills the sphere where their centers of mass lie. We should recall that both the low and the high
number of clumps limits are those of an initial homogeneous condition, unfavorable to violent
relaxation.
At the level of phase space, we have performed two types of comparison, one involving the
-1 0 1 2
0
0.5
1
1.5
2
-1 0 1 2
-6
-4
-2
0
2
4
-60 -40 -20 0
0
0.01
0.02
0.03
0.04
E
-1 0 1 2
-1
0
1
10−3
10−2
10−1
−60
−40
−20
0
0.5
1
J2E 10
−3
10−2
10−1
−60
−40
−20
0
0.5
1
J2E
10−3 10−2 10−1
−70
−60
−50
−40
−30
−20
−10
0
J2
E
10−3 10−2 10−1
−70
−60
−50
−40
−30
−20
−10
0
J2
E
Figure 2: Comparison between the results of one simulation with u = 0.23, called C2 4 (starting from 20 cold
clumps), and the best-fit f (ν) model (ν = 5/8, Ψ = 5.4). Left set of four panels. Density as measured from the
simulation (error bars) and model profile (top left). Residuals from the best-fit profile (top right). Anisotropy
profile of the simulation (error bars) and profile (bottom left). Energy density distribution N(E) (bottom right).
Right set of four panels. Final phase space density N(E, J2) (left column), compared with that of the best fitting
f (ν) model (right column). The model curve for N(E) and surface for N(E, J2) have been computed by a Monte
Carlo sampling of phase space.
energy density distribution N(E) and the other based on N(E, J2). The chosen normalization
factors are such that: M =
∫
N(E)dE =
∫
N(E, J2)dEdJ2. In Fig. 2 we plot the final energy
density distribution for a simulation run called C2 4 with respect to the predictions of the best-
fit model identified from the study of the density and pressure anisotropy distributions. The
agreement is good, especially for strongly bound particles. At the deeper level of N(E, J2),
simulations and models also agree quite well, as illustrated in the right set of four panels of
Fig. 2. For the case shown, the distribution contour lines basically match in the range from
Emin to E ≈ −30 (in the relevant units); however, the theoretical model shows a peak located
near the origin, not observed in the simulations.
References
1. Barnes, J., Hut, P., Nature 324, 466 (1986)
2. Bertin, G., Ciotti, L., Del Principe, M., Astron. Astrophys. 386, 149 (2002)
3. Bertin, G., Stiavelli, M., Astron. Astrophys. 137, 26 (1984)
4. Bertin, G., Stiavelli, M., Rep. Prog. Phys. 56, 493 (1993)
5. Bertin, G., Trenti, M., Astrophys. J. 584, 729 (2003)
6. Bertin, G., Trenti, M., AIP Conf. Proc. 703, 318 (2004)
7. Dehnen, W., Astrophys. J. Letters 536, L39 (2000)
8. Londrillo, P., Messina, A., Stiavelli, M., Mon. Not. Roy. Astron. Soc. 250, 54 (1991)
9. Lynden-Bell, D., Mon. Not. Roy. Astron. Soc. 136, 101 (1967)
10. Stiavelli, M., Bertin, G., Mon. Not. Roy. Astron. Soc. 229, 61 (1987)
11. Trenti, M., Ph.D. Thesis, Scuola Normale Superiore, Pisa, in preparation
12. Trenti, M., Bertin, G., submitted (2004)
13. van Albada, T.S., Mon. Not. Roy. Astron. Soc. 201, 939 (1982)


Dynamical models and numerical simulations of incomplete violent relaxation

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.315.1178

Dynamical models and numerical simulations of incomplete violent relaxation

Abstract

Similar works

Full text

Available Versions

AIR Universita degli studi di Milano