Heisenberg picture master equation time evolution

src/QuantumOptics.jl

@@ -22,6 +22,7 @@ module timeevolution
     include("timeevolution_base.jl")
     include("time_dependent_operators.jl")
     include("master.jl")
+    include("master_heisenberg.jl")
     include("schroedinger.jl")
     include("mcwf.jl")
     include("bloch_redfield_master.jl")

src/master.jl

@@ -4,6 +4,7 @@
 Integrate the master equation with dmaster_h as derivative function.
 
 Further information can be found at [`master`](@ref).
+See [`master_h_heisenberg`](@ref) for the Heisenberg picture.
 """
 function master_h(tspan, rho0::Operator, H::AbstractOperator, J;
                 rates=nothing,

src/master_heisenberg.jl

@@ -0,0 +1,73 @@
+
+
+"""
+    timeevolution.master_h_heisenberg(tspan, A0, H, J; <keyword arguments>)
+
+Integrate the Heisenberg picture master equation with `dmaster_h_heisenberg!` as derivative function.
+See [`master_h`](@ref) for the Schrödinger picture.
+"""
+function master_h_heisenberg(tspan, A0::Operator, H::AbstractOperator, J;
+                           rates=nothing,
+                           Jdagger=dagger.(J),
+                           fout=nothing,
+                           kwargs...)
+    _check_const(H)
+    _check_const.(J)
+    _check_const.(Jdagger)
+    tspan, A0 = _promote_time_and_state(A0, H, J, tspan)
+    tmp = copy(A0)
+    dmaster_heisenberg_(t, A, dA) = dmaster_h_heisenberg!(dA, H, J, Jdagger, rates, A, tmp)
+    integrate_master(tspan, dmaster_heisenberg_, A0, fout; kwargs...)
+end
+
+
+
+"""
+    dmaster_h_heisenberg!(dA, H, J, Jdagger, rates::Nothing, A, dA_cache)
+
+Update `drho` according to a heisenberg picture Lindblad equation.
+See [`dmaster_h`](@ref).
+"""
+function dmaster_h_heisenberg!(dA, H, J, Jdagger, rates::Nothing, A, dA_cache)
+    QuantumOpticsBase.mul!(dA,H,A,eltype(A)(im),zero(eltype(A)))
+    QuantumOpticsBase.mul!(dA,A,H,-eltype(A)(im),one(eltype(A)))
+    for i=1:length(J)
+        QuantumOpticsBase.mul!(dA_cache,Jdagger[i],A)
+        QuantumOpticsBase.mul!(dA,dA_cache,J[i],true,true)
+
+        QuantumOpticsBase.mul!(dA_cache,Jdagger[i],J[i],eltype(A)(-0.5),zero(eltype(A)))
+        QuantumOpticsBase.mul!(dA,dA_cache,A,true,true)
+        QuantumOpticsBase.mul!(dA,A,dA_cache,true,true)
+    end
+    return dA
+end
+
+
+function dmaster_h_heisenberg!(dA, H, J, Jdagger, rates::AbstractVector, A, dA_cache)
+    QuantumOpticsBase.mul!(dA,H,A,eltype(A)(im),zero(eltype(A)))
+    QuantumOpticsBase.mul!(dA,A,H,-eltype(A)(im),one(eltype(A)))
+    for i=1:length(J)
+        QuantumOpticsBase.mul!(dA_cache,Jdagger[i],A, eltype(A)(rates[i]),zero(eltype(A)))
+        QuantumOpticsBase.mul!(dA,dA_cache,J[i],true,true)
+
+        QuantumOpticsBase.mul!(dA_cache,Jdagger[i],J[i],eltype(A)(-0.5)*rates[i],zero(eltype(A)))
+        QuantumOpticsBase.mul!(dA,dA_cache,A,true,true)
+        QuantumOpticsBase.mul!(dA,A,dA_cache,true,true)
+    end
+    return dA
+end
+
+
+function dmaster_h_heisenberg!(dA, H, J, Jdagger, rates::AbstractMatrix, A, dA_cache)
+    QuantumOpticsBase.mul!(dA,H,A,eltype(A)(im),zero(eltype(A)))
+    QuantumOpticsBase.mul!(dA,A,H,-eltype(A)(im),one(eltype(A)))
+    for j=1:length(J), i=1:length(J)
+        QuantumOpticsBase.mul!(dA_cache,Jdagger[i],A, eltype(A)(rates[i,j]),zero(eltype(A)))
+        QuantumOpticsBase.mul!(dA,dA_cache,J[j],true,true)
+
+        QuantumOpticsBase.mul!(dA_cache,Jdagger[i],J[j],eltype(A)(-0.5)*rates[i,j],zero(eltype(A)))
+        QuantumOpticsBase.mul!(dA,dA_cache,A,true,true)
+        QuantumOpticsBase.mul!(dA,A,dA_cache,true,true)
+    end
+    return dA
+end

test/test_timeevolution_master_heisenberg.jl

@@ -0,0 +1,72 @@
+using QuantumOptics
+using LinearAlgebra
+using Test
+
+
+# Tests for the master_h_heisenberg function
+# ==========================================
+# We check that time evolution of <Ztot> is the same in the Schrödinger and Heisenberg
+# pictures for a time reversal symmetry breaking Hamiltonian with non hermitian jump operators.
+
+
+@testset "master_h_heisenberg" begin
+
+N = 4
+tmax = 10.0
+dt = 0.1
+b = SpinBasis(1//2)
+basis = tensor([b for _=1:N]...)
+sx(i) = embed(basis, i, sigmax(b))
+sy(i) = embed(basis, i, sigmay(b))
+sz(i) = embed(basis, i, sigmaz(b))
+
+
+# time reversal symmetry breaking Hamiltonian
+function ising(N)
+    H = sum(sz(i) * sz(i+1) for i in 1:N-1)
+    H += -sum(sx(i) for i in 1:N)
+    H += -sum(sz(i) * sx(i+1)* sy(i+2) for i in 1:N-2)
+    return H
+end
+
+Ztot(N) = sum(sz(i) for i in 1:N-1)
+
+H = ising(N)
+
+jump_ops = [sz(i)-im*sy(i) for i in 1:N]
+
+ψ0 = tensor([spinup(b) for _=1:N]...)
+ρ0 = ψ0 ⊗ dagger(ψ0)
+times = 0:dt:tmax
+
+A0 = Ztot(N)
+
+fout_s(t, ρ) = real(expect(A0, ρ)) # shrödinger picture observable
+fout_h(t, A) = real(expect(A, ρ0)) # heisenberg picture observable
+
+n = length(jump_ops)
+
+# test with no rates
+tout_s, exp_val_s = timeevolution.master_h(times, ρ0, H, jump_ops; fout=fout_s)
+tout_h, exp_val_h = timeevolution.master_h_heisenberg(times, A0, H, jump_ops; fout=fout_h)
+@test norm(exp_val_s - exp_val_h) < 1e-5
+
+# test with diagonal rates
+rates = rand(n)*0.05
+tout_s, exp_val_s = timeevolution.master_h(times, ρ0, H, jump_ops; fout=fout_s, rates=rates)
+tout_h, exp_val_h = timeevolution.master_h_heisenberg(times, A0, H, jump_ops; fout=fout_h, rates=rates)
+@test norm(exp_val_s - exp_val_h) < 1e-5
+
+# test with rates matrix
+rates = rand(n, n)*0.05
+tout_s, exp_val_s = timeevolution.master_h(times, ρ0, H, jump_ops; fout=fout_s, rates=rates)
+tout_h, exp_val_h = timeevolution.master_h_heisenberg(times, A0, H, jump_ops; fout=fout_h, rates=rates)
+@test norm(exp_val_s - exp_val_h) < 1e-5
+
+# uncomment to plot the results
+# using Plots
+# p = plot(tout_s,exp_val_s)
+# plot!(p, tout_h, exp_val_h)
+# display(p)
+# sleep(100)
+end