add entanglement entropy (#6)

* add entanglement entropy

* tests for entanglement_entropy

* Improve docstring, dispatch for entanglement_entropy

* Bump version v0.1.1 -> v0.1.2

Author: PhilipVinc

Project.toml

@@ -1,6 +1,6 @@
 name = "QuantumOpticsBase"
 uuid = "4f57444f-1401-5e15-980d-4471b28d5678"
-version = "v0.1.1"
+version = "v0.1.2"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

src/QuantumOpticsBase.jl

@@ -28,7 +28,7 @@ export bases, Basis, GenericBasis, CompositeBasis, basis,
         metrics, tracenorm, tracenorm_h, tracenorm_nh,
                 tracedistance, tracedistance_h, tracedistance_nh,
                 entropy_vn, fidelity, ptranspose, PPT,
-                negativity, logarithmic_negativity,
+                negativity, logarithmic_negativity, entanglement_entropy,
         PauliBasis, PauliTransferMatrix, DensePauliTransferMatrix,
                 ChiMatrix, DenseChiMatrix, avg_gate_fidelity
 

src/metrics.jl

@@ -236,3 +236,46 @@ function avg_gate_fidelity(x::T, y::T) where T <: Union{PauliTransferMatrix{B, B
     dim = 2 ^ length(x.basis_l)
     return (tr(transpose(x.data) * y.data) + dim) / (dim^2 + dim)
 end
+
+"""
+    entanglement_entropy(state, partition, [entropy_fun=entropy_vn])
+
+Computes the entanglement entropy of `state` between the list of sites `partition`
+and the rest of the system. The state must be defined in a composite basis.
+
+If `state isa AbstractOperator` the operator-space entanglement entropy is
+computed, which has the property
+```julia
+entanglement_entropy(dm(ket)) = 2 * entanglement_entropy(ket)
+```
+
+By default the computed entropy is the Von-Neumann entropy, but a different
+function can be provided (for example to compute the entanglement-renyi entropy).
+"""
+function entanglement_entropy(psi::Ket{B}, partition::Vector{Int}, 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
+
+function entanglement_entropy(rho::DenseOperator{B,B}, partition::Vector{Int}, args...) where {B<:CompositeBasis}
+    # check that sites is within the range
+    hilb = rho.basis_l
+    all(partition .<= length(hilb.bases)) || throw(ArgumentError("Indices in partition must be within the bounds of the composite basis."))
+    length(partition) <= length(hilb.bases) || throw(ArgumentError("Partition cannot include the whole system."))
+
+    # build the doubled hilbert space for the vectorised dm, normalized like a Ket.
+    b_doubled = hilb^2
+    rho_vec = normalize!(Ket(b_doubled, vec(rho.data)))
+
+    return entanglement_entropy(rho_vec,
+                                vcat(partition, partition.+length(hilb.bases)),
+                                args...)
+end
+
+entanglement_entropy(state, partition::Number, args...) =
+    entanglement_entropy(state, [partition], args...)
+entanglement_entropy(state, partition, args...) =
+    entanglement_entropy(state, collect(partition), args...)

test/test_metrics.jl

@@ -90,6 +90,17 @@ rho = rho ⊗ dm(g)
 @test logarithmic_negativity(rho, 1) == logarithmic_negativity(rho, 2) ≈ 1.0
 @test isapprox(logarithmic_negativity(rho, 3), 0.0, atol=1e-15)
 
+# Entanglement Entropy
+b1 = SpinBasis(1//2)
+psi = 1/sqrt(2)*(spinup(b1)⊗spindown(b1) + spindown(b1)⊗spinup(b1))
+rho_ent = dm(psi)
+rho_mix = DenseOperator(rho_ent.basis_l, diagm(ComplexF64[1.0,1.0,1.0,1.0]))
+@test entanglement_entropy(rho_mix, 1) ≈ 0.0
+@test entanglement_entropy(rho_ent, 1, entropy_vn) ≈ 2 * log(2)
+@test entanglement_entropy(rho_ent, 1) ≈ 2 * entanglement_entropy(psi, 1)
+@test_throws ArgumentError entanglement_entropy(rho_mix, (1,2))
+@test_throws ArgumentError entanglement_entropy(rho_mix, 3)
+
 q2 = PauliBasis(2)
 CNOT = DenseOperator(q2, q2, diagm(0 => [1,1,0,0], 1 => [0,0,1], -1 => [0,0,1]))
 CNOT_sop = SuperOperator(CNOT)