Initial conditions of nbody problem
Note: I have edited the original question to try to be clearer.
I'm trying to simulate a gravitational nbody problem with a <PERSON>-Hut algorithm, in which the bodies start in a homogeneous sphere with random velocities. Because of this, the initial angular momentum is almost zero. I have tested the code with things like the solar system and compared its results to a brute force algorithm and it seems to work fine.
Despite using very low time steps (I'm using a velocity verlet integrator, if that matters) the angular momentum increases a lot, not being conserved at all. This image shows a typical increase of the momentum (later I will show more detailled data)
You can see that the energy is being conserved perfectly but the angular momentum keeps increasing with the bodies being accelerated. 617
All the angular momentums are computed from the actual center of mass of the cluster.
When a body has a close encounter and gets expelled from the cluster, both the angular momentum and the energy increase a lot:
I have tried different timesteps (from $\delta t= 10^{-8}t_0$ to $\delta t=10^{-3}t_0$) and gravitational softening parameters $\epsilon$ (from $F_{ij} \propto r_{ij}(\epsilon + r_{ij}^3)^{-1}$) from $1$ light year to $10^{-5}$ light years (the initial separation between bodies is $6$ ly), and they all have the same general results (there are 8 series in each line but the deviation is minimal):
The sudden rise of the lines is caused by a body being expelled by a close encounter. My principal concern is the initial increase of the angular momentum, because it doesn't come with a variation in the total energy. Because real galaxies have angular momentum, I thought that maybe I should start the simulation with a non zero angular momentum.
My questions are:
* What are the standard initial conditions in this kind of problems? They involve a small initial angular momentum or they use a "spinning" cluster to begin the simulation?
* Which variation in angular momentum and energy can be considered reasonable in nbody simulations?
Edit: I wasn't aplying gravitational softening, but it does not change the results.
Edit: This is the relevant code of the angular momentum, I can see a flaw but who knows...
double acc = 0;
for(int i=0;i<nstars;i++){
double r[3] = { x[i]-cmx , y[i]-cmy , z[i]-cmz };
double v[3] = { vx[i] , vy[i] , vz[i] };
double li[3] = { r[1]*v[2]-r[2]*v[1], -(r[0]*v[2] - r[2]*v[0]), r[0]*v[1] - r[1]*v[0] };
// Li = Rᵢ×Vᵢ
acc += m[i]*norm(li);
// acc += m⋅R×V
}
return acc;
cmx, cmy, cmz are the coordinates of the mass (just checked, it doesn't move a bit and it's where it should be). Using a fixed reference frame makes no difference.
Edit: Just tried gravitational softening (it was badly implemented). It helps a lot to avoid bodies fly away, but makes no difference in the constant increase of the angular momentum. 512