From 1db1848c2d25e0a98aea5d66d025b16da21c776a Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Tue, 8 Oct 2024 15:49:05 +0900
Subject: [PATCH 1/6] add documentation for `sesolve`

---
 .../src/users_guide/time_evolution/sesolve.md | 110 +++++++++++++++++-
 1 file changed, 109 insertions(+), 1 deletion(-)

diff --git a/docs/src/users_guide/time_evolution/sesolve.md b/docs/src/users_guide/time_evolution/sesolve.md
index 7df270a22..e3a06d754 100644
--- a/docs/src/users_guide/time_evolution/sesolve.md
+++ b/docs/src/users_guide/time_evolution/sesolve.md
@@ -2,4 +2,112 @@
 
 ## [Schrödinger Equation Solver](@id doc-TE:Schrödinger-Equation-Solver)
 
-This page is still under construction, please visit [API](@ref doc-API) first.
+### Unitary evolution
+
+The dynamics of a closed (pure) quantum system is governed by the Schrödinger equation
+
+```math
+i\hbar\frac{\partial}{\partial t}\Psi(\vec{x}, t) = \hat{H}\Psi(\vec{x}, t),
+```
+
+where ``\Psi(\vec{x}, t)`` is the wave function, ``\hat{H}`` is the Hamiltonian, and ``\hbar`` is reduced Planck constant. In general, the Schrödinger equation is a partial differential equation (PDE) where both 
+``\Psi`` and ``\hat{H}`` are functions of space ``\vec{x}`` and time ``t``. For computational purposes it is useful to expand the PDE in a set of basis functions that span the Hilbert space of the Hamiltonian, and to write the equation in matrix and vector form, namely
+
+```math
+i\hbar\frac{d}{dt}|\psi(t)\rangle = \hat{H}|\psi(t)\rangle,
+```
+
+where ``|\psi(t)\rangle`` is the state vector, and the Hamiltonian ``\hat{H}`` is now under matrix representation. This matrix equation can, in principle, be solved by diagonalizing the Hamiltonian matrix ``\hat{H}``. In practice, however, it is difficult to perform this diagonalization unless the size of the Hilbert space (dimension of the matrix ``\hat{H}``) is small. Analytically, it is a formidable task to calculate the dynamics for systems with more than two states. If, in addition, we consider dissipation due to the inevitable interaction with a surrounding environment, the computational complexity grows even larger, and we have to resort to numerical calculations in all realistic situations. This illustrates the importance of numerical calculations in describing the dynamics of open quantum systems, and the need for efficient and accessible tools for this task.
+
+The Schrödinger equation, which governs the time-evolution of closed quantum systems, is defined by its Hamiltonian and state vector. In the previous sections, [Manipulating States and Operators](@ref doc:Manipulating-States-and-Operators) and [Tensor Products and Partial Traces](@ref doc:Tensor-products-and-Partial-Traces), we showed how Hamiltonians and state vectors are constructed in `QuantumToolbox.jl`. Given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state vector ``|\psi(0)\rangle`` using the `QuantumToolbox` time evolution problem [`sesolveProblem`](@ref) or directly call the function [`sesolve`](@ref). It evolves the state vector ``|\psi(t)\rangle`` and evaluates the expectation values for a set of operators `e_ops` at each given time points, using an ordinary differential equation solver provided by the powerful julia package[`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/).
+
+### Example
+
+```@setup sesolve
+using QuantumToolbox
+```
+
+For example, the time evolution of a quantum spin-``\frac{1}{2}`` system (initialized in spin-``\uparrow``) with tunneling rate ``0.1`` is calculated, and the expectation values of the Pauli-Z operator ``\hat{\sigma}_z`` is also evaluated, with the following code
+
+```@example sesolve
+H = 2 * π * 0.1 * sigmax()
+ψ0 = basis(2, 0) # spin-up
+tlist = LinRange(0.0, 10.0, 20)
+
+prob = sesolveProblem(H, ψ0, tlist, e_ops = [sigmaz()])
+sol = sesolve(prob)
+```
+
+!!! note "Note"
+    Here, we generate the time evolution problem by [`sesolveProblem`](@ref) first and then put it into the function [`sesolve`](@ref). One can also directly call [`sesolve`](@ref), which we also demonstrates in below.
+
+The function returns an instance of [`TimeEvolutionSol`](@ref), as described in the previous section [Time Evolution Solutions](@ref doc-TE:Time-Evolution-Solutions). The attribute `expect` in `sol`ution is a list of expectation values for the operator(s) that are passed to the `e_ops` keyword argument. 
+
+```@example sesolve
+sol.expect
+```
+
+Passing multiple operators to `e_ops` as a `Vector` results in the expectation values for each operators at each time points.
+
+```@example sesolve
+tlist = LinRange(0.0, 10.0, 100)
+sol = sesolve(H, ψ0, tlist, e_ops = [sigmaz(), sigmay()])
+```
+
+!!! note "Note"
+    Here, we call [`sesolve`](@ref) directly instead of pre-defining [`sesolveProblem`](@ref) first (as shown previously).
+
+```@example sesolve
+times = sol.times
+print(size(times))
+```
+
+```@example sesolve
+expt = sol.expect
+print(size(expt))
+```
+
+We can therefore plot the expectation values:
+
+```@example sesolve
+using CairoMakie
+CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
+
+expt_z = real(expt[1,:])
+expt_y = real(expt[2,:])
+
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
+lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
+lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
+
+axislegend(ax, position = :rb)
+
+fig
+```
+
+If the keyword argument `e_ops` is not specified (or given as an empty `Vector`), the [`sesolve`](@ref) and functions return a [`TimeEvolutionSol`](@ref) that contains a list of state vectors which corresponds to the time points specified in `tlist`:
+
+```@example sesolve
+tlist = [0, 10]
+sol = sesolve(H, ψ0, tlist) # or specify: e_ops = []
+
+sol.states
+```
+
+This is because the `states` will be saved depend on the keyword argument `saveat` in `kwargs`. If `e_ops` is specified, the default value of `saveat=[tlist[end]]` (only save the final state), otherwise, `saveat=tlist` (saving the states corresponding to `tlist`). 
+
+One can also specify `e_ops` and `saveat` separately:
+
+```@example sesolve
+tlist = [0, 5, 10]
+sol = sesolve(H, ψ0, tlist, e_ops = [sigmay()], saveat = tlist)
+```
+
+```@example sesolve
+sol.expect
+```
+
+```@example sesolve
+sol.states
+```

From bf00a2895b2df3041dd48ae375648f46bb139482 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Tue, 8 Oct 2024 17:03:19 +0900
Subject: [PATCH 2/6] fix typo in `sesolve`

---
 docs/src/users_guide/time_evolution/sesolve.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/users_guide/time_evolution/sesolve.md b/docs/src/users_guide/time_evolution/sesolve.md
index e3a06d754..febe15d4c 100644
--- a/docs/src/users_guide/time_evolution/sesolve.md
+++ b/docs/src/users_guide/time_evolution/sesolve.md
@@ -19,7 +19,7 @@ i\hbar\frac{d}{dt}|\psi(t)\rangle = \hat{H}|\psi(t)\rangle,
 
 where ``|\psi(t)\rangle`` is the state vector, and the Hamiltonian ``\hat{H}`` is now under matrix representation. This matrix equation can, in principle, be solved by diagonalizing the Hamiltonian matrix ``\hat{H}``. In practice, however, it is difficult to perform this diagonalization unless the size of the Hilbert space (dimension of the matrix ``\hat{H}``) is small. Analytically, it is a formidable task to calculate the dynamics for systems with more than two states. If, in addition, we consider dissipation due to the inevitable interaction with a surrounding environment, the computational complexity grows even larger, and we have to resort to numerical calculations in all realistic situations. This illustrates the importance of numerical calculations in describing the dynamics of open quantum systems, and the need for efficient and accessible tools for this task.
 
-The Schrödinger equation, which governs the time-evolution of closed quantum systems, is defined by its Hamiltonian and state vector. In the previous sections, [Manipulating States and Operators](@ref doc:Manipulating-States-and-Operators) and [Tensor Products and Partial Traces](@ref doc:Tensor-products-and-Partial-Traces), we showed how Hamiltonians and state vectors are constructed in `QuantumToolbox.jl`. Given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state vector ``|\psi(0)\rangle`` using the `QuantumToolbox` time evolution problem [`sesolveProblem`](@ref) or directly call the function [`sesolve`](@ref). It evolves the state vector ``|\psi(t)\rangle`` and evaluates the expectation values for a set of operators `e_ops` at each given time points, using an ordinary differential equation solver provided by the powerful julia package[`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/).
+The Schrödinger equation, which governs the time-evolution of closed quantum systems, is defined by its Hamiltonian and state vector. In the previous sections, [Manipulating States and Operators](@ref doc:Manipulating-States-and-Operators) and [Tensor Products and Partial Traces](@ref doc:Tensor-products-and-Partial-Traces), we showed how Hamiltonians and state vectors are constructed in `QuantumToolbox.jl`. Given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state vector ``|\psi(0)\rangle`` using the `QuantumToolbox` time evolution problem [`sesolveProblem`](@ref) or directly call the function [`sesolve`](@ref). It evolves the state vector ``|\psi(t)\rangle`` and evaluates the expectation values for a set of operators `e_ops` at each given time points, using an ordinary differential equation solver provided by the powerful julia package [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/).
 
 ### Example
 
