Low rank master equation
In this tutorial, we will show how to solve the master equation using the low-rank method. For a detailed explanation of the method, we recommend to read the article (Gravina and Savona, 2024).
As a test, we will consider the dissipative Ising model with a transverse field. The Hamiltonian is given by
where the sums are over nearest neighbors, and the collapse operators are given by
We start by importing the packages
julia
using QuantumToolbox
using CairoMakie
CairoMakie.enable_only_mime!(MIME"image/svg+xml"())Define lattice
julia
Nx, Ny = 2, 3
latt = Lattice(Nx = Nx, Ny = Ny)Lattice{Int64, LinearIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}}(2, 3, 6, [1 3 5; 2 4 6], CartesianIndices((2, 3)))Define lr-space dimensions
julia
N_cut = 2 # Number of states of each mode
N_modes = latt.N # Number of modes
N = N_cut^N_modes # Total number of states
M = latt.N + 1 # Number of states in the LR basis7Define lr states. Take as initial state all spins up. All other N states are taken as those with miniman Hamming distance to the initial state.
julia
ϕ = Vector{QuantumObject{Vector{ComplexF64},KetQuantumObject,M-1}}(undef, M)
ϕ[1] = kron(fill(basis(2, 1), N_modes)...)
i = 1
for j in 1:N_modes
global i += 1
i <= M && (ϕ[i] = MultiSiteOperator(latt, j=>sigmap()) * ϕ[1])
end
for k in 1:N_modes-1
for l in k+1:N_modes
global i += 1
i <= M && (ϕ[i] = MultiSiteOperator(latt, k=>sigmap(), l=>sigmap()) * ϕ[1])
end
end
for i in i+1:M
ϕ[i] = QuantumObject(rand(ComplexF64, size(ϕ[1])[1]), dims = ϕ[1].dims)
normalize!(ϕ[i])
endDefine the initial state
julia
z = hcat(get_data.(ϕ)...)
B = Matrix(Diagonal([1 + 0im; zeros(M - 1)]))
S = z' * z # Overlap matrix
B = B / tr(S * B) # Normalize B
ρ = QuantumObject(z * B * z', dims = ntuple(i->N_cut, Val(N_modes))); # Full density matrixQuantum Object: type=Operator dims=[2, 2, 2, 2, 2, 2] size=(64, 64) ishermitian=true
64×64 Matrix{ComplexF64}:
0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
⋮ ⋱
0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0imDefine the Hamiltonian and collapse operators
julia
# Define Hamiltonian and collapse operators
Jx = 0.9
Jy = 1.04
Jz = 1.0
hx = 0.0
hy = 0.0
hz = 0.0
γ = 1
Sx = mapreduce(i->MultiSiteOperator(latt, i=>sigmax()), +, 1:latt.N)
Sy = mapreduce(i->MultiSiteOperator(latt, i=>sigmay()), +, 1:latt.N)
Sz = mapreduce(i->MultiSiteOperator(latt, i=>sigmaz()), +, 1:latt.N)
H, c_ops = DissipativeIsing(Jx, Jy, Jz, hx, hy, hz, γ, latt; boundary_condition = Val(:periodic_bc), order = 1)
e_ops = (Sx, Sy, Sz)
tl = range(0, 10, 100)Full evolution
julia
sol_me = mesolve(H, ρ, tl, c_ops; e_ops = [e_ops...]);
Strue = entropy_vn(sol_me.states[end], base=2) / latt.N0.29324728062546046Low Rank evolution
Define functions to be evaluated during the low-rank evolution
julia
function f_purity(p, z, B)
N = p.N
M = p.M
S = p.S
T = p.temp_MM
mul!(T, S, B)
return tr(T^2)
end
function f_trace(p, z, B)
N = p.N
M = p.M
S = p.S
T = p.temp_MM
mul!(T, S, B)
return tr(T)
end
function f_entropy(p, z, B)
C = p.A0
σ = p.Bi
mul!(C, z, sqrt(B))
mul!(σ, C', C)
return entropy_vn(Qobj(Hermitian(σ), type=Operator), base=2)
endf_entropy (generic function with 1 method)Define the options for the low-rank evolution
julia
opt = (err_max = 1e-3, p0 = 0.0, atol_inv = 1e-6, adj_condition = "variational", Δt = 0.0);
sol_lr = lr_mesolve(H, z, B, tl, c_ops; e_ops = e_ops, f_ops = (f_purity, f_entropy, f_trace), opt = opt);TimeEvolutionLRSol{Vector{Float64}, Vector{QuantumObject{Matrix{ComplexF64}, OperatorQuantumObject, 6}}, Matrix{ComplexF64}, SciMLBase.ReturnCode.T, OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Float64, Float64, Vector{Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}}, Vector{Int64}}([10.0], QuantumObject{Matrix{ComplexF64}, OperatorQuantumObject, 6}[Quantum Object: type=Operator dims=[2, 2, 2, 2, 2, 2] size=(64, 64) ishermitian=false
64×64 Matrix{ComplexF64}:
0.00024987+7.83654e-20im … 0.00606777-0.00919129im
-2.24542e-9+1.82222e-9im 1.60118e-8+1.29703e-7im
-2.0359e-9+1.22481e-9im -1.03291e-8+1.04889e-7im
1.71251e-5-0.000450592im -0.0175574-0.0116535im
-1.99883e-9+1.25438e-9im -8.35726e-9+1.04033e-7im
-6.74301e-6-0.000441917im … -0.0177936-0.0103975im
-1.00213e-5-0.000439887im -0.0177892-0.0102042im
2.98276e-9+5.04937e-9im 2.77433e-7+4.19569e-9im
-2.24268e-9+1.4597e-9im -2.0363e-10+1.23572e-7im
-1.00215e-5-0.000439887im -0.0177892-0.0102042im
⋮ ⋱
9.46533e-9-1.70044e-8im … -4.9143e-7-8.35857e-7im
3.84988e-9+3.56795e-9im 2.49017e-7-8.75575e-8im
-0.00116913+0.000295455im -0.0176962+0.059073im
-0.00119314+0.000260925im -0.0199029+0.0591602im
1.20231e-8-1.84941e-8im -6.62026e-7-1.04809e-6im
-0.0013331+0.000124241im … -0.0298358+0.061477im
1.25907e-9-1.58089e-8im -6.96845e-7-4.31998e-7im
1.90703e-9-5.534e-9im -8.03311e-8-2.25274e-7im
0.00606777+0.00919129im 0.583662+8.67362e-18im], ComplexF64[0.0 + 0.0im 0.0 + 0.0im … -1.3598090128702001e-5 - 4.336808689942018e-19im -1.3360101014645888e-5 - 2.0599841277224584e-18im; 0.0 + 0.0im 0.0 + 0.0im … 2.1364363166417194e-5 - 1.3417001884508117e-18im 2.0895616506237766e-5 + 5.692061405548898e-19im; -6.0 + 0.0im -5.987982183309875 + 1.4836204593182117e-20im … -4.564252066985382 - 1.5178830414797062e-17im -4.5643347215734815 + 2.1250362580715887e-17im], ComplexF64[1.0 + 0.0im 0.9995913405802064 - 4.591774807899561e-41im … 0.5162555177271383 - 4.042666582026977e-18im 0.5162788335256082 - 8.13819405667546e-18im; -0.0 + 0.0im 0.0031973769163587424 + 0.0im … 1.7339220302154035 + 0.0im 1.733838373076687 + 0.0im; 1.0 + 0.0im 1.0000000000000002 + 0.0im … 0.9999999999999999 + 4.0389678347315804e-28im 0.9999999999999999 - 4.480729941655347e-28im], SciMLBase.ReturnCode.Success, OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-8, 1.0e-6, Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.0024566430410880675 - 0.005902277956175032im -9.36276574483052e-6 + 4.5088367031826844e-6im … -3.843286776774304e-5 - 0.00019439199885410747im -2.0065796786539022e-5 + 1.722148370682773e-6im; 6.45440555759243e-5 + 5.5170447711710806e-5im -0.012179068121601846 - 0.005169411608331753im … -0.0213052553738994 - 0.004824570153473949im 0.004775610170541168 + 0.038045509659198284im; … ; -0.00016236411363600132 - 0.00010172578020124223im 0.0016796356375573698 - 0.010460241152385651im … 0.006427489256206535 + 0.02854169179756393im 0.038808337153689275 - 0.012020737452856704im; 0.9758709758035521 - 0.006272785828394867im 4.406375431201654e-6 - 2.1061301158014263e-6im … 6.811993690306314e-5 + 3.62475200437046e-5im 3.899992551180602e-6 - 2.1600589372443332e-5im]], Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.6675295551781188 + 1.9209929114941938e-27im -4.132793205296643e-7 + 2.114314610411267e-6im … 1.4491452892522766e-5 - 1.168014098213419e-5im 2.219279108776238e-6 - 2.05814051598504e-6im; -4.132793205296643e-7 - 2.114314610411267e-6im 0.04057260672868191 + 1.1675841065320728e-28im … 8.131899538213679e-6 + 3.652762630015909e-7im -7.569673608243792e-6 + 4.7468018446205837e-7im; … ; 1.4491452892522766e-5 + 1.168014098213419e-5im 8.131899538213679e-6 - 3.652762630015909e-7im … 0.0003199637192733957 + 9.207802589287645e-31im 4.773126500374891e-6 + 6.5021573957207665e-6im; 2.219279108776238e-6 + 2.05814051598504e-6im -7.569673608243792e-6 - 4.7468018446205837e-7im … 4.773126500374891e-6 - 6.5021573957207665e-6im 0.0003010755324602631 + 8.664245038950806e-31im]], [36])Plot the results
julia
m_me = real(sol_me.expect[3, :]) / Nx / Ny
m_lr = real(sol_lr.expect[3, :]) / Nx / Ny
fig = Figure(size = (500, 350), fontsize = 15)
ax = Axis(fig[1, 1], xlabel = L"\gamma t", ylabel = L"M_{z}", xlabelsize = 20, ylabelsize = 20)
lines!(ax, tl, m_lr, label = L"LR $[M=M(t)]$", linewidth = 2)
lines!(ax, tl, m_me, label = "Fock", linewidth = 2, linestyle = :dash)
axislegend(ax, position = :rb)
ax2 = Axis(fig[1, 2], xlabel = L"\gamma t", ylabel = "Value", xlabelsize = 20, ylabelsize = 20)
lines!(ax2, tl, 1 .- real(sol_lr.fexpect[1, :]), label = L"$1-P$", linewidth = 2)
lines!(
ax2,
tl,
1 .- real(sol_lr.fexpect[3, :]),
label = L"$1-\mathrm{Tr}(\rho)$",
linewidth = 2,
linestyle = :dash,
color = :orange,
)
lines!(ax2, tl, real(sol_lr.fexpect[2, :]) / Nx / Ny, color = :blue, label = L"S", linewidth = 2)
hlines!(ax2, [Strue], color = :blue, linestyle = :dash, linewidth = 2, label = L"S^{\,\mathrm{true}}_{\mathrm{ss}}")
axislegend(ax2, position = :rb)
fig