Add more tutorials for HierarchicalEOM.jl (qutip#15)

Author: ytdHuang

.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.

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()
+```

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

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()
+```
+