@@ -41,7 +41,7 @@ sol = sesolve(prob)
 !!! note "Note"
     Here, we generate the time evolution problem by [`sesolveProblem`](@ref) first and then put it into the function [`sesolve`](@ref). One can also directly call [`sesolve`](@ref), which we also demonstrates in below.
 
-The function returns an instance of [`TimeEvolutionSol`](@ref), as described in the previous section [Time Evolution Solutions](@ref doc-TE:Time-Evolution-Solutions). The attribute `expect` in `sol`ution is a list of expectation values for the operator(s) that are passed to the `e_ops` keyword argument. 
+The function returns [`TimeEvolutionSol`](@ref), as described in the previous section [Time Evolution Solutions](@ref doc-TE:Time-Evolution-Solutions). The attribute `expect` in `sol`ution is a list of expectation values for the operator(s) that are passed to the `e_ops` keyword argument. 
 
 ```@example sesolve
 sol.expect
@@ -95,7 +95,7 @@ sol = sesolve(H, ψ0, tlist) # or specify: e_ops = []
 sol.states
 ```
 
-This is because the `states` will be saved depend on the keyword argument `saveat` in `kwargs`. If `e_ops` is specified, the default value of `saveat=[tlist[end]]` (only save the final state), otherwise, `saveat=tlist` (saving the states corresponding to `tlist`). 
+This is because the `states` will be saved depend on the keyword argument `saveat` in `kwargs`. If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). 
 
 One can also specify `e_ops` and `saveat` separately:
 

