Fix entanglement and introduce concurrence

CHANGELOG.md

@@ -11,11 +11,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Change save callbacks from `PresetTimeCallback` to `FunctionCallingCallback`. ([#410])
 - Align `eigenstates` and `eigenenergies` to QuTiP. ([#411])
 - Introduce `vector_to_operator` and `operator_to_vector`. ([#413])
-- Introduce some entropy related functions. ([#416])
+- Introduce some entropy related functions. ([#414], [#416])
   - `entropy_linear`
   - `entropy_mutual`
   - `entropy_conditional`
   - `entropy_relative`
+- Fix `entanglement` and introduce `concurrence`. ([#414], [#418], [#419])
 
 ## [v0.27.0]
 Release date: 2025-02-14
@@ -154,4 +155,7 @@ Release date: 2024-11-13
 [#410]: https://github.com/qutip/QuantumToolbox.jl/issues/410
 [#411]: https://github.com/qutip/QuantumToolbox.jl/issues/411
 [#413]: https://github.com/qutip/QuantumToolbox.jl/issues/413
+[#414]: https://github.com/qutip/QuantumToolbox.jl/issues/414
 [#416]: https://github.com/qutip/QuantumToolbox.jl/issues/416
+[#418]: https://github.com/qutip/QuantumToolbox.jl/issues/418
+[#419]: https://github.com/qutip/QuantumToolbox.jl/issues/419

docs/src/resources/api.md

@@ -259,6 +259,8 @@ entropy_linear
 entropy_mutual
 entropy_conditional
 entanglement
+concurrence
+negativity
 tracedist
 fidelity
 ```
@@ -282,7 +284,6 @@ BlockDiagonalForm
 
 ```@docs
 wigner
-negativity
 ```
 
 ## [Linear Maps](@id doc-API:Linear-Maps)

src/entropy.jl

@@ -1,9 +1,9 @@
 #=
-Entropy related functions.
+Entropy related functions and some entanglement measures.
 =#
 
 export entropy_vn, entropy_relative, entropy_linear, entropy_mutual, entropy_conditional
-export entanglement
+export entanglement, concurrence
 
 @doc raw"""
     entropy_vn(ρ::QuantumObject; base::Int=0, tol::Real=1e-15)
@@ -121,7 +121,7 @@ function entropy_relative(
 
     # the relative entropy is guaranteed to be ≥ 0
     # so we calculate the value to 0 to avoid small violations of the lower bound.
-    return max(0.0, dot(p_vals, log_p) - dot(p, P, log_q))
+    return max(0.0, dot(p_vals, log_p) - dot(p, P, log_q)) # use 0.0 to make sure it always return value in Float-type
 end
 
 @doc raw"""
@@ -195,7 +195,7 @@ Calculates the [entanglement entropy](https://en.wikipedia.org/wiki/Entropy_of_e
 
 # Notes
 
-- `ρ` can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+- `ρ` can be either a [`Ket`](@ref) or an [`Operator`](@ref). But should be a pure state.
 - `sel` specifies the indices of the remaining sub-system. See also [`ptrace`](@ref).
 - `kwargs` are the keyword arguments for calculating Von Neumann entropy. See also [`entropy_vn`](@ref).
 """
@@ -204,8 +204,43 @@ function entanglement(
     sel::Union{Int,AbstractVector{Int},Tuple},
     kwargs...,
 ) where {OpType<:Union{KetQuantumObject,OperatorQuantumObject}}
-    _ρ = normalize(ρ)
-    ρ_tr = ptrace(_ρ, sel)
+    p = purity(ρ)
+    isapprox(p, 1; atol = 1e-2) || throw(
+        ArgumentError(
+            "The entanglement entropy only works for normalized pure state, the purity of the given state: $(p) ≉ 1",
+        ),
+    )
+
+    ρ_tr = ptrace(ρ, sel)
     val = entropy_vn(ρ_tr; kwargs...)
-    return (val > 0) * val
+    return max(0.0, val)  # use 0.0 to make sure it always return value in Float-type
+end
+
+@doc raw"""
+    concurrence(ρ::QuantumObject)
+
+Calculate the [concurrence](https://en.wikipedia.org/wiki/Concurrence_(quantum_computing)) for a two-qubit state.
+
+# Notes
+
+- `ρ` can be either a [`Ket`](@ref) or an [`Operator`](@ref).
+"""
+function concurrence(ρ::QuantumObject{OpType}) where {OpType<:Union{KetQuantumObject,OperatorQuantumObject}}
+    (ρ.dimensions == Dimensions((Space(2), Space(2)))) || throw(
+        ArgumentError(
+            "The `concurrence` only works for a two-qubit state, invalid dims = $(_get_dims_string(ρ.dimensions)).",
+        ),
+    )
+
+    _ρ = ket2dm(ρ).data
+    σy = sigmay()
+    σyσy = kron(σy, σy).data
+    ρ_tilde = σyσy * conj(_ρ) * σyσy
+
+    # we use the alternative way to calculate concurrence (more efficient)
+    # calculate the square root of each eigenvalues (in decreasing order) of the non-Hermitian matrix: ρ * ρ_tilde
+    # note that we add abs here to avoid problems with sqrt for very small negative numbers
+    λ = sqrt.(abs.(real(eigvals(_ρ * ρ_tilde; sortby = x -> -real(x)))))
+
+    return max(0.0, λ[1] - λ[2] - λ[3] - λ[4]) # use 0.0 to make sure it always return value in Float-type
 end

test/core-test/entropy_and_metric.jl

@@ -51,17 +51,39 @@
     end
 end
 
-@testset "Entanglement" begin
-    g = fock(2, 1)
-    e = fock(2, 0)
-    state = normalize(kron(g, e) + kron(e, g))
-    rho = state * state'
-    @test entanglement(state, 1) / log(2) ≈ 1
-    @test entanglement(rho, 1) / log(2) ≈ 1
+@testset "entanglement and concurrence" begin
+    # bell state
+    ψb = bell_state(Val(1), Val(0))
+    ρb = ket2dm(ψb)
+    @test entanglement(ψb, 1) / log(2) ≈ 1
+    @test entanglement(ρb, 1) / log(2) ≈ 1
+    @test concurrence(ψb) ≈ 1
+    @test concurrence(ρb) ≈ 1
+
+    # separable pure state
+    ψs = kron(rand_ket(2), rand_ket(2))
+    @test entanglement(ψs, 1) + 1 ≈ 1
+    @test entanglement(ψs, 2) + 1 ≈ 1
+    @test concurrence(ψs) + 1 ≈ 1
+
+    # this only works for "pure" two-qubit states
+    ψr = rand_ket((2, 2)) # might be an entangled two-qubit state
+    val = concurrence(ψr)
+    @test isapprox(val, sqrt(2 * entropy_linear(ptrace(ψr, 1))); atol = 1e-5) # √(2 * (1 - Tr(ρA^2)))
+    @test isapprox(val, sqrt(2 * entropy_linear(ptrace(ψr, 2))); atol = 1e-5) # √(2 * (1 - Tr(ρB^2)))
+
+    @test_throws ArgumentError entanglement(rand_dm((2, 2)), 1)
+    @test_throws ArgumentError concurrence(rand_dm((2, 3)))
+    @test_throws ArgumentError concurrence(rand_dm(4))
 
     @testset "Type Stability (entanglement)" begin
-        @inferred entanglement(state, 1)
-        @inferred entanglement(rho, 1)
+        @inferred entanglement(ψb, 1)
+        @inferred entanglement(ρb, 1)
+    end
+
+    @testset "Type Stability (concurrence)" begin
+        @inferred concurrence(ψb)
+        @inferred concurrence(ρb)
     end
 end