Quick Start

Content

• Import HierarchicalEOM.jl
• System and Bath
• HEOM Liouvillian superoperator
• Time Evolution
• Stationary State
• Reduced Density Operator
• Expectation Value
• Multiple Baths

Import HierarchicalEOM.jl

Here are the functions in HierarchicalEOM.jl that we will use in this tutorial (Quick Start):

import HierarchicalEOM
import HierarchicalEOM: Boson_DrudeLorentz_Pade, M_Boson, evolution, SteadyState, getRho, BosonBath, Expect

Note that you can also type using HierarchicalEOM to import everything you need in HierarchicalEOM.jl. To check the versions of dependencies of HierarchicalEOM.jl , run the following function

HierarchicalEOM.versioninfo()

                                   __
                                  /  \
 __     __                     __ \__/ __
|  |   |  |                   /  \    /  \
|  |   |  | ______   ______   \__/_  _\__/
|  |___|  |/  __  \ /  __  \ / '   \/     \
|   ___   |  |__)  |  /  \  |    _     _   |
|  |   |  |   ____/| (    ) |   / \   / \  |
|  |   |  |  |____ |  \__/  |  |   | |   | |
|__|   |__|\______) \______/|__|   |_|   |_|

Julia framework for Hierarchical Equations of Motion
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
Copyright © NCKU-QFORT 2022 and later.
Lead  developer : Yi-Te Huang
Other developers:
    Simon Cross, Neill Lambert, Po-Chen Kuo and Shen-Liang Yang

Package information:
====================================
HierarchicalEOM Ver. 1.6.0
LinearSolve     Ver. 2.22.1
OrdinaryDiffEq  Ver. 6.80.1
FastExpm        Ver. 1.1.0
JLD2            Ver. 0.4.48

System information:
====================================
Julia Version: 1.10.4
OS       : Linux (x86_64-linux-gnu)
CPU      : 4 × AMD EPYC 7763 64-Core Processor
Memory   : 15.606 GB
WORD_SIZE: 64
LIBM     : libopenlibm
LLVM     : libLLVM-15.0.7 (ORCJIT, znver3)
BLAS     : libopenblas64_.so (ilp64)

System and Bath

Let us consider a simple two-level system coupled to a Drude-Lorentz bosonic bath. The system Hamiltonian, $H_{sys}$, and the bath spectral density, $J_D$, are

\[H_{sys}=\frac{\epsilon \sigma_z}{2} + \frac{\Delta \sigma_x}{2} ~~\text{and}\]

\[J_{D}(\omega)=\frac{2\lambda W\omega}{W^2+\omega^2},\]

System Hamiltonian and initial state

You can construct system hamiltonian, initial state, and coupling operators by standard julia bulit-in types: Vector, SparseVector, Matrix, SparseMatrix.

Moreover, it is also convenient to use QuantumToolbox or QuantumOptics framework. They both provide many useful functions to create arbitrary quantum states and operators which can be combined in all the expected ways.

We demonstrate this tutorial by QuantumToolbox:

import QuantumToolbox: sigmaz, sigmax, basis, ket2dm

# The system Hamiltonian
ϵ = 0.5 # energy of 2-level system
Δ = 1.0 # tunneling term

Hsys = 0.5 * ϵ * sigmaz() + 0.5 * Δ * sigmax()

# System initial state
ρ0 = ket2dm(basis(2, 0));

Bath Properties

Now, we demonstrate how to describe the bath using the built-in implementation of $J_D(\omega)$ under Pade expansion by calling Boson_DrudeLorentz_Pade

λ = 0.1  # coupling strength
W = 0.5  # band-width (cut-off frequency)
kT = 0.5  # the product of the Boltzmann constant k and the absolute temperature T

Q = sigmaz() # system-bath coupling operator

N = 2 # Number of expansion terms to retain:

# Padé expansion:
bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N)

BosonBath object with (system) dim = 2 and 3 exponential-expansion terms

For other different expansions of the different spectral density correlation functions, please refer to Bosonic Bath and Fermionic Bath.

HEOM Liouvillian superoperator

For bosonic bath, we can construct the HEOM Liouvillian superoperator matrix by calling M_Boson

tier = 5 # maximum tier of hierarchy
L = M_Boson(Hsys, tier, bath)

Boson type HEOMLS matrix acting on even-parity ADOs
system dim = 2
number of ADOs N = 56
data =
224×224 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 1160 stored entries:
⎡⠿⣧⣄⠀⠘⠤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠃⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠙⠻⣦⡄⠀⠑⠠⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠒⡄⠀⠉⢱⣶⡄⠀⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠈⠑⡀⠀⠉⢻⣶⠀⠈⠒⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠈⠢⡀⠀⢿⣷⣀⠱⣀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠘⠀⢄⡘⠻⣦⣀⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠠⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠀⢜⠛⣤⡤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⢤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠋⠛⣤⡄⠀⠠⠄⠀⠀⠀⠀⠀⠠⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠑⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢻⣶⡄⠈⠢⡀⠀⠀⠀⠀⠈⠒⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠉⠰⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠆⡀⠉⢱⣶⡀⠰⣀⠀⠀⠀⠀⠀⠰⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠳⠄⡀⠀⠀⠀⠀⠀⠀⠀⠈⠢⢀⡈⠿⢇⣀⢁⠀⠀⠀⠀⠀⠀⢁⡀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⢀⡀⠀⠀⠀⠀⠀⠀⠀⠘⠄⢘⠿⢇⣘⡀⠀⠀⠀⠀⠀⠘⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠂⢄⠀⡀⠀⠀⠀⠀⠀⠀⠒⠸⠛⣤⣀⠀⢀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢢⡀⠀⠀⠀⠀⠀⠀⠀⠘⢻⣶⡌⠒⡄⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠐⢦⠀⠀⠀⠀⠀⠐⢢⠉⢱⣶⡔⢢⠀⠀⠀⠢⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠰⣄⠀⠀⠀⠀⠉⠰⣉⢻⣶⠆⠀⠀⠀⠰⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠂⠰⣄⠀⠀⠀⠈⠁⠻⣦⣀⢄⡀⢄⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⠰⡄⠀⠀⠀⢜⠻⣦⡤⠀⠠⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠘⠂⠠⣌⠀⠋⠛⣤⡤⠤⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⠆⠀⣏⠱⣦⎦

