Time Evolution

Introduction

HierarchicalEOM.jl implements various methods and solvers to simulate the open quantum system dynamics. The HEOM Liouvillian superoperator (HEOMLS) matrix $\hat{\mathcal{M}}$ characterizes the dynamics of the reduce state and in the full extended space of all auxiliary density operators (ADOs) $\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t)$, namely

\[\begin{equation} \partial_{t}\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t)=\hat{\mathcal{M}}\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t) \end{equation}\]

Output of `evolution`

To solve the dynamics of the reduced state and also all the ADOs, you only need to call evolution. Different methods are implemented with different input parameters of the function which makes it easy to switch between different methods. The output of the function evolution for each methods will always be in the type Vector{ADOs}, which contains a list of ADOs corresponds to the given time steps.

Expectation Values

Given an observable $A$ and the ADOs $\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t)$, one can calculate the expectation value by

\[\langle A(t) \rangle = \textrm{Tr}\left[A \rho^{(0,0,p)}_{ \vert }(t)\right],\]

where, $m=n=0$ represents the reduced density operator, see ADOs for more details.

One can directly calculate the expectation values using the function Expect together with the output of evolution:

A::AbstractMatrix # observable
ados_list = evolution(...) # the input parameters depend on the different methods you choose.

Elist = Expect(A, ados_list)

Here, Elist contains the expectation values corresponding to the ados_list (i.e., the reduced density operator in each time step).

Common and optional parameters

Furthermore, there are two common optional parameters for all the methods provided below:

verbose::Bool : To display verbose output and progress bar during the process or not. Defaults to true.
filename::String : If filename was specified, the ADOs at each time point will be saved into the JLD2 file during the solving process. Default to Empty String: "".

If the filename is specified, the function will automatically save (update) the ADOs to the file (with ".jld2" behind the filename) once it obtains the solution of the specified time step. The saving method is based on the package JLD2.jl, which saves and loads Julia data structures in a format comprising a subset of HDF5.

tlist = 0:0.5:5
ados_list = evolution(..., tlist, ...; filename="test", ...)

The solution of the ADOs for each time step in tlist is saved in the file named "test.jld2".

To retrieve the solution (ADOs) from a previously saved file "text.jld2", just read the file with the methods provided by JLD2.jl. The solution for a specific time step can be extract by using the string of the time step as the "key". For example, if you want to obtain the solution at time 1.5, which is one of the time steps in tlist:

using HierarchicalEOM, JLD2 # rembember to import these before retrieving the solution

t = 1.5
filename = "test.jld2"
jldopen(filename, "r") do file
    ados = file[string(t)]
end

Ordinary Differential Equation Method

