Add SU(2) phase space tools (#223)

* I implemented coherent spin states

* add SU(2) phase space tools

* Some efficiency and style changes

* Fix mistake in wignersu2

* tidy up SU(2) phase space tools

* further tidy up

* even further tidy up

* even further tidy up with the test

* even further tidy up with the test XD

* Small style changes

* Fix required packages to v0.6 compatible releases

* optimise tests

REQUIRE

@@ -4,3 +4,5 @@ OrdinaryDiffEq 3.19.1
 DiffEqCallbacks 1.1
 StochasticDiffEq 4.4.5
 RecursiveArrayTools
+GSL 0.3.6
+WignerSymbols 0.1.1

src/QuantumOptics.jl

@@ -24,7 +24,7 @@ export bases, Basis, GenericBasis, CompositeBasis, basis,
         nlevel, NLevelBasis, transition, nlevelstate,
         manybody, ManyBodyBasis, fermionstates, bosonstates,
                 manybodyoperator, onebodyexpect, occupation,
-        phasespace, qfunc, wigner, coherentspinstate,
+        phasespace, qfunc, wigner, coherentspinstate, qfuncsu2, wignersu2,
         metrics, tracenorm, tracenorm_h, tracenorm_nh,
                 tracedistance, tracedistance_h, tracedistance_nh,
                 entropy_vn, fidelity, ptranspose, PPT,

src/phasespace.jl

@@ -1,9 +1,11 @@
 module phasespace
 
-export qfunc, wigner, coherentspinstate
+export qfunc, wigner, coherentspinstate, qfuncsu2, wignersu2
 
 using ..bases, ..states, ..operators, ..operators_dense, ..fock, ..spin
 
+import WignerSymbols: clebschgordan
+import GSL: sf_legendre_sphPlm
 
 """
     qfunc(a, α)
@@ -204,21 +206,216 @@ function _clenshaw(L::Int, abs2_2α::Float64, ρ::Matrix{Complex128})
     end
 end
 
-############################ coherent spin state ##############################
-function coherentspinstate(b::SpinBasis, theta::Number, phi::Number,
-        result = Ket(b, Vector{Complex128}(length(b))))
-    #theta = real(theta)
-    #phi = real(phi)
+"""
+    coherentspinstate(b::SpinBasis, θ::Real, ϕ::Real)
+
+A coherent spin state |θ,ϕ⟩ is analogous to the coherent state of the linear harmonic
+oscillator. Coherent spin states represent a collection of identical two-level
+systems and can be described by two angles θ and ϕ (although this
+parametrization is not unique), similarly to a qubit on the
+Bloch sphere.
+"""
+function coherentspinstate(b::SpinBasis, theta::Real, phi::Real,
+    result = Ket(b, Vector{Complex128}(length(b))))
     data = result.data
-    N = length(b)-1
+    N = BigInt(length(b)-1)
+    sinth = sin(0.5theta)
+    costh = cos(0.5theta)
+    expphi = exp(0.5im*phi)
+    expphi_con = conj(expphi)
     @inbounds for n=0:N
-        data[n+1] =
-        sqrt(factorial(BigInt(N))/(factorial(BigInt(n)) *
-        factorial(BigInt(N-n)))) *
-        (sin(theta/2.)*exp(1im*phi/2.))^n *
-        (cos(theta/2.)*exp(-1im*phi/2.))^(N-n)
+        data[n+1] = sqrt(binomial(N, n)) * (sinth*expphi)^n * (costh*expphi_con)^(N-n)
+    end
+    return result
+end
+
+"""
+    qfuncsu2(ket,Ntheta;Nphi=2Ntheta)
+    qfuncsu2(rho,Ntheta;Nphi=2Ntheta)
+
+Husimi Q SU(2) representation ``⟨θ,ϕ|ρ|θ,ϕ⟩/π`` for the given state.
+
+The function calculates the SU(2) Husimi representation of a state on the
+generalised bloch sphere (0 < θ < π and 0 < ϕ < 2 π) with a given
+resolution `(Ntheta, Nphi)`.
+
+    qfuncsu2(rho,θ,ϕ)
+    qfuncsu2(ket,θ,ϕ)
+
+This version calculates the Husimi Q SU(2) function at a position given by θ and ϕ.
+"""
+function qfuncsu2(psi::Ket, Ntheta::Int; Nphi::Int=2Ntheta)
+    b = psi.basis
+    psi_bra_data = psi.data'
+    lb = float(b.spinnumber)
+    @assert isa(b, SpinBasis)
+    result = Array{Float64}(Ntheta,Nphi)
+    @inbounds  for i = 0:Ntheta-1, j = 0:Nphi-1
+        result[i+1,j+1] = (2*lb+1)/(4pi)*abs2(psi_bra_data*coherentspinstate(b,pi-i*pi/(Ntheta-1),j*2pi/(Nphi-1)-pi).data)
+    end
+    return result
+end
+
+function qfuncsu2(rho::DenseOperator, Ntheta::Int; Nphi::Int=2Ntheta)
+    b = basis(rho)
+    lb = float(b.spinnumber)
+    @assert isa(b, SpinBasis)
+    result = Array{Float64}(Ntheta,Nphi)
+    @inbounds  for i = 0:Ntheta-1, j = 0:Nphi-1
+        c = coherentspinstate(b,pi-i*1pi/(Ntheta-1),j*2pi/(Nphi-1)-pi)
+        result[i+1,j+1] = abs((2*lb+1)/(4pi)*c.data'*rho.data*c.data)
     end
     return result
 end
 
+function qfuncsu2(psi::Ket, theta::Real, phi::Real)
+    b = basis(psi)
+    psi_bra_data = psi.data'
+    lb = float(b.spinnumber)
+    @assert isa(b, SpinBasis)
+    result = (2*lb+1)/(4pi)*abs2(psi_bra_data*coherentspinstate(b,theta,phi).data)
+    return result
+end
+
+function qfuncsu2(rho::DenseOperator, theta::Real, phi::Real)
+    b = basis(rho)
+    lb = float(b.spinnumber)
+    @assert isa(b, SpinBasis)
+    c = coherentspinstate(b,theta,phi)
+    result = abs((2*lb+1)/(4pi)*c.data'*rho.data*c.data)
+    return result
+end
+
+"""
+    wignersu2(ket,Ntheta;Nphi=2Ntheta)
+    wignersu2(rho,Ntheta;Nphi=2Ntheta)
+
+Wigner SU(2) representation for the given state with a resolution `(Ntheta, Nphi)`.
+
+The function calculates the SU(2) Wigner representation of a state on the
+generalised bloch sphere (0 < θ < π and 0 < ϕ < 2 π) with a given resolution by
+decomposing the state into the basis of spherical harmonics.
+
+    wignersu2(rho,θ,ϕ)
+    wignersu2(ket,θ,ϕ)
+
+This version calculates the Wigner SU(2) function at a position given by θ and ϕ
+"""
+function wignersu2(rho::DenseOperator, theta::Real, phi::Real)
+
+    N = length(basis(rho))-1
+
+    ### Tensor generation ###
+    BandT = Array{Vector{Float64}}(N,N+1)
+    BandT[1,1] = collect(linspace(-N/2, N/2, N+1))
+    BandT[1,2] = -collect(sqrt.(linspace(1, N, N)).*sqrt.(linspace((N)/2, 1/2, N)))
+    BandT[2,1] = clebschgordan(1,0,1,0,2,0)*BandT[1,1].*BandT[1,1] -
+        clebschgordan(1,-1,1,1,2,0)*[zeros(N+1-length(BandT[1,2])); BandT[1,2].*BandT[1,2]] -
+        clebschgordan(1,1,1,-1,2,0)*[BandT[1,2].*BandT[1,2]; zeros(N+1-length(BandT[1,2]))]
+    BandT[2,2] = clebschgordan(1,0,1,1,2,1)BandT[1,1][1:N].*BandT[1,2]+
+        clebschgordan(1,1,1,0,2,1)*BandT[1,2][1:N].*BandT[1,1][2:end]
+    BandT[2,3] = BandT[1,2][1:N+1-(2)].*BandT[1,2][2:end]
+
+    @inbounds for S=2:N-1
+        BandT[S+1,1] = clebschgordan(1,0,S,0,S+1,0)*BandT[1,1].*BandT[S,1] -
+            [zeros(N+1-length(BandT[1,2])); clebschgordan(1,-1,S,1,S+1,0)*BandT[1,2].*BandT[S,2]] -
+            clebschgordan(1,1,S,-1,S+1,0)*[BandT[1,2].*BandT[S,2]; zeros(N+1-length(BandT[1,2]))]
+        BandT[S+1,S+1] = clebschgordan(1,0,S,S,S+1,S)BandT[1,1][1:N+1-S].*BandT[S,S+1]+
+            clebschgordan(1,1,S,S-1,S+1,S)*BandT[1,2][1:N+1-S].*BandT[S,S][2:end]
+        BandT[S+1,S+2] = BandT[1,2][1:N+1-(S+1)].*BandT[S,S+1][2:end]
+        for  M=1:S-1
+            BandT[S+1,M+1] = clebschgordan(1, 0, S, M, S+1,M)*BandT[1,1][1:N+1-M].*BandT[S,M+1] +
+                clebschgordan(1,1,S,M-1,S+1,M)*BandT[1,2][1:N+1-M].*BandT[S,M][2:end] -
+                clebschgordan(1,-1,S,M+1,S+1,M)*[zeros(1); BandT[1,2][1:N-M].*BandT[S,M+2][1:N-M]]
+        end
+
+    end
+
+    NormT = zeros(N)
+    @inbounds for S = 1:N
+        NormT[S] = sum(BandT[S,1].^2)
+    end
+
+    @inbounds for S = 1:N, M = 0:S
+        BandT[S, M + 1] =  BandT[S, M + 1]/sqrt(NormT[S])
+    end
+
+    ### State decomposition ###
+    c = rho.data;
+    EVT = Array{Complex128}(N,N+1)
+    @inbounds for S = 1:N, M = 0:S
+        EVT[S,M+1] = conj(sum(BandT[S,M+1].*diag(c,M)))
+    end
+
+    wignermap = _wignersu2int(N,theta,phi, EVT)
+    return wignermap*sqrt((N+1)/(4pi))
+end
+
+function wignersu2(rho::DenseOperator, Ntheta::Int; Nphi::Int=2Ntheta)
+
+    N = length(basis(rho))-1
+
+    ### Tensor generation ###
+    BandT = Array{Vector{Float64}}(N,N+1)
+    BandT[1,1] = collect(linspace(-N/2, N/2, N+1))
+    BandT[1,2] = -collect(sqrt.(linspace(1, N, N)).*sqrt.(linspace((N)/2, 1/2, N)))
+    BandT[2,1] = clebschgordan(1,0,1,0,2,0)*BandT[1,1].*BandT[1,1] -
+        clebschgordan(1,-1,1,1,2,0)*[zeros(N+1-length(BandT[1,2])); BandT[1,2].*BandT[1,2]] -
+        clebschgordan(1,1,1,-1,2,0)*[BandT[1,2].*BandT[1,2]; zeros(N+1-length(BandT[1,2]))]
+    BandT[2,2] = clebschgordan(1,0,1,1,2,1)BandT[1,1][1:N].*BandT[1,2]+
+        clebschgordan(1,1,1,0,2,1)*BandT[1,2][1:N].*BandT[1,1][2:end]
+    BandT[2,3] = BandT[1,2][1:N+1-(2)].*BandT[1,2][2:end]
+
+    @inbounds for S=2:N-1
+        BandT[S+1,1] = clebschgordan(1,0,S,0,S+1,0)*BandT[1,1].*BandT[S,1] -
+            [zeros(N+1-length(BandT[1,2])); clebschgordan(1,-1,S,1,S+1,0)*BandT[1,2].*BandT[S,2]] -
+            clebschgordan(1,1,S,-1,S+1,0)*[BandT[1,2].*BandT[S,2]; zeros(N+1-length(BandT[1,2]))]
+        BandT[S+1,S+1] = clebschgordan(1,0,S,S,S+1,S)BandT[1,1][1:N+1-S].*BandT[S,S+1]+
+            clebschgordan(1,1,S,S-1,S+1,S)*BandT[1,2][1:N+1-S].*BandT[S,S][2:end]
+        BandT[S+1,S+2] = BandT[1,2][1:N+1-(S+1)].*BandT[S,S+1][2:end]
+        for  M=1:S-1
+            BandT[S+1,M+1] = clebschgordan(1, 0, S, M, S+1,M)*BandT[1,1][1:N+1-M].*BandT[S,M+1] +
+                clebschgordan(1,1,S,M-1,S+1,M)*BandT[1,2][1:N+1-M].*BandT[S,M][2:end] -
+                clebschgordan(1,-1,S,M+1,S+1,M)*[0.0im; BandT[1,2][1:N-M].*BandT[S,M+2][1:N-M]]
+        end
+    end
+
+    NormT = zeros(N)
+    @inbounds for S = 1:N
+        NormT[S] = sum(BandT[S,1].^2)
+    end
+    @inbounds for S = 1:N, M = 0:S
+        BandT[S, M + 1] = BandT[S, M + 1]/sqrt(NormT[S])
+    end
+
+    ### State decomposition ###
+    c = rho.data;
+    EVT = Array{Complex128}(N,N+1)
+    @inbounds for S = 1:N, M = 0:S
+        EVT[S,M+1] = conj(sum(BandT[S,M+1].*diag(c,M)))
+    end
+
+    wignermap = Array{Float64}(Ntheta,Nphi)
+    @inbounds for i = 1:Ntheta, j = 1:Nphi
+        wignermap[i,j] = _wignersu2int(N,i*1pi/(Ntheta-1),j*2pi/(Nphi-1)-pi, EVT)
+    end
+    return wignermap*sqrt((N+1)/(4pi))
+end
+
+function _wignersu2int(N::Integer, theta::Real, phi::Real, EVT::Array{Complex128, 2})
+    UberBand = sqrt(1/(1+N))*ylm(0,0,theta,phi)
+    @inbounds for S = 1:N
+        @inbounds for M = 1:S
+            UberBand += 2*real(EVT[S,M+1]*conj(ylm(S,M,theta,phi)))
+        end
+        UberBand += EVT[S,1]*ylm(S,0,theta,phi)
+    end
+    UberBand
+end
+wignersu2(psi::Ket, args...) = wignersu2(dm(psi), args...)
+
+function ylm(l::Integer, m::Integer, theta::Real, phi::Real)
+    sf_legendre_sphPlm(l,m,cos(theta))*e^(1im*m*phi)
+end
+
 end #module

test/test_phasespace.jl

@@ -86,14 +86,42 @@ for (i,x)=enumerate(X), (j,y)=enumerate(Y)
 end
 
 # Test SU(2) phasespace
-b = SpinBasis(50)
+b = SpinBasis(5)
 theta = π*rand()
 phi =2π*rand()
 css = coherentspinstate(b,theta,phi)
+dmcss = dm(css)
+csssx = coherentspinstate(b,π/2,0)
+dmcsssx = dm(csssx)
+rs = randstate(b)
+dmrs = dm(rs)
 sx = expect(sigmax(b)/2,css); # eigenstate of jx operator
 sy = expect(sigmay(b)/2,css); # eigenstate of jy
 sz = expect(sigmaz(b)/2,css); # eigenstate of jz operator
 ssq = sx^2 + sy^2 + sz^2
 
+qsu2sx = qfuncsu2(csssx,theta,phi)
+qsu2sxdm = qfuncsu2(dmcsssx,theta,phi)
+res = 250
+costhetam = Array{Float64}(res,2*res)
+for i = 1:res, j = 1:2*res
+    costhetam[i,j] = sin(i*1pi/(res-1))
+end
+wsu2 = sum(wignersu2(rs,res).*costhetam)*(π/res)^2
+wsu2dm = sum(wignersu2(dmrs,res).*costhetam)*(π/res)^2
+qsu2 = sum(qfuncsu2(rs,res).*costhetam)*(π/res)^2
+qsu2dm = sum(qfuncsu2(dmrs,res).*costhetam)*(π/res)^2
+
 @test ssq ≈ float(b.spinnumber)^2
+
+@test isapprox(qsu2sxdm, (2*float(b.spinnumber)+1)/(4pi)*(0.5*(sin(theta)cos(phi)+1))^(2*float(b.spinnumber)),atol=1e-2)
+@test isapprox(qsu2sx, (2*float(b.spinnumber)+1)/(4pi)*(0.5*(sin(theta)cos(phi)+1))^(2*float(b.spinnumber)),atol=1e-2)
+@test isapprox(qsu2, 1.0, atol=1e-2)
+@test isapprox(qsu2dm, 1.0, atol=1e-2)
+
+@test isapprox(wignersu2(csssx,π/2,0), (4*float(b.spinnumber)+1)/(4pi), atol=1e-2)
+@test isapprox(wignersu2(dmcsssx,π/2,0), (4*float(b.spinnumber)+1)/(4pi),atol=1e-2)
+@test isapprox(wsu2, 1.0, atol=1e-2)
+@test isapprox(wsu2dm, 1.0, atol=1e-2)
+
 end # testset