To learn more about the HEOM Liouvillian superoperator matrix (including other types: M_Fermion, M_Boson_Fermion), please refer to HEOMLS Matrices.

Time Evolution

Next, we can calculate the time evolution for the entire auxiliary density operators (ADOs) by calling evolution

tlist = 0:0.2:50
ados_list = evolution(L, ρ0, tlist);

Solving time evolution for ADOs by Ordinary Differential Equations method...
[DONE]

To learn more about evolution, please refer to Time Evolution.

Stationary State

We can also solve the stationary state of the auxiliary density operators (ADOs) by calling SteadyState.

ados_steady = SteadyState(L)

56 Auxiliary Density Operators with even-parity and (system) dim = 2

To learn more about SteadyState, please refer to Stationary State.

Reduced Density Operator

To obtain the reduced density operator, one can either access the first element of auxiliary density operator (ADOs) or call getRho:

# reduce density operator in third time step of the evolution
ρ = ados_list[3][1]
ρ = getRho(ados_list[3])

# reduce density operator in stationary state
ρ = ados_steady[1]
ρ = getRho(ados_steady);

One of the great features of HierarchicalEOM.jl is that we allow users to not only considering the density operator of the reduced state but also easily take high-order terms into account without struggling in finding the indices (see Auxiliary Density Operators and Hierarchy Dictionary for more details).

Expectation Value

We can now compare the results obtained from evolution and SteadyState:

# Define the operators that measure the populations of the two
# system states:
P00 = ket2dm(basis(2, 0))
P11 = ket2dm(basis(2, 1))

# Define the operator that measures the 0, 1 element of density matrix
# (corresponding to coherence):
P01 = basis(2, 0) * basis(2, 1)'

# for steady state
p00_s = Expect(P00, ados_steady)
p01_s = Expect(P01, ados_steady)

# for time evolution
p00_e = Expect(P00, ados_list)
p01_e = Expect(P01, ados_list);

Plot the results

using Plots, LaTeXStrings

plot(tlist, p00_e, label = L"\textrm{P}_{00}", linecolor = :blue, linestyle = :solid, linewidth = 3, grid = false)
plot!(tlist, p01_e, label = L"\textrm{P}_{01}", linecolor = :red, linestyle = :solid, linewidth = 3)
plot!(
    tlist,
    ones(length(tlist)) .* p00_s,
    label = L"\textrm{P}_{00} \textrm{(Steady State)}",
    linecolor = :blue,
    linestyle = :dash,
    linewidth = 3,
)
plot!(
    tlist,
    ones(length(tlist)) .* p01_s,
    label = L"\textrm{P}_{01} \textrm{(Steady State)}",
    linecolor = :red,
    linestyle = :dash,
    linewidth = 3,
)

xlabel!("time")
ylabel!("Population")

Multiple Baths

HierarchicalEOM.jl also supports for system to interact with multiple baths.

All you need to do is to provide a list of baths instead of a single bath

Note that, for the following, we use the built-in linear algebra in Julia (instead of QuantumToolbox.jl) to construct the operators

# The system Hamiltonian
Hsys = [
    0.25 1.50 2.50
    1.50 0.75 3.50
    2.50 3.50 1.25
]

# System initial state
ρ0 = [
    1 0 0
    0 0 0
    0 0 0
];

# Construct one bath for each system state:
# note that `BosonBath[]` make the list created in type: Vector{BosonBath}
baths = BosonBath[]
for i in 1:3
    # system-bath coupling operator: |i><i|
    Q = zeros(3, 3)
    Q[i, i] = 1

    push!(baths, Boson_DrudeLorentz_Pade(Q, λ, W, kT, N))
end

L = M_Boson(Hsys, tier, baths)

tlist = 0:0.025:5
ados_list = evolution(L, ρ0, tlist)

# Projector for each system state:
P00 = [1 0 0; 0 0 0; 0 0 0]
P11 = [0 0 0; 0 1 0; 0 0 0]
P22 = [0 0 0; 0 0 0; 0 0 1]

# calculate population for each system state:
p0 = Expect(P00, ados_list)
p1 = Expect(P11, ados_list)
p2 = Expect(P22, ados_list)

plot(tlist, p0, linewidth = 3, linecolor = "blue", label = L"P_0", grid = false)
plot!(tlist, p1, linewidth = 3, linecolor = "orange", label = L"P_1")
plot!(tlist, p2, linewidth = 3, linecolor = :green, label = L"P_2")
xlabel!("time")
ylabel!("Population")

Note that this example can also be found in qutip documentation.

This page was generated using Literate.jl.