what can i use to plot point with density in 3d?

Question

i have a set of point in a 3d space (with x,y,z value) that represent the position of physical object and i want to see the movement of these object. if a have 10 physical object and 100 point i want to plot the first 0-10 point then 10-20 and erase 0-10 then plot 20-30 and erase 10-20 ,so i can see the movement of the object. the number of frame will be (number of point/number of physical body). i already tried matplotlib with the animation function but with 10000 physical body and 3000000 point i get barely 10 frame/second. my question is what can i use to see the movement of these object with point that represent there position overtime. the requirement that is 1 it need to be fast to get more then 10 frame/second 2 be able to see the density of point 3 if i can see the density of point it can be ether in 3d or in 2d but i will prefer 3d

i already tried matplotlib and seaborn

it can be in python or c++ but i prefere python

this is the code that i was using with matplotlib the code to compute these point is write in c++ the whole thing is for a N-body simulation if it can help

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
import matplotlib.pyplot as plt
import numpy as np
import math
import time
from pathlib import Path

data = Path("C:/Users/JF/source/repos/NBodyC++/NBodyC++/")
dataopen=data/"gravity.cvs"
gravity=open(dataopen)
w= np.loadtxt(gravity,delimiter=",")
nombreplanete=1000
l=0
n=0
n2=nombreplanete
inter=1
fig = plt.figure("black")
fig.patch.set_facecolor('black')

ax = fig.add_subplot(111, projection='3d',)
ax.tick_params( colors='white')
ax.set_xlabel('X',color="white")
ax.set_ylabel('Y',color="white")
ax.set_zlabel('Z',color="white")
ax.set_xlim(-500,500)
ax.set_ylim(-500,500)
ax.set_zlim(-500,500)

ax.set_facecolor("black")


start_time = time.time()
def draw_graph(i):
    global etoile
    global nombreplanete
    global l
    global n 
    global n2
    
    if n2>nombreplanete:
        ax.collections[0].remove()

    ax.scatter(w[n:n2,0],w[n:n2,1],w[n:n2,2] ,s=5, c='white', marker='.')
    n+=nombreplanete
    n2+=nombreplanete
    l+=1
    
    if l%10==0:
        print("nombre d'image par seconde:",l/(time.time() - start_time),"iteration:",l)
                   
anima = animation.FuncAnimation(plt.gcf(),draw_graph,interval=1)
plt.show()