### Build Documentation with Make

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/README.md

This command initiates the documentation build process using the 'make' utility. It assumes 'make' is installed and configured for the project. The output will be generated in the 'docs/build/' directory.

```shell
make docs
```

--------------------------------

### Simulate System Interacting with Multiple Baths

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

This example shows how to set up and run a simulation where a quantum system interacts with multiple baths. It defines system Hamiltonian, initial state, projectors, and then constructs a list of baths. The HEOMsolve function is used to obtain the time evolution of populations, which are then plotted. Dependencies include Plots, LaTeXStrings, and HierarchicalEOM.jl specific functions like BosonBath, Boson_DrudeLorentz_Pade, M_Boson, and HEOMsolve.

```julia
# The system Hamiltonian
Hsys = Qobj([
    0.25 1.50 2.50
    1.50 0.75 3.50
    2.50 3.50 1.25
])

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

# Projector for each system state:
P00 = ket2dm(basis(3, 0))
P11 = ket2dm(basis(3, 1))
P22 = ket2dm(basis(3, 2));

# Construct one bath for each system state:
# note that `BosonBath[]` make the list created in type: Vector{BosonBath}
baths = BosonBath[]
for i in 0:2
    # system-bath coupling operator: |i><i|
    Q = ket2dm(basis(3, i))
    push!(baths, Boson_DrudeLorentz_Pade(Q, λ, W, kT, N))
end

L = M_Boson(Hsys, tier, baths)

tlist = 0:0.025:5
sol = HEOMsolve(L, ρ0, tlist; e_ops = [P00, P11, P22])

# calculate population for each system state:
p0 = real(sol.expect[1, :])
p1 = real(sol.expect[2, :])
p2 = real(sol.expect[3, :])

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")
```

--------------------------------

### Install and Load HierarchicalEOM.jl Package

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Installs the HierarchicalEOM.jl package using Julia's package manager and demonstrates how to load it and check its version information.

```julia
using Pkg
Pkg.add("HierarchicalEOM")

# Load the package and check version information
using HierarchicalEOM
HierarchicalEOM.versioninfo()
```

--------------------------------

### Import HierarchicalEOM.jl Functions

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

Imports core functions from the HierarchicalEOM.jl library for use in simulations. It also shows how to import all functions using 'using HierarchicalEOM' and check dependency versions.

```julia
import HierarchicalEOM
import HierarchicalEOM: Boson_DrudeLorentz_Pade, M_Boson, HEOMsolve, getRho, BosonBath

# To import everything:
# using HierarchicalEOM

# To check dependency versions:
HierarchicalEOM.versioninfo()
```

--------------------------------

### Install HierarchicalEOM.jl Package

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/README.md

Installs the HierarchicalEOM.jl package using Julia's package manager. This can be done directly in the Julia REPL or within the Pkg REPL.

```julia
using Pkg
Pkg.add("HierarchicalEOM")
```

```julia-REPL
(1.10) pkg> add HierarchicalEOM
```

--------------------------------

### Define Quantum System and Operators with QuantumToolbox.jl

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

Defines a two-level system Hamiltonian, initial state, and measurement operators using the QuantumToolbox.jl framework. This includes creating Qobj objects for operators like sigma_z, sigma_x, and density matrices.

```julia
import QuantumToolbox: Qobj, sigmaz, sigmax, basis, ket2dm, expect, steadystate

# System parameters
ϵ = 0.5 # energy of 2-level system
Δ = 1.0 # tunneling term

# System Hamiltonian
Hsys = 0.5 * ϵ * sigmaz() + 0.5 * Δ * sigmax()

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

# Measurement operators
P00 = ket2dm(basis(2, 0))
P11 = ket2dm(basis(2, 1))
P01 = basis(2, 0) * basis(2, 1)'
```

--------------------------------

### Develop and Instantiate Project Dependencies with Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/README.md

This Julia command sets up the project environment for documentation building. It adds the local HierarchicalEOM package as a development dependency and then instantiates all necessary packages. This ensures the documentation is built against the local code.

```julia
julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
```

--------------------------------

### Build Documentation Manually with Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/README.md

This command manually triggers the documentation build process using Julia. It executes the 'docs/make.jl' script, which is responsible for generating the documentation files. This method is an alternative to using 'make'.

```julia
julia --project=docs docs/make.jl
```

--------------------------------

### Load and Get Version Info for HierarchicalEOM.jl

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/README.md

Loads the HierarchicalEOM.jl package into the current Julia session and displays version information. This is useful for verifying the installation and understanding the package's current state.

```julia
using HierarchicalEOM
HierarchicalEOM.versioninfo()
HierarchicalEOM.about()
```

