N-Body Simulation of Globular Cluster

Question

I am trying to code an n-body simulation using the leap frog scheme to simulate a globular cluster however am running into the issue of particles being flung out of the system from (I think) approaching to close to the other particles as well as this leading to a large increase in potential energy. I have tried to write some code to account for collisions, but am unsure what collision radius to use if I am using a softening factor as well.

The initial positions are randomly generated within a sphere and all have 0 initial velocities.

If the softening factor is set to 0.1 * Radius of the sphere of generation and the collision radius set to the radius of the sun (leading to no collisions) with N bodies of solar mass I get the following results:

From the graphs of energy, we can see there a huge spike in potential energy when the bodies originally get close to one another. And over time my energy trails off to 0, it think just due to numerical calculations.

Is there a way to stop these stars being flung out. I am trying to investigate how long it takes the initial cluster to form an equilibrium state.

sphere generation:

sun_mass = 1.989e30


N = 100
mass = sun_mass
M = np.full([N],mass)
R = 1e13
epsilon = 0.1 * R
collision_radius = 7e8
V = np.zeros([N,3])
M = np.full([N],mass)
P = np.zeros([N,3])


t = 50 * 365 * 24 * 60 * 60
dt = 30 * 24 * 60 * 60



print(t/dt)

np.random.seed(54321)

for i in range(N):
    
    phi = np.random.uniform(0,2*np.pi)
    costheta = np.random.uniform(-1,1)
    u = np.random.uniform(0,1)
    
    theta = np.arccos( costheta )
    r = R * (u) **(1/3)
    
    
    x = r * np.sin( theta) * np.cos( phi )
    y = r * np.sin( theta) * np.sin( phi )
    z = r * np.cos( theta )

    P[i] = (x,y,z)

Programe:

def programe(position, mass, velocity, softening, time, dt, R, collision_radius):
    
    no_of_time_steps = (round(time/dt))

    all_positions = []
    all_velocities = []
    #print(all_positions)
    #print(len(all_positions[0]))
    kinetic_energy = []
    potential_energy = []
    total_energy = []
    
        
    for i in range(no_of_time_steps):
        
        position, mass, velocity = detect_collisions(position, mass, velocity, collision_radius)
        
        all_positions.append(position)
        all_velocities.append(velocity)
    
        'graph'
        plots = np.round(np.linspace(0,no_of_time_steps,num=500))
        for k in range(len(plots)):
            if i == plots[k]:
                print("test")
                #print(i)
                graph(position, R, k)
        
        'energies'
        kinetic_energy.append(calc_kinetic_energy(velocity, mass))
        potential_energy.append(calc_potential_energy(position, mass))
        total_energy.append(calc_total_energy(position, velocity, mass))
        
        'leap frog'
        velocity = calc_next_v_half(position, mass, velocity, softening, dt)    
        position = calc_next_position(position, mass, velocity, dt)    
        velocity = calc_next_v_half(position, mass, velocity, softening, dt)
        

    
    all_positions = np.array(all_positions)
    graphing(all_positions, time, dt, kinetic_energy, potential_energy, total_energy, no_of_time_steps, R)
    #print(len(mass))

    return(all_positions, all_velocities, kinetic_energy, potential_energy, total_energy)

collision detection:

def indexOf(Array, item):
    for i in range(len(Array)):
        if (Array[i] == item).all():
            return i
        
def detect_collisions(position, mass, velocity, collision_radius):
    i = 0
    
    newP = position
    newM = mass
    newV = velocity
    #print(len(position), len(newM))
    while i < len(newM):
        if newM[i] == 0:
            i += 1
            continue
        j = i + 1
        while j < len(newP):
           #print(j)
            if newM[j] == 0:
                j += 1
                continue
            p1 = position[i]
            p2 = position[j]
            
            if calc_seperation(p1,p2) < collision_radius:
                index1 = indexOf(position, p1)
                index2 = indexOf(position, p2)
                print('collision', index1, index2)
                newM, newV, newP = handle_collision(newM, newV, newP, [index1, index2])
            j += 1
        i += 1
    
    return(newP, newM, newV)

def handle_collision(M, V, P, indexes):
    if M[indexes[0]] > M[indexes[1]]:
        primary = indexes[0]
        secondary = indexes[1]
    else:
        primary = indexes[1]
        secondary = indexes[0]

    primaryMomentum = M[primary] * V[primary]
    secondaryMomentum = M[secondary] * V[secondary]
    newMomentum = primaryMomentum + secondaryMomentum
    
    newMass = M[primary] + M[secondary]
    newVelocity = newMomentum / newMass
    M[primary] = newMass
    V[primary] = newVelocity
    M[secondary] = 0
    V[secondary] = 0
    P[secondary] = 0

    newM = np.delete(M, secondary, axis=0)
    newV = np.delete(V, secondary, axis=0)
    newP = np.delete(P, secondary, axis=0)

    return (newM, newV, newP)