From 21aed4f1107d201a65be2e549ac242eeb50a39bd Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Tue, 8 Oct 2024 17:03:46 +0900
Subject: [PATCH 3/6] add documentation for von Neumann eq.

---
 .../src/users_guide/time_evolution/mesolve.md | 55 ++++++++++++++++++-
 1 file changed, 54 insertions(+), 1 deletion(-)

diff --git a/docs/src/users_guide/time_evolution/mesolve.md b/docs/src/users_guide/time_evolution/mesolve.md
index b63c85d74..ae95ef3a3 100644
--- a/docs/src/users_guide/time_evolution/mesolve.md
+++ b/docs/src/users_guide/time_evolution/mesolve.md
@@ -2,8 +2,61 @@
 
 ## [Lindblad Master Equation Solver](@id doc-TE:Lindblad-Master-Equation-Solver)
 
-This page is still under construction, please visit [API](@ref doc-API) first.
+```@setup mesolve
+using QuantumToolbox
+```
 
 ### Von Neumann equation
 
+While the evolution of the state vector in a closed quantum system is deterministic (as we discussed in the previous section: [Schrödinger Equation Solver](@ref doc-TE:Schrödinger-Equation-Solver)), open quantum systems are stochastic in nature. The effect of an environment on the system of interest is to induce stochastic transitions between energy levels, and to introduce uncertainty in the phase difference between states of the system. The state of an open quantum system is therefore described in terms of ensemble averaged states using the density matrix formalism. A density matrix ``\hat{\rho}`` describes a probability distribution of quantum states ``|\psi_n\rangle`` in a matrix representation, namely
+
+```math
+\hat{\rho} = \sum_n p_n |\psi_n\rangle\langle\psi_n|,
+```
+
+where ``p_n`` is the classical probability that the system is in the quantum state ``|\psi_n\rangle``. The time evolution of a density matrix ``\hat{\rho}`` is the topic of the remaining portions of this section.
+
+The time evolution of the density matrix ``\hat{\rho}(t)`` under closed system dynamics is governed by the von Neumann equation: 
+
+```math
+\frac{d\hat{\rho}(t)}{dt} = -\frac{i}{\hbar}\left[\hat{H}, \hat{\rho}(t)\right],
+```
+
+where ``[\cdot, \cdot]`` represents the commutator. In `QuantumToolbox`, given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state using the `QuantumToolbox` time evolution problem [`mesolveProblem`](@ref) or directly call the function [`mesolve`](@ref). It evolves the density matrix ``\hat{\rho}(t)`` and evaluates the expectation values for a set of operators `e_ops` at each given time points, using an ordinary differential equation solver provided by the powerful julia package [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/).
+
+```@example mesolve
+H = 0.5 * sigmax()
+state0 = basis(2, 0)         # state vector
+state0 = ket2dm(basis(2, 0)) # density matrix
+tlist = LinRange(0.0, 10.0, 20)
+
+sol = mesolve(H, state0, tlist, e_ops = [sigmaz()])
+```
+
+!!! note "Type of initial state"
+    The initial state `state0` here can be given as a state vector ``|\psi(0)\rangle`` (in the type of [`Ket`](@ref)) or a density matrix ``\hat{\rho}(0)`` (in the type of [`Operator`](@ref)). If it is given as a [`Ket`](@ref), it will be transformed to density matrix ``\hat{\rho}(0) = |\psi(0)\rangle\langle\psi(0)|`` internally in [`mesolve`](@ref).
+
+The function returns [`TimeEvolutionSol`](@ref), as described in the previous section [Time Evolution Solutions](@ref doc-TE:Time-Evolution-Solutions). The stored `states` will always be in the type of [`Operator`](@ref) (density matrix).
+
+```@example mesolve
+sol.states
+```
+
+Here, only the final state is stored because the `states` will be saved depend on the keyword argument `saveat` in `kwargs`. If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). 
+
+One can also specify `e_ops` and `saveat` separately:
+
+```@example mesolve
+tlist = [0, 5, 10]
+sol = mesolve(H, state0, tlist, e_ops = [sigmay()], saveat = tlist)
+```
+
+```@example mesolve
+sol.expect
+```
+
+```@example mesolve
+sol.states
+```
+
 ### The Lindblad master equation
\ No newline at end of file

From 98e11d6b307121a7731171eb78dcc4a388d5a4c8 Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Tue, 8 Oct 2024 17:19:29 +0900
Subject: [PATCH 4/6] remove title from each time evolution page to fix typo in
 sidebar

---
 docs/src/users_guide/time_evolution/intro.md          |  4 ++--
 docs/src/users_guide/time_evolution/mcsolve.md        |  4 +---
 docs/src/users_guide/time_evolution/mesolve.md        |  8 +++-----
 docs/src/users_guide/time_evolution/sesolve.md        |  8 +++-----
 docs/src/users_guide/time_evolution/solution.md       | 10 ++++------
 docs/src/users_guide/time_evolution/stochastic.md     |  6 ++----
 docs/src/users_guide/time_evolution/time_dependent.md |  4 +---
 7 files changed, 16 insertions(+), 28 deletions(-)

