Low rank master equation

Author

Luca Gravina

Last Update

2025-01-13

In this tutorial, we will show how to solve the master equation using the low-rank method. For a detailed explanation of the method, we recommend to read the Ref. (Gravina and Savona 2024).

As a test, we will consider the dissipative Ising model with a transverse field. The Hamiltonian is given by

\[ \hat{H} = \frac{J_x}{2} \sum_{\langle i,j\rangle} \sigma_i^x \sigma_j^x + \frac{J_y}{2} \sum_{\langle i,j\rangle} \sigma_i^y \sigma_j^y + \frac{J_z}{2} \sum_{\langle i,j\rangle} \sigma_i^z \sigma_j^z - \sum_i h_i \sigma_i^z + h_x \sum_i \sigma_i^x + h_y \sum_i \sigma_i^y + h_z \sum_i \sigma_i^z \] K where the sums are over nearest neighbors, and the collapse operators are given by

\[ c_i = \sqrt{\gamma} \sigma_i^{-} \]

We start by importing the packages

using QuantumToolbox
using CairoMakie

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
N_modes = latt.N
N = N_cut^N_modes
M = latt.N + 1
7

Define lr states. Take as initial state all spins up. All other N states are taken as those with minimum Hamming distance to the initial state.

ϕ = Vector{QuantumObject{Vector{ComplexF64},KetQuantumObject,M-1}}(undef, M)
ϕ[1] = kron(fill(basis(2, 1), N_modes)...)

i = 1
for j in 1:N_modes
    global i += 1
    i <= M && (ϕ[i] = MultiSiteOperator(latt, j=>sigmap()) * ϕ[1])
end
for k in 1:N_modes-1
    for l in k+1:N_modes
        global i += 1
        i <= M && (ϕ[i] = MultiSiteOperator(latt, k=>sigmap(), l=>sigmap()) * ϕ[1])
    end
end
for i in i+1:M
    ϕ[i] = QuantumObject(rand(ComplexF64, size(ϕ[1])[1]), dims = ϕ[1].dims)
    normalize!(ϕ[i])
end

Define the initial state

z = hcat(get_data.(ϕ)...)
B = Matrix(Diagonal([1 + 0im; zeros(M - 1)]))
S = z' * z # Overlap matrix
B = B / tr(S * B) # Normalize B

ρ = QuantumObject(z * B * z', dims = ntuple(i->N_cut, Val(N_modes))); # Full density matrix

Define the Hamiltonian and collapse operators

# Define Hamiltonian and collapse operators
Jx = 0.9
Jy = 1.04
Jz = 1.0
hx = 0.0
hy = 0.0
hz = 0.0
γ = 1

Sx = mapreduce(i->MultiSiteOperator(latt, i=>sigmax()), +, 1:latt.N)
Sy = mapreduce(i->MultiSiteOperator(latt, i=>sigmay()), +, 1:latt.N)
Sz = mapreduce(i->MultiSiteOperator(latt, i=>sigmaz()), +, 1:latt.N)

H, c_ops = DissipativeIsing(Jx, Jy, Jz, hx, hy, hz, γ, latt; boundary_condition = Val(:periodic_bc), order = 1)
e_ops = (Sx, Sy, Sz)

tl = range(0, 10, 100)
0.0:0.10101010101010101:10.0

Full evolution

sol_me = mesolve(H, ρ, tl, c_ops; e_ops = [e_ops...]);
Strue = entropy_vn(sol_me.states[end], base=2) / latt.N
Progress: [==============================] 100.0% --- Elapsed Time: 0h 00m 00s (ETA: 0h 00m 00s)
0.29324728062546046

Low 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)
    return entropy_vn(Qobj(Hermitian(σ), type=Operator), base=2)
end
f_entropy (generic function with 1 method)

Define the options for the low-rank evolution

opt = (err_max = 1e-3, p0 = 0.0, atol_inv = 1e-6, adj_condition = "variational", Δt = 0.0);

sol_lr = lr_mesolve(H, z, B, tl, c_ops; e_ops = e_ops, f_ops = (f_purity, f_entropy, f_trace), opt = opt);
Progress: 2%Progress: 3%Progress: 4%Progress: 5%Progress: 6%Progress: 7%Progress: 8%Progress: 9%Progress: 10%Progress: 11%Progress: 12%Progress: 13%Progress: 14%Progress: 15%Progress: 16%Progress: 17%Progress: 18%Progress: 19%Progress: 20%Progress: 21%Progress: 22%Progress: 23%Progress: 24%Progress: 25%Progress: 26%Progress: 27%Progress: 28%Progress: 29%Progress: 30%Progress: 31%Progress: 32%Progress: 33%Progress: 34%Progress: 35%Progress: 36%Progress: 37%Progress: 38%Progress: 39%Progress: 40%Progress: 41%Progress: 42%Progress: 43%Progress: 44%Progress: 45%Progress: 46%Progress: 47%Progress: 48%Progress: 49%Progress: 50%Progress: 51%Progress: 52%Progress: 53%Progress: 54%Progress: 55%Progress: 56%Progress: 57%Progress: 58%Progress: 59%Progress: 60%Progress: 61%Progress: 62%Progress: 63%Progress: 64%Progress: 65%Progress: 66%Progress: 67%Progress: 68%Progress: 69%Progress: 70%Progress: 71%Progress: 72%Progress: 73%Progress: 74%Progress: 75%Progress: 76%Progress: 77%Progress: 78%Progress: 79%Progress: 80%Progress: 81%Progress: 82%Progress: 83%Progress: 84%Progress: 85%Progress: 86%Progress: 87%Progress: 88%Progress: 89%Progress: 90%Progress: 91%Progress: 92%Progress: 93%Progress: 94%Progress: 95%Progress: 96%Progress: 97%Progress: 98%Progress: 99%Progress: 100%

Plot the results

m_me = real(sol_me.expect[3, :]) / Nx / Ny
m_lr = real(sol_lr.expect[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(sol_lr.fexpect[1, :]), label = L"$1-P$", linewidth = 2)
lines!(
    ax2,
    tl,
    1 .- real(sol_lr.fexpect[3, :]),
    label = L"$1-\mathrm{Tr}(\rho)$",
    linewidth = 2,
    linestyle = :dash,
    color = :orange,
)
lines!(ax2, tl, real(sol_lr.fexpect[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

Version Information

QuantumToolbox.versioninfo()

 QuantumToolbox.jl: Quantum Toolbox in Julia
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
Copyright © QuTiP team 2022 and later.
Current admin team:
    Alberto Mercurio and Yi-Te Huang

Package information:
====================================
Julia              Ver. 1.11.2
QuantumToolbox     Ver. 0.24.0
SciMLOperators     Ver. 0.3.12
LinearSolve        Ver. 2.38.0
OrdinaryDiffEqCore Ver. 1.14.1

System information:
====================================
OS       : Linux (x86_64-linux-gnu)
CPU      : 4 × AMD EPYC 7763 64-Core Processor
Memory   : 15.615 GB
WORD_SIZE: 64
LIBM     : libopenlibm
LLVM     : libLLVM-16.0.6 (ORCJIT, znver3)
BLAS     : libopenblas64_.so (ilp64)
Threads  : 4 (on 4 virtual cores)

References

Gravina, Luca, and Vincenzo Savona. 2024. “Adaptive variational low-rank dynamics for open quantum systems.” Phys. Rev. Res. 6 (April): 023072. https://doi.org/10.1103/PhysRevResearch.6.023072.