Low rank master equation
In this tutorial, we will show how to solve the master equation using the low-rank method. For a detailed explaination 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
\[\hat{H} = \frac{J_x}{2} \sum_{\langle i,j \rangle} \sigma_i^x \sigma_j^x + \frac{J_y}{2} \sum_{\langle i,j \rangle} \sigma_i^y \sigma_j^y + \frac{J_z}{2} \sum_{\langle i,j \rangle} \sigma_i^z \sigma_j^z - \sum_i h_i \sigma_i^z + h_x \sum_i \sigma_i^x + h_y \sum_i \sigma_i^y + h_z \sum_i \sigma_i^z,\]
where the sums are over nearest neighbors, and the collapse operators are given by
\[c_i = \sqrt{\gamma} \sigma_i^-.\]
We start by importing the packages
using QuantumToolbox
using CairoMakie
CairoMakie.enable_only_mime!(MIME"image/svg+xml"())Define lattice
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
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.
ϕ = 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] = SingleSiteOperator(sigmap(), j, latt) * ϕ[1])
end
for k in 1:N_modes-1
for l in k+1:N_modes
global i += 1
i <= M && (ϕ[i] = SingleSiteOperator(sigmap(), k, latt) * SingleSiteOperator(sigmap(), l, latt) * ϕ[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
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
# 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->SingleSiteOperator(sigmax(), i, latt), +, 1:latt.N)
Sy = mapreduce(i->SingleSiteOperator(sigmay(), i, latt), +, 1:latt.N)
Sz = mapreduce(i->SingleSiteOperator(sigmaz(), i, latt), +, 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
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
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
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
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