The first method is implemented by solving the ordinary differential equation (ODE) as shown above. HierarchicalEOM.jl wraps some of the functions in DifferentialEquations.jl, which is a very rich numerical library for solving the differential equations and provides many ODE solvers. It offers quite a few options for the user to tailor the solver to their specific needs. The default solver (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to benchmarks for DifferentialEquations solvers and also the documentation of DifferentialEquations.jl.

Extension for CUDA.jl

HierarchicalEOM.jl provides an extension to support GPU (CUDA.jl) acceleration for solving the time evolution (only for ODE method with time-independent system Hamiltonian), but this feature requires Julia 1.9+ and HierarchicalEOM 1.1+. See here for more details.

Given the initial state as Density Operator (`AbstractMatrix` type)

See the docstring of this method:

HierarchicalEOM.HeomAPI.evolution — Method

evolution(M, ρ0, tlist; solver, reltol, abstol, maxiters, save_everystep, verbose, filename, SOLVEROptions...)

Solve the time evolution for auxiliary density operators based on ordinary differential equations with initial state is given in the type of density-matrix (ρ0).

Parameters

M::AbstractHEOMLSMatrix : the matrix given from HEOM model
ρ0 : system initial state (density matrix)
tlist::AbstractVector : Denote the specific time points to save the solution at, during the solving process.
solver : solver in package DifferentialEquations.jl. Default to DP5().
reltol::Real : Relative tolerance in adaptive timestepping. Default to 1.0e-6.
abstol::Real : Absolute tolerance in adaptive timestepping. Default to 1.0e-8.
maxiters::Real : Maximum number of iterations before stopping. Default to 1e5.
save_everystep::Bool : Saves the result at every step. Defaults to false.
verbose::Bool : To display verbose output and progress bar during the process or not. Defaults to true.
filename::String : If filename was specified, the ADOs at each time point will be saved into the JLD2 file "filename.jld2" during the solving process.
SOLVEROptions : extra options for solver

For more details about solvers and extra options, please refer to DifferentialEquations.jl

Returns

ADOs_list : The auxiliary density operators in each time point.

source

# the time-independent HEOMLS matrix
M::AbstractHEOMLSMatrix  

# the initial state of the system density operator
ρ0::AbstractMatrix

# specific time points to save the solution during the solving process.  
tlist = 0:0.5:2 # [0.0, 0.5, 1.0, 1.5, 2.0]

ados_list = evolution(M, ρ0, tlist)

Given the initial state as Auxiliary Density Operators

Note

This method is usually used when you want to solve the time evolution again with the initial state are given from the last time point of the previous result.

See the docstring of this method:

HierarchicalEOM.HeomAPI.evolution — Method

evolution(M, ados, tlist; solver, reltol, abstol, maxiters, save_everystep, verbose, filename, SOLVEROptions...)

Solve the time evolution for auxiliary density operators based on ordinary differential equations with initial state is given in the type of ADOs.

Parameters

M::AbstractHEOMLSMatrix : the matrix given from HEOM model
ados::ADOs : initial auxiliary density operators
tlist::AbstractVector : Denote the specific time points to save the solution at, during the solving process.
solver : solver in package DifferentialEquations.jl. Default to DP5().
reltol::Real : Relative tolerance in adaptive timestepping. Default to 1.0e-6.
abstol::Real : Absolute tolerance in adaptive timestepping. Default to 1.0e-8.
maxiters::Real : Maximum number of iterations before stopping. Default to 1e5.
save_everystep::Bool : Saves the result at every step. Defaults to false.
verbose::Bool : To display verbose output and progress bar during the process or not. Defaults to true.
filename::String : If filename was specified, the ADOs at each time point will be saved into the JLD2 file "filename.jld2" during the solving process.
SOLVEROptions : extra options for solver

For more details about solvers and extra options, please refer to DifferentialEquations.jl

Returns

ADOs_list : The auxiliary density operators in each time point.

source

# the time-independent HEOMLS matrix
M::AbstractHEOMLSMatrix  

# the initial state of the ADOs (usually obtianed from previous solving result)
ados::ADOs      

# specific time points to save the solution during the solving process. 
tlist = 0:0.5:2 # [0.0, 0.5, 1.0, 1.5, 2.0]

ados_list = evolution(M, ados, tlist)

Propagator Method

The second method is implemented by directly construct the propagator of a given HEOMLS matrix $\hat{\mathcal{M}}$. Because $\hat{\mathcal{M}}$ is time-independent, the equation above can be solved analytically as

\[\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t)=\hat{\mathcal{G}}(t)\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(0),\]

where $\hat{\mathcal{G}}(t)\equiv \exp(\hat{\mathcal{M}}t)$ is the propagator for all ADOs corresponding to $\hat{\mathcal{M}}$.

To construct the propagator, we wrap the function in the package fastExpm.jl, which is optimized for the exponentiation of either large-dense or sparse matrices.

Given the initial state as Density Operator (`AbstractMatrix` type)

See the docstring of this method:

HierarchicalEOM.HeomAPI.evolution — Method

evolution(M, ρ0, Δt, steps; threshold, nonzero_tol, verbose, filename)

Solve the time evolution for auxiliary density operators based on propagator (generated by FastExpm.jl) with initial state is given in the type of density-matrix (ρ0).

This method will return the time evolution of ADOs corresponds to tlist = 0 : Δt : (Δt * steps)

Parameters

M::AbstractHEOMLSMatrix : the matrix given from HEOM model
ρ0 : system initial state (density matrix)
Δt::Real : A specific time step (time interval).
steps::Int : The number of time steps
threshold::Real : Determines the threshold for the Taylor series. Defaults to 1.0e-6.
nonzero_tol::Real : Strips elements smaller than nonzero_tol at each computation step to preserve sparsity. Defaults to 1.0e-14.
verbose::Bool : To display verbose output and progress bar during the process or not. Defaults to true.
filename::String : If filename was specified, the ADOs at each time point will be saved into the JLD2 file "filename.jld2" during the solving process.

