# Stable disk for $N$-body simulation

I want to simulate a Galaxy of $N$ particles which I have to generate first. What I have done so fare does not lead to the result I want to see. I probably did an error and don't see it. Here comes what I did:

1) Generate uniformly distributed particles inside a unit circle with mass $m_i = 1/N$:

For every point I do:

   R = 1;
phi = rand(0,2*pi);
m = 1/N;


Doing that I get the following distribution:

2) Generate velocities such that the circle does not explode or implode:

I want to generate something where $F_{\text{Gravity}} = F_{\text{Centrifugal}}$

Since $F_{\text{Centrifugal}} = \frac{mv^2}{r}$ and $F_{\text{Gravity}} = G\frac{M(r) m}{r^2}$ I can equate both equations and get for the velocity $v$

$M(r)$ is the mass inside radius $r$. I can write $M_{\text{tot}} = \pi R^2$ and $M=\pi r^2$ (density $\rho=1$). From those equations I get $M(r) = \frac{r^2}{R^2} M_{\text{tot}}$. I can now put it back to $v$ and get:

$v = \sqrt{\frac{Gr^2M_{\text{tot}}}{rR^2}}$

Since $M_{\text{tot}} = G = R = 1$ I get:

$v = \sqrt{r}$

So, for every point I do:

   velocity = sqrt(radius);
v_x = -velocity*sin(phi);
v_y = velocity*cos(phi);


Then I get the following result:

Ok, it seems that the circle explodes somehow. Particles from the center move outwards. Is my calculation for $v$ wrong, or does it simply not work for my purpose? Do I need another velocity distribution? In the end, for large $t$, I can't see any spiral arms forming or any other structure. I tried $N=4000$ and $N=50000$. Do I need more particles? Did I something wrong?

• the phi for the velocities should be chosen randomly (uniform in $[0,2\pi)$). – AccidentalFourierTransform Feb 25 '17 at 11:54
• @AccidentalFourierTransform I do that. I take the same phi I already take for the position of the particle. Is that wrong? Do I have to generate phi a second time for the velocity? – Samuel Feb 25 '17 at 12:10
• yes, that is wrong, and yes, you need to generate phi a second time for the velocity (try to think why!). – AccidentalFourierTransform Feb 25 '17 at 12:12
• By taking the same phi both for position and velocities, you are making all velocities to be radial (i.e., no tangential component). This is the reason for the initial "outward explosion". – AccidentalFourierTransform Feb 25 '17 at 12:34
• If the distribution of $\phi_v$ is uniform, then the initial cluster has no angular momentum. You may want to try with non-uniform $\phi_v$. – AccidentalFourierTransform Feb 25 '17 at 13:17

Your logic is correct, including the choice of tangential velocity as

 x = r*cos(phi);
y = r*sin(phi);
vx=-v*sin(phi);
vy= v*cos(phi);


Where you (tried to) set v in centrifugal balance. However, the gravitational attraction of a disc is different from that of a sphere, in particular your assumption that the acceleration equals GM(r)/r^2 is wrong. The best way to find the correct form is to obtain the forces from your N-body model and take the azimuthal average of the radial acceleration in radial bins. Your initial model should satisfy the virial theorem, I.e. the kinetic energy should be -0.5 times the gravitational energy.

But how did you compute the forces and how did you integrate the equations of motion? You should avoid the exact N-body forces, because these diverge at close encounters and instead use softened forces. You essentially add for each particle-particle interaction a constant vertical offset h, the softening length. Set this to about 0.1 times the radius.

The time integration is best done using the leapfrog integrator. For example

for(I:particles)
I.vel += 0.5*tau*I.acc
I.pos += tau*I.vel
time += tau
compute_accelerations(particles);
for(I:particles)
I.vel += 0.5*tau*I.acc


to integrate one time step of length tau. Of course, you must compute the initial accelerations before the first step. The step size tau should be small enough for the total energy to be conserved by one part in 1000 or better.

• I used an tree algorithm with an Runge Kutta integrator. Is my approach correct if I use a sphere (3D)? – Samuel Feb 26 '17 at 1:27
• You can assign tangential velocities in equilibrium to particles in a sphere, but the sphere will not rotate and is not a disc. Did you use softening? How precisely did you decide when to open a cell with the tree code? -- perhaps you were not careful with that and ran into the "exploding galaxies bug"? – Walter Feb 26 '17 at 8:13
• I already simulated a Plummer sphere and get good results. So I think there is an problem with the initial condition of the velocity. I tried your idea but get worse results. The sphere explodes. But at the very beginning there are some spiral arms visible! Maybe I did something wrong. Could you please carry out your idea and show how you would derive the right velocity for such a system? That would be great. – Samuel Feb 26 '17 at 17:37
• Unfortunately, you failed to provide sufficient information. How/what are the accelerations? How well are total energy, linear and angular momentum conserved???? – Walter Feb 26 '17 at 19:31