Thancks for @ascendants 's mention, I found a other way in sciPy.

from scipy.spatial import cKDTree

data: list[] = ...
kdtree = cKDTree(data)

# `query_ball_point` return `list[list[int]]`, it is a list of neighbors indeces
neighbors = kdtree.query_ball_point(data, r=Radius, workers=-1) # set workders to infinity many be quikly

# settle a for loop to fix it
num: list[int] = []
for i in neighbors:
    num.append(len(i)-1)  # look like it will count itself, so decrease 1

I had run it with a lammps result, the performance is good with 4300 atom