4. Populating Velocity Space
Alice: Now that we completely understand how you choose the radial distance
for a particle, there is only one thing left to do: to choose the magnitude of
its velocity. In your code you had a more complicated construction than the
one-liner you used for the position. For the velocity you wrote:
You are throwing dice a number of times, until some criterion is satisfied,
and then you move on to a new one-liner that gives you the value for the
variable velocity, which is the magnitude of the velocity.
Bob: Here is the recipe for that part. Aarseth et al. start
with the observation that the maximum velocity allowed at a radius r
is the escape velocity
Alice: I see that you are using again the choice of units given by
Aarseth et al., where
Bob: Yes, they are by far the most convenient, they save time when
writing down and manipulating the equations, and they make it also less
likely to make mistakes.
Alice: I'm not sure about the last part of what you said. The advantage
of keeping the full physical quantities is that you always have a few extra
checks you can make, at the end: if the physical dimensions of length, mass
and time are not exactly the same, at the left and right hand side of an
equation, the equation must be wrong. If you work only with dimensionless
quantities, you loose that advantage.
Bob: I don't consider it an advantage to have to do so much more work that
you are likely to make more mistakes, so that you can then happily catch them.
But clearly we are talking about matters of taste, and we have already decided
we will present our results either way, now that we know exactly how
to transform in both directions, between physical and dimensionless variables.
Alice: Yes, that is what we decided. How did you find the probability
distribution for the velocities?
Bob: I started with the distribution function for the energy of the
particles:
Bob: It's just what it is, for Plummer's model. I found it in a useful
table in The Gravitational Million-Body Problem, by Douglas Heggie
and Piet Hut [Cambridge University Press, 2003]. It is table 8.1 on p. 73,
a whole page full of useful properties of Plummer's model. By the
way,
Alice: It is a power law in energy, and it looks like a polytrope.
Ah, yes, now I remember: Plummer's model is nothing else but a polytrope
of index five. Polytropes are defined in general for index
But it is not fair to use such remembered knowledge, or equations that
you pluck from a book. Our whole approach is aimed at spelling out
everything, both to help us in our research, to see new aspects we
might otherwise have overlooked, and to help us in our teaching, to
make things crystal clear to the students.
Bob: How do you expect to get new insight into the fact that the distribution
function of Plummer's model is a seven-halves power of the energy? That fact
will not change, no matter how long you stare at it.
Alice: That's not the point. Once we spell out in complete detail how
you derive and verify all aspects of one recipe for one star cluster model,
you can then follow that example approach to construct any other recipe for
any other cluster model. In other words, we are using Plummer's model as
a case study in order to present a meta-recipe for constructing recipes for
constructing models for star clusters.
Bob: I had no idea we were doing something that fancy. But whatever words
you want to hang on it, I cannot deny that it is a good thing to check things
for yourself, and most importantly, I'll have to explain at least some
of these things in class pretty soon, so okay, let's go through the derivation.
However, we would probably lose the thread of our argument if we would
go into deriving eq. (37) right now. Let's postpone
that a bit, and first see whether we can reconstruct what I have written
in my program. As you can expect, here too I just followed the recipe from
Aarseth et al.. Let us first see whether we can at least derive that
recipe, assuming the validity of eq. (37).
Alice: Sounds like a good plan.
Bob: Here is what I have understood, so far, of the recipe. Given the
distribution function for the energy, you have to transform that into an
equation for the magnitude for the velocity. What makes life simpler, is
the fact that we are comparing particles with different velocity choices at
a given point, so we know that their potential energies are all the same.
This means that the probability
So I guess the next step is to invert
Bob: To say the least. Therefore, for the velocities, they
choose a different approach. If you plot the function
You can imagine that you can obtain a distribution of the required
weighting by throwing darts at that graph. If you hit a point
somewhere above the graph, you throw a new dart, and you keep throwing
new darts until you hit a point below the graph, anywhere between
Alice: That is a clever solution. It is called a rejection technique.
Didn't John von Neumann first apply that? You start by allowing more
solutions than the minimal set of correct ones, and then you weed out
the incorrect ones, by rejecting them.
Bob: Indeed. And to make the procedure efficient, you don't want to
throw darts way above the graph, so you limit yourself to the maximum
value that the graph attains in the interval of interest, or perhaps
a slightly higher value. The authors of the paper choose a value of 0.1.
Alice: Is that a safe value? Let's check for ourselves. The derivative
of
Bob: You can now see what I did when I wrote:
Alice: Ah, yes. So x stands for
Bob: The only thing left to do now is to derive the form of the
distribution function.
Alice: Yes. Let us see how far we can get. Clearly we need a little help
from our friends: here is the classic reference Galactic Dynamics,
by James Binney and Scott Tremaine [Princeton University Press, 1987].
I have found it to be a very useful book, whenever I had to look up some
properties of particular models in stellar dynamics. It also has helped
me a number of times to refresh my memory about Jeans equations, the
tensor virial theorem and those sort of things.
Bob: I see that your browsing was rather uneven: there is a piece, about
one third along the way, which has a gray strip. Let me open it up there.
Aha, chapter 4: Equilibria of Collisionless Systems. How come those pages
are so well-used? I thought you were mainly interested in collisional
systems?
Alice: I didn't realize my copy of Binney and Tremaine betrayed my past
browsing so well. You're right, that's the chapter I tend to consult most.
And precisely because we are interested in collisional systems, we have
to find a way to start our simulations with a collisionless system.
In other words, we run a numerical simulation of a collisional N-body
system in order to see the effects of two-body relaxation and all that.
But we need to have a well-defined starting point. A formal way to define
this is to say that we ignore collisions from time
I'm sure you don't like such a formal way of speaking about it. But
you can't complain: you started asking philosophical questions about
well-read chapters in someone else's book! So there.
Bob: Okay, let's go get our distribution function. Does this book
tell us how to do that?
Alice: Yes, but let us first see how far we can get under our own steam.
Let us go back all the way to how we got started, namely with the
Plummer potential:
Let us define
In general,
Each star moving with velocity
Now you see why I started with that funny tilde over the distribution
function: since we will be dealing primarily with
Bob: I still would be happy to not use tildes at all, but we'll each
follow our own way.
Alice: There is one more thing we definitely have to stress here.
Even though we have switched to
This is something students always get confused about. It is tempting
to write
Bob: Noted!
Alice: The most fundamental stage for a star system in stellar dynamics is
phase space. Because our physical eyes see what happens in
configuration space, we think about the density as a rather basic function.
But really, what we call density is already a type of shadow: it is a
projection down to 3D from the distribution function in 6D.
In a star cluster with spherical symmetry in configuration space, we
can show this projection effect as follows. The mass density in stars
at a distance
Also, we know that the distribution function has to be zero for positive
energies, since unbound particles will escape, and their density in an
(almost) infinite universe will become (almost) zero:
For a given radial distance
This means that we can rewrite eq. (46) as:
Now
So we get for the density:
Bob: I'm glad you didn't!
Alice: I decided to take my hat off for you, or at least f's hat.
Even for me, too much formality gets bothersome.
Bob: Looking over your shoulder, I see that your eq. (49)
is exactly eq. (4-137) in Binney and Tremaine, if you pull out their
Alice: Yes. Let us follow their next few steps, for our particular
Plummer case. They remark that
Bob: That makes sense. If
Alice: Right. But here everything is monotonic: the density keeps dropping
off when you move away from the center, the radius keeps increasing,
and the potential keeps increasing as well, or equivalently,
Bob: So how do we get from
And the conversion factor is nothing less than the inverse of the
gravitational acceleration, which I derived early on, with the
prophetic remark that it would come in handy later on.
Alice: No.
Bob: No?
Alice: No, there is no conversion factor. Think about what density
means, in this context. The mass density
This is the same thing as what we had seen with the phase space
distribution function
Bob: You're right! How tricky. Or how stupid of me.
Alice: Let's call it tricky. Believe me, I've made these mistakes
often enough myself.
Bob: So eq. (49) becomes simply:
Alice: Indeed. If we divide both sides of the equation by
Alice: In order to get a simple looking right-hand side in their next
equation (4-139), which they obtain by differentiating eq. (50)
with respect to
Alice: Yes, and they give the solution of inverting Eq.
(51), as their eq. (4-140a):
Alice: Not so quick: how do we know that they inverted eq.
(51) correctly?
Bob: There you go again! You don't take anything on authority, do you?
Alice: Well, no, preferably not. The whole reason to embark on this
exercise is to prove things for ourselves, from scratch. I won't call
it starting from scratch if we simply accept their inversion. We
might as well have accepted their expression for the distribution function
in the first place.
Bob: Which would have been fine with me, to be honest. But okay, we got
started, we may as well finish. Though I really begin to doubt that our
students will have the stomach to go through all of this.
Alice: The best ones will. And those are the ones we'd like to stay in
touch with. So this may a good way to find good students, to see who is
interested in figuring this out, all the way to the bottom.
Bob: I'm not so sure. You might just get formal and pedantic students.
I prefer to work with students who get codes to run and who get
results with codes, rather than students who can solve this inversion
equation -- what do they call it? -- an Abel integral equation.
Alice: The very best students will be able to do both.
Bob: As able as Abel.
Alice: Such a silly pun would work better if you had an accent. Come,
let's take our usual approach: first let's assume that
eqs. (52) and (53) are indeed correct, and let
us check whether we can then derive the correct distribution function
for Plummer's model. Having done that, we will come back, and check
for ourselves whether we can prove this Abel integral equation
inversion.
Bob: Something tells me it will be a long day.
Alice: Perhaps. Whatever it takes! But I agree with you, this may
well take longer than we thought. Let's continue tomorrow.
Bob: That sounds better already. See you then! 4.1. Choosing a Velocity
, itself a function that
depends on the radius. The escape velocity can be determined by
requiring that a particle at radius r has exactly zero total energy,
i.e. its kinetic energy is just enough for a parabolic escape
to infinity. Since the potential energy for a test particle moving in
Plummer's model is given as
.
.
4.2. A Meta-Recipe
here is the mass of a single star, assuming that all
stars have equal mass. If you multiply both sides of the equation with
, you get 
, the mass density of stars in phase
space.
through a distribution function:
for 
, otherwise you would get
the whole universe filled with escapers.
4.3. Following the Recipe
to have an absolute
value for the velocity 
at radial position
is given by
can be written in terms of the escape velocity 
as
. If we introduce the
variable
, for a given fixed
. The distribution function for 
then
becomes, in terms of 
, proportional to the following function:
.
This looks exactly like the problem we had for determining the radial
positions. There we knew the density, i.e. the probability function
to find a particle at a given position. By integrating the density we
obtained the cumulative mass function 
, and then we
inverted that to obtain 
.
. However, that
doesn't look so easy.
, then the height of that curve, for each
value, gives you the relative probability that
would lie in a region of small fixed width around that
value.
and 
. If you follow that procedure,
you are automatically guaranteed that you score more hits at places
where the graph is higher, and exactly so in proportion to the height
of a graph.
4.4. A Rejection Technique
is
, we set the derivative to zero,
and solve for 
:
and y stands for
. You keep throwing darts until you find a y
value under the graph. That gives you the corresponding x value.
Since this value is equal to 
, you have to
multiply x with the escape velocity 
, which we found
in eq. (36) to be:
4.5. Distribution Function
all the way to 
. Then, at 
, we switch on
the effects of close encounters, and through our simulations we can see
in what way our star clusters then begin to deviate from the dynamical
equilibrium we started with.
as the density of stars
in phase space. This means that in a six-dimensional volume element
, you will find a number of particles equal to
. We will restrict
ourselves to equal-mass systems, where each star has a mass
. This means that the volume element 
contains on average an amount of mass equal to
.
will be time dependent, i.e.
of the form 
, but here we will restrict
ourselves to dynamical equilibrium situations, where 
is constant, at least on a dynamical time scale, on the order of the crossing
time, the time it will take for a typical particle to cross the system.
at radial position
has a kinetic energy 
and a
potential energy 
. It is most convenient to talk
about specific energies, namely the energy per unit mass, 
,
which is given as:
where the energy
is the sum of the kinetic and potential energy
per unit mass of the particles that reside in the volume element
.
,
I preferred to use the 
notation for a distribution
function with an 
dependence, leaving the 
notation for the dependence of Cartesian phase space coordinates.
dependence, we still
have to write 
, in order to find the number
of stars in the volume element 
.
instead, but that would be wrong: the
element 
would include all particles with specific
energy between 
and 
, globally in the
system, which will amount to far more mass than is present in the
local volume element 
around the position
in phase space.
4.6. The Density as a Projection
from the center is
. In that
case only the magnitude 
can come into play.
We can thus substitute:
, the escape velocity is given
by definition as the velocity for which a particle can just reach infinity,
which means that the total energy is zero for that particle, and therefore
also the energy per unit mass. As you already showed right at the beginning
of our discussion, this implies:
and 
are negative.
It is easier to introduce quantities that are positive. Following the
notation used by Binney and Tremaine, we can define:
:
again since the functional form is different: 
.
4.7. A Change of Variables
in front of the integral.
is a monotonic function
of 
, and therefore it is possible to use 
as the independent variable in the expression for 
.
would not be monotonic,
there would be some values 
and 
with
for which 
.
In that case, if you would write 
, you would not
know whether you would be referring to 
or
.
keeps decreasing.
to 
?
That can't be difficult. We have
is the
amount of mass in stars within the volume 
, and not
within the whole shell of the cluster between 
and
. The latter amount of mass would be
.
: the physical meaning does not change,
where we write it as 
or as 
.
In both cases we multiply 
with the volume element
.
4.8. Deprojecting the Density
, we get Binney and Tremaine's eq. (4-138):
:
,
when you give it 
. What we want is to obtain
, and we can certainly figure out how to compute
. I see the logic now.