diff --git a/docs/src/users_guide/time_evolution/intro.md b/docs/src/users_guide/time_evolution/intro.md
index 6a7f2f5ab..37488c956 100644
--- a/docs/src/users_guide/time_evolution/intro.md
+++ b/docs/src/users_guide/time_evolution/intro.md
@@ -12,10 +12,10 @@ Pages = [
     "stochastic.md",
     "time_dependent.md",
 ]
-Depth = 2:3
+Depth = 1:2
 ```
 
-## [Introduction](@id doc-TE:Introduction)
+# [Introduction](@id doc-TE:Introduction)
 
 Although in some cases, we want to find the stationary states of a quantum system, often we are interested in the dynamics: how the state of a system or an ensemble of systems evolves with time. `QuantumToolbox` provides many ways to model dynamics.
 
diff --git a/docs/src/users_guide/time_evolution/mcsolve.md b/docs/src/users_guide/time_evolution/mcsolve.md
index b6504e394..b0cfa5591 100644
--- a/docs/src/users_guide/time_evolution/mcsolve.md
+++ b/docs/src/users_guide/time_evolution/mcsolve.md
@@ -1,5 +1,3 @@
-# Time Evolution and Quantum System Dynamics
-
-## [Monte-Carlo Solver](@id doc-TE:Monte-Carlo-Solver)
+# [Monte-Carlo Solver](@id doc-TE:Monte-Carlo-Solver)
 
 This page is still under construction, please visit [API](@ref doc-API) first.
\ No newline at end of file
diff --git a/docs/src/users_guide/time_evolution/mesolve.md b/docs/src/users_guide/time_evolution/mesolve.md
index ae95ef3a3..d408870a8 100644
--- a/docs/src/users_guide/time_evolution/mesolve.md
+++ b/docs/src/users_guide/time_evolution/mesolve.md
@@ -1,12 +1,10 @@
-# Time Evolution and Quantum System Dynamics
-
-## [Lindblad Master Equation Solver](@id doc-TE:Lindblad-Master-Equation-Solver)
+# [Lindblad Master Equation Solver](@id doc-TE:Lindblad-Master-Equation-Solver)
 
 ```@setup mesolve
 using QuantumToolbox
 ```
 
-### Von Neumann equation
+## Von Neumann equation
 
 While the evolution of the state vector in a closed quantum system is deterministic (as we discussed in the previous section: [Schrödinger Equation Solver](@ref doc-TE:Schrödinger-Equation-Solver)), open quantum systems are stochastic in nature. The effect of an environment on the system of interest is to induce stochastic transitions between energy levels, and to introduce uncertainty in the phase difference between states of the system. The state of an open quantum system is therefore described in terms of ensemble averaged states using the density matrix formalism. A density matrix ``\hat{\rho}`` describes a probability distribution of quantum states ``|\psi_n\rangle`` in a matrix representation, namely
 
@@ -59,4 +57,4 @@ sol.expect
 sol.states
 ```
 
-### The Lindblad master equation
\ No newline at end of file
+## The Lindblad master equation
\ No newline at end of file
diff --git a/docs/src/users_guide/time_evolution/sesolve.md b/docs/src/users_guide/time_evolution/sesolve.md
index febe15d4c..f58d0ff09 100644
--- a/docs/src/users_guide/time_evolution/sesolve.md
+++ b/docs/src/users_guide/time_evolution/sesolve.md
@@ -1,8 +1,6 @@
-# Time Evolution and Quantum System Dynamics
+# [Schrödinger Equation Solver](@id doc-TE:Schrödinger-Equation-Solver)
 
-## [Schrödinger Equation Solver](@id doc-TE:Schrödinger-Equation-Solver)
-
-### Unitary evolution
+## Unitary evolution
 
 The dynamics of a closed (pure) quantum system is governed by the Schrödinger equation
 
@@ -21,7 +19,7 @@ where ``|\psi(t)\rangle`` is the state vector, and the Hamiltonian ``\hat{H}`` i
 
 The Schrödinger equation, which governs the time-evolution of closed quantum systems, is defined by its Hamiltonian and state vector. In the previous sections, [Manipulating States and Operators](@ref doc:Manipulating-States-and-Operators) and [Tensor Products and Partial Traces](@ref doc:Tensor-products-and-Partial-Traces), we showed how Hamiltonians and state vectors are constructed in `QuantumToolbox.jl`. Given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state vector ``|\psi(0)\rangle`` using the `QuantumToolbox` time evolution problem [`sesolveProblem`](@ref) or directly call the function [`sesolve`](@ref). It evolves the state vector ``|\psi(t)\rangle`` and evaluates the expectation values for a set of operators `e_ops` at each given time points, using an ordinary differential equation solver provided by the powerful julia package [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/).
 