--------------------------------

### Setup for Demonstrating Memory Savings (Julia)

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/heom_matrix/lazy_operators.md

Sets up the necessary environment for demonstrating memory savings in HierarchicalEOM.jl. It uses the SparseArrays package and defines a spin-boson model with specific parameters for the system and baths.

```julia
using HierarchicalEOM
using SparseArrays

# spin-boson model
Hsys = 1.0 * sigmaz() + 0.25 * sigmax()

Δ  = 0.1 # coupling strength
W  = 0.2 # band-width  of the environment
kT = 0.5 # the product of the Boltzmann constant k and the absolute temperature T
N  = 10  # Number of exponential terms
baths = [
    Boson_DrudeLorentz_Pade(sigmax(), 0.01, 0.5, 0.5, 3),
    Boson_DrudeLorentz_Pade(sigmax(), 0.01, 0.5, 0.5, 3),
]
```

--------------------------------

### Construct HEOM Liouvillian Superoperator Matrix for Bosonic Baths

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

This snippet demonstrates the setup for constructing the HEOM Liouvillian superoperator matrix (HEOMLS) for a system coupled to bosonic baths. It initializes a system Hamiltonian using Pauli matrices and standard operators.

```julia
using HierarchicalEOM
using QuantumToolbox

# System setup
ϵ = 0.5
Δ = 1.0
Hsys = 0.5 * ϵ * sigmaz() + 0.5 * Δ * sigmax()
```

--------------------------------

### Compare Memory Usage of Assembly Modes (Julia)

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/heom_matrix/lazy_operators.md

This example uses `Base.summarysize` to calculate and compare the memory usage in megabytes (MB) of matrices created with different assembly modes (`L_full`, `L_combine`, `L_none`). It prints detailed information about the system dimension, number of ADOs, matrix dimensions, non-zero elements, and the percentage of memory reduction achieved by using `:combine` and `:none` modes compared to `:full`.

```julia
mem_full = Base.summarysize(L_full.data) / 1024^2
mem_combine = Base.summarysize(L_combine.data) / 1024^2
mem_none = Base.summarysize(L_none.data) / 1024^2

print(
    "System dimension (d)     : $(size(Hsys, 1))\n",
    "Number of ADOs   (N_ADO) : $(L_full.N)\n",
    "L_full matrix size       : $(size(L_full, 1)) × $(size(L_full, 2))\n",
    "L_full non-zero elements : $(nnz(L_full.data.A))\n",
    "\n",
    "Memory usage:\n",
    "  full    : $(round(mem_full, digits=3)) MB\n",
    "  combine : $(round(mem_combine, digits=3)) MB\n",
    "  none    : $(round(mem_none, digits=3)) MB\n",
    "Memory reduction:\n",
    "  combine vs full : $(round(100 * (1 - mem_combine/mem_full), digits=1))%\n",
    "   none   vs full : $(round(100 * (1 - mem_none/mem_full), digits=1))%\n",
)
```

--------------------------------

### Define Bosonic Bath with Pade Expansion

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

Constructs a bosonic bath spectral density using the Drude-Lorentz model with Pade expansion. This function requires system-bath coupling operator, coupling strength, bandwidth, temperature, and the number of expansion terms.

```julia
# Bath parameters
λ = 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 for the bath
bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N)
```

--------------------------------

### Calculate and Plot Time Evolution and Steady State Populations

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

This snippet demonstrates how to extract population data from a simulation result (sol.expect) for both time evolution and steady-state. It then uses the Plots.jl library to visualize these populations over time, including dashed lines for steady-state values. Dependencies include Plots and LaTeXStrings.

```julia
# for time evolution
p00_e = real(sol.expect[1, :]) # P00 is the 1st element in e_ops
p01_e = real(sol.expect[3, :]); # P01 is the 3rd element in e_ops

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


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")
```

--------------------------------

### Setting up the Single-Impurity Anderson Model

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/extensions/CUDA.md

Provides the Julia code for setting up the parameters and Hamiltonian for the single-impurity Anderson model, a common example used to demonstrate the HierarchicalEOM.jl CUDA extension.

```julia
ϵ  = -5
U  = 10
Γ  = 2
μ  = 0
W  = 10
kT = 0.5
N  = 5
tier = 3

tlist = 0:0.1:10
ωlist = -10:1:10

σm = sigmam()
σz = sigmaz()
II = qeye(2)
d_up = tensor(     σm, II)
d_dn = tensor(-1 * σz, σm)
ψ0   = tensor(basis(2, 0), basis(2, 0))
Hsys = ϵ * (d_up' * d_up + d_dn' * d_dn) + U * (d_up' * d_up * d_dn' * d_dn)

bath_up = Fermion_Lorentz_Pade(d_up, Γ, μ, W, kT, N)
bath_dn = Fermion_Lorentz_Pade(d_dn, Γ, μ, W, kT, N)
bath_list = [bath_up, bath_dn]

# even HEOMLS matrix
M_even_cpu = M_Fermion(Hsys, tier, bath_list)

# odd HEOMLS matrix
M_odd_cpu  = M_Fermion(Hsys, tier, bath_list, ODD)
```

