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.2932472806254766Low 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}}([10.0], Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.9758685522789969 - 0.006273695106632848im 2.492860914637139e-6 - 7.156869837890943e-7im … 3.775797443628507e-5 + 2.9466380992070495e-5im 3.268728560347496e-5 - 1.709317159884015e-5im; 0.0003981574745855075 + 0.00013310747075959934im 0.00257174247123672 - 0.010061680639511747im … 0.054923538663834016 - 0.02137275699750448im 0.005841112052344932 - 0.041740957520038806im; … ; -0.0002493314214458491 - 7.203926799350146e-5im -0.01241673945732015 - 0.00568813699603282im … 0.007031452325110198 + 0.025002371274166428im 0.019476740386133236 + 0.003706407258705394im; 0.0024681064343023674 - 0.005906389264660715im -6.161371270946958e-6 - 2.950353589551846e-6im … 8.349409961913042e-5 - 0.00016726765453363302im -4.002673859618563e-5 - 1.2590844796523399e-5im]], Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.6675357410773601 + 1.6008422608799466e-27im -8.240146084828846e-7 - 3.890252226983361e-8im … 1.36033770612411e-5 - 1.1377093821294141e-5im 2.546882477495092e-6 - 1.907521227709345e-6im; -8.240146084828846e-7 + 3.890252226983361e-8im 0.04054650574801683 + 9.723608181290597e-29im … -1.0598317074863365e-5 - 3.783592954232962e-5im -7.341015160578552e-6 - 4.427779211125318e-6im; … ; 1.36033770612411e-5 + 1.1377093821294141e-5im -1.0598317074863365e-5 + 3.783592954232962e-5im … 0.000314711721328375 + 7.5472187104712665e-31im -2.7653378343184173e-6 + 5.894423047904932e-6im; 2.546882477495092e-6 + 1.907521227709345e-6im -7.341015160578552e-6 + 4.427779211125318e-6im … -2.7653378343184173e-6 - 5.894423047904932e-6im 0.00031572068582301694 + 7.57141506286647e-31im]], ComplexF64[0.0 + 0.0im 0.0 + 0.0im … 1.3478110667420468e-6 + 1.0842021724855044e-19im 1.3273272125564981e-6 - 2.7647155398380363e-18im; 0.0 + 0.0im 0.0 + 0.0im … -2.3337804213768073e-6 - 1.3010426069826053e-18im -2.271404987073803e-6 - 4.0657581468206416e-19im; -6.0 + 0.0im -5.987982183309876 - 1.9413917363256665e-19im … -4.564265490409348 + 7.112366251504909e-17im -4.564347160978902 + 2.1250362580715887e-17im; 6.0 + 0.0im 6.013514592692537 + 6.053762417912302e-18im … 8.148099395377217 + 5.828670879282072e-16im 8.14799372925165 - 4.163336342344337e-16im], ComplexF64[1.0 + 0.0im 0.9995913405802064 + 3.545424123549448e-38im … 0.5162593322096489 - 5.2348535528430034e-18im 0.5162824748057558 - 3.33048019326886e-18im; -0.0 + 0.0im 0.0031973769163571412 + 0.0im … 1.7339106629513463 + 0.0im 1.733827512551668 + 0.0im; 1.0 + 0.0im 1.0000000000000002 + 2.350988701644575e-38im … 1.0000000000000002 + 3.988480736797436e-27im 1.0 + 3.5088533064230605e-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