Previous | ToC | Up | Next |

** 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:

x = 0.0 y = 0.1 while y > x*x*(1.0-x*x)**3.5 x = frand(0,1) y = frand(0,0.1) end velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25)

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 , 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

** 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, 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.

** 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
through a distribution function:

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 to have an absolute value for the velocity at radial position is given by

So I guess the next step is to invert . However, that doesn't look so easy.

** Bob**: To say the least. Therefore, for the velocities, they
choose a different approach. If you plot the function
, 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.

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 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.

** 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 is

** Bob**: You can now see what I did when I wrote:

x = 0.0 y = 0.1 while y > x*x*(1.0-x*x)**3.5 x = frand(0,1) y = frand(0,0.1) end velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25)

** Alice**: Ah, yes. So ` x` stands 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:

** 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 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.

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 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 .

In general, 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.

Each star moving with velocity 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:

Now you see why I started with that funny tilde over the distribution function: since we will be dealing primarily with , I preferred to use the notation for a distribution function with an dependence, leaving the notation for the dependence of Cartesian phase space coordinates.

** 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 dependence, we still
have to write , in order to find the number
of stars in the volume element .

This is something students always get confused about. It is tempting to write 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.

** 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 from the center is

where the integral is over all of velocity space. If in addition to spherical symmetry, the distribution function is isotropic, we know that there is no direction dependence in . In that case only the magnitude can come into play. We can thus substitute:

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 , 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:

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
in front of the integral.

** Alice**: Yes. Let us follow their next few steps, for our particular
Plummer case. They remark that is a monotonic function
of , and therefore it is possible to use
as the independent variable in the expression for .

** Bob**: That makes sense. If 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
.

** 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,
keeps decreasing.

** Bob**: So how do we get from to ?
That can't be difficult. We have

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 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
.

This is the same thing as what we had seen with the phase space distribution function : the physical meaning does not change, where we write it as or as . In both cases we multiply with the volume element .

** 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
, we get Binney and Tremaine's eq. (4-138):

** 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!

Previous | ToC | Up | Next |