-### Example
+## Example
 
 ```@setup sesolve
 using QuantumToolbox
diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md
index efaea742f..efd31d0d2 100644
--- a/docs/src/users_guide/time_evolution/solution.md
+++ b/docs/src/users_guide/time_evolution/solution.md
@@ -1,12 +1,10 @@
-# Time Evolution and Quantum System Dynamics
-
-## [Time Evolution Solutions](@id doc-TE:Time-Evolution-Solutions)
+# [Time Evolution Solutions](@id doc-TE:Time-Evolution-Solutions)
 
 ```@setup TE-solution
 using QuantumToolbox
 ```
 
-### Solution
+## Solution
 `QuantumToolbox` utilizes the powerful [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/) to simulate different kinds of quantum system dynamics. Thus, we will first look at the data structure used for returning the solution (`sol`) from [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/). The solution stores all the crucial data needed for analyzing and plotting the results of a simulation. A generic structure [`TimeEvolutionSol`](@ref) contains the following properties for storing simulation data:
 
 | **Fields (Attributes)** | **Description** |
@@ -19,7 +17,7 @@ using QuantumToolbox
 | `sol.reltol` | The relative tolerance which is used during the solving process. |
 | `sol.retcode` (or `sol.converged`) | The returned status from the solver. |
 
-### Accessing data in solutions
+## Accessing data in solutions
 
 To understand how to access the data in solution, we will use an example as a guide, although we do not worry about the simulation details at this stage. The Schrödinger equation solver ([`sesolve`](@ref)) used in this example returns [`TimeEvolutionSol`](@ref):
 
@@ -88,6 +86,6 @@ Here, the solution contains only one (final) state. Because the `states` will be
 
 Some other solvers can have other output.
 
-### Multiple trajectories solution
+## Multiple trajectories solution
 
 This part is still under construction, please visit [API](@ref doc-API) first.
\ No newline at end of file
diff --git a/docs/src/users_guide/time_evolution/stochastic.md b/docs/src/users_guide/time_evolution/stochastic.md
index 49259e78b..2006195ab 100644
--- a/docs/src/users_guide/time_evolution/stochastic.md
+++ b/docs/src/users_guide/time_evolution/stochastic.md
@@ -1,7 +1,5 @@
-# Time Evolution and Quantum System Dynamics
-
-## [Stochastic Solver](@id doc-TE:Stochastic-Solver)
+# [Stochastic Solver](@id doc-TE:Stochastic-Solver)
 
 This page is still under construction, please visit [API](@ref doc-API) first.
 
-### Stochastic Schrodinger equation
\ No newline at end of file
+## Stochastic Schrodinger equation
\ No newline at end of file
diff --git a/docs/src/users_guide/time_evolution/time_dependent.md b/docs/src/users_guide/time_evolution/time_dependent.md
index 37274c343..f9346e884 100644
--- a/docs/src/users_guide/time_evolution/time_dependent.md
+++ b/docs/src/users_guide/time_evolution/time_dependent.md
@@ -1,5 +1,3 @@
-# Time Evolution and Quantum System Dynamics
-
-## [Solving Problems with Time-dependent Hamiltonians](@id doc-TE:Solving-Problems-with-Time-dependent-Hamiltonians)
+# [Solving Problems with Time-dependent Hamiltonians](@id doc-TE:Solving-Problems-with-Time-dependent-Hamiltonians)
 
 This page is still under construction, please visit [API](@ref doc-API) first.
\ No newline at end of file

From b034fd0a60bc12ae42a8c4343dbfbd71e40a272c Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Wed, 9 Oct 2024 13:37:26 +0900
Subject: [PATCH 5/6] minor changes

---
 docs/src/users_guide/time_evolution/sesolve.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/users_guide/time_evolution/sesolve.md b/docs/src/users_guide/time_evolution/sesolve.md
index f58d0ff09..fce63f7ac 100644
--- a/docs/src/users_guide/time_evolution/sesolve.md
+++ b/docs/src/users_guide/time_evolution/sesolve.md
@@ -19,7 +19,7 @@ where ``|\psi(t)\rangle`` is the state vector, and the Hamiltonian ``\hat{H}`` i
 
 The Schrödinger equation, which governs the time-evolution of closed quantum systems, is defined by its Hamiltonian and state vector. In the previous sections, [Manipulating States and Operators](@ref doc:Manipulating-States-and-Operators) and [Tensor Products and Partial Traces](@ref doc:Tensor-products-and-Partial-Traces), we showed how Hamiltonians and state vectors are constructed in `QuantumToolbox.jl`. Given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state vector ``|\psi(0)\rangle`` using the `QuantumToolbox` time evolution problem [`sesolveProblem`](@ref) or directly call the function [`sesolve`](@ref). It evolves the state vector ``|\psi(t)\rangle`` and evaluates the expectation values for a set of operators `e_ops` at each given time points, using an ordinary differential equation solver provided by the powerful julia package [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/).
 
-## Example
+## Example: Spin dynamics
 
 ```@setup sesolve
 using QuantumToolbox

From 88fd4690cf9fd5302ced591d31b68078b4956fde Mon Sep 17 00:00:00 2001
From: Yi-Te Huang <y.t.d.huang@phys.ncku.edu.tw>
Date: Wed, 9 Oct 2024 13:37:45 +0900
Subject: [PATCH 6/6] update documentation for `mesolve`

---
 docs/src/users_guide/states_and_operators.md  |   2 +-
 docs/src/users_guide/steadystate.md           |  14 +-
 .../src/users_guide/time_evolution/mesolve.md | 173 +++++++++++++++++-
 3 files changed, 179 insertions(+), 10 deletions(-)

