Previous ToC Up Next

4. Populating Velocity Space

4.1. Choosing a Velocity

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

per unit mass of the test particle, we can equate that to the kinetic energy, also per unit mass, of a particle moving at the escape velocity:

This is the maximum velocity allowed at radius r, and we also know that the minimum velocity at radius r is zero. The question is: what is the probability distribution for .

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.

4.2. A Meta-Recipe

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:

Alice: Where did you get this expression from?

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:

for negative energy, in other words for bound particles, while for , otherwise you would get the whole universe filled with escapers.

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.

4.3. Following the Recipe

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

where the energy per unit mass can be written in terms of the escape velocity as . If we introduce the variable

we can write , for a given fixed . The distribution function for then becomes, in terms of , proportional to the following function:

Alice: And the range of admissible q values is . 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 .

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.

4.4. A Rejection Technique

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

To find the extrema for , we set the derivative to zero, and solve for :

Indeed: 0.1 is a rather tight upper limit.

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

Okay, I understand the procedure now, and it looks correct.

4.5. Distribution Function

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:

In the case of spherical symmetry and isotropy, we have where the energy is the sum of the kinetic and potential energy per unit mass of the particles that reside in the volume element .

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!

4.6. The Density as a Projection

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:


So we get for the density:

Here both and are negative. It is easier to introduce quantities that are positive. Following the notation used by Binney and Tremaine, we can define:

We can write again explicitly how the density in configuration space is a projection down from phase space, through an integration over :

Strictly speaking I should have changed the symbol for again since the functional form is different: .

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.

4.7. A Change of Variables

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.

4.8. Deprojecting the Density

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

Bob: Why would they divide by that funny factor?

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 :

Bob: Ah, and here they do their deprojecting: they claim that this equation can be inverted. Eq. (51) gives you , when you give it . What we want is to obtain , and we can certainly figure out how to compute . I see the logic now.

Alice: Yes, and they give the solution of inverting Eq. (51), as their eq. (4-140a):

and in alternative form as their (4-140b):

Bob: Great! That will make it finally possible for us to determine the distribution function: with a little luck we can solve that integral.

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