Operator entanglement of an MPO

Hi,

I would like to compute the operator entanglement of an MPO, defined as S_{op} = - \sum_{\alpha} \lambda_\alpha^2 \log_2{\lambda_\alpha^2} where \lambda_\alpha are the Schmidt values across a cut of the system (see this paper). This is the same definition as the entanglement entropy for MPS states but extended to MPOs. Unlike the MPS case, the operator entanglement is not a genuine measure of entanglement. For example, the operator entanglement across bond b (with N_A on the left and N_B sites on the right) of the maximally mixed state \rho_{mm} = \frac{1}{2^N} \mathbb{I}^{\otimes N} is

S_{op}( \rho_{mm}, b ) = N_A \log_2(2)

since the reduced density matrix on \rho_{N_A} is

\rho_{N_A} = \frac{1}{2^N_A} \mathbb{I}^{\otimes N_A} = \sum_{\alpha_A} \frac{1}{2^{N_A}} | \alpha_A \rangle \langle \alpha_A | \\

where | \alpha_A \rangle is an orthonormal basis, making the schmidt values \lambda_\alpha ^2 = \frac{1}{2^{N_A}} the Schmidt values.

I thought I could compute the operator entanglement the same was as computing the von Neumann entanglement for the MPS, but I’m not getting the expected result for the maximally mixed state. Do you know what I’m doing wrong?

using ITensors, ITensorMPS

function operator_entanglement(state::Union{MPS, MPO}, bond::Int)
    state = orthogonalize(state, bond)

    U,S,V = svd(state[bond], (linkinds(state, bond-1)..., siteinds(state, bond)...))
    SvN = 0.0
    for n=1:dim(S, 1)
        p = S[n,n]^2
        SvN -= p * log2(p)
    end
    return SvN
end

let 
    N = 4
    s = siteinds("Qubit", N)

    # Operator entanglement of the maximally mixed state
    max_mixed = MPO(s, "Id")/2^N
    Na = 3
    Sop_mix = operator_entanglement(max_mixed, Na+1)
    @show maxlinkdim(max_mixed)
    println("Operator Entanglement = ", Sop_mix)
    println("Exact Operator Entanglement = ", Na*log2(2))
    println("Does it match the expected value?", Sop_mix == Na * log2(2) )

end

gives the following output:

maxlinkdim(max_mixed) = 1
Operator Entanglement = 0.25
Exact Operator Entanglement = 3.0
Does it match the expected value?false