diff --git a/docs/src/users_guide/states_and_operators.md b/docs/src/users_guide/states_and_operators.md
index b836c3489..3cd714820 100644
--- a/docs/src/users_guide/states_and_operators.md
+++ b/docs/src/users_guide/states_and_operators.md
@@ -404,4 +404,4 @@ t = 0.8
 exp(L * t)
 ```
 
-See the section [Time Evolution and Quantum System Dynamics](@ref doc:Time-Evolution-and-Quantum-System-Dynamics) for more details.
+See the section [Lindblad Master Equation Solver](@ref doc-TE:Lindblad-Master-Equation-Solver) for more details.
diff --git a/docs/src/users_guide/steadystate.md b/docs/src/users_guide/steadystate.md
index a69409f5c..2c53efb03 100644
--- a/docs/src/users_guide/steadystate.md
+++ b/docs/src/users_guide/steadystate.md
@@ -64,7 +64,9 @@ solver = SteadyStateLinearSolver(alg = MKLLUFactorization())
 
 See [`LinearSolve.jl` Solvers](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/) for more details.
 
-## Example: Harmonic Oscillator in Thermal Bath
+## Example: Harmonic oscillator in thermal bath
+
+Here, we demonstrate [`steadystate`](@ref) by using the example with the harmonic oscillator in thermal bath from the previous section ([Lindblad Master Equation Solver](@ref doc-TE:Lindblad-Master-Equation-Solver)).
 
 ```@example steady_state_example
 using QuantumToolbox
@@ -77,16 +79,16 @@ a = destroy(N)
 H = a' * a
 ψ0 = basis(N, 10)  # initial state
 κ = 0.1  # coupling to oscillator
-n_th = 2  # temperature with average of 2 excitations
+n_th = 2 # temperature with average of 2 excitations
 
 # collapse operators 
-c_op_list = [
+c_ops = [
     sqrt(κ * (n_th + 1)) * a, # emission
     sqrt(κ *  n_th     ) * a' # absorption
 ]
 
 # find steady-state solution
-ρ_ss = steadystate(H, c_op_list)
+ρ_ss = steadystate(H, c_ops)
 
 # find expectation value for particle number in steady state
 e_ops = [a' * a]
@@ -95,11 +97,11 @@ exp_ss = real(expect(e_ops[1], ρ_ss))
 tlist = LinRange(0, 50, 100)
 
 # monte-carlo
-sol_mc = mcsolve(H, ψ0, tlist, c_op_list, e_ops=e_ops, ntraj=100, progress_bar=false)
+sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=100, progress_bar=false)
 exp_mc = real(sol_mc.expect[1, :])
 
 # master eq.
-sol_me = mesolve(H, ψ0, tlist, c_op_list, e_ops=e_ops, progress_bar=false)
+sol_me = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops, progress_bar=false)
 exp_me = real(sol_me.expect[1, :])
 
 # plot the results
diff --git a/docs/src/users_guide/time_evolution/mesolve.md b/docs/src/users_guide/time_evolution/mesolve.md
index d408870a8..fdda3e9bb 100644
--- a/docs/src/users_guide/time_evolution/mesolve.md
+++ b/docs/src/users_guide/time_evolution/mesolve.md
@@ -17,10 +17,15 @@ where ``p_n`` is the classical probability that the system is in the quantum sta
 The time evolution of the density matrix ``\hat{\rho}(t)`` under closed system dynamics is governed by the von Neumann equation: 
 
 ```math
-\frac{d\hat{\rho}(t)}{dt} = -\frac{i}{\hbar}\left[\hat{H}, \hat{\rho}(t)\right],
+\begin{equation}
+\label{von-Neumann-Eq}
+\frac{d}{dt}\hat{\rho}(t) = -\frac{i}{\hbar}\left[\hat{H}, \hat{\rho}(t)\right],
+\end{equation}
 ```
 
-where ``[\cdot, \cdot]`` represents the commutator. In `QuantumToolbox`, given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state using the `QuantumToolbox` time evolution problem [`mesolveProblem`](@ref) or directly call the function [`mesolve`](@ref). It evolves the density matrix ``\hat{\rho}(t)`` and evaluates the expectation values for a set of operators `e_ops` at each given time points, using an ordinary differential equation solver provided by the powerful julia package [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/).
+where ``[\cdot, \cdot]`` represents the commutator. The above equation is equivalent to the Schrödinger equation described in the [previous section](@ref doc-TE:Schrödinger-Equation-Solver) under the density matrix formalism.
+
+In `QuantumToolbox`, given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state using the `QuantumToolbox` time evolution problem [`mesolveProblem`](@ref) or directly call the function [`mesolve`](@ref). It evolves the density matrix ``\hat{\rho}(t)`` and evaluates the expectation values for a set of operators `e_ops` at each given time points, using an ordinary differential equation solver provided by the powerful julia package [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/).
 
 ```@example mesolve
 H = 0.5 * sigmax()
@@ -57,4 +62,166 @@ sol.expect
 sol.states
 ```
 
