diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md
index 0f1cf32..03491df 100644
--- a/.github/pull_request_template.md
+++ b/.github/pull_request_template.md
@@ -3,6 +3,7 @@ Thank you for contributing to Tutorials for Quantum Toolbox in `Julia`! Please m
 
 - [ ] Please read [Contributing to QuantumToolbox.jl](https://qutip.org/QuantumToolbox.jl/stable/resources/contributing).
 - [ ] The (last update) `date` were modified for new or updated tutorials.
+- [ ] For new tutorials, a `Version Information` section was added at the end and displays the output of `versioninfo()`.
 - [ ] All tutorials were able to render locally by running: `make render`.
 
 Request for a review after you have completed all the tasks. If you have not finished them all, you can also open a [Draft Pull Request](https://github.blog/2019-02-14-introducing-draft-pull-requests/) to let the others know this on-going work.
diff --git a/HierarchicalEOM.jl/SIAM.qmd b/HierarchicalEOM.jl/SIAM.qmd
new file mode 100644
index 0000000..89b079a
--- /dev/null
+++ b/HierarchicalEOM.jl/SIAM.qmd
@@ -0,0 +1,103 @@
+---
+title: "Single-impurity Anderson model"
+author: Yi-Te Huang
+date: 2025-01-14  # last update (keep this comment as a reminder)
+
+engine: julia
+---
+
+## Introduction
+
+The investigation of the Kondo effect in single-impurity Anderson model is crucial as it serves both as a valuable testing ground for the theories of the Kondo effect and has the potential to lead to a better understanding of this intrinsic many-body phenomena. For further detailed discussions of this model (under different parameters) using [`HierarchicalEOM.jl`](https://github.com/qutip/HierarchicalEOM.jl), we recommend to read the article [@HierarchicalEOM-jl2023].
+
+## Hamiltonian
+
+We consider a single-level electronic system [which can be populated by a spin-up ($\uparrow$) or spin-down ($\downarrow$) electron] coupled to a fermionic reservoir ($\textrm{f}$). The total Hamiltonian is given by $H_{\textrm{T}}=H_\textrm{s}+H_\textrm{f}+H_\textrm{sf}$, where each terms takes the form
+
+$$
+\begin{aligned}
+H_{\textrm{s}}  &= \epsilon \left(d^\dagger_\uparrow d_\uparrow + d^\dagger_\downarrow d_\downarrow \right) + U\left(d^\dagger_\uparrow d_\uparrow d^\dagger_\downarrow d_\downarrow\right),\\
+H_{\textrm{f}}  &=\sum_{\sigma=\uparrow,\downarrow}\sum_{k}\epsilon_{\sigma,k}c_{\sigma,k}^{\dagger}c_{\sigma,k},\\
+H_{\textrm{sf}} &=\sum_{\sigma=\uparrow,\downarrow}\sum_{k}g_{k}c_{\sigma,k}^{\dagger}d_{\sigma} + g_{k}^*d_{\sigma}^{\dagger}c_{\sigma,k}.
+\end{aligned}
+$$
+
+Here, $d_\uparrow$ $(d_\downarrow)$ annihilates a spin-up (spin-down) electron in the system, $\epsilon$ is the energy of the electron, and $U$ is the Coulomb repulsion energy for double occupation. Furthermore, $c_{\sigma,k}$ $(c_{\sigma,k}^{\dagger})$ annihilates (creates) an electron in the state $k$ (with energy $\epsilon_{\sigma,k}$) of the reservoir.
+
+Now, we need to build the system Hamiltonian and initial state with the package [`QuantumToolbox.jl`](https://github.com/qutip/QuantumToolbox.jl) to construct the operators.
+
+```{julia}
+using HierarchicalEOM  # this automatically loads `QuantumToolbox`
+using CairoMakie       # for plotting results
+```
+
+```{julia}
+ϵ = -5
+U = 10
+σm = sigmam() # σ-
+σz = sigmaz() # σz
+II = qeye(2)  # identity matrix
+
+# construct the annihilation operator for both spin-up and spin-down
+# (utilize Jordan–Wigner transformation)
+d_up = tensor(σm, II)
+d_dn = tensor(-1 * σz, σm)
+Hsys = ϵ * (d_up' * d_up + d_dn' * d_dn) + U * (d_up' * d_up * d_dn' * d_dn)
+```
+
+## Construct bath objects
+
+We assume the fermionic reservoir to have a [Lorentzian-shaped spectral density](https://qutip.org/HierarchicalEOM.jl/stable/bath_fermion/Fermion_Lorentz/#doc-Fermion-Lorentz), and we utilize the Padé decomposition. Furthermore, the spectral densities depend on the following physical parameters:
+
+- the coupling strength $\Gamma$ between system and reservoirs
+- the band-width $W$
+- the product of the Boltzmann constant $k$ and the absolute temperature $T$ : $kT$
+- the chemical potential $\mu$
+- the total number of exponentials for the reservoir $2(N + 1)$
+
+```{julia}
+Γ = 2
+μ = 0
+W = 10
+kT = 0.5
+N = 5
+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]
+```
+
+## Construct HEOMLS matrix
+
+```{julia}
+tier = 3
+M_even = M_Fermion(Hsys, tier, bath_list)
+M_odd = M_Fermion(Hsys, tier, bath_list, ODD)
+```
+
+## Solve stationary state of ADOs
+
+```{julia}
+ados_s = steadystate(M_even)
+```
+
+## Calculate density of states (DOS)
+
+```{julia}
+ωlist = -10:1:10
+dos = DensityOfStates(M_odd, ados_s, d_up, ωlist)
+```
+
+plot the results
+
+```{julia}
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = L"\omega")
+lines!(ax, ωlist, dos)
+
+fig
+```
+
+## Version Information
+```{julia}
+HierarchicalEOM.versioninfo()
+```
diff --git a/HierarchicalEOM.jl/cavityQED.qmd b/HierarchicalEOM.jl/cavityQED.qmd
index a68634d..3e2a481 100644
--- a/HierarchicalEOM.jl/cavityQED.qmd
+++ b/HierarchicalEOM.jl/cavityQED.qmd
@@ -41,12 +41,12 @@ Here, $H_{\textrm{b}}$ describes a bosonic reservoir where $b_{k}$ $(b_{k}^{\dag
 Now, we need to build the system Hamiltonian and initial state with the package [`QuantumToolbox.jl`](https://github.com/qutip/QuantumToolbox.jl) to construct the operators.
 
 ```{julia}
-using HierarchicalEOM
-using CairoMakie
+using HierarchicalEOM  # this automatically loads `QuantumToolbox`
+using CairoMakie       # for plotting results
 ```
 
 ```{julia}
-N = 3 ## system cavity Hilbert space cutoff
+N = 3 # system cavity Hilbert space cutoff
 ωA = 2
 ωc = 2
 g = 0.1
diff --git a/HierarchicalEOM.jl/dynamical_decoupling.qmd b/HierarchicalEOM.jl/dynamical_decoupling.qmd
new file mode 100644
index 0000000..cbf4953
--- /dev/null
+++ b/HierarchicalEOM.jl/dynamical_decoupling.qmd
@@ -0,0 +1,172 @@
+---
+title: "Driven systems and dynamical decoupling"
+author: Yi-Te Huang
+date: 2025-01-14  # last update (keep this comment as a reminder)
+
+engine: julia
+---
+
+Inspirations taken from an example in QuTiP-BoFiN article [@QuTiP-BoFiN2023].
+
+## Introduction
+
+In this example, we show how to solve the time evolution with time-dependent Hamiltonian problems in hierarchical equations of motion approach.
+
+Here, we study dynamical decoupling which is a common tool used to undo the dephasing effect from the environment even for finite pulse duration.  
+
+## Hamiltonian
+
+We consider a two-level system coupled to a bosonic reservoir ($\textrm{b}$). The total Hamiltonian is given by $H_{\textrm{T}}=H_\textrm{s}+H_\textrm{b}+H_\textrm{sb}$, where each terms takes the form
+
+$$
+\begin{aligned}
+H_{\textrm{s}}(t) &= H_0 + H_{\textrm{D}}(t),\\
+H_0               &= \frac{\omega_0}{2} \sigma_z,\\
+H_{\textrm{b}}    &=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k},\\
+H_{\textrm{sb}}   &=\sigma_z\sum_{k}g_{\alpha,k}(b_k + b_k^{\dagger}).
+\end{aligned}
+$$
+
+Here, $H_{\textrm{b}}$ describes a bosonic reservoir where $b_{k}$ $(b_{k}^{\dagger})$ is the bosonic annihilation (creation) operator associated to the $k$th mode (with frequency $\omega_{k}$).
+
+Furthermore, to observe the time evolution of the coherence, we consider the initial state to be
+$$
+ψ(t=0)=\frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right)
+$$
+
+Now, we need to build the system Hamiltonian and initial state with the package [`QuantumToolbox.jl`](https://github.com/qutip/QuantumToolbox.jl) to construct the operators.
+
+```{julia}
+using HierarchicalEOM  # this automatically loads `QuantumToolbox`
+using CairoMakie       # for plotting results
+```
+
+```{julia}
+ω0 = 0.0
+σz = sigmaz()
+σx = sigmax()
+H0 = 0.5 * ω0 * σz
+
+# Define the operator that measures the 0, 1 element of density matrix
+ρ01 = Qobj([0 1; 0 0])
+
+ψ0 = (basis(2, 0) + basis(2, 1)) / √2
+```
+
+The time-dependent driving term $H_{\textrm{D}}(t)$ has the form
+
+$$
+H_{\textrm{D}}(t) = \sum_{n=1}^N f_n(t) \sigma_x 
+$$
+
+where the pulse is chosen to have duration $\tau$ together with a delay $\Delta$ between each pulses, namely
+
+$$
+f_n(t)
+= \begin{cases}
+  V & \textrm{if}~~(n-1)\tau + n\Delta \leq t \leq n (\tau + \Delta),\\
+  0 & \textrm{otherwise}.
+  \end{cases} 
+$$
+
+Here, we set the period of the pulses to be $\tau V = \pi/2$. Therefore, we consider two scenarios with fast and slow pulses:
+
+```{julia}
+# a function which returns the amplitude of the pulse at time t
+function pulse(p, t)
+    τ = 0.5 * π / p.V
+    period = τ + p.Δ
+
+    if (t % period) < τ
+        return p.V
+    else
+        return 0
+    end
+end
+
+tlist = 0:0.4:400
+amp_fast = 0.50
+amp_slow = 0.01
+delay = 20
+
+fastTuple = (V = amp_fast, Δ = delay)
+slowTuple = (V = amp_slow, Δ = delay)
+
+# plot
+fig = Figure(size = (600, 350))
+ax = Axis(fig[1, 1], xlabel = L"t")
+lines!(ax, tlist, [pulse(fastTuple, t) for t in tlist], label = "Fast Pulse", linestyle = :solid)
+lines!(ax, tlist, [pulse(slowTuple, t) for t in tlist], label = "Slow Pulse", linestyle = :dash)
+
+axislegend(ax, position = :rt)
+
+fig
+```
+
+## Construct bath objects
+
+We assume the bosonic reservoir to have a [Drude-Lorentz Spectral Density](https://qutip.org/HierarchicalEOM.jl/stable/bath_boson/Boson_Drude_Lorentz/#Boson-Drude-Lorentz), and we utilize the Padé decomposition. Furthermore, the spectral densities depend on the following physical parameters: 
+
+- the coupling strength $\Gamma$ between system and reservoir
+- the band-width $W$
+- the product of the Boltzmann constant $k$ and the absolute temperature $T$ : $kT$
+- the total number of exponentials for the reservoir $(N + 1)$
+
+```{julia}
+Γ = 0.0005
+W = 0.005
+kT = 0.05
+N = 3
+bath = Boson_DrudeLorentz_Pade(σz, Γ, W, kT, N)
+```
+
+## Construct HEOMLS matrix
+
+::: {.callout-warning}
+Only provide the **time-independent** part of system Hamiltonian when constructing HEOMLS matrices (the time-dependent part `H_t` should be given when solving the time evolution).
+:::
+
+```{julia}
+tier = 6
+M = M_Boson(H0, tier, bath)
+```
+
+## time evolution with time-independent Hamiltonian
+```{julia}
+noPulseSol = HEOMsolve(M, ψ0, tlist; e_ops = [ρ01]);
+```
+
+## Solve time evolution with time-dependent Hamiltonian
+
+We need to provide a `QuantumToolbox.QuantumObjectEvolution` (named as `H_D` in this case)
+
+```{julia}
+H_D = QobjEvo(σx, pulse)
+```
+
+The keyword argument `params` in `HEOMsolve` will be passed to the argument `p` in user-defined function (`pulse` in this case) directly:
+
+```{julia}
+fastPulseSol = HEOMsolve(M, ψ0, tlist; e_ops = [ρ01], H_t = H_D, params = fastTuple)
+slowPulseSol = HEOMsolve(M, ψ0, tlist; e_ops = [ρ01], H_t = H_D, params = slowTuple)
+```
+
+## Plot the coherence
+
+```{julia}
+fig = Figure(size = (600, 350))
+ax = Axis(fig[1, 1], xlabel = L"t", ylabel = L"\rho_{01}")
+lines!(ax, tlist, real(fastPulseSol.expect[1, :]), label = "Fast Pulse", linestyle = :solid)
+lines!(ax, tlist, real(slowPulseSol.expect[1, :]), label = "Slow Pulse", linestyle = :dot)
+lines!(ax, tlist, real(noPulseSol.expect[1, :]),   label = "no Pulse",   linestyle = :dash)
+
+axislegend(ax, position = :lb)
+
+fig
+```
+
+## Version Information
+```{julia}
+HierarchicalEOM.versioninfo()
+```
+
diff --git a/HierarchicalEOM.jl/electronic_current.qmd b/HierarchicalEOM.jl/electronic_current.qmd
new file mode 100644
index 0000000..053868a
--- /dev/null
+++ b/HierarchicalEOM.jl/electronic_current.qmd
@@ -0,0 +1,169 @@
+---
+title: "Electronic Current"
+author: Yi-Te Huang
+date: 2025-01-14  # last update (keep this comment as a reminder)
+
+engine: julia
+---
+
+Inspirations taken from [qutip documentation](https://qutip.org/documentation.html).
+
+## Introduction
+
+In this example, we demonstrate how to compute an environmental observable: the electronic current. For further detailed discussions of calculating the electronic current using [`HierarchicalEOM.jl`](https://github.com/qutip/HierarchicalEOM.jl), we recommend to read the article [@HierarchicalEOM-jl2023].
+
+## Hamiltonian
+We consider a single-level charge system coupled to two [left (L) and right (R)] fermionic reservoirs ($\textrm{f}$). The total Hamiltonian is given by $H_{\textrm{T}}=H_\textrm{s}+H_\textrm{f}+H_\textrm{sf}$, where each terms takes the form
+
+$$
+\begin{aligned}
+H_{\textrm{s}}  &= \epsilon d^\dagger d,\\
+H_{\textrm{f}}  &=\sum_{\alpha=\textrm{L},\textrm{R}}\sum_{k}\epsilon_{\alpha,k}c_{\alpha,k}^{\dagger}c_{\alpha,k},\\
+H_{\textrm{sf}} &=\sum_{\alpha=\textrm{L},\textrm{R}}\sum_{k}g_{\alpha,k}c_{\alpha,k}^{\dagger}d + g_{\alpha,k}^* d^{\dagger}c_{\alpha,k}.
+\end{aligned}
+$$
+
+Here, $d$ $(d^\dagger)$ annihilates (creates) an electron in the system and $\epsilon$ is the energy of the electron. Furthermore, $c_{\alpha,k}$ $(c_{\alpha,k}^{\dagger})$ annihilates (creates) an electron in the state $k$ (with energy $\epsilon_{\alpha,k}$) of the $\alpha$-th reservoir.
+
+Now, we need to build the system Hamiltonian and initial state with the package [`QuantumToolbox.jl`](https://github.com/qutip/QuantumToolbox.jl) to construct the operators.
+
+```{julia}
+using HierarchicalEOM  # this automatically loads `QuantumToolbox`
+using CairoMakie       # for plotting results
+```
+
+```{julia}
+d = destroy(2) # annihilation operator of the system electron
+
+# The system Hamiltonian
+ϵ = 1.0 # site energy
+Hsys = ϵ * d' * d
+
+# System initial state
+ψ0 = basis(2, 0);
+```
+
+## Construct bath objects
+
+We assume the fermionic reservoir to have a [Lorentzian-shaped spectral density](https://qutip.org/HierarchicalEOM.jl/stable/bath_fermion/Fermion_Lorentz/#doc-Fermion-Lorentz), and we utilize the Padé decomposition. Furthermore, the spectral densities depend on the following physical parameters: 
+
+- the coupling strength $\Gamma$ between system and reservoirs
+- the band-width $W$
+- the product of the Boltzmann constant $k$ and the absolute temperature $T$ : $kT$
+- the chemical potential $\mu$
+- the total number of exponentials for the reservoir $2(N + 1)$
+
+```{julia}
+Γ = 0.01
+W = 1
+kT = 0.025851991
+μL =  1.0 # Left  bath
+μR = -1.0 # Right bath
+N = 2
+bath_L = Fermion_Lorentz_Pade(d, Γ, μL, W, kT, N)
+bath_R = Fermion_Lorentz_Pade(d, Γ, μR, W, kT, N)
+baths = [bath_L, bath_R]
+```
+
+## Construct HEOMLS matrix
+
+```{julia}
+tier = 5
+M = M_Fermion(Hsys, tier, baths)
+```
+
+## Solve time evolution of ADOs
+
+```{julia}
+tlist = 0:0.5:100
+ados_evolution = HEOMsolve(M, ψ0, tlist).ados
+```
+
+## Solve stationary state of ADOs
+```{julia}
+ados_steady = steadystate(M)
+```
+
+## Calculate current
+
+Within the influence functional approach, the expectation value of the electronic current from the $\alpha$-fermionic bath into the system can be written in terms of the first-level-fermionic ($n=1$) auxiliary density operators, namely
+
+$$
+\langle I_\alpha(t) \rangle =(-e) \frac{d\langle \mathcal{N}_\alpha\rangle}{dt}=i e \sum_{q\in\textbf{q}}(-1)^{\delta_{\nu,-}} ~\textrm{Tr}\left[d^{\bar{\nu}}\rho^{(0,1,+)}_{\vert \textbf{q}}(t)\right],
+$$
+
+where $e$ represents the value of the elementary charge, and $\mathcal{N}_\alpha=\sum_k c^\dagger_{\alpha,k}c_{\alpha,k}$ is the occupation number operator for the $\alpha$-fermionic bath.
+
+Given an ADOs, we provide a function which calculates the current from the $\alpha$-fermionic bath into the system with the help of [Hierarchy Dictionary](https://qutip.org/HierarchicalEOM.jl/stable/hierarchy_dictionary/#doc-Hierarchy-Dictionary).
+
+```{julia}
+# `bathIdx`:  
+# - 1 means 1st bath (bath_L)
+# - 2 means 2nd bath (bath_R)
+function Ic(ados, M::M_Fermion, bathIdx::Int)
+    # the hierarchy dictionary
+    HDict = M.hierarchy
+
+    # we need all the indices of ADOs for the first level
+    idx_list = HDict.lvl2idx[1]
+    I = 0.0im
+    for idx in idx_list
+        ρ1 = ados[idx]  # 1st-level ADO
+
+        # find the corresponding bath index (α) and exponent term index (k)
+        nvec = HDict.idx2nvec[idx]
+        for (α, k, _) in getIndexEnsemble(nvec, HDict.bathPtr)
+            if α == bathIdx
+                exponent = M.bath[α][k]
+                if exponent.types == "fA"     # fermion-absorption
+                    I += tr(exponent.op' * ρ1)
+                elseif exponent.types == "fE" # fermion-emission
+                    I -= tr(exponent.op' * ρ1)
+                end
+                break
+            end
+        end
+    end
+
+    eV_to_Joule = 1.60218E-19  # unit conversion
+
+    # (e / ħ) * I  [change unit to μA] 
+    return 1.519270639695384E15 * real(1im * I) * eV_to_Joule * 1E6
+end
+```
+
+## steady current
+```{julia}
+Is_L = ones(length(tlist)) .* Ic(ados_steady, M, 1)
+Is_R = ones(length(tlist)) .* Ic(ados_steady, M, 2)
+```
+
+## time evolution current
+```{julia}
+Ie_L = []
+Ie_R = []
+for ados in ados_evolution
+    push!(Ie_L, Ic(ados, M, 1))
+    push!(Ie_R, Ic(ados, M, 2))
+end
+```
+
+plot the result
+
+```{julia}
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = "time", ylabel = "Current")
+lines!(ax, tlist, Ie_L, label = "Bath L", color = :blue, linestyle = :solid)
+lines!(ax, tlist, Ie_R, label = "Bath R", color = :red, linestyle = :solid)
+lines!(ax, tlist, Is_L, label = "Bath L (Steady State)", color = :blue, linestyle = :dash)
+lines!(ax, tlist, Is_R, label = "Bath R (Steady State)", color = :red, linestyle = :dash)
+
+axislegend(ax, position = :rt)
+
+fig
+```
+
+## Version Information
+```{julia}
+HierarchicalEOM.versioninfo()
+```
diff --git a/bibliography.bib b/bibliography.bib
index 0ec9ddd..29654cf 100644
--- a/bibliography.bib
+++ b/bibliography.bib
@@ -11,3 +11,32 @@ @article{gravina2024adaptive
   publisher = {American Physical Society},
   doi = {10.1103/PhysRevResearch.6.023072}
 }
+
+@article{HierarchicalEOM-jl2023,
+  doi = {10.1038/s42005-023-01427-2},
+  url = {https://doi.org/10.1038/s42005-023-01427-2},
+  year = {2023},
+  month = {Oct},
+  publisher = {Nature Portfolio},
+  volume = {6},
+  number = {1},
+  pages = {313},
+  author = {Huang, Yi-Te and Kuo, Po-Chen and Lambert, Neill and Cirio, Mauro and Cross, Simon and Yang, Shen-Liang and Nori, Franco and Chen, Yueh-Nan},
+  title = {An efficient {J}ulia framework for hierarchical equations of motion in open quantum systems},
+  journal = {Communications Physics}
+}
+
+@article{QuTiP-BoFiN2023,
+  title = {QuTiP-BoFiN: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics},
+  author = {Lambert, Neill and Raheja, Tarun and Cross, Simon and Menczel, Paul and Ahmed, Shahnawaz and Pitchford, Alexander and Burgarth, Daniel and Nori, Franco},
+  journal = {Phys. Rev. Res.},
+  volume = {5},
+  issue = {1},
+  pages = {013181},
+  numpages = {18},
+  year = {2023},
+  month = {Mar},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevResearch.5.013181},
+  url = {https://link.aps.org/doi/10.1103/PhysRevResearch.5.013181}
+}