--------------------------------

### Perform Time Evolution with HEOMsolve

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

Calculates the time evolution of the system's auxiliary density operators (ADOs) using the HEOMsolve function. It requires the Liouvillian superoperator, initial state, time list, and optionally, expectation operators.

```julia
tlist = 0:0.2:50
sol = HEOMsolve(L, ρ0, tlist; e_ops = [P00, P11, P01])
```

--------------------------------

### Solve for Stationary State

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

Solves for the stationary state of the auxiliary density operators (ADOs) using the steadystate function. This function takes the Liouvillian superoperator as input.

```julia
ados_steady = steadystate(L)
```

--------------------------------

### Construct HEOM Liouvillian Superoperator

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

Builds the HEOM Liouvillian superoperator matrix for a bosonic bath. This function takes the system Hamiltonian, the maximum tier of the hierarchy, and the defined bath object as input.

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

--------------------------------

### Get Total Number of ADOs in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/ADOs.md

Demonstrates how to get the total number of auxiliary density operators stored within an ADOs object using the length function in Julia. This is equivalent to accessing the ADOs.N field.

```julia
ados::ADOs
length(ados)
```

--------------------------------

### Construct Fermion Lorentz Padé Bath in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/bath_fermion/Fermion_Lorentz.md

Constructs a `Fermion_Lorentz_Pade` bath object for simulating fermionic environments using the Padé expansion. Requires system coupling operator, environmental parameters (coupling strength, chemical potential, band-width, temperature), and the number of terms for the expansion.

```julia
ds # coupling operator
Γ  # coupling strength
μ  # chemical potential of the environment
W  # band-width  of the environment
kT # the product of the Boltzmann constant k and the absolute temperature T
N  # Number of exponential terms for each correlation functions (C^{+} and C^{-})
bath = Fermion_Lorentz_Pade(ds, Γ, μ, W, kT, N - 1)
```

--------------------------------

### Get HierarchicalEOM Version Information

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/libraryAPI.md

Retrieves version information for the HierarchicalEOM library. This is useful for tracking dependencies and ensuring compatibility.

```julia
HierarchicalEOM.versioninfo()
```

--------------------------------

### Get HierarchicalEOM Citation Information

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/libraryAPI.md

Provides citation information for the HierarchicalEOM library, allowing users to properly reference the software in their publications.

```julia
HierarchicalEOM.cite()
```

--------------------------------

### Get HierarchicalEOM About Information

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/libraryAPI.md

Displays 'about' information for the HierarchicalEOM library, potentially including author details, license, or other relevant metadata.

```julia
HierarchicalEOM.about()
```

--------------------------------

### Get Length and Element Type of ADOs

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/libraryAPI.md

Returns the number of auxiliary density operators in an ADOs object and the data type of its elements.

```julia
Base.length(A::ADOs)
Base.eltype(A::ADOs)
```

--------------------------------

### Get Size of HEOM Superoperator Matrix

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/libraryAPI.md

Returns the dimensions of an HEOMSuperOp or AbstractHEOMLSMatrix object. Can query the total size or the size along a specific dimension.

```julia
Base.size(M::HEOMSuperOp)
Base.size(M::HEOMSuperOp, dim::Int)
Base.size(M::AbstractHEOMLSMatrix)
Base.size(M::AbstractHEOMLSMatrix, dim::Int)
```

--------------------------------

### Get Index Ensemble from Hierarchy Dictionary

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/libraryAPI.md

Retrieves an ensemble of indices from a HierarchyDict or MixHierarchyDict. This is useful for iterating through or selecting specific parts of the hierarchical state space.

```julia
getIndexEnsemble
```

--------------------------------

### Construct Fermionic Bath with Lorentz Spectral Density (Pade/Matsubara)

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Constructs a fermionic bath using the Lorentz spectral density with Pade or Matsubara expansion. Parameters include the system annihilation operator, coupling strength, chemical potential, band-width, temperature, and the number of exponential terms per correlation function.

