newtonian gravity - Velocity distribution in Plummer's models and others mass distributions

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?

Answer

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.