Setup structure for tutorials website

HierarchicalEOM.jl/cavityQED.qmd

@@ -0,0 +1,264 @@
+---
+title: "Cavity QED system"
+author: Shen-Liang Yang, Yi-Te Huang
+date: last-modified
+date-format: iso
+
+engine: julia
+---
+
+## Package Imports
+```{julia}
+using HierarchicalEOM
+using CairoMakie
+```
+
+## Introduction
+
+Cavity quantum electrodynamics (cavity QED) is an important topic for studying the interaction between atoms (or other particles) and light confined in a reflective cavity, under conditions where the quantum nature of photons is significant.
+
+## Hamiltonian
+
+The Jaynes-Cummings model is a standard model in the realm of cavity QED. It illustrates the interaction between a two-level atom ($\textrm{A}$) and a quantized single-mode within a cavity ($\textrm{c}$).
+
+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.
+
+$$
+\begin{aligned}
+H_{\textrm{s}}&=H_{\textrm{A}}+H_{\textrm{c}}+H_{\textrm{int}},\\
+H_{\textrm{A}}&=\frac{\omega_A}{2}\sigma_z,\\
+H_{\textrm{c}}&=\omega_{\textrm{c}} a^\dagger a,\\
+H_{\textrm{int}}&=g (a^\dagger\sigma^-+a\sigma^+),
+\end{aligned}
+$$
+
+where $\sigma^-$ ($\sigma^+$) is the annihilation (creation) operator of the atom, and $a$ ($a^\dagger$) is the annihilation (creation) operator of the cavity.
+
+Furthermore, we consider the system is coupled to a bosonic reservoir ($\textrm{b}$). The total Hamiltonian is given by $H_{\textrm{Total}}=H_\textrm{s}+H_\textrm{b}+H_\textrm{sb}$, where $H_\textrm{b}$ and $H_\textrm{sb}$ takes the form
+
+$$
+\begin{aligned}
+H_{\textrm{b}}    &=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k},\\
+H_{\textrm{sb}}   &=(a+a^\dagger)\sum_{k}g_{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}$). Also, $H_{\textrm{sb}}$ illustrates the interaction between the cavity and the bosonic 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}
+N = 3 ## system cavity Hilbert space cutoff
+ωA = 2
+ωc = 2
+g = 0.1
+
+# operators
+a_c = destroy(N)
+I_c = qeye(N)
+σz_A = sigmaz()
+σm_A = sigmam()
+I_A = qeye(2)
+
+# operators in tensor-space
+a = tensor(a_c, I_A)
+σz = tensor(I_c, σz_A)
+σm = tensor(I_c, σm_A)
+
+# Hamiltonian
+H_A = 0.5 * ωA * σz
+H_c = ωc * a' * a
+H_int = g * (a' * σm + a * σm')
+H_s = H_A + H_c + H_int
+
+# initial state
+ψ0 = tensor(basis(N, 0), basis(2, 0));
+```
+
+## 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.01
+W = 1
+kT = 0.025
+N = 20
+Bath = Boson_DrudeLorentz_Pade(a + a', Γ, W, kT, N)
+```
+
+Before incorporating the correlation function into the HEOMLS matrix, it is essential to verify (by using `correlation_function`) if the total number of exponentials for the reservoir sufficiently describes the practical situation.
+
+```{julia}
+tlist_test = 0:0.1:10;
+Bath_test = Boson_DrudeLorentz_Pade(a + a', Γ, W, kT, 1000);
+Ct = correlation_function(Bath, tlist_test);
+Ct2 = correlation_function(Bath_test, tlist_test)
+
+# plot
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = L"t", ylabel = L"C(t)")
+lines!(ax, tlist_test, real(Ct2), label = L"$N=1000$ (real part)", linestyle = :solid)
+lines!(ax, tlist_test, real(Ct),  label = L"$N=20$ (real part)", linestyle = :dash)
+lines!(ax, tlist_test, imag(Ct2), label = L"$N=1000$ (imag part)", linestyle = :solid)
+lines!(ax, tlist_test, imag(Ct),  label = L"$N=20$ (imag part)", linestyle = :dash)
+
+axislegend(ax, position = :rt)
+
+fig
+```
+
+## Construct HEOMLS matrix
+
+Here, we consider an incoherent pumping to the atom, which can be described by an Lindblad dissipator (see [here](https://qutip.org/HierarchicalEOM.jl/stable/heom_matrix/master_eq/) for more details).
+
+Furthermore, we set the [important threshold](https://qutip.org/HierarchicalEOM.jl/stable/heom_matrix/HEOMLS_intro/#doc-Importance-Value-and-Threshold) to be `1e-6`.
+
+```{julia}
+pump = 0.01
+J_pump = sqrt(pump) * σm'
+tier = 2
+M_Heom = M_Boson(H_s, tier, threshold = 1e-6, Bath)
+M_Heom = addBosonDissipator(M_Heom, J_pump)
+```
+
+## Solve time evolution of ADOs
+
+```{julia}
+t_list = 0:1:500
+sol_H = HEOMsolve(M_Heom, ψ0, t_list; e_ops = [σz, a' * a])
+```
+
+## Solve stationary state of ADOs
+
+```{julia}
+steady_H = steadystate(M_Heom);
+```
+
+## Expectation values
+
+observable of atom: $\sigma_z$
+
+```{julia}
+σz_evo_H = real(sol_H.expect[1, :])
+σz_steady_H = expect(σz, steady_H)
+```
+
+observable of cavity: $a^\dagger a$ (average photon number)
+
+```{julia}
+np_evo_H = real(sol_H.expect[2, :])
+np_steady_H = expect(a' * a, steady_H)
+```
+
+plot results
+
+```{julia}
+fig = Figure(size = (600, 350))
+
+ax1 = Axis(fig[1, 1], xlabel = L"t")
+lines!(ax1, t_list, σz_evo_H, label = L"\langle \sigma_z \rangle", linestyle = :solid)
+lines!(ax1, t_list, ones(length(t_list)) .* σz_steady_H, label = L"\langle \sigma_z \rangle ~~(\textrm{steady})", linestyle = :dash)
+axislegend(ax1, position = :rt)
+
+ax2 = Axis(fig[2, 1], xlabel = L"t")
+lines!(ax2, t_list, np_evo_H, label = L"\langle a^\dagger a \rangle", linestyle = :solid)
+lines!(ax2, t_list, ones(length(t_list)) .* np_steady_H, label = L"\langle a^\dagger a \rangle ~~(\textrm{steady})", linestyle = :dash)
+axislegend(ax2, position = :rt)
+
+fig
+```
+
+## Power spectrum
+
+```{julia}
+ω_list = 1:0.01:3
+psd_H = PowerSpectrum(M_Heom, steady_H, a, ω_list)
+
+# plot
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = L"\omega")
+lines!(ax, ω_list, psd_H)
+
+fig
+```
+
+## Compare with Master Eq. approach
+
+The Lindblad master equations which describes the cavity couples to an extra bosonic reservoir with [Drude-Lorentzian spectral density](https://qutip.org/HierarchicalEOM.jl/stable/bath_boson/Boson_Drude_Lorentz/#Boson-Drude-Lorentz) is given by
+
+```{julia}
+# Drude_Lorentzian spectral density
+Drude_Lorentz(ω, Γ, W) = 4 * Γ * W * ω / ((ω)^2 + (W)^2)
+
+# Bose-Einstein distribution
+n_b(ω, kT) = 1 / (exp(ω / kT) - 1)
+
+# build the jump operators
+jump_op = [
+    sqrt(Drude_Lorentz(ωc, Γ, W) * (n_b(ωc, kT) + 1)) * a,
+    sqrt(Drude_Lorentz(ωc, Γ, W) * (n_b(ωc, kT)))     * a',
+    J_pump
+];
+
+# construct the HEOMLS matrix for master equation
+M_master = M_S(H_s)
+M_master = addBosonDissipator(M_master, jump_op)
+
+# time evolution
+sol_M = HEOMsolve(M_master, ψ0, t_list; e_ops = [σz, a' * a]);
+
+# steady state
+steady_M = steadystate(M_master);
+
+# expectation value of σz
+σz_evo_M = real(sol_M.expect[1, :])
+σz_steady_M = expect(σz, steady_M)
+
+# average photon number
+np_evo_M = real(sol_M.expect[2, :])
+np_steady_M = expect(a' * a, steady_M);
+```
+
+plot results
+
+```{julia}
+fig = Figure(size = (600, 350))
+
+ax1 = Axis(fig[1, 1], xlabel = L"t")
+lines!(ax1, t_list, σz_evo_M, label = L"\langle \sigma_z \rangle", linestyle = :solid)
+lines!(ax1, t_list, ones(length(t_list)) .* σz_steady_M, label = L"\langle \sigma_z \rangle ~~(\textrm{steady})", linestyle = :dash)
+axislegend(ax1, position = :rt)
+
+ax2 = Axis(fig[2, 1], xlabel = L"t")
+lines!(ax2, t_list, np_evo_M, label = L"\langle a^\dagger a \rangle", linestyle = :solid)
+lines!(ax2, t_list, ones(length(t_list)) .* np_steady_M, label = L"\langle a^\dagger a \rangle ~~(\textrm{steady})", linestyle = :dash)
+axislegend(ax2, position = :rt)
+
+fig
+```
+
+We can also calculate the power spectrum
+
+```{julia}
+ω_list = 1:0.01:3
+psd_M = PowerSpectrum(M_master, steady_M, a, ω_list)
+
+# plot
+fig = Figure(size = (500, 350))
+ax = Axis(fig[1, 1], xlabel = L"\omega")
+lines!(ax, ω_list, psd_M)
+
+fig
+```
+
+Due to the weak coupling between the system and an extra bosonic environment, the Master equation's outcome is expected to be similar to the results obtained from the HEOM method.
+
+## Version Information
+```{julia}
+HierarchicalEOM.versioninfo()
+```

HierarchicalEOM.jl/toc.md

@@ -0,0 +1,17 @@
+---
+title: "Tutorials for `HierarchicalEOM.jl`"
+listing:
+  id: HierarchicalEOM-listings
+  type: table
+  date-format: iso
+  sort: false
+  sort-ui: false
+  fields: [date, title, author]
+  contents:
+    - "cavityQED.qmd"
+---
+
+The following tutorials demonstrate and introduce specific functionality of `HierarchicalEOM.jl`.
+
+::: {#HierarchicalEOM-listings}
+:::

Makefile

@@ -6,7 +6,7 @@ include _environment
 default: help
 
 render:
-	${JULIA} --project=@. -e 'import Pkg; Pkg.instantiate(); Pkg.resolve(); Pkg.precompile(); using QuantumToolbox, HierarchicalEOM;'
+	${JULIA} --project=@. -e 'import Pkg; Pkg.resolve(); Pkg.instantiate(); Pkg.precompile(); using QuantumToolbox, HierarchicalEOM;'
 	${JULIA} --project=@. -e 'using QuantumToolbox, HierarchicalEOM;'
 	${QUARTO} render
 

Project.toml

@@ -2,3 +2,4 @@
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 HierarchicalEOM = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
+QuartoNotebookRunner = "4c0109c6-14e9-4c88-93f0-2b974d3468f4"