```julia
using HierarchicalEOM
using QuantumToolbox

# Fermionic bath parameters
ds = sigmam()  # system annihilation operator
Γ  = 2.0       # coupling strength
μ  = 0.0       # chemical potential
W  = 10.0      # band-width
kT = 0.5       # temperature
N  = 5         # number of exponential terms per correlation function

# Create bath with Pade expansion
bath = Fermion_Lorentz_Pade(ds, Γ, μ, W, kT, N)

# Alternative: Matsubara expansion
bath_matsubara = Fermion_Lorentz_Matsubara(ds, Γ, μ, W, kT, N)

print(bath)  # Shows number of exponential terms (always even for fermionic baths)
```

--------------------------------

### Get Element Type of HEOM Matrix

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/libraryAPI.md

Returns the data type of the elements within an HEOMSuperOp or AbstractHEOMLSMatrix. This is useful for understanding the numerical precision and type of calculations performed.

```julia
Base.eltype(M::HEOMSuperOp)
Base.eltype(M::AbstractHEOMLSMatrix)
```

--------------------------------

### Cite QuantumToolbox.jl

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/README.md

BibTeX entry for citing the QuantumToolbox.jl package, a dependency of HierarchicalEOM.jl. Citing this package is recommended as HierarchicalEOM.jl is built upon its functionalities for simulating open quantum systems.

```bibtex
@article{QuantumToolbox.jl2025,
  title = {Quantum{T}oolbox.jl: {A}n efficient {J}ulia framework for simulating open quantum systems},
  author = {Mercurio, Alberto and Huang, Yi-Te and Cai, Li-Xun and Chen, Yueh-Nan and Savona, Vincenzo and Nori, Franco},
  journal = {{Quantum}},
  issn = {2521-327X},
  publisher = {{Verein zur F{"{o}}rderung des Open Access Publizierens in den Quantenwissenschaften}},
  volume = {9},
  pages = {1866},
  month = sep,
  year = {2025},
  doi = {10.22331/q-2025-09-29-1866},
  url = {https://doi.org/10.22331/q-2025-09-29-1866}
}
```

--------------------------------

### Construct Bosonic Bath with Drude-Lorentz Spectral Density (Pade/Matsubara)

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Constructs a bosonic bath using the Drude-Lorentz spectral density, employing either Pade or Matsubara expansion for the correlation function. Requires system coupling operator, coupling strength, band-width, temperature, and the number of exponential terms.

```julia
using HierarchicalEOM
using QuantumToolbox

# Bath parameters
Vs = sigmaz()  # system-bath coupling operator
Δ  = 0.1       # coupling strength
W  = 0.5       # band-width (cut-off frequency)
kT = 0.5       # temperature (k_B * T)
N  = 4         # number of exponential terms

# Create bath with Pade expansion
bath = Boson_DrudeLorentz_Pade(Vs, Δ, W, kT, N)

# Alternative: Matsubara expansion
bath_matsubara = Boson_DrudeLorentz_Matsubara(Vs, Δ, W, kT, N)

# Verify correlation function
tlist = 0:0.01:5
c_values = correlation_function(bath, tlist)
```

--------------------------------

### Get Rho and ADO from ADOs Object

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/libraryAPI.md

Retrieves the density matrix or a specific auxiliary density operator from an ADOs object. These functions are essential for accessing and manipulating the state information within the HEOM framework.

```julia
getRho
getADO
```

--------------------------------

### Construct Underdamped Bosonic Bath using Matsubara Expansion (Julia)

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/bath_boson/Boson_Underdamped.md

This code snippet demonstrates how to construct an underdamped bosonic bath using the `Boson_Underdamped_Matsubara` function in Julia. It requires parameters defining the system-bath coupling, bath properties, and temperature. The function returns a bath object representing the spectral density.

```julia
Vs # coupling operator
Δ  # coupling strength
W  # band-width of the environment
ω0 # resonance frequency of the environment
kT # the product of the Boltzmann constant k and the absolute temperature T
N  # Number of exponential terms
bath = Boson_Underdamped_Matsubara(Vs, Δ, W, ω0, kT, N - 2)
```

--------------------------------

### Extract Reduced Density Operator

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/quick_start.md

Retrieves the reduced density operator from the results of time evolution or stationary state calculation. The reduced density operator is typically the first element of the auxiliary density operator list.

```julia
# From time evolution results
ados_list = sol.ados
ρ_final = ados_list[end][1]  # index `1` represents the reduced density operator
ρ_final_alt = getRho(ados_list[end])

# From stationary state results
ρ_steady = ados_steady[1]
ρ_steady_alt = getRho(ados_steady)
```

--------------------------------

### Pardiso Solvers (Intel MKL)

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/LS_solvers.md

Provides solvers based on Intel's MKL Pardiso library. Requires the Pardiso.jl package. Offers both factorization (MKLPardisoFactorize) and iterative (MKLPardisoIterate) approaches.

