贝尔对角态纠缠熵计算结果异常排查求助

我最近在尝试计算贝尔对角态的纠缠熵，试了好几种Python实现，但所有结果都显示这些态的纠缠熵全为1，这显然是错误的——根据相关研究，只有贝尔对角态构成的四面体顶点处的态才是纠缠态，其他态的纠缠熵应该有所不同才对。我完全找不到问题出在哪，麻烦各位帮忙看看！

以下是我的代码：

import numpy as np
import qutip as q


def tetrahedron(a1, a2, a3):

    is_bellD_state = False

    plano1 = a1 - a2 + a3
    plano2 = a1 + a2 - a3
    plano3 = a1 - a2 - a3
    plano4 = a1 + a2 + a3
    
    if plano1 >= -1 and plano2 >= -1 and plano3 <= 1 and plano4 <= 1:
        is_bellD_state = True

    return is_bellD_state


def Bell_Diagonal_Density_Matrix(e):

    bell_density_matrices = []

    for b in ['00', '01', '10', '11']:
        bell_density_matrices.append( q.bell_state(b) * q.bell_state(b).dag() )
    
    rho = q.Qobj(np.zeros((4,4)), dims=[[2,2],[2,2]])

    for i in range(4):
        rho = rho + e[i]*bell_density_matrices[i]
    
    return rho


def Entanglement_Entropy(rho_A):

    eigvals = np.linalg.eigvalsh(rho_A)
    eigvals = eigvals[eigvals > 0]

    S = 0

    for v in eigvals:
        S = S - v*np.log(v)/np.log(2)

    return S



## MAIN ##

a = np.arange(-1, 1, 0.1)
for a1 in a:
    for a2 in a:
        for a3 in a:

            if tetrahedron(a1, a2, a3):
            
                M = [[1,1,-1,1], [1,1,1,-1], [1,-1,1,1], [1,-1,-1,-1]]
            
                e = np.dot(M, [1, a1, a2, a3])/4
            
                if np.isclose(sum(e), 1):
                
                    rho = Bell_Diagonal_Density_Matrix(e)
                
                    rho_A = rho.ptrace(0)
                
                    S = Entanglement_Entropy(rho_A)
                    print(S)

我之前怀疑是自己手动计算偏迹出错，所以改用了QuTiP库，但结果还是不对，我也不确定是不是QuTiP的使用方式出了问题。

备注：内容来源于stack exchange，提问作者C. R. L.