Mysterious Python Segmentation fault

Question

Problem

I have received Segmentation fault from a pure Python package : pytreegrav.

What is mysterious is that larger size of arrays don't lead to Segmentation fault, while smaller size of arrays do lead to.

What I did & received

The code outputs

Kinetic complete.  
Segmentation fault (core dumped)

To identify the place producing Segmentation fault, I put print("Kinetic complete.") and print("Potential complete."). I received the first but not the latter, so I think ptgPotential causes Segmentation fault.
I added del and gc.collect() and it doesn't make any change.
I tried sys.settrace and gdb python refering to What causes a Python segmentation fault?. However I couldn't use it well.
I just put sys.settrace(None) in the first block of my code, and (gdb) backtrace returned

#0  0x00007fffbd5e6bd0 in pytreegrav::octree::Octree::BuildTree$2415(instance::jitclass::Octree$237fffd6b84c10$3cSizes$3aarray$28float64$2c$201d$2c$20A$29$2cDeltas$3aarray$28float64$2c$201d$2c$20A$29$2cCoordinates$3aarray$28float64$2c$202d$2c$20A$29$2cMasses$3aarray$28float64$2c$201d$2c$20A$29$2cQuadrupoles$3aarray$28float64$2c$203d$2c$20A$29$2cHasQuads$3abool$2cNumParticles$3aint64$2cNumNodes$3aint64$2cSoftenings$3aarray$28float64$2c$201d$2c$20A$29$2cNextBranch$3aarray$28int64$2c$201d$2c$20A$29$2cFirstSubnode$3aarray$28int64$2c$201d$2c$20A$29$2cTreewalkIndices$3aarray$28int64$2c$201d$2c$20A$29$3e, Array, Array, Array) ()

... (long backtrace)  

#46 0x000055555573afe5 in _start () at ../sysdeps/x86_64/elf/start.S:103

(gdb) info locals returned No symbol table info available.

Script

It's long, so I cut off a large part of the code.

from pytreegrav import Potential as ptgPotential

...(cut)

for j,i in enumerate(snaplist):
    with h5py.File(datafile, 'r') as f:
        pos_DM = np.asarray(f["PartType1/Coordinates"])
        vel_DM = np.asarray(f["PartType1/Velocities"], dtype="float64") * 1e-2  # kpc/Myr
        mas_DM = np.asarray(f["PartType1/Masses"]) * 1e8  # M_sun
        IDs_DM = np.asarray(f["PartType1/ParticleIDs"], dtype="int32")
    
    # Withdraw satellite particles having been bound in previous snapshot
#------------------------------#
#    pos_DM_sate_bnd = pos_DM
#    vel_DM_sate_bnd = vel_DM
#    mas_DM_sate_bnd = mas_DM
#    IDs_DM_sate_bnd = IDs_DM
#    num_DM_sate_bnd = len(IDs_DM_sate_bnd)
#------------------------------#
    idxs_DM_sate_bnd = np.empty_like(IDs_DM_sate_bnd)
    libc1.get_bound_idxs(IDs_DM, IDs_DM_sate_bnd, idxs_DM_sate_bnd)
    pos_DM_sate_bnd = pos_DM[idxs_DM_sate_bnd,:]
    vel_DM_sate_bnd = vel_DM[idxs_DM_sate_bnd,:]
    mas_DM_sate_bnd = mas_DM[idxs_DM_sate_bnd]
    IDs_DM_sate_bnd = IDs_DM[idxs_DM_sate_bnd]
    num_DM_sate_bnd = len(idxs_DM_sate_bnd)
#------------------------------#
    del pos_DM; del vel_DM; del mas_DM; del IDs_DM; gc.collect()


# Not important section from here...  
    # Calculate BcV
    BcV_DM_sate_bnd = calc_CoM(vel_DM_sate_bnd, num_DM_sate_bnd)

    # Calculate kinetic energies
    num_DM_sate_bnd_ctp = ctypes.c_int(num_DM_sate_bnd)
    mass_sate_ctp = ctypes.c_double(mas_DM_sate_bnd[0])
    kin_ene = np.empty(num_DM_sate_bnd, dtype="float64")
    libc1.calc_kin_ene( vel_DM_sate_bnd - BcV_DM_sate_bnd, mass_sate_ctp, num_DM_sate_bnd_ctp, kin_ene)
    del vel_DM_sate_bnd; gc.collect()
    print("Kinetic complete.")
# to here.
    

    # Calculate potential energies
    eps_DM = np.repeat(eps, num_DM_sate_bnd)
    pot_ene = ptgPotential(pos_DM_sate_bnd, mas_DM_sate_bnd, softening=eps_DM, G=4.493e-12, theta=.6, parallel=True, method="tree")
    print("Potential complete.")

...

Helps for you to read the script

The function ptgPotential is a calculator of potential energy, and is also the causer of Segmentation fault.
The flow is following:
1. Read the particle data.
2. Pick out a part of the particles using idxs_DM_sate_bnd. If I use the comment-out part (switch the section #-----#), the code uses all particles.
  2.5. (Not important) Calculate the baryocentric velocity (BcV).
3. (Not important) Calculate kinetic energies of each particles.
4. Calculate potential energies of each particles.
pytreegrav consists of only python, but it uses jit compiler.

Context

This is an analysis code for N-body gravitational simulation.
The number of particles (N) = 23075000 (~10^7).

I want to calculate bounding energy (= kinetic energy - potential energy).

I use C library to calculate kinetic energy, and it works with no problem.

For potential energy, I use a python package "pytreegrav". And it seems to cause Segmentation fault.

What is mysterious, is that it doesn't raise Segmentation fault with all (N~10^7) particles while it raises if I pick out and use N=8075000 (~10^6) particles.

I suspect that it's not pytreegrav's fault, so I put this question here just in case.

Environment

The machine is a cluster of an institute. I don't have access to the details of OS.

Linux
Python 3.8.8 (Anaconda3-2022.05-Linux-x86_64)
Numpy 1.20.1
Numba 0.53.1