```julia
using Pardiso
using LinearSolve
MKLPardisoFactorize()
MKLPardisoIterate()
```

--------------------------------

### Import LinearSolve.jl

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/LS_solvers.md

Imports the LinearSolve.jl package, which provides algorithms for solving linear systems. This is a prerequisite for using any of the listed solvers.

```julia
using LinearSolve
```

--------------------------------

### Accessing TimeEvolutionHEOMSol Fields in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/time_evolution.md

Demonstrates how to access various fields within the `TimeEvolutionHEOMSol` object, which stores the results of a quantum system dynamics simulation. These fields include simulation parameters, time points, auxiliary density operators (ADOs), and calculated expectation values.

```julia
sol::TimeEvolutionHEOMSol

sol.Btier      # the tier (cutoff level) for bosonic hierarchy
sol.Ftier      # the tier (cutoff level) for fermionic hierarchy
sol.times      # The list of time points at which the expectation values are calculated during the evolution.
sol.times_ados # The list of time points at which the ADOs are stored during the evolution.
sol.ados       # The list of result ADOs corresponding to each time point in `times_ados`.
sol.expect     # The expectation values corresponding to each time point in `times`.
sol.retcode    # The return code from the solver.
sol.alg        # The algorithm which is used during the solving process.
sol.abstol     # The absolute tolerance which is used during the solving process.
sol.reltol     # The relative tolerance which is used during the solving process.
```

--------------------------------

### Construct Fermion Lorentz Matsubara Bath in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/bath_fermion/Fermion_Lorentz.md

Constructs a `Fermion_Lorentz_Matsubara` bath object for simulating fermionic environments using the Matsubara expansion. Requires system coupling operator, environmental parameters (coupling strength, chemical potential, band-width, temperature), and the number of terms for the expansion.

```julia
ds # coupling operator
Γ  # coupling strength
μ  # chemical potential of the environment
W  # band-width  of the environment
kT # the product of the Boltzmann constant k and the absolute temperature T
N  # Number of exponential terms for each correlation functions (C^{+} and C^{-})
bath = Fermion_Lorentz_Matsubara(ds, Γ, μ, W, kT, N - 1)
```

--------------------------------

### UMFPACKFactorization Solver

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/LS_solvers.md

Utilizes the UMFPACK factorization algorithm. This solver is recommended when the sparsity pattern of the system has significant structure, making it efficient for complex systems and baths.

```julia
UMFPACKFactorization()
```

--------------------------------

### Construct M_Fermion Matrix for Fermionic Baths (Julia)

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/heom_matrix/M_Fermion.md

Demonstrates how to construct the HEOM Liouvillian superoperator matrix for fermionic baths using the M_Fermion function. It shows examples for both EVEN and ODD parity.

```julia
Hs::QuantumObject # system Hamiltonian
tier = 3
Bath::FermionBath

# create HEOMLS matrix in both EVEN and ODD parity
M_even = M_Fermion(Hs, tier, Bath)
M_odd  = M_Fermion(Hs, tier, Bath, ODD)
```

--------------------------------

### Calculating Expectation Values using QuantumToolbox.expect in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/time_evolution.md

Illustrates an alternative method for calculating expectation values using the `QuantumToolbox.expect` function. This approach takes an observable and a list of ADOs (obtained from the simulation results) as input.

```julia
A::QuantumObject # observable
sol = HEOMsolve(...) # the input parameters depend on the different methods you choose.
ados_list = sol.ados

Elist = expect(A, ados_list)
```

--------------------------------

### Constructing HEOM Matrices with Parity - Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/Parity.md

Demonstrates how to construct HEOM matrices (M) for different bath types (bosonic, fermionic, or combined) while specifying the parity (EVEN or ODD) of the auxiliary density operators (ADOs). This is crucial for systems involving fermionic components where operator commutation relations depend on parity.

```julia
Hs::QuantumObject   # system Hamiltonian
Bbath::BosonBath    # bosonic   bath object
Fbath::FermionBath  # fermionic bath object
Btier::Int          # bosonic   truncation level 
Ftier::Int          # fermionic truncation level 

# create HEOMLS matrix in EVEN or ODD parity
M_even = M_S(Hs, EVEN)
M_odd  = M_S(Hs, ODD)

M_even = M_Boson(Hs, Btier, Bbath, EVEN) 
M_odd  = M_Boson(Hs, Btier, Bbath, ODD) 

M_even = M_Fermion(Hs, Ftier, Fbath, EVEN) 
M_odd  = M_Fermion(Hs, Ftier, Fbath, ODD)

M_even = M_Boson_Fermion(Hs, Btier, Ftier, Bbath, Fbath, EVEN) 
M_odd  = M_Boson_Fermion(Hs, Btier, Ftier, Bbath, Fbath, ODD)
```

