Add generalized pauli operators (aka clock and shift matrices)

Author: akirakyle

src/QuantumOpticsBase.jl

@@ -55,7 +55,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,40 @@ 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.
+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.
+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)
+
+Pauli Y operator. Only defined for a two level system.
+"""
+function pauliy(::Type{T}, b::NLevelBasis) where T
+    length(b) == 2 || throw(ArgumentError("pauliy only defined for two level system"))
+    SparseOperator(b, spdiagm(-1 => T[-1im], 1 => T[1im]))
+end
+pauliy(b::NLevelBasis) = pauliy(ComplexF64,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