The Plummer's sphere is an model for the mass density in a globular cluster of stars. For an $N$-body simulation I have initialized the position of $N$ masses with a Monte-Carlo technique but cannot find a way of initializing the velocity initial conditions.

Is there a simple function that given a position in a Plummers sphere assigns a velocity to a given mass? Lots of sites list the velocity for a circular orbit but is this a good approximation to a globular cluster and how should it be treated off the $x$-$y$ plane?


The Plummer model has a potential of the form

$$ \Phi(r)=-\frac{1}{\sqrt{r^2+1}} $$ (obviously ignoring all constants). Equating the above with the kinetic energy, you get $$ \frac12v_e^2+\Phi(r)=0\to v_e(r)=\sqrt{2}\left(r+1\right)^{-1/4} $$ This velocity is the maximum velocity you can have at a radius $r$, so we must have that $0\leq v\leq v_{e}$. In order to get the velocity, we need to know its distribution function; this can be found in this book. It reduces down to $$ g(v)dv\propto\left(-E\right)^{7/2}v^2dv $$ where $E=-v_e^2+\frac12v^2$. If we let $q=v/v_e$, then the above becomes $$ g(q)=\left(1-q\right)^{7/2}q^2 $$ with $0\leq q\leq1$. It then becomes a matter of using either rejection sampling (which I generally do not recommend) or the Newton root-finding method (which I usually do recommend) to find a fit. Then you just randomize the angles $\theta,\,\phi$ and get your x,y,z velocities:

vel = x*sqrt(2.0)*(1.0+r(i)*r(i))**(-0.25)
theta = acos(rand(-1.0,1.0))
phi = rand(0.0, 2.0*pi)
vx(i) = vel*sin(theta)*cos(phi)
vy(i) = vel*sin(theta)*sin(phi)
vz(i) = vel*cos(theta)

You can find more information in this great pdf book.

  • 1
    $\begingroup$ +1 for the nice summary. I'll mention that an optimally chosen range will get rejection sampling to work a bit over 46% of the time in this case, so I don't think it's too bad. It's the approach used by Aarseth, Henon, and Wielen, whose half-page appendix covers the initialization problem quite well. $\endgroup$ – user10851 Apr 7 '14 at 20:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.