--------------------------------

### Get Reduced Density Operator from ADOs in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/ADOs.md

Shows how to retrieve the system's reduced density operator from an ADOs object using the getRho function in Julia. The reduced density operator is a specific type of QuantumObject.

```julia
ados::ADOs
ρ = getRho(ados)
```

--------------------------------

### KLUFactorization Solver

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/LS_solvers.md

Employs the KLU factorization algorithm. This solver is suitable for systems with less structured sparsity patterns, offering better performance in such cases.

```julia
KLUFactorization()
```

--------------------------------

### Querying Mixed Bath Nvecs and ADO Indices using MixHierarchyDict in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/hierarchy_dictionary.md

Demonstrates how to use MixHierarchyDict in Julia to get bosonic and fermionic Nvecs from an ADO index, and vice versa. It also shows how to retrieve lists of ADO indices for specific bosonic or fermionic hierarchy levels.

```julia
# obtain the nvec(s) correspond to 10-th ADO
nvec_b, nvec_f = HDict.idx2nvec[10]

# obtain the index of the ADO corresponds to the given nvec
nvec_b::Nvec
nvec_f::Nvec
idx = HDict.nvec2idx[(nvec_b, nvec_f)]

# obtain a list of indices which corresponds to all ADOs in 3rd-bosonic-level of hierarchy
idx_list = HDict.Blvl2idx[3] 

# obtain a list of indices which corresponds to all ADOs in 4rd-fermionic-level of hierarchy
idx_list = HDict.Flvl2idx[4]
```

--------------------------------

### Construct Fermionic Bath with Exponential Decomposition

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Creates a FermionBath object for fermionic environments. It requires a system annihilation operator (must have odd parity) and exponential terms for absorption (η_absorb, γ_absorb) and emission (η_emit, γ_emit) processes. The code demonstrates bath creation and iteration through its terms, along with calculating correlation functions.

```julia
using HierarchicalEOM
using QuantumToolbox

# System annihilation operator (must have odd parity)
ds = sigmam()

# Exponential terms for absorption and emission processes
η_absorb = [0.5 + 0.1im, 0.3 + 0.0im]
γ_absorb = [0.1 - 0.05im, 0.5 - 0.05im]
η_emit   = [0.5 + 0.1im, 0.3 + 0.0im]
γ_emit   = [0.1 + 0.05im, 0.5 + 0.05im]  # complex conjugate of γ_absorb

bath = FermionBath(ds, η_absorb, γ_absorb, η_emit, γ_emit)

# Check bath information
print(bath)  # FermionBath object with 4 exponential-expansion terms

# Access exponential terms
for exp in bath
    println(exp)  # Types: "fA" (absorption) or "fE" (emission)
end

# Calculate correlation functions
tlist = 0:0.1:10
cp_list, cm_list = correlation_function(bath, tlist)  # C^+ and C^- correlation functions
```

--------------------------------

### Solving Steady State on CPU (ODE Method)

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/extensions/CUDA.md

Calculates the steady state of a quantum system using the `steadystate` function with the ODE method and a CPU-based HEOMLS matrix and initial state.

```julia
ados_ss = steadystate(M_even_cpu, ψ0);
```

--------------------------------

### Querying Nvec and ADO Indices using HierarchyDict in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/hierarchy_dictionary.md

Shows how to use the HierarchyDict in Julia to retrieve the Nvec object corresponding to a specific ADO index, and conversely, find the ADO index for a given Nvec. It also demonstrates how to get a list of ADO indices for a particular hierarchy level.

```julia
# obtain the nvec corresponds to 10-th ADO
nvec = HDict.idx2nvec[10]

# obtain the index of the ADO corresponds to the given nvec
nvec::Nvec
idx = HDict.nvec2idx[nvec]

# obtain a list of indices which corresponds to all ADOs in 3rd-level of hierarchy
idx_list = HDict.lvl2idx[3]
```

--------------------------------

### Construct FermionBath Object in Julia

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/bath_fermion/fermionic_bath_intro.md

Constructs a FermionBath object to describe the interaction between a quantum system and a fermionic environment. It requires the system's annihilation operator and lists of coefficients (eta) and decay rates (gamma) for absorption and emission processes. The lengths of these lists must be equal, and gamma lists must be complex conjugates.

```julia
bath = FermionBath(ds, η_absorb, γ_absorb, η_emit, γ_emit)
```

--------------------------------

### Construct BosonBath (Combined Real/Imaginary Parts)

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/bath_boson/bosonic_bath_intro.md

