Low Rank Master Equation
We start by importing the packages
using Plots
using LaTeXStrings
using QuantumToolbox;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=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. # 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.
hx = 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.29324728071737893Low 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)
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)
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.,
atol_inv = 1e-6,
adj_condition="variational",
Δt = 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.9765077144286615 - 0.006219079334131151im 3.2738471569113803e-7 - 1.2433471897187135e-7im … -6.041837417576852e-5 - 5.104087727021705e-5im -7.765622853291598e-5 - 9.71957521734838e-6im; -3.2129191790072773e-6 + 8.087303783668407e-8im 0.0022678070883443095 - 0.010329981916013795im … -0.05404631687980407 - 0.005849192984818809im -0.034612086469765393 + 0.00614939201344269im; … ; -5.118646376613913e-6 - 1.2877734070820306e-5im -0.012499998833705214 - 0.005168867520479114im … 0.009743669962713142 - 0.030474128946095323im -0.011403424060781445 - 0.013910409417949476im; 0.002555201012396967 - 0.006054926255348682im -5.153994446765919e-6 + 8.471326158858638e-6im … 0.00020920502670602656 - 0.00011539567778750413im 4.6442622956752305e-5 - 9.348971972053967e-6im]], Base.ReshapedArray{ComplexF64, 2, SubArray{ComplexF64, 1, Vector{ComplexF64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}[[0.6667701299939592 - 2.4264498604930655e-25im -2.905934722054598e-7 - 1.6201572430339277e-7im … 1.573099961485082e-5 - 1.272074949756447e-5im 1.7111520631743565e-6 - 1.8000758891826222e-6im; -2.905934722054598e-7 + 1.6201572430339277e-7im 0.040572988483578896 - 1.4764956889514336e-26im … -1.0317461348245335e-6 + 1.8644462488532083e-5im 3.30123653071663e-7 - 9.08827931179431e-6im; … ; 1.573099961485082e-5 + 1.272074949756447e-5im -1.0317461348245335e-6 - 1.8644462488532083e-5im … 0.0003126801411520382 - 1.1378774344378459e-28im -1.72807072659889e-5 + 4.262386051411162e-6im; 1.7111520631743565e-6 + 1.8000758891826222e-6im 3.30123653071663e-7 + 9.08827931179431e-6im … -1.72807072659889e-5 - 4.262386051411162e-6im 0.00029207375150039913 - 1.0628885153992282e-28im]], ComplexF64[0.0 + 0.0im 0.0 + 0.0im … -1.0560591333953762e-6 + 5.421010862427522e-19im -1.0652403984604333e-6 - 1.8431436932253575e-18im; 0.0 + 0.0im 0.0 + 0.0im … 3.823589705409525e-6 + 8.673617379884035e-19im 3.637521248456418e-6 + 3.577867169202165e-18im; -6.0 + 0.0im -5.987982183320124 - 5.811416879613635e-20im … -4.5642437343926625 - 6.245004513516506e-17im -4.56432675976526 + 8.673617379884035e-19im; 6.0 + 0.0im 6.0135145926944364 - 7.497557741935917e-17im … 8.148126910955343 - 5.134781488891349e-16im 8.148018212334872 - 1.0547118733938987e-15im], ComplexF64[1.0 + 0.0im 0.9995913405809307 + 1.456507381741672e-38im … 0.516253382141987 - 4.1298876929083915e-18im 0.5162767010000036 - 3.658063043294913e-18im; -0.0 + 0.0im 0.0031973769106007264 + 0.0im … 1.7339345602122123 + 0.0im 1.7338510328083667 + 0.0im; 1.0 + 0.0im 1.0000000000000007 - 8.816207631167156e-39im … 1.0 - 3.0776934900654643e-25im 1.0 - 3.1503949110906327e-25im], [7, 7, 7, 7, 7, 7, 7, 8, 15, 17 … 36, 36, 36, 36, 36, 36, 36, 36, 36, 36])Plot the results
fig = plot(layout=(1,2), size=(800,400), legend=:topleft, xlabel=L"\gamma t")
m_me = real(mesol.expect[3,:])/Nx/Ny
m_lr = real(lrsol.expvals[3,:])/Nx/Ny
plot!(fig[1], tl, m_lr, label=raw"LR $[M=M(t)]$", lw=2)
plot!(fig[1], tl, m_me, ls=:dash, label="Fock", lw=2)
ylabel!(fig[1], L"M_{z}")
plot!(fig[2], tl, 1 .-real(lrsol.funvals[1,:]), label=L"$1-P$", lw=2)
plot!(fig[2], tl, 1 .-real(lrsol.funvals[3,:]), c=:orange, label=L"$1-\rm{Tr}(\rho)$", lw=2, ls=:dash)
plot!(fig[2], tl, real(lrsol.funvals[2,:])/Nx/Ny, c=:blue, label=L"S", lw=2)
hline!(fig[2], [Strue], c=:blue, ls=:dash, lw=2, label=L"S^{\rm \,true}_{\rm ss}")
ylabel!(fig[2], "value")
xlabel!(fig[2], L"\gamma t")