introduce concurrence

Author: ytdHuang

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

@@ -3,7 +3,7 @@ Entropy related functions.
 =#
 
 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"""
@@ -209,3 +209,28 @@ function entanglement(
     val = entropy_vn(ρ_tr; kwargs...)
     return (val > 0) * val
 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 support 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,38 @@
     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 concurrence(rand_dm(4))
+    @test_throws ArgumentError concurrence(rand_dm((2, 3)))
 
     @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