contextaul revisions of rabi.qmd

TendonFFF · TendonFFF · commit 4ded5e97e70a · 2025-01-16T18:23:18.000+08:00
diff --git a/QuantumToolbox.jl/time_evolution/rabi.qmd b/QuantumToolbox.jl/time_evolution/rabi.qmd
@@ -8,10 +8,10 @@ engine: julia
 
 Inspirations taken from this [QuTiP tutorial](https://nbviewer.org/urls/qutip.org/qutip-tutorials/tutorials-v5/time-evolution/004_rabi-oscillations.ipynb) by J.R. Johansson, P.D. Nation, and C. Staufenbiel
 
-In this notebook, the usage of `QuantumToolbox.sesolve()` and `QuantumToolbox.mesolve()` will be demonstrated with Jaynes-Cummings model to observe Rabi Oscillation in the isolated case and the dissipative case. In dissipative case, a bosonic environment would interact with JC model's cavity and two-level atom.
+In this notebook, the usage of `QuantumToolbox.sesolve` and `QuantumToolbox.mesolve` will be demonstrated with Jaynes-Cummings model to observe Rabi oscillations in the isolated case and the dissipative case. In dissipative case, a bosonic environment would interact with the cavity and two-level atom in JC model.
 
 ## Introduction to Jaynes-Cumming model
-Jaynes-Cummings (JC) model, a simplest quantum mechanical model for light-matter interaction, describes an atom interacting with an external electromagnetic field. To simplify the interaction at the time of the model proposal, JC model considered a two-level atom interacting with a single bosonic mode (or you can consider it a single-mode cavity), i.e., a two-dimensional Hilbert space interacting with a Fock space. <!-- (with finite truncation $N$) -->
+Jaynes-Cummings (JC) model, a simplest quantum mechanical model for light-matter interaction, describes an atom interacting with an external electromagnetic field. To simplify the interaction, JC model considered a two-level atom interacting with a single bosonic mode (or you can consider it a single-mode cavity).
 
 The Hamiltonian of JC model is given by
 
@@ -27,116 +27,91 @@ where
 
 with
 
-- $\omega_\text{a}$: Frequency gap of the two-level atom
+- $\omega_\text{a}$: Frequency of the two-level atom
 - $\omega_\text{c}$: Frequency of the cavity's EM mode (This is not specified whether to be in resonance with the atom or not.)
 - $\Omega$         : Coupling strength between the atom and the cavity
-- $\hat{\sigma}_z$ : Pauli's $z$ matrix. Equivalent to $|e\rangle\langle e| - |g\rangle\langle g|$
-- $\hat{a}$        : Annihilation operator in Fock space <!-- $N$-truncated -->
+- $\hat{\sigma}_z$ : Pauli-$Z$ matrix. Equivalent to $|e\rangle\langle e| - |g\rangle\langle g|$
+- $\hat{a}$        : Annihilation operator of single-mode cavity <!-- $N$-truncated -->
 - $\hat{\sigma}$   : Lowering operator of atom. Equivalent to $|g\rangle\langle e|$
 
-By applying rotating wave approximation (RWA), the counter rotating terms ($\hat{\sigma} \cdot \hat{a}$ and its hermition conjugate), which rotate considerably faster than the others in the interaction picture, are ignored, yielding
+By applying [rotating wave approximation (RWA)](https://en.wikipedia.org/wiki/Rotating-wave_approximation), the counter rotating terms ($\hat{\sigma} \cdot \hat{a}$ and its Hermitian conjugate), which rotate considerably faster than the others in the interaction picture, are ignored, yielding
 
 $$
 \hat{H}_\text{tot} \approx \hat{H}_{\text{a}} + \hat{H}_{\text{c}} + \Omega \left( \hat{\sigma} \cdot \hat{a}^\dagger + \hat{\sigma}^\dagger \cdot \hat{a} \right)
 $$
 
 ### Usage of `QuantumToolbox` for JC model in general
 
+##### import:
 ```{julia}
 using QuantumToolbox
-using CairoMakie     # For high-quality plotting
+import QuantumToolbox: ⊗, sigmaz, sigmam, destroy, qeye, basis, fock, n_thermal, mesolve, sesolve, versioninfo
+using CairoMakie
+import CairoMakie: Figure, Axis, lines!, axislegend, xlims!, display
 ```
 
 ```{julia}
-N = 15 # Fock space truncated dimension
+N = 2 # Fock space truncated dimension
 
 ωa = 1
-ωc = 1 * ωa    # this frequency considers cavity and atom are IN resonance
-ωc_ = 0.8 * ωa # this frequency considers cavity and atom are OFF resonance
+ωc = 1 * ωa    # considering cavity and atom are in resonance
 σz = sigmaz() ⊗ qeye(N) # order of tensor product should be consistent throughout
 a  = qeye(2)  ⊗ destroy(N)  
 Ω  = 0.05
 σ  = sigmam() ⊗ qeye(N)
 
 Ha = ωa / 2 * σz
 Hc = ωc * a' * a # the symbol `'` after a `QuantumObject` act as adjoint
-Hc_ = ωc_ * a' * a
 Hint = Ω * (σ * a' + σ' * a)
 
-Htot  = Ha + Hc  + Hint
-Htot_ = Ha + Hc_ + Hint
+Htot  = Ha + Hc + Hint
 
 print(Htot)
-println("\n####################\n")
-print(Htot_)
 ```
 
 ## Isolated case
-For the case of JC system being isolated, i.e., with no interaction with the surrounding environment, the system's time-evolution is governed solely by Schrödinger equation $\hat{H}|\psi\rangle = \partial_t|\psi(t)\rangle$. Using `QuantumToolbox.sesolve()` is ideal for pure state evolution.
+For the case of JC model being isolated, i.e., no interaction with the surrounding environment, the time-evolution is governed solely by Schrödinger equation $\hat{H}|\psi(t)\rangle = \partial_t|\psi(t)\rangle$. Using `QuantumToolbox.sesolve` is ideal for pure state evolution.
 
+For the context of [Rabi problem](https://en.wikipedia.org/wiki/Rabi_problem), we set the initial state $\psi_0 = |e\rangle \otimes |0\rangle$ where $|e\rangle$ is the excited state of atom and $|0\rangle$ is the vacuum state of cavity.
 
 ```{julia}
-# define initial state ψ0
-ψ0 = basis(2,0) ⊗ fock(N, 1) 
-ψ0 = ψ0 / norm(ψ0)
+e_ket = basis(2,0) 
+ψ0 = e_ket ⊗ fock(N, 0)
 
-tlist = 0:2:200 # a list of time points of interest
+tlist = 0:2.5:1000 # a list of time points of interest
 
 # define a list of operators whose expectation value dynamics exhibit Rabi oscillation
 eop_ls = [
-    a' * a, # number operator of cavity Fock space
-    σ' * σ  # excited state population in atom
+    a' * a,                      # number operator of cavity
+    (e_ket * e_ket') ⊗ qeye(N), # excited state population in atom
 ]
 
-sol  = sesolve(Htot , ψ0, tlist; e_ops = eop_ls)
-sol_ = sesolve(Htot_, ψ0, tlist; e_ops = eop_ls)
+sol = sesolve(Htot , ψ0, tlist; e_ops = eop_ls)
 print(sol)
-println("\n####################\n")
-print(sol_)
 ```
 
-Create the first figure for expectation value of photon number
+Compare the dynamics of $| e \rangle\langle e|$ alongside $a^\dagger a$
 ```{julia}
-n  = real.(sol.expect[1, :])
-n_ = real.(sol_.expect[1, :])
-fign = Figure(size = (600, 350))
-axn = Axis(
-    fign[1, 1],
-    title = L"$\langle n\rangle$",
+n = real.(sol.expect[1, :])
+e = real.(sol.expect[2, :])
+fig_se = Figure(size = (600, 350))
+ax_se = Axis(
+    fig_se[1, 1],
     xlabel = L"time $[1/\omega_a]$", 
     ylabel = "expectation value", 
     xlabelsize = 15, 
     ylabelsize = 15,
     width = 400,
     height = 220
 )
-lines!(axn, tlist, n,  label = "resonant")
-lines!(axn, tlist, n_, label = "off resonant")
-axislegend(axn; position = :rt, labelsize = 12)
-display(fign);
+xlims!(ax_se, 0, 400)
+lines!(ax_se, tlist, n, label = L"$\langle a^\dagger a \rangle$")
+lines!(ax_se, tlist, e, label = L"$P_e$")
+axislegend(ax_se; position = :rt, labelsize = 15)
+display(fig_se);
 ```
 
-Create the second figure for excited state population
-```{julia}
-e  = real.(sol.expect[2, :])
-e_ = real.(sol_.expect[2, :]) 
-fige = Figure(size = (600, 350))
-axe = Axis(
-    fige[1, 1], 
-    title = L"$\langle e \rangle$",
-    xlabel = L"time $[1/\omega_a]$", 
-    ylabel = "probability", 
-    xlabelsize = 15, 
-    ylabelsize = 15,
-    width = 400,
-    height = 220
-)
-ylims!(axe, 0, 1)
-lines!(axe, tlist, e,  label = "resonant")
-lines!(axe, tlist, e_, label = "off resonant")
-axislegend(axe; position = :rb, labelsize = 12)
-display(fige);
-```
+In the above plot, the behaviour of the energy exchange between the atom and the cavity is clearly visible, addressing the Rabi problem.
 
 ## Dissipative case
 
@@ -157,16 +132,13 @@ where for the $l$-th mode
 - $\beta_l$ is the coupling strength with the cavity
 - $\hat{b}_l$ is the annihilation operator
 
-
-Following the RWA approach previously mentioned and the standard procedure of Born-Markovian master equation, we obtain $\kappa$, the cavity dissipation rate, and $\gamma$, the atom dissipation rate.
-
-Therefore, the time evolution of the dissipative JC model can be described by the master equation
+Following the RWA approach previously mentioned and the standard procedure of [Born-Markovian approximation](https://en.wikiversity.org/wiki/Open_Quantum_Systems/The_Quantum_Optical_Master_Equation), we obtain $\kappa$, the cavity dissipation rate, and $\gamma$, the atom dissipation rate. Therefore, the time evolution of the dissipative JC model can be described by the [Lindblad master equation](https://en.wikipedia.org/wiki/Lindbladian)
 
 $$
 \dot{\hat{\rho}} = -\frac{i}{\hbar} [\hat{H}, \hat{\rho}] + \sum_{i = 1}^4 \mathcal{D}[\sqrt{\Gamma_i} \hat{S}_i]\left(\hat{\rho}\right)
 $$
 
-where the $(\Gamma_i, \hat{S}_i)$ pairs are
+where $\sqrt{\Gamma_i} \hat{S}_i$ are the collapse operators, given by 
 
 |$i$| $\Gamma_i$ | $\hat{S}_i$ |
 |---| :--------: | :---------: |
@@ -175,18 +147,21 @@ where the $(\Gamma_i, \hat{S}_i)$ pairs are
 |3|$\gamma \cdot n(\omega_a, T)$|$\hat{\sigma}^\dagger$|
 |4|$\gamma \cdot [1 + n(\omega_a, T)]$|$\hat{\sigma}$|
 
-with $n(\omega, T)$ being the Bose-Einstein distribution for the EM environment and $\mathcal{D}[\hat{\mathcal{O}}]\left(\cdot\right) = \hat{\mathcal{O}} \left(\cdot\right) \hat{\mathcal{O}}^\dagger - \frac{1}{2} \{ \hat{\mathcal{O}}^\dagger \hat{\mathcal{O}}, \cdot \}$ being the Lindblad dissipator.
+with $n(\omega, T)$ being the Bose-Einstein distribution for the EM environment and 
+$$
+\mathcal{D}[\hat{\mathcal{O}}]\left(\cdot\right) = \hat{\mathcal{O}} \left(\cdot\right) \hat{\mathcal{O}}^\dagger - \frac{1}{2} \{ \hat{\mathcal{O}}^\dagger \hat{\mathcal{O}}, \cdot \}
+$$ 
+being the Lindblad dissipator.
 
 ### Solve for evolutions in dissipative case
 
 We can now define variables in `julia` and solve the evolution of dissipative JC model 
 ```{julia}
-# Collapse operators for interaction with the environment with variable 
-# dissipation rates, cavity frequency, and thermal energy of the environment
-# `n_thermal()` gives Bose-Einstein distribution
-cop_ls(_γ, _κ, _ωc, _KT) = (
-    √(_κ * n_thermal(_ωc, _KT)) * a', 
-    √(_κ * (1 + n_thermal(_ωc, _KT))) * a, 
+# Collapse operators for interaction with the environment with variable dissipation rates 
+# and thermal energy of the environment. `n_thermal()` gives Bose-Einstein distribution
+cop_ls(_γ, _κ, _KT) = (
+    √(_κ * n_thermal(ωc, _KT)) * a', 
+    √(_κ * (1 + n_thermal(ωc, _KT))) * a, 
     √(_γ * n_thermal(ωa, _KT)) * σ', 
     √(_γ * (1 + n_thermal(ωa, _KT))) * σ, 
 )
@@ -195,68 +170,58 @@ cop_ls(_γ, _κ, _ωc, _KT) = (
 ```{julia}
 # use the same ψ0, tlist, and eop_ls from isolated case
 γ = 4e-3
-κ = 0.005
+κ = 7e-3
 KT = 0 # for theoretical vacuum EM field
 
 # `mesolve()` only has one additional keyword argument `c_ops` from `sesolve()`
-sol_me  = mesolve(Htot,  ψ0, tlist, cop_ls(γ, κ, ωc,  KT), e_ops = eop_ls)
-sol_me_ = mesolve(Htot_, ψ0, tlist, cop_ls(γ, κ, ωc_, KT), e_ops = eop_ls)
+sol_me  = mesolve(Htot,  ψ0, tlist, cop_ls(γ, κ, KT), e_ops = eop_ls)
 
 print(sol_me)
-println("\n####################\n")
-print(sol_me_)
 ```
 
-Create the first figure for the expectation value of photon number
 ```{julia}
-n_me  = real.(sol_me.expect[1, :])
-n_me_ = real.(sol_me_.expect[1, :])
-fign_me = Figure(size = (600, 350))
-axn_me = Axis(
-    fign_me[1, 1],
-    title = L"$\langle n\rangle$",
+n_me = real.(sol_me.expect[1, :])
+e_me = real.(sol_me.expect[2, :])
+
+fig_me = Figure(size = (600, 350))
+ax_me = Axis(
+    fig_me[1, 1],
     xlabel = L"time $[1/\omega_a]$", 
     ylabel = "expectation value", 
     xlabelsize = 15, 
     ylabelsize = 15,
     width = 400,
     height = 220
-    )
-lines!(axn_me, tlist, n_me,  label = "resonant")
-lines!(axn_me, tlist, n_me_, label = "off resonant")
-axislegend(axn_me; position = :rt, labelsize = 12)
-display(fign_me);
+)
+lines!(ax_me, tlist, n_me, label = L"\langle a^\dagger a \rangle")
+lines!(ax_me, tlist, e_me, label = L"$P_e$")
+axislegend(ax_me; position = :rt, labelsize = 15)
+display(fig_me);
 ```
 
-Create the second figure for atom excited state population
+From the above example, one can see that the dissipative system is losing energy over time and asymptoting to zero. We can further consdider the near-vacuum environment with finite temperature.
+
 ```{julia}
-e_me  = real.(sol_me.expect[2, :])
-e_me_ = real.(sol_me_.expect[2, :]) 
-fige_me = Figure(size = (600, 350))
-axe_me = Axis(
-    fige_me[1, 1], 
-    title = L"$\langle e \rangle$",
+sol_me_  = mesolve(Htot,  ψ0, tlist, cop_ls(γ, κ, 0.3 * ωa), e_ops = eop_ls) # replace KT with finite temperature
+
+n_me_ = real.(sol_me_.expect[1, :])
+e_me_ = real.(sol_me_.expect[2, :])
+fig_me_ = Figure(size = (600, 350))
+ax_me_ = Axis(
+    fig_me_[1, 1],
     xlabel = L"time $[1/\omega_a]$", 
-    ylabel = "probability", 
+    ylabel = "expectation value", 
     xlabelsize = 15, 
     ylabelsize = 15,
     width = 400,
     height = 220
-    )
-ylims!(axe_me, 0, 1)
-lines!(axe_me, tlist, e_me,  label = "resonant")
-lines!(axe_me, tlist, e_me_, label = "off resonant")
-axislegend(axe_me; position = :rt, labelsize = 12)
-display(fige_me);
+)
+lines!(ax_me_, tlist, n_me_, label = L"\langle a^\dagger a \rangle")
+lines!(ax_me_, tlist, e_me_, label = L"$P_e$")
+axislegend(ax_me_; position = :rt, labelsize = 15)
+display(fig_me_);
 ```
-
-<!-- 
-### Note for RWA
-
-Without applying RWA, things are more interesting in the \emph{Ultrastrong Coupling} (USC) regime, where $\Omega \sim \omega_c$. In that regime, the \emph{localized} master equation (the one we used for this tutorial) fails. As that phenomenon is worth itself as an independent tutorial, we will pass that on to [this tutorial for `liouvillian_generalized()`]().
--->
-
-
+Despite the asymptotic behaviour persisting, one can see that they no longer approach zero and instead find a steady condition above zero. That is, the system eventually becomes thermalized by the environment.
 
 ## Version Information
 ```{julia}
