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(alg = Tsit5(), 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.9758680491596291 - 0.006273233699143627im 2.6082250601415474e-6 - 8.042205231015986e-7im … 2.2545949418629325e-5 + 6.825395504179448e-6im 1.1688617296036082e-5 - 1.4503086192035739e-5im; 0.00037110212527516773 + 5.468475025042419e-5im 0.0026248465372509682 - 0.0100687744191106im … -0.043981864757323935 + 0.05873275332355478im 0.0197805009767266 + 0.005746599995520802im; … ; -0.0001864036386563782 - 0.00010751309345776612im -0.012445363273406515 - 0.0056479415824247136im … -0.011742010622634939 - 0.029612745632784054im -0.0014978945579874872 + 0.010529711609283829im; 0.0024683161830855345 - 0.005906451660800724im -5.936728049495166e-6 - 2.9425356160346778e-6im … 9.161076019912218e-5 - 0.00020853743565589714im 3.015882507657663e-6 + 1.2599287249351672e-5im]], Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.6675354158987169 + 2.148497777208022e-27im -3.449187286945449e-7 + 4.254277809164585e-7im … 1.5867891284263173e-5 - 1.4281803225370318e-5im 2.0392845209154294e-6 - 1.4486694094205966e-6im; -3.449187286945449e-7 - 4.254277809164585e-7im 0.04054554136975638 + 1.3049795326789782e-28im … -9.38590727749871e-6 + 4.983618692124325e-6im 1.0778217411080972e-5 - 1.020650110529565e-5im; … ; 1.5867891284263173e-5 + 1.4281803225370318e-5im -9.38590727749871e-6 - 4.983618692124325e-6im … 0.0003102744731612362 + 9.986346792009773e-31im -4.7511991643280847e-7 + 1.2413213174753718e-6im; 2.0392845209154294e-6 + 1.4486694094205966e-6im 1.0778217411080972e-5 + 1.020650110529565e-5im … -4.7511991643280847e-7 - 1.2413213174753718e-6im 0.0003203952732559365 + 1.0312090055799833e-30im]], ComplexF64[0.0 + 0.0im 0.0 + 0.0im … -4.608715993708137e-6 - 2.6020852139652106e-18im -4.491893472506668e-6 - 1.0842021724855044e-19im; 0.0 + 0.0im 0.0 + 0.0im … 5.60411592229212e-6 - 8.944667923005412e-19im 5.5359662690587514e-6 - 1.3552527156068805e-18im; -6.0 + 0.0im -5.987982183309876 - 1.9413917363256665e-19im … -4.5642644153688074 + 7.28583859910259e-17im -4.564346140456851 - 1.3444106938820255e-17im; 6.0 + 0.0im 6.013514592692537 + 6.053762417912302e-18im … 8.148096744356959 + 6.38378239159465e-16im 8.147990994054492 - 5.551115123125783e-16im], ComplexF64[1.0 + 0.0im 0.9995913405802064 + 3.545424123549448e-38im … 0.5162587835918441 - 8.665285209690444e-18im 0.5162819619102399 - 2.88565590616302e-18im; -0.0 + 0.0im 0.0031973769163571412 + 0.0im … 1.7339139277465936 + 0.0im 1.73383069377884 + 0.0im; 1.0 + 0.0im 1.0000000000000002 + 2.350988701644575e-38im … 1.0 + 1.5714109232002555e-27im 1.0 - 1.8869552852886602e-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