Surge in energy at start of N-body simulation

Question

I am trying to use an n-body simulation to simulate the collapse of a globular cluster, starting from a point where the bodies all have 0 initial velocities generated in random positions of a sphere of a fixed radius and I am investigating the time taken to form an equilibrium system. However, when running my code when the particles initially come close together there is a large surge in potential energy. Is there a way to avoid this? I am using the leapfrog approximation to calculate subsequent velocities and positions of the particles.

initial conditions

sun_mass = 1.989e30


N = 50
mass = 25 * sun_mass
M = np.full([N],mass)    
R = 3.086e+16          # 1pc
epsilon = 0.1 * R
collision_radius = 7e8
V = np.zeros([N,3])
P = np.zeros([N,3])

t = 50000000 * 365 * 24 * 60 * 60
dt = 18000000 * 24 * 60 * 60



print(t/dt)

np.random.seed(54321)



for i in range(N):
    
    #M[i] = np.random.uniform(sun_mass, 100*sun_mass, 1)
    
    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)

code to run programme:

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

    #all_positions = np.full((no_of_time_steps, len(mass), 3), 0.0)
    all_positions = []
    all_velocities = []
    #print(all_positions)
    #print(len(all_positions[0]))
    kinetic_energy = []
    potential_energy = []
    total_energy = []
    
    #previous_velocity = calc_previous_half_velocity(velocity, calc_acceleration(position, mass, softening), dt)
        
    for i in range(no_of_time_steps):
        
        position, mass, velocity = detect_collisions(position, mass, velocity, collision_radius)
        
        #all_positions[i] = position
        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, mass, R, i, dt, 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)
        

    #columns_to_delete = len(all_velocities)
    #print(no_of_time_steps)
    #print(len(kinetic_energy), len(potential_energy), len(total_energy))
    
    
    all_positions = np.array(all_positions)
    all_velocities = np.array(all_velocities)
    #print(all_positions)
    graphing(all_positions, all_velocities, mass, 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)
        

leapfrog functions:

'LeapFrog functions'

'Acceleration Calculation'
def calc_acceleration(position, mass, softening):
    
    G = 6.67 * 10**(-11)
    
    N = position.shape[0] # N = number of rows in particle_positions array
    acceleration = np.zeros([N,3])
    
    #print(N)
    for i in range(N):
        #print(i)
        for j in range(N):
            if i != j:
                #print("j", j)
                dx = position[i,0] - position[j,0]
                dy = position[i,1] - position[j,1]
                dz = position[i,2] - position[j,2]
                
                #print(dx,dy,dz)
                
                inv_r3 = ((dx**2 + dy**2 + dz**2 + softening**2)**(-1.5))
                
                acceleration[i,0] += - G * mass[j] * dx * inv_r3
                acceleration[i,1] += - G * mass[j] * dy * inv_r3
                acceleration[i,2] += - G * mass[j] * dz * inv_r3

    return(acceleration) #, print(acceleration))

def calc_previous_half_velocity(velocity, acceleration, dt):
    previous_velocity = np.zero_like(velocity)
    
    previous_velocity = velocity - acceleration * dt/2
    return(previous_velocity)

def calc_next_v_half(position, mass, velocity, softening, dt):
    half_velocity = np.zeros_like(velocity)
    half_velocity = velocity + calc_acceleration(position, mass, softening) * dt/2
    return(half_velocity)
    
def calc_next_velocity(position, mass, velocity, softening, dt):
    next_velocity = np.zeros_like(velocity)
    
    next_velocity = velocity + calc_acceleration(position, mass, softening) * dt
    
    return(next_velocity)


def calc_next_position(position, mass, velocity, dt):
    next_position = np.zeros_like(position)
    
    next_position = position + velocity * dt
    
    return(next_position)

other functions used:

def calc_CoM(position, mass):
        sumMR = np.zeros(3)
        sumM = 0.0
        position = np.array(position)
        for i in range(len(mass)):
            sumM += mass[i]
            sumMR += mass[i] * position[i]
        return(sumMR / sumM)


def calc_total_mass2(mass):
    total_mass = 0.0
    print((mass[0]))
    for i in range(len(mass)):
        total_mass += mass[i]
    return(total_mass)

def calc_total_mass(mass):
    total_mass = sum(mass)
    return(total_mass)

def calc_CoM_seperation(position, mass, i):
    new_mass = np.array(mass) 
    new_pos = np.array(position)
    new_mass[i] = 0.0
    new_pos[i] = 0.0
    position = np.array(position)
    r = np.linalg.norm(calc_CoM(new_pos, new_mass) - position[i])
    return(r)

def calc_seperation(p1, p2):
    return np.linalg.norm(p1-p2)