For more details, please refer to FastExpm.jl

Returns

ADOs_list : The auxiliary density operators of each time step.

source

# the time-independent HEOMLS matrix
M::AbstractHEOMLSMatrix  

# the initial state of the system density operator
ρ0::AbstractMatrix

# A specific time interval (time step)
Δt = 0.5

# The number of time steps for the propagator to apply
steps = 4

# equivalent to tlist = 0 : Δt : (Δt * steps)
ados_list = evolution(M, ρ0, Δt, steps)

Given the initial state as Auxiliary Density Operators

Note

This method is usually used when you want to solve the time evolution again with the initial state are given from the last time point of the previous result.

See the docstring of this method:

HierarchicalEOM.HeomAPI.evolution — Method

evolution(M, ados, Δt, steps; threshold, nonzero_tol, verbose, filename)

Solve the time evolution for auxiliary density operators based on propagator (generated by FastExpm.jl) with initial state is given in the type of ADOs.

This method will return the time evolution of ADOs corresponds to tlist = 0 : Δt : (Δt * steps)

Parameters

M::AbstractHEOMLSMatrix : the matrix given from HEOM model
ados::ADOs : initial auxiliary density operators
Δt::Real : A specific time step (time interval).
steps::Int : The number of time steps
threshold::Real : Determines the threshold for the Taylor series. Defaults to 1.0e-6.
nonzero_tol::Real : Strips elements smaller than nonzero_tol at each computation step to preserve sparsity. Defaults to 1.0e-14.
verbose::Bool : To display verbose output and progress bar during the process or not. Defaults to true.
filename::String : If filename was specified, the ADOs at each time point will be saved into the JLD2 file "filename.jld2" during the solving process.

For more details, please refer to FastExpm.jl

Returns

ADOs_list : The auxiliary density operators of each time step.

source

# the time-independent HEOMLS matrix
M::AbstractHEOMLSMatrix  

# the initial state of the ADOs (usually obtianed from previous solving result)
ados::ADOs

# A specific time interval (time step)
Δt = 0.5

# The number of time steps for the propagator to apply
steps = 4

# equivalent to tlist = 0 : Δt : (Δt * steps)
ados_list = evolution(M, ados, Δt, steps)

Time Dependent Problems

In general, the time-dependent system Hamiltonian can be separated into the time-independent and time-dependent parts, namely

\[H_s (t) = H_0 + H_1(t).\]

We again wrap some of the functions in DifferentialEquations.jl to solve the time-dependent problems here.

To deal with the time-dependent system Hamiltonian problem in HierarchicalEOM.jl, we first construct the HEOMLS matrices $\hat{\mathcal{M}}$ with time-independent Hamiltonian $H_0$:

M = M_S(H0, ...)
M = M_Boson(H0, ...)
M = M_Fermion(H0, ...)
M = M_BosonFermion(H0, ...)

To solve the dynamics characterized by $\hat{\mathcal{M}}$ together with the time-dependent part of system Hamiltonian $H_1(t)$, you can call either of the following two functions (one takes the type of initial state as density matrix and the other one takes ADOs):

evolution(M::AbstractHEOMLSMatrix, ρ0, tlist::AbstractVector, H::Function, param::Tuple = ())
evolution(M::AbstractHEOMLSMatrix, ados::ADOs, tlist::AbstractVector, H::Function, param::Tuple = ()).

Here, the definition of user-defined function H must be in the form H(p::Tuple, t) and returns the time-dependent part of system Hamiltonian (in AbstractMatrix type) at any given time point t. The parameter p should be a Tuple which contains all the extra parameters you need for the function H. For example:

function H_pump(p, t)  
    p0, p1, p2 = p 
    # in this case, p should be passed in as a tuple: (p0, p1, p2) 
    
    σx = [0 1; 1 0] # Pauli-X matrix
    return (sin(p0 * t) + sin(p1 * t) + sin(p2 * t)) * σx
end

The parameter tuple p will be passed to your function H directly from one of the parameter in evolution called param:

M::AbstractHEOMLSMatrix
ρ0::AbstractMatrix
tlist = 0:0.1:10
p = (0.1, 1, 10)

