Merge pull request #158 from ytdHuang/dev/spin-j

Constructing general spin-j operators

docs/src/api.md

@@ -93,6 +93,13 @@ sigmam
 sigmax
 sigmay
 sigmaz
+jmat
+spin_Jx
+spin_Jy
+spin_Jz
+spin_Jm
+spin_Jp
+spin_J_set
 destroy
 create
 fdestroy

src/qobj/operators.jl

@@ -2,6 +2,7 @@
 Functions for generating (common) quantum operators.
 =#
 
+export jmat, spin_Jx, spin_Jy, spin_Jz, spin_Jm, spin_Jp, spin_J_set
 export sigmam, sigmap, sigmax, sigmay, sigmaz
 export destroy, create, eye, qeye, projection
 export fdestroy, fcreate
@@ -139,40 +140,181 @@ julia> fock(20, 4)' * a_d * fock(20, 3)
 """
 create(N::Int) = QuantumObject(spdiagm(-1 => Array{ComplexF64}(sqrt.(1:N-1))), Operator, [N])
 
+@doc raw"""
+    jmat(j::Real, which::Symbol)
+
+Generate higher-order Spin-`j` operators, where `j` is the spin quantum number and can be a non-negative integer or half-integer
+
+The parameter `which` specifies which of the following operator to return.
+- `:x`: ``S_x``
+- `:y`: ``S_y``
+- `:z`: ``S_z``
+- `:+`: ``S_+``
+- `:-`: ``S_-``
+
+Note that if the parameter `which` is not specified, returns a set of Spin-`j` operators: ``(S_x, S_y, S_z)``
+
+# Examples
+```
+julia> jmat(0.5, :x)
+Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
+2×2 SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
+     ⋅      0.5+0.0im
+ 0.5+0.0im      ⋅
+
+julia> jmat(0.5, :-)
+Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=false
+2×2 SparseMatrixCSC{ComplexF64, Int64} with 1 stored entry:
+     ⋅          ⋅    
+ 1.0+0.0im      ⋅
+```
+"""
+jmat(j::Real, which::Symbol) = jmat(j, Val(which))
+jmat(j::Real) = (jmat(j, Val(:x)), jmat(j, Val(:y)), jmat(j, Val(:z)))
+function jmat(j::Real, ::Val{:x})
+    N = 2 * j
+    ((floor(N) != (N)) || (j < 0)) &&
+        throw(ArgumentError("The spin quantum number (j) j must be a non-negative integer or half-integer."))
+
+    σ = _jm(j)
+    return QuantumObject((σ' + σ) / 2, Operator, [Int(N + 1)])
+end
+function jmat(j::Real, ::Val{:y})
+    N = 2 * j
+    ((floor(N) != (N)) || (j < 0)) &&
+        throw(ArgumentError("The spin quantum number (j) j must be a non-negative integer or half-integer."))
+
+    σ = _jm(j)
+    return QuantumObject((σ' - σ) / 2im, Operator, [Int(N + 1)])
+end
+function jmat(j::Real, ::Val{:z})
+    N = 2 * j
+    ((floor(N) != (N)) || (j < 0)) &&
+        throw(ArgumentError("The spin quantum number (j) j must be a non-negative integer or half-integer."))
+
+    return QuantumObject(_jz(j), Operator, [Int(N + 1)])
+end
+function jmat(j::Real, ::Val{:+})
+    N = 2 * j
+    ((floor(N) != (N)) || (j < 0)) &&
+        throw(ArgumentError("The spin quantum number (j) j must be a non-negative integer or half-integer."))
+
+    return QuantumObject(adjoint(_jm(j)), Operator, [Int(N + 1)])
+end
+function jmat(j::Real, ::Val{:-})
+    N = 2 * j
+    ((floor(N) != (N)) || (j < 0)) &&
+        throw(ArgumentError("The spin quantum number (j) j must be a non-negative integer or half-integer."))
+
+    return QuantumObject(_jm(j), Operator, [Int(N + 1)])
+end
+jmat(j::Real, ::Val{T}) where {T} = throw(ArgumentError("Invalid spin operator: $(T)"))
+
+function _jm(j::Real)
+    m = j:(-1):-j
+    data = sqrt.(j * (j + 1) .- m .* (m .- 1))[1:end-1]
+    return spdiagm(-1 => Array{ComplexF64}(data))
+end
+_jz(j::Real) = spdiagm(0 => Array{ComplexF64}(j .- (0:Int(2 * j))))
+
+@doc raw"""
+    spin_Jx(j::Real)
+
+``S_x`` operator for Spin-`j`, where `j` is the spin quantum number and can be a non-negative integer or half-integer
+
+See also [`jmat`](@ref).
+"""
+spin_Jx(j::Real) = jmat(j, Val(:x))
+
+@doc raw"""
+    spin_Jy(j::Real)
+
+``S_y`` operator for Spin-`j`, where `j` is the spin quantum number and can be a non-negative integer or half-integer
+
+See also [`jmat`](@ref).
+"""
+spin_Jy(j::Real) = jmat(j, Val(:y))
+
+@doc raw"""
+    spin_Jz(j::Real)
+
+``S_z`` operator for Spin-`j`, where `j` is the spin quantum number and can be a non-negative integer or half-integer
+
+See also [`jmat`](@ref).
+"""
+spin_Jz(j::Real) = jmat(j, Val(:z))
+
+@doc raw"""
+    spin_Jm(j::Real)
+
+``S_-`` operator for Spin-`j`, where `j` is the spin quantum number and can be a non-negative integer or half-integer
+
+See also [`jmat`](@ref).
+"""
+spin_Jm(j::Real) = jmat(j, Val(:-))
+
+@doc raw"""
+    spin_Jp(j::Real)
+
+``S_+`` operator for Spin-`j`, where `j` is the spin quantum number and can be a non-negative integer or half-integer
+
+See also [`jmat`](@ref).
+"""
+spin_Jp(j::Real) = jmat(j, Val(:+))
+
+@doc raw"""
+    spin_J_set(j::Real)
+
+A set of Spin-`j` operators ``(S_x, S_y, S_z)``, where `j` is the spin quantum number and can be a non-negative integer or half-integer
+
+Note that this functions is same as `jmat(j)`. See also [`jmat`](@ref).
+"""
+spin_J_set(j::Real) = jmat(j)
+
 @doc raw"""
     sigmap()
 
 Pauli ladder operator ``\hat{\sigma}_+ = \hat{\sigma}_x + i \hat{\sigma}_y``.
