Low rank master equation
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 = Nx * Ny + 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}}(undef, M)
ϕ[1] = kron(repeat([basis(2, 0)], N_modes)...)
global i = 1
for j in 1:N_modes
global i += 1
i <= M && (ϕ[i] = mb(sp, j, latt) * ϕ[1])
end
for k in 1:N_modes-1
for l in k+1:N_modes
global i += 1
i <= M && (ϕ[i] = mb(sp, k, latt) * mb(sp, 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(broadcast(x -> x.data, ϕ)...)
p0 = 0.0 # Population of the lr states other than the initial state
B = Matrix(Diagonal([1 + 0im; p0 * ones(M - 1)]))
S = z' * z # Overlap matrix
B = B / tr(S * B) # Normalize B
ρ = QuantumObject(z * B * z', dims = ones(Int, N_modes) * N_cut); # Full density matrixQuantum Object: type=Operator dims=[2, 2, 2, 2, 2, 2] size=(64, 64) ishermitian=true
64×64 Matrix{ComplexF64}:
1.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 0.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
γ = 1
Sx = sum([mb(sx, i, latt) for i in 1:latt.N])
Sy = sum([mb(sy, i, latt) for i in 1:latt.N])
Sz = sum([mb(sz, i, latt) for i in 1:latt.N])
SFxx = sum([mb(sx, i, latt) * mb(sx, j, latt) for i in 1:latt.N for j in 1:latt.N])
H, c_ops = TFIM(Jx, Jy, Jz, hx, γ, latt; bc = pbc, order = 1)
e_ops = (Sx, Sy, Sz, SFxx)
tl = LinRange(0, 10, 100);100-element LinRange{Float64, Int64}:
0.0, 0.10101, 0.20202, 0.30303, 0.40404, …, 9.69697, 9.79798, 9.89899, 10.0Full evolution
@time mesol = mesolve(H, ρ, tl, c_ops; e_ops = [e_ops...]);
A = Matrix(mesol.states[end].data)
λ = eigvals(Hermitian(A))
Strue = -sum(λ .* log2.(λ)) / latt.N;0.29324728062546734Low 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)
λ = eigvals(Hermitian(σ))
λ = λ[λ.>1e-10]
return -sum(λ .* log2.(λ))
end;f_entropy (generic function with 1 method)Define the options for the low-rank evolution
opt =
LRMesolveOptions(err_max = 1e-3, p0 = 0.0, atol_inv = 1e-6, adj_condition = "variational", Δt = 0.0);
@time lrsol = lr_mesolve(H, z, B, tl, c_ops; e_ops = e_ops, f_ops = (f_purity, f_entropy, f_trace), opt = opt);LRTimeEvolutionSol{Vector{Float64}, Vector{Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}}, Matrix{ComplexF64}, Vector{Int64}}([0.0, 0.10101010101010102, 0.20202020202020204, 0.30303030303030304, 0.4040404040404041, 0.5050505050505051, 0.6060606060606061, 0.7070707070707071, 0.8080808080808082, 0.9090909090909092 … 9.09090909090909, 9.191919191919192, 9.292929292929292, 9.393939393939394, 9.494949494949495, 9.595959595959595, 9.696969696969697, 9.797979797979798, 9.8989898989899, 10.0], Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.975870158217177 - 0.0062733283402879175im 4.519295912638144e-6 - 1.2961438268448496e-6im … 7.34314745047746e-5 - 5.9312268904361754e-5im 3.061697633799469e-5 - 3.6096001741519286e-6im; 7.440221117122334e-5 + 8.91425323513208e-5im 0.002205479036235476 - 0.00997657005222155im … -0.029412714507766516 - 0.01956655648028395im -0.020774395504225152 - 0.0006418548176660775im; … ; -6.17718235284902e-6 - 2.926106253183959e-5im -0.012542139986330317 - 0.005517545235872339im … 0.01848130781325862 - 0.007646337552937583im 0.003367023624732646 - 0.01026696828884103im; 0.002457949115040511 - 0.005900048882616366im -1.7996252057768343e-5 + 9.633825072978356e-6im … -6.466384722065548e-5 - 2.470745261961446e-5im -1.674559278467308e-5 - 7.10287099739902e-6im]], Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.6675323028396343 - 2.641376125598628e-27im -1.2217747336520244e-6 + 5.886939567354092e-7im … 1.3947751258207069e-5 - 1.0170577640932473e-5im 2.261022004402645e-6 - 1.5858908674238402e-6im; -1.2217747336520244e-6 - 5.886939567354092e-7im 0.04057558309908829 - 1.6055459192650106e-28im … -6.71835641136147e-8 - 2.4150263662724667e-5im -1.140093117912524e-5 + 1.8113420283170253e-5im; … ; 1.3947751258207069e-5 + 1.0170577640932473e-5im -6.71835641136147e-8 + 2.4150263662724667e-5im … 0.0003176522199324979 - 1.2569264234912577e-30im 3.815730845202987e-6 + 3.1966037982961808e-6im; 2.261022004402645e-6 + 1.5858908674238402e-6im -1.140093117912524e-5 - 1.8113420283170253e-5im … 3.815730845202987e-6 - 3.1966037982961808e-6im 0.00031494411896017727 - 1.2462106675292326e-30im]], ComplexF64[0.0 + 0.0im 0.0 + 0.0im … 3.3786313230450642e-6 - 8.673617379884035e-19im 3.3140913457884932e-6 - 3.604972223514302e-18im; 0.0 + 0.0im 0.0 + 0.0im … -5.325574989928426e-6 - 5.421010862427522e-19im -5.195696904924929e-6 - 1.81603863891322e-18im; -6.0 + 0.0im -5.987982183309875 - 1.8243808117344763e-19im … -4.564263714975137 - 7.806255641895632e-18im -4.564345654448981 + 5.811323644522304e-17im; 6.0 + 0.0im 6.013514592692536 - 3.8194631541223296e-17im … 8.148098529295455 - 2.2898349882893854e-16im 8.147992896631594 + 6.245004513516506e-16im], ComplexF64[1.0 + 0.0im 0.9995913405802062 + 1.1922112809822968e-38im … 0.5162581912415639 - 3.815963511559993e-18im 0.5162814480901279 - 1.6125745767007885e-17im; -0.0 + 0.0im 0.0031973769163574704 + 0.0im … 1.733915564776959 + 0.0im 1.7338318986613197 + 0.0im; 1.0 + 0.0im 1.0 + 5.877471754111438e-39im … 1.0 + 1.672385119068545e-27im 1.0 - 4.020035173006276e-27im], [7, 7, 7, 7, 7, 7, 7, 8, 15, 17 … 36, 36, 36, 36, 36, 36, 36, 36, 36, 36])Plot the results
m_me = real(mesol.expect[3, :]) / Nx / Ny
m_lr = real(lrsol.expvals[3, :]) / Nx / Ny
fig = Figure(size = (800, 400), 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(lrsol.funvals[1, :]), label = L"$1-P$", linewidth = 2)
lines!(
ax2,
tl,
1 .- real(lrsol.funvals[3, :]),
label = L"$1-\mathrm{Tr}(\rho)$",
linewidth = 2,
linestyle = :dash,
color = :orange,
)
lines!(ax2, tl, real(lrsol.funvals[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