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.29324728062546573Low 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.9758698989252109 - 0.006271649796513063im -1.2826871522200105e-6 + 5.014789266833318e-7im … 3.499774840193694e-5 - 2.7866848914528426e-5im 7.297286139252118e-6 - 4.0726757493804194e-6im; 0.00013433505858931414 + 7.229565812173492e-5im 0.0019421023536524333 - 0.010229021806554705im … 0.020339230022600806 - 0.021251931738810438im -0.015765883630352848 + 0.0016034781375374515im; … ; -3.617588342560502e-5 - 3.440597576763145e-6im -0.012444032592489822 - 0.005644742476205572im … 0.0073651636498951656 + 0.012562141979339454im 0.0018218465863280272 - 0.010334722682188935im; 0.002458256333708311 - 0.005903488091690724im 7.167860001713305e-6 + 6.980269778432505e-6im … 1.711223415598168e-5 - 8.887746304136125e-5im -1.9645469690389544e-5 - 1.5339185914557983e-5im]], Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.6675326204079556 - 2.1505951413230616e-27im -6.904907776002374e-7 + 9.788047232776618e-7im … 1.6340577434937217e-5 - 1.3563304811454035e-5im 2.0751438904462665e-6 - 1.273654853260319e-6im; -6.904907776002374e-7 - 9.788047232776618e-7im 0.04057680269368733 - 1.3072660729319542e-28im … 5.9200096720292375e-6 - 4.6522179382207906e-5im 1.88274345211871e-6 + 1.613261389613719e-5im; … ; 1.6340577434937217e-5 + 1.3563304811454035e-5im 5.9200096720292375e-6 + 4.6522179382207906e-5im … 0.0003138680389573182 - 1.0111911522551284e-30im 5.765386223043992e-6 - 4.521982973014431e-7im; 2.0751438904462665e-6 + 1.273654853260319e-6im 1.88274345211871e-6 - 1.613261389613719e-5im … 5.765386223043992e-6 + 4.521982973014431e-7im 0.00032085473399777036 - 1.0337002431198112e-30im]], ComplexF64[0.0 + 0.0im 0.0 + 0.0im … -4.244564344886993e-6 + 1.0299920638612292e-18im -4.162866407233736e-6 - 7.995991022080595e-19im; 0.0 + 0.0im 0.0 + 0.0im … 6.415569687269213e-6 + 2.981555974335137e-19im 6.2784328292197e-6 - 4.336808689942018e-19im; -6.0 + 0.0im -5.987982183309876 + 1.5277690545447886e-20im … -4.564264299147102 + 7.806255641895632e-17im -4.56434599605206 + 2.2551405187698492e-17im; 6.0 + 0.0im 6.013514592692536 + 5.709650890972888e-17im … 8.148091613908333 + 8.743006318923108e-16im 8.14798607815379 - 2.636779683484747e-16im], ComplexF64[1.0 + 0.0im 0.9995913405802064 + 1.188953990727943e-37im … 0.5162587083042479 + 3.208628836298161e-19im 0.5162818763239461 + 6.449038555012288e-18im; -0.0 + 0.0im 0.0031973769163571525 + 0.0im … 1.7339159387457204 + 0.0im 1.7338327048837956 + 0.0im; 1.0 + 0.0im 1.0000000000000002 + 1.4179400606793843e-37im … 1.0000000000000002 - 3.1775317262302356e-27im 1.0000000000000002 + 2.4928004604983973e-28im], [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 = (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(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