ados_list = evolution(M, ρ0, tlist, H_pump, p)

Warning

If you don't need any extra param in your case, you still need to put a redundant one in the definition of H, for example:

function H_pump(p, t)
    σx = [0 1; 1 0] # Pauli-X matrix
    return sin(0.1 * t) * σx
end

M::AbstractHEOMLSMatrix
ρ0::AbstractMatrix
tlist = 0:0.1:10

ados_list = evolution(M, ρ0, tlist, H_pump)

Note

The default value for param in evolution is an empty tuple ().

Given the initial state as Density Operator (`AbstractMatrix` type)

See the docstring of this method:

HierarchicalEOM.HeomAPI.evolution — Function

evolution(M, ρ0, tlist, H, param; solver, reltol, abstol, maxiters, save_everystep, verbose, filename, SOLVEROptions...)

Solve the time evolution for auxiliary density operators with time-dependent system Hamiltonian based on ordinary differential equations with initial state is given in the type of density-matrix (ρ0).

Parameters

M::AbstractHEOMLSMatrix : the matrix given from HEOM model (with time-independent system Hamiltonian)
ρ0 : system initial state (density matrix)
tlist::AbstractVector : Denote the specific time points to save the solution at, during the solving process.
H::Function : a function for time-dependent part of system Hamiltonian. The function will be called by H(param, t) and should return the time-dependent part system Hamiltonian matrix at time t with AbstractMatrix type.
param::Tuple: the tuple of parameters which is used to call H(param, t) for the time-dependent system Hamiltonian. Default to empty tuple ().
solver : solver in package DifferentialEquations.jl. Default to DP5().
reltol::Real : Relative tolerance in adaptive timestepping. Default to 1.0e-6.
abstol::Real : Absolute tolerance in adaptive timestepping. Default to 1.0e-8.
maxiters::Real : Maximum number of iterations before stopping. Default to 1e5.
save_everystep::Bool : Saves the result at every step. Defaults to false.
verbose::Bool : To display verbose output and progress bar during the process or not. Defaults to true.
filename::String : If filename was specified, the ADOs at each time point will be saved into the JLD2 file "filename.jld2" during the solving process.
SOLVEROptions : extra options for solver

For more details about solvers and extra options, please refer to DifferentialEquations.jl

Returns

ADOs_list : The auxiliary density operators in each time point.

source

Given the initial state as Auxiliary Density Operators

Note

This method is usually used when you want to solve the time evolution again with the initial state are given from the last time point of the previous result.

See the docstring of this method:

HierarchicalEOM.HeomAPI.evolution — Function

evolution(M, ados, tlist, H, param; solver, reltol, abstol, maxiters, save_everystep, verbose, filename, SOLVEROptions...)

Solve the time evolution for auxiliary density operators with time-dependent system Hamiltonian based on ordinary differential equations with initial state is given in the type of ADOs.

Parameters

M::AbstractHEOMLSMatrix : the matrix given from HEOM model (with time-independent system Hamiltonian)
ados::ADOs : initial auxiliary density operators
tlist::AbstractVector : Denote the specific time points to save the solution at, during the solving process.
H::Function : a function for time-dependent part of system Hamiltonian. The function will be called by H(param, t) and should return the time-dependent part system Hamiltonian matrix at time t with AbstractMatrix type.
param::Tuple: the tuple of parameters which is used to call H(param, t) for the time-dependent system Hamiltonian. Default to empty tuple ().
solver : solver in package DifferentialEquations.jl. Default to DP5().
reltol::Real : Relative tolerance in adaptive timestepping. Default to 1.0e-6.
abstol::Real : Absolute tolerance in adaptive timestepping. Default to 1.0e-8.
maxiters::Real : Maximum number of iterations before stopping. Default to 1e5.
save_everystep::Bool : Saves the result at every step. Defaults to false.
verbose::Bool : To display verbose output and progress bar during the process or not. Defaults to true.
filename::String : If filename was specified, the ADOs at each time point will be saved into the JLD2 file "filename.jld2" during the solving process.
SOLVEROptions : extra options for solver

For more details about solvers and extra options, please refer to DifferentialEquations.jl

Returns

ADOs_list : The auxiliary density operators in each time point.

source