Add generalized pauli operators (aka clock and shift matrices) (#193)

Author: akirakyle

src/QuantumOpticsBase.jl

@@ -58,7 +58,7 @@ export Basis, GenericBasis, CompositeBasis, basis,
                 PositionBasis, MomentumBasis, samplepoints, spacing, gaussianstate,
                 position, momentum, potentialoperator, transform,
         #nlevel
-                NLevelBasis, transition, nlevelstate,
+                NLevelBasis, transition, nlevelstate, paulix, pauliz, pauliy,
         #manybody
                 ManyBodyBasis, FermionBitstring, fermionstates, bosonstates,
                 manybodyoperator, onebodyexpect,

src/nlevel.jl

@@ -31,3 +31,65 @@ function nlevelstate(::Type{T}, b::NLevelBasis, n::Integer) where T
     basisstate(T, b, n)
 end
 nlevelstate(b::NLevelBasis, n::Integer) = nlevelstate(ComplexF64, b, n)
+
+"""
+    paulix([T=ComplexF64,] b::NLevelBasis, pow=1)
+
+Generalized Pauli ``X`` operator for the given `N` level system. This is also
+called the shift matrix, and generalizes the qubit pauli matrices to qudits,
+while preserving the operators being unitary. A different generalization is
+given by the angular momentum operators which preserves the operators being
+Hermitian and is implemented for the `SpinBasis` (see [`sigmax`](@ref)).
+
+Powers of `paulix` together with [`pauliz`](@ref) form a complete, orthornormal
+(under Hilbert–Schmidt norm) operator basis.
+
+Returns `X^pow`.
+"""
+function paulix(::Type{T}, b::NLevelBasis, pow=1) where T
+    N = length(b)
+    SparseOperator(b, spdiagm(pow => fill(one(T), N-pow),
+                              pow-N => fill(one(T), pow)))
+end
+paulix(b::NLevelBasis, pow=1) = paulix(ComplexF64, b, pow)
+
+"""
+    pauliz([T=ComplexF64,] b::NLevelBasis, pow=1)
+
+Generalized Pauli ``Z`` operator for the given `N` level system. This is also
+called the clock matrix, and generalizes the qubit pauli matrices to qudits,
+while preserving the operators being unitary. A different generalization is
+given by the angular momentum operators which preserves the operators being
+Hermitian and is implemented for the `SpinBasis` (see [`sigmaz`](@ref)).
+
+Powers of `pauliz` together with [`paulix`](@ref) form a complete, orthornormal
+(under Hilbert–Schmidt norm) operator basis.
+
+Returns `Z^pow`.
+"""
+function pauliz(::Type{T}, b::NLevelBasis, pow=1) where T
+    N = length(b)
+    ω = exp(2π*1im*pow/N)
+    SparseOperator(b, spdiagm(0 => T[ω^n for n=1:N]))
+end
+pauliz(b::NLevelBasis, pow=1) = pauliz(ComplexF64, b, pow)
+
+"""
+    pauliy([T=ComplexF64,] b::NLevelBasis)
+
+Generalized Pauli ``Y`` operator for the given `N` level system.
+
+Returns `Y^pow = ω^pow X^pow Z^pow` where `ω = ω = exp(2π*1im*pow/N)` and
+`N = length(b)` when odd or `N = 2*length(b)` when even. This is due to the
+ order of the generalized Pauli group in even versus odd dimensions.
+
+See [`paulix`](@ref) and [`pauliz`](@ref) for more details.
+"""
+function pauliy(::Type{T}, b::NLevelBasis, pow=1) where T
+    N = length(b)
+    N = N%2 == 0 ? 2N : N
+    ω = exp(2π*1im*pow/N)
+    exp(2π*1im*pow/N)*paulix(T,b,pow)*pauliz(T,b,pow)
+end
+
+pauliy(b::NLevelBasis, pow=1) = pauliy(ComplexF64,b)

src/spin.jl

@@ -3,7 +3,13 @@ import QuantumInterface: SpinBasis
 """
     sigmax([T=ComplexF64,] b::SpinBasis)
 
-Pauli ``σ_x`` operator for the given Spin basis.
+Angular momentum ``σ_x`` operator for the given `SpinBasis`. For
+`SpinBasis(1//2)` this is the Pauli ``X`` operator while for higher spins this
+is the Hermitian operator that gives the quantized angular momentum observable
+along the ``x`` direction.
+
+See [`paulix`](@ref) for the generalization of the qubit pauli operators to
+qudits while preserving them being unitary instead of Hermitian.
 """
 function sigmax(::Type{T}, b::SpinBasis) where T
     N = length(b)
@@ -16,7 +22,13 @@ sigmax(b::SpinBasis) = sigmax(ComplexF64,b)
 """
     sigmay([T=ComplexF64,] b::SpinBasis)
 
-Pauli ``σ_y`` operator for the given Spin basis.
+Angular momentum ``σ_y`` operator for the given `SpinBasis`. For
+`SpinBasis(1//2)` this is the Pauli ``Y`` operator while for higher spins this
+is the Hermitian operator that gives the quantized angular momentum observable
+along the ``y`` direction.
+
+See [`pauliy`](@ref) for the generalization of the qubit pauli operators to
+qudits while preserving them being unitary instead of Hermitian.
 """
 function sigmay(::Type{T}, b::SpinBasis) where T
     N = length(b)
@@ -29,7 +41,13 @@ sigmay(b::SpinBasis) = sigmay(ComplexF64,b)
 """
     sigmaz([T=ComplexF64,] b::SpinBasis)
 
-Pauli ``σ_z`` operator for the given Spin basis.
+Angular momentum ``σ_z`` operator for the given `SpinBasis`. For
+`SpinBasis(1//2)` this is the Pauli ``Z`` operator while for higher spins this
+is the Hermitian operator that gives the quantized angular momentum observable
+along the ``z`` direction.
+
+See [`pauliz`](@ref) for the generalization of the qubit pauli operators to
+qudits while preserving them being unitary instead of Hermitian.
 """
 function sigmaz(::Type{T}, b::SpinBasis) where T
     N = length(b)

test/test_nlevel.jl

@@ -4,6 +4,8 @@ using LinearAlgebra
 
 @testset "nlevel" begin
 
+D(op1::AbstractOperator, op2::AbstractOperator) = abs(tracedistance_nh(dense(op1), dense(op2)))
+
 N = 3
 b = NLevelBasis(N)
 
@@ -26,4 +28,52 @@ b = NLevelBasis(N)
 @test norm(dagger(nlevelstate(b, 1))*nlevelstate(b, 2)) == 0.
 @test norm(dagger(nlevelstate(b, 1))*transition(b, 1, 2)*nlevelstate(b, 2)) == 1.
 
+for N=[2, 3, 4, 5]
+    b = NLevelBasis(N)
+    I = identityoperator(b)
+    Zero = SparseOperator(b)
+    px = paulix(b)
+    pz = pauliz(b)
+
+    @test 1e-14 > abs(tr(px))
+    @test 1e-14 > abs(tr(pz))
+    @test 1e-14 > D(px^N, I)
+    @test 1e-14 > D(pz^N, I)
+    for m=2:N
+        @test 1e-14 > D(px^m, paulix(b,m))
+        @test 1e-14 > D(pz^m, pauliz(b,m))
+    end
+
+    for m=1:N,n=1:N
+        @test 1e-14 > D(px^m * pz^n, exp(2π*1im*m*n/N) * pz^n * px^m)
+    end
+end
+
+
+# Test special relations for qubit pauli
+b = NLevelBasis(2)
+I = identityoperator(b)
+Zero = SparseOperator(b)
+px = paulix(b)
+py = pauliy(b)
+pz = pauliz(b)
+
+antikommutator(x, y) = x*y + y*x
+
+@test 1e-14 > D(antikommutator(px, px), 2*I)
+@test 1e-14 > D(antikommutator(px, py), Zero)
+@test 1e-14 > D(antikommutator(px, pz), Zero)
+@test 1e-14 > D(antikommutator(py, px), Zero)
+@test 1e-14 > D(antikommutator(py, py), 2*I)
+@test 1e-14 > D(antikommutator(py, pz), Zero)
+@test 1e-14 > D(antikommutator(pz, px), Zero)
+@test 1e-14 > D(antikommutator(pz, py), Zero)
+@test 1e-14 > D(antikommutator(pz, pz), 2*I)
+
+# Test if involutory for spin 1/2
+@test 1e-14 > D(px*px, I)
+@test 1e-14 > D(py*py, I)
+@test 1e-14 > D(pz*pz, I)
+@test 1e-14 > D(-1im*px*py*pz, I)
+
 end # testset