`entanglement_entropy` can be a footgun (it has different meaning for Kets and Operators) -- consider deprecating in favor of separate functions

[acroscarrillo]

Edit (by @Krastanov): This is expected behavior. See first comment below for possible ways forward.

Original post below:

There seems to be a bug in the entanglement entropy calculation, it seems like it might be a factor of 2 off?

Here's the code to reproduce the error:

I create the following state and check is normalised:

julia> ψ_0
Ket(dim=8)
  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]
 0.5773502691896257 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
 0.5773502691896257 + 0.0im
                0.0 + 0.0im
 0.5773502691896257 + 0.0im

julia> ψ_0'*ψ_0
1.0 + 0.0im

Now I calculate the entanglement between the A={1,2} and B={3} as follows:

julia> entanglement_entropy(dm(ψ_0),[3])
1.2730283365896256 + 0.0im

However, this is wrong. I have done the calculation analytically (two methods, SVD and partial tracing) and it should be $-\frac{2}{3}\log(\frac{2}{3})-\frac{1}{3}\log(\frac{1}{3}) \approx 0.6365471166931516$, which could look like $1.2730283365896256 / 2$ but not exactly?

---

It's interesting cause your ptrace agrees with my calculations since

julia> ptrace(dm(ψ_0),[3]).data
4×4 Matrix{ComplexF64}:
 0.333333+0.0im       0.0+0.0im  0.0+0.0im       0.0+0.0im
      0.0+0.0im  0.333333+0.0im  0.0+0.0im  0.333333+0.0im
      0.0+0.0im       0.0+0.0im  0.0+0.0im       0.0+0.0im
      0.0+0.0im  0.333333+0.0im  0.0+0.0im  0.333333+0.0im

julia> eigvals(ptrace(dm(ψ_0),[3]).data)
4-element Vector{Float64}:
 0.0
 0.0
 0.3333333333333333
 0.6666666666666666

and by inspection this spectrum does indeed agree with my calculations. I am a bit puzzled cause in the source code you guys have

function entropy_vn(rho::DenseOpType{B,B}; tol=1e-15) where B
    evals::Vector{ComplexF64} = eigvals(rho.data)
    entr = zero(eltype(rho))
    for d ∈ evals
        if !(abs(d) < tol)
            entr -= d*log(d)
        end
    end
    return entr
end

so it should be doing the same provided the spectrum is the same, which it should be since it is being passed the same density matrix.

function entanglement_entropy(psi::Ket{B}, partition, entropy_fun=entropy_vn) where B<:CompositeBasis
    # check that sites are within the range
    @assert all(partition .<= length(psi.basis.bases))

    rho = ptrace(psi, partition)
    return entropy_fun(rho)
end

To check I am not crazy:

julia> -0.3333333333333333*log(0.3333333333333333)-0.6666666666666666*log(0.6666666666666666)
0.6365141682948128

I guess it could be that I am miss using the partition argument, I didn't find the documentation very clear.