-## The Lindblad master equation
\ No newline at end of file
+## The Lindblad master equation
+
+The standard approach for deriving the equations of motion for a system interacting with its environment is to expand the scope of the system to include the environment. The combined quantum system is then closed, and its evolution is governed by the von Neumann equation given in Eq. \eqref{von-Neumann-Eq}
+
+```math
+\begin{equation}
+\label{tot-von-Neumann-Eq}
+\frac{d}{dt}\hat{\rho}_{\textrm{tot}}(t) = -\frac{i}{\hbar}\left[\hat{H}_{\textrm{tot}}, \hat{\rho}_{\textrm{tot}}(t)\right].
+\end{equation}
+```
+
+Here, the total Hamiltonian
+
+```math
+\hat{H}_{\textrm{tot}} = \hat{H}_{\textrm{sys}} + \hat{H}_{\textrm{env}} + \hat{H}_{\textrm{int}},
+```
+
+includes the original system Hamiltonian ``\hat{H}_{\textrm{sys}}``, the Hamiltonian for the environment ``\hat{H}_{\textrm{env}}``, and a term representing the interaction between the system and its environment ``\hat{H}_{\textrm{int}}``. Since we are only interested in the dynamics of the system, we can, perform a partial trace over the environmental degrees of freedom in Eq. \eqref{tot-von-Neumann-Eq}, and thereby obtain a master equation for the motion of the original system density matrix ``\hat{\rho}_{\textrm{sys}}(t)=\textrm{Tr}_{\textrm{env}}[\hat{\rho}_{\textrm{tot}}(t)]``. The most general trace-preserving and completely positive form of this evolution is the Lindblad master equation for the reduced density matrix, namely
+
+```math
+\begin{equation}
+\label{Lindblad-master-Eq}
+\frac{d}{dt}\hat{\rho}_{\textrm{sys}}(t) = -\frac{i}{\hbar}\left[\hat{H}_{\textrm{sys}}, \hat{\rho}_{\textrm{sys}}(t)\right] + \sum_n \hat{C}_n \hat{\rho}_{\textrm{sys}}(t) \hat{C}_n^\dagger - \frac{1}{2} \hat{C}_n^\dagger \hat{C}_n \hat{\rho}_{\textrm{sys}}(t) - \frac{1}{2} \hat{\rho}_{\textrm{sys}}(t) \hat{C}_n^\dagger \hat{C}_n
+\end{equation}
+```
+
+where ``\hat{C}_n \equiv \sqrt{\gamma_n}\hat{A}_n`` are the collapse operators, ``\hat{A}_n`` are the operators acting on the system in ``\hat{H}_{\textrm{int}}`` which describes the system-environment interaction, and ``\gamma_n`` are the corresponding rates. The derivation of Eq. \eqref{Lindblad-master-Eq} may be found in several sources, and will not be reproduced here. Instead, we emphasize the approximations that are required to arrive at the master equation in the form of Eq. \eqref{Lindblad-master-Eq} from physical arguments, and hence perform a calculation in `QuantumToolbox`:
+
+- **Separability:** At ``t = 0``, there are no correlations between the system and environment, such that the total density matrix can be written as a tensor product, namely ``\hat{\rho}_{\textrm{tot}}(0)=\hat{\rho}_{\textrm{sys}}(0)\otimes\hat{\rho}_{\textrm{env}}(0)``.
+- **Born approximation:** Requires: (i) the state of the environment does not significantly change as a result of the interaction with the system; (ii) the system and the environment remain separable throughout the evolution. These assumptions are justified if the interaction is weak, and if the environment is much larger than the system. In summary, ``\hat{\rho}_{\textrm{tot}}(t)\approx\hat{\rho}_{\textrm{sys}}(t)\otimes\hat{\rho}_{\textrm{env}}(0)``.
+- **Markov approximation:** The time-scale of decay for the environment ``\tau_{\textrm{env}}`` is much shorter than the smallest time-scale of the system dynamics, i.e., ``\tau_{\textrm{sys}}\gg\tau_{\textrm{env}}``. This approximation is often deemed a “short-memory environment” as it requires the environmental correlation functions decay in a fast time-scale compared to those of the system.
+- **Secular approximation:** Stipulates that elements in the master equation corresponding to transition frequencies satisfy ``|\omega_{ab}-\omega_{cd}| \ll 1/\tau_{\textrm{sys}}``, i.e., all fast rotating terms in the interaction picture can be neglected. It also ignores terms that lead to a small renormalization of the system energy levels. This approximation is not strictly necessary for all master-equation formalisms (e.g., the Block-Redfield master equation), but it is required for arriving at the Lindblad form in Eq. \eqref{Lindblad-master-Eq} which is used in [`mesolve`](@ref).
+
+For systems with environments satisfying the conditions outlined above, the Lindblad master equation in Eq. \eqref{Lindblad-master-Eq} governs the time-evolution of the system density matrix, giving an ensemble average of the system dynamics. In order to ensure that these approximations are not violated, it is important that the decay rates ``\gamma_n`` be smaller than the minimum energy splitting in the system Hamiltonian. Situations that demand special attention therefore include, for example, systems strongly coupled to their environment, and systems with degenerate or nearly degenerate energy levels.
+
+What is new in the master equation compared to the Schrödinger equation (or von Neumann equation) are processes that describe dissipation in the quantum system due to its interaction with an environment. For example, evolution that includes incoherent processes such as relaxation and dephasing. These environmental interactions are defined by the operators ``\hat{A}_n`` through which the system couples to the environment, and rates ``\gamma_n`` that describe the strength of the processes.
+
+In `QuantumToolbox`, the function [`mesolve`](@ref) can also be used for solving the master equation. This is done by passing a list of collapse operators (`c_ops`) as the fourth argument of the [`mesolve`](@ref) function in order to define the dissipation processes of the master equation in Eq. \eqref{Lindblad-master-Eq}. As we mentioned above, each collapse operator ``\hat{C}_n`` is the product of ``\sqrt{\gamma_n}`` (the square root of the rate) and ``\hat{A}_n`` (operator which describes the dissipation process). 
+
+Furthermore, `QuantumToolbox` solves the master equation in the [`SuperOperator`](@ref) formalism. That is, a Liouvillian [`SuperOperator`](@ref) will be generated internally in [`mesolve`](@ref) by the input system Hamiltonian ``\hat{H}_{\textrm{sys}}`` and the collapse operators ``\hat{C}_n``. One can also generate the Liouvillian [`SuperOperator`](@ref) manually for special purposes, and pass it as the first argument of the [`mesolve`](@ref) function. To do so, it is useful to read the section [Superoperators and Vectorized Operators](@ref doc:Superoperators-and-Vectorized-Operators), and also the docstrings of the following functions:
+- [`spre`](@ref)
+- [`spost`](@ref)
+- [`sprepost`](@ref)
+- [`liouvillian`](@ref)
+- [`lindblad_dissipator`](@ref)
+
+## Example: Spin dynamics
+
+Using the example with the dynamics of spin-``\frac{1}{2}`` from the previous section ([Schrödinger Equation Solver](@ref doc-TE:Schrödinger-Equation-Solver)), we can easily add a relaxation process (describing the dissipation of energy from the spin to the environment), by adding `[sqrt(γ) * sigmax()]` in the fourth parameter of the [`mesolve`](@ref) function.
+
+```@example mesolve
+H = 2 * π * 0.1 * sigmax()
+ψ0 = basis(2, 0) # spin-up
+tlist = LinRange(0.0, 10.0, 100)
+
+γ = 0.05
+c_ops = [sqrt(γ) * sigmax()]
+
+sol = mesolve(H, ψ0, tlist, c_ops, e_ops = [sigmaz(), sigmay()])
+```
+
+We can therefore plot the expectation values:
+
+```@example mesolve
+using CairoMakie
+CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
+
+times = sol.times
+expt_z = real(sol.expect[1,:])
+expt_y = real(sol.expect[2,:])
+
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
+lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
+lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
+
+axislegend(ax, position = :rt)
+
+fig
+```
+
+## Example: Harmonic oscillator in thermal bath
+
+Consider a harmonic oscillator (single-mode cavity) couples to a thermal bath. If the single-mode cavity initially is in a `10`-photon [`fock`](@ref) state, the dynamics is calculated with the following code:
+
+```@example mesolve
+# Define parameters
+N = 20  # number of basis states to consider
+a = destroy(N)
+H = a' * a
+ψ0 = fock(N, 10)  # initial state
+κ = 0.1  # coupling to oscillator
+n_th = 2 # temperature with average of 2 excitations
+tlist = LinRange(0, 50, 100)
+
+# collapse operators 
+c_ops = [
+    sqrt(κ * (n_th + 1)) * a, # emission
+    sqrt(κ *  n_th     ) * a' # absorption
+]
+
+# find expectation value for particle number
+e_ops = [a' * a]
+
+sol = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops)
+Num = real(sol.expect[1, :])
+
+# plot the results
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], 
+    title = L"Decay of Fock state $|10\rangle$ in a thermal environment with $\langle n\rangle=2$",
+    xlabel = "Time", 
+    ylabel = "Number of excitations",
+)
+lines!(ax, tlist, Num)
+
+fig
+```
+
+## Example: Two-level atom coupled to dissipative single-mode cavity
+
+Consider a two-level atom coupled to a dissipative single-mode cavity through a dipole-type interaction, which supports a coherent exchange of quanta between the two systems. If the atom initially is in its ground state and the cavity in a `5`-photon [`fock`](@ref) state, the dynamics is calculated with the following code:
+
+!!! note "Generate Liouviilian superoperator manually"
+    In this example, we demonstrate how to generate the Liouvillian [`SuperOperator`](@ref) manually and pass it as the first argument of the [`mesolve`](@ref) function.
+
+```@example mesolve
+# two-level atom
+σm = tensor(destroy(2), qeye(10))
+H_a = 2 * π * σm' * σm
+
+# dissipative single-mode cavity
+γ = 0.1 # dissipation rate
+a  = tensor(qeye(2), destroy(10))
+H_c = 2 * π * a' * a
+c_ops = [sqrt(γ) * a]
+
+# coupling between two-level atom and single-mode cavity
+g = 0.25 # coupling strength
+H_I = 2 * π * g * (σm * a' + σm' * a)
+
+ψ0 = tensor(basis(2,0), fock(10, 5)) # initial state
+tlist = LinRange(0.0, 10.0, 200)
+
+# generate Liouvillian superoperator manually
+L = liouvillian(H_a + H_c + H_I, c_ops)
+sol = mesolve(L, ψ0, tlist, e_ops=[σm' * σm, a' * a])
+
+times = sol.times
+
+# expectation value of Number operator
+N_atom   = real(sol.expect[1,:])
+N_cavity = real(sol.expect[2,:])
+
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
+lines!(ax, times, N_atom, label = "atom excitation probability", linestyle = :solid)
+lines!(ax, times, N_cavity, label = "cavity photon number", linestyle = :dash)
+
+axislegend(ax, position = :rt)
+
+fig
+```