+
+See also [`jmat`](@ref).
 """
-sigmap() = destroy(2)
+sigmap() = jmat(0.5, Val(:+))
 
 @doc raw"""
     sigmam()
 
 Pauli ladder operator ``\hat{\sigma}_- = \hat{\sigma}_x - i \hat{\sigma}_y``.
+
+See also [`jmat`](@ref).
 """
-sigmam() = create(2)
+sigmam() = jmat(0.5, Val(:-))
 
 @doc raw"""
     sigmax()
 
 Pauli operator ``\hat{\sigma}_x = \hat{\sigma}_- + \hat{\sigma}_+``.
+
+See also [`jmat`](@ref).
 """
-sigmax() = sigmam() + sigmap()
+sigmax() = rmul!(jmat(0.5, Val(:x)), 2)
 
 @doc raw"""
     sigmay()
 
 Pauli operator ``\hat{\sigma}_y = i \left( \hat{\sigma}_- - \hat{\sigma}_+ \right)``.
+
+See also [`jmat`](@ref).
 """
-sigmay() = 1im * (sigmam() - sigmap())
+sigmay() = rmul!(jmat(0.5, Val(:y)), 2)
 
 @doc raw"""
     sigmaz()
 
 Pauli operator ``\hat{\sigma}_z = \comm{\hat{\sigma}_+}{\hat{\sigma}_-}``.
+
+See also [`jmat`](@ref).
 """
-sigmaz() = sigmap() * sigmam() - sigmam() * sigmap()
+sigmaz() = rmul!(jmat(0.5, Val(:z)), 2)
 
 @doc raw"""
     eye(N::Int; type=OperatorQuantumObject, dims=[N])

test/states_and_operators.jl

@@ -18,6 +18,35 @@
     @test all(eigenenergies(ρ_B) .>= 0)
     @test_throws DimensionMismatch rand_dm(4, dims = [2])
 
+    # Pauli matrices and general Spin-j operators
+    J0 = Qobj(spdiagm(0 => [0.0im]))
+    Jx, Jy, Jz = spin_J_set(0.5)
+    @test spin_Jx(0) == spin_Jy(0) == spin_Jz(0) == spin_Jp(0) == spin_Jm(0) == J0
+    @test sigmax() ≈ 2 * spin_Jx(0.5) ≈ 2 * Jx ≈ Qobj(sparse([0.0im 1; 1 0]))
+    @test sigmay() ≈ 2 * spin_Jy(0.5) ≈ 2 * Jy ≈ Qobj(sparse([0.0im -1im; 1im 0]))
+    @test sigmaz() ≈ 2 * spin_Jz(0.5) ≈ 2 * Jz ≈ Qobj(sparse([1 0.0im; 0 -1]))
+    @test sigmap() ≈ spin_Jp(0.5) ≈ Qobj(sparse([0.0im 1; 0 0]))
+    @test sigmam() ≈ spin_Jm(0.5) ≈ Qobj(sparse([0.0im 0; 1 0]))
+    for which in [:x, :y, :z, :+, :-]
+        @test jmat(2.5, which).dims == [6] # 2.5 * 2 + 1 = 6
+
+        for j_wrong in [-1, 8.7]  # Invalid j
+            @test_throws ArgumentError jmat(j_wrong, which)
+        end
+    end
+    @test_throws ArgumentError jmat(0.5, :wrong)
+
+    # test commutation relations for general Spin-j operators
+    # the negative signs follow the rule of Levi-Civita symbol (ϵ_ijk)
+    j = 2
+    Jx, Jy, Jz = spin_J_set(j)
+    @test commutator(Jx, Jy) ≈ 1im * Jz
+    @test commutator(Jy, Jz) ≈ 1im * Jx
+    @test commutator(Jz, Jx) ≈ 1im * Jy
+    @test commutator(Jx, Jz) ≈ -1im * Jy
+    @test commutator(Jy, Jx) ≈ -1im * Jz
+    @test commutator(Jz, Jy) ≈ -1im * Jx
+
     # test commutation relations for fermionic creation and annihilation operators
     sites = 4
     SIZE = 2^sites