Constructs a BosonBath object when the real and imaginary parts of the correlation function are combined. Requires the system interaction operator `Vs` and lists of coefficients `η` and decay rates `γ` for the exponential terms. The lengths of `η` and `γ` must be identical.

```julia
bath = BosonBath(Vs, η, γ)
```

--------------------------------

### Work with Auxiliary Density Operators (ADOs)

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Demonstrates how to work with Auxiliary Density Operators (ADOs) which encode environmental memory effects in the HierarchicalEOM framework. ADOs can be obtained from steady-state calculations or time evolution. The code shows how to access ADO properties, extract the reduced density operator, and calculate expectation values.

```julia
using HierarchicalEOM
using QuantumToolbox

# Obtain ADOs from steady state or time evolution
Hsys = 0.5 * sigmaz() + 0.25 * sigmax()
bath = Boson_DrudeLorentz_Pade(sigmaz(), 0.1, 0.5, 0.5, 2)
M = M_Boson(Hsys, 5, bath)

ados = steadystate(M)

# Access ADOs properties
println("Number of ADOs: ", ados.N)
println("Dimensions: ", ados.dimensions)
println("Parity: ", ados.parity)

# Get reduced density operator
ρ = getRho(ados)  # Equivalent to ados[1]

# Access specific ADO by index
ado_1 = getADO(ados, 1)   # Reduced density operator
ado_10 = getADO(ados, 10) # 10th ADO

# Bracket access
ρ = ados[1]           # First ADO (reduced density operator)
ado = ados[10]        # 10th ADO
ado_range = ados[3:10] # ADOs from index 3 to 10
ado_last = ados[end]  # Last ADO

# Iteration
for (idx, ado) in enumerate(ados)
    if idx == 1
        println("Reduced density operator trace: ", tr(ado))
    end
end

# Calculate expectation value
A = sigmaz()
exp_val = expect(A, ados)
println("Expectation value of σz: ", real(exp_val))

# Expectation values for a list of ADOs
ados_list = [ados, ados]  # Example list
exp_list = expect(A, ados_list)
```

--------------------------------

### Construct Bosonic Bath with Exponential Decomposition

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Creates a BosonBath object by defining system coupling operators and exponential coefficients (η) and decay rates (γ) for correlation functions. Supports combined or separated real/imaginary parts for coefficients and provides methods to access and iterate through bath terms, as well as calculate the correlation function.

```julia
using HierarchicalEOM
using QuantumToolbox

# System coupling operator
Vs = sigmaz()

# Method 1: Combined real and imaginary parts
η = [0.5 + 0.1im, 0.3 + 0.0im, 0.2 + 0.0im]  # exponential coefficients
γ = [0.1 + 0.0im, 0.5 + 0.0im, 1.0 + 0.0im]  # decay rates
bath = BosonBath(Vs, η, γ)

# Method 2: Separated real and imaginary parts
η_real = [0.5, 0.3, 0.2]
γ_real = [0.1, 0.5, 1.0]
η_imag = [0.1, 0.05]
γ_imag = [0.2, 0.8]
bath_sep = BosonBath(Vs, η_real, γ_real, η_imag, γ_imag)

# Access exponential terms
print(bath)  # BosonBath object with 3 exponential-expansion terms
e = bath[1]  # Access first exponential term
for exp in bath
    println(exp)  # Iterate through all exponential terms
end

# Calculate correlation function
tlist = 0:0.1:10
c_list = correlation_function(bath, tlist)
```

--------------------------------

### Memory Optimization with Lazy Operators

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Demonstrates memory optimization using lazy evaluation for HEOMLS matrix representation. It compares the memory usage of full sparse matrices versus combined and individual lazy operators, showing significant memory reduction.

```julia
using HierarchicalEOM
using QuantumToolbox
using LinearSolve

# System setup
Hsys = 1.0 * sigmaz() + 0.25 * sigmax()
baths = [
    Boson_DrudeLorentz_Pade(sigmax(), 0.01, 0.5, 0.5, 3),
    Boson_DrudeLorentz_Pade(sigmax(), 0.01, 0.5, 0.5, 3),
]

tier = 10

# Full sparse matrix (default) - high memory usage
L_full = M_Boson(Hsys, tier, baths; assemble=Val(:full))

# Combined lazy operators (recommended) - low memory usage
L_lazy = M_Boson(Hsys, tier, baths; assemble=Val(:combine))

# Individual lazy operators (for analysis)
L_none = M_Boson(Hsys, tier, baths; assemble=Val(:none))

# Compare memory usage
mem_full = Base.summarysize(L_full.data) / 1024^2
mem_lazy = Base.summarysize(L_lazy.data) / 1024^2
println("Full matrix: $(round(mem_full, digits=2)) MB")
println("Lazy operators: $(round(mem_lazy, digits=2)) MB")
println("Memory reduction: $(round(100 * (1 - mem_lazy/mem_full), digits=1))%")

# Lazy operators work with all solvers
ρ0 = ket2dm(basis(2, 0))
tlist = 0:0.5:10
sol = HEOMsolve(L_lazy, ρ0, tlist)

# Steady state with lazy operators (requires matrix-free solver)
ados_ss = steadystate(L_lazy; solver=KrylovJL_GMRES())
```

--------------------------------

### Manage Hierarchy with HierarchyDict

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Explains the use of `HierarchyDict` for managing the mapping between Auxiliary Density Operator (ADO) indices and their corresponding hierarchy levels and bath exponential terms. It allows conversion between ADO indices and N-vectors, retrieval of indices at specific hierarchy levels, and analysis of ADOs based on bath information. This is crucial for understanding the structure of the HEOM hierarchy.

```julia
using HierarchicalEOM
using QuantumToolbox

# Setup
Hsys = 0.5 * sigmaz() + 0.25 * sigmax()
bath = Boson_DrudeLorentz_Pade(sigmaz(), 0.1, 0.5, 0.5, 3)
M = M_Boson(Hsys, 5, bath)

# Get hierarchy dictionary
HDict = M.hierarchy

# Convert between index and Nvec
nvec = HDict.idx2nvec[10]    # Get Nvec for 10th ADO
idx = HDict.nvec2idx[nvec]   # Convert back to index

# Get all ADO indices at a specific hierarchy level
level_3_indices = HDict.lvl2idx[3]  # All ADOs at level 3

# Analyze ADOs with bath information
ados = steadystate(M)

for (idx, ado) in enumerate(ados)
    nvec = HDict.idx2nvec[idx]

    # Get bath and exponential term indices
    for (α, k, n) in getIndexEnsemble(nvec, HDict.bathPtr)
        # α: bath index
        # k: exponential term index in α-th bath
        # n: repetition count
        exponent = M.bath[α][k]  # Access the exponential term
    end
end

# For mixed bosonic-fermionic systems
M_mixed = M_Boson_Fermion(Hsys, 3, 3, bath, Fermion_Lorentz_Pade(sigmam(), 1.0, 0.0, 5.0, 0.5, 2))
HDict_mixed = M_mixed.hierarchy

# Get both bosonic and fermionic Nvec
nvec_b, nvec_f = HDict_mixed.idx2nvec[10]

# Access level indices separately
boson_level_2 = HDict_mixed.Blvl2idx[2]   # Bosonic level 2
fermion_level_2 = HDict_mixed.Flvl2idx[2] # Fermionic level 2
```

--------------------------------

### Solving Steady State on CPU (Linear-Solve)

Source: https://github.com/qutip/hierarchicaleom.jl/blob/main/docs/src/extensions/CUDA.md

Calculates the steady state of a quantum system using the `steadystate` function with the linear-solve method and a CPU-based HEOMLS matrix.

```julia
ados_ss = steadystate(M_even_cpu);
```

--------------------------------

### Calculate Steady State using steadystate

Source: https://context7.com/qutip/hierarchicaleom.jl/llms.txt

Calculates the stationary state of a quantum system using the `steadystate` function. This can be done without an initial condition (linear solve) or with an initial density operator or auxiliary density operators (ADOs). The function returns the steady-state ADOs, from which the reduced density operator and expectation values can be extracted.

```julia
using HierarchicalEOM
using QuantumToolbox

# System setup
Hsys = 0.5 * 0.5 * sigmaz() + 0.5 * 1.0 * sigmax()
bath = Boson_DrudeLorentz_Pade(sigmaz(), 0.1, 0.5, 0.5, 2)
M = M_Boson(Hsys, 5, bath)

# Method 1: Linear solve (no initial condition required)
ados_ss = steadystate(M)

# Method 2: ODE evolution with density operator initial state
ρ0 = ket2dm(basis(2, 0))
ados_ss = steadystate(M, ρ0)

# Method 3: ODE evolution with ADOs initial state
initial_ados = ADOs(zeros(ComplexF64, M.N * M.sup_dim), M.N, M.dimensions, M.parity)
ados_ss = steadystate(M, initial_ados)

# Extract reduced density operator
ρ_steady = getRho(ados_ss)
println("Steady state trace: ", tr(ρ_steady))

# Calculate expectation values
P00 = ket2dm(basis(2, 0))
p00_ss = expect(P00, ados_ss)
println("P00 at steady state: ", real(p00_ss))
```