[Python] Entanglement entropy of Bell Diagonal States returns always one

I've been trying to compute the entanglement entropy for Bell Diagonal states. I tried many implementations in python and all of them gave me that all possible states have this entropy equal to 1.

This is completely wrong, because the only entangled states should be the points of the tetrahedron formed by the Bell Diagonal states, as it is shown in Figure 7 of this article.

I have no idea of what I'm doing wrong.

This is my code:

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*bell_density_matrices

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)


I decided to use qutip because I thought I was computing the partial trace wrong, but I'm also not sure if I'm using qutip correctly.

Continue reading...