|
| 1 | +--- |
| 2 | +title: "Vacuum Rabi oscillation" |
| 3 | +author: Li-Xun Cai |
| 4 | +date: 2025-01-14 # last update (keep this comment as a reminder) |
| 5 | + |
| 6 | +engine: julia |
| 7 | +--- |
| 8 | + |
| 9 | +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 |
| 10 | + |
| 11 | +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. |
| 12 | + |
| 13 | +## `using` Packages |
| 14 | + |
| 15 | +```{julia} |
| 16 | +using QuantumToolbox |
| 17 | +using CairoMakie # For high-quality plotting |
| 18 | +``` |
| 19 | + |
| 20 | +## Introduction to Jaynes-Cumming model |
| 21 | +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$) --> |
| 22 | + |
| 23 | +The Hamiltonian of JC model is given by |
| 24 | + |
| 25 | +$$ |
| 26 | +\hat{H}_\text{tot} = \hat{H}_{\text{a}} + \hat{H}_{\text{c}} + \hat{H}_{\text{int}} |
| 27 | +$$ |
| 28 | + |
| 29 | +where |
| 30 | + |
| 31 | +- $\hat{H}_{\text{a}} = \frac{\omega_\text{a}}{2} \hat{\sigma}_z$: Hamiltonian for atom alone |
| 32 | +- $\hat{H}_{\text{c}} = \omega_\text{c} \hat{a}^\dagger \hat{a}$: Hamiltonian for cavity alone |
| 33 | +- $\hat{H}_{\text{int}} = \Omega \left( \hat{\sigma}^\dagger + \hat{\sigma} \right) \cdot \left( \hat{a}^\dagger + \hat{a} \right)$: Interaction Hamiltonian for coherent interaction |
| 34 | + |
| 35 | +with |
| 36 | + |
| 37 | +- $\omega_\text{a}$: Frequency gap of the two-level atom |
| 38 | +- $\omega_\text{c}$: Frequency of the cavity's EM mode (This is not specified whether to be in resonance with the atom or not.) |
| 39 | +- $\Omega$ : Coupling strength between the atom and the cavity |
| 40 | +- $\hat{\sigma}_z$ : Pauli's $z$ matrix. Equivalent to $|e\rangle\langle e| - |g\rangle\langle g|$ |
| 41 | +- $\hat{a}$ : Annihilation operator in Fock space <!-- $N$-truncated --> |
| 42 | +- $\hat{\sigma}$ : Lowering operator of atom. Equivalent to $|g\rangle\langle e|$ |
| 43 | + |
| 44 | +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 |
| 45 | + |
| 46 | +$$ |
| 47 | +\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) |
| 48 | +$$ |
| 49 | + |
| 50 | +### Usage of `QuantumToolbox` for JC model in general |
| 51 | + |
| 52 | +Since the packages are already imported with `using` at the beginning, we hereby start by defining some parameters and system Hamiltonian. |
| 53 | + |
| 54 | +```{julia} |
| 55 | +N = 15 # Fock space truncated dimension |
| 56 | +
|
| 57 | +ωa = 1 |
| 58 | +ωc = 1 * ωa # this frequency considers cavity and atom are IN resonance |
| 59 | +ωc_ = 0.8 * ωa # this frequency considers cavity and atom are OFF resonance |
| 60 | +σz = sigmaz() ⊗ qeye(N) # order of tensor product should be consistent throughout |
| 61 | +a = qeye(2) ⊗ destroy(N) |
| 62 | +Ω = 0.05 |
| 63 | +σ = sigmam() ⊗ qeye(N) |
| 64 | +
|
| 65 | +Ha = ωa / 2 * σz |
| 66 | +Hc = ωc * a' * a # the symbol `'` after a `QuantumObject` act as † |
| 67 | +Hc_ = ωc_ * a' * a |
| 68 | +Hint = Ω * (σ * a' + σ' * a) |
| 69 | +
|
| 70 | +Htot = Ha + Hc + Hint |
| 71 | +Htot_ = Ha + Hc_ + Hint |
| 72 | +
|
| 73 | +print(Htot) |
| 74 | +println("\n####################\n") |
| 75 | +print(Htot_) |
| 76 | +``` |
| 77 | + |
| 78 | +## Isolated case |
| 79 | +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. |
| 80 | + |
| 81 | + |
| 82 | +```{julia} |
| 83 | +# define initial state ψ0 |
| 84 | +ψ0 = basis(2,0) ⊗ fock(N, 1) |
| 85 | +ψ0 = ψ0 / norm(ψ0) |
| 86 | +
|
| 87 | +tlist = 0:2:200 # a list of time points of interest |
| 88 | +
|
| 89 | +# define a list of operators whose expectation value dynamics exhibit Rabi oscillation |
| 90 | +eop_ls = [ |
| 91 | + a' * a, # number operator of cavity Fock space |
| 92 | + σ' * σ # excited state population in atom |
| 93 | +] |
| 94 | +
|
| 95 | +sol = sesolve(Htot , ψ0, tlist; e_ops = eop_ls) |
| 96 | +sol_ = sesolve(Htot_, ψ0, tlist; e_ops = eop_ls) |
| 97 | +print(sol) |
| 98 | +println("\n####################\n") |
| 99 | +print(sol_) |
| 100 | +``` |
| 101 | + |
| 102 | +```{julia} |
| 103 | +# Create the first figure for expectation value of photon number |
| 104 | +n = real.(sol.expect[1, :]) |
| 105 | +n_ = real.(sol_.expect[1, :]) |
| 106 | +fign = Figure(size = (600, 350)) |
| 107 | +axn = Axis( |
| 108 | + fign[1, 1], |
| 109 | + title = "⟨n⟩", |
| 110 | + xlabel = "time [1/ωa]", |
| 111 | + ylabel = "expectation value", |
| 112 | + xlabelsize = 15, |
| 113 | + ylabelsize = 15, |
| 114 | + width = 400, |
| 115 | + height = 220 |
| 116 | +) |
| 117 | +lines!(axn, tlist, n, label = "resonant") |
| 118 | +lines!(axn, tlist, n_, label = "off resonant") |
| 119 | +axislegend(axn; position = :rt, labelsize = 12) |
| 120 | +display(fign); |
| 121 | +``` |
| 122 | + |
| 123 | +```{julia} |
| 124 | +# Create the second figure for excited state population |
| 125 | +e = real.(sol.expect[2, :]) |
| 126 | +e_ = real.(sol_.expect[2, :]) |
| 127 | +fige = Figure(size = (600, 350)) |
| 128 | +axe = Axis( |
| 129 | + fige[1, 1], |
| 130 | + title = "⟨e⟩", |
| 131 | + xlabel = "time [1/ωa]", |
| 132 | + ylabel = "probability", |
| 133 | + xlabelsize = 15, |
| 134 | + ylabelsize = 15, |
| 135 | + width = 400, |
| 136 | + height = 220 |
| 137 | +) |
| 138 | +ylims!(axe, 0, 1) |
| 139 | +lines!(axe, tlist, e, label = "resonant") |
| 140 | +lines!(axe, tlist, e_, label = "off resonant") |
| 141 | +axislegend(axe; position = :rb, labelsize = 12) |
| 142 | +display(fige); |
| 143 | +``` |
| 144 | + |
| 145 | +## Dissipative case |
| 146 | + |
| 147 | +In contrast to isolated evolution, a factual system interacts with its surrounding environments, resulting in energy/particle exchange. We are currently interested in observing JC model's Rabi oscillation with the addition of interaction with the external EM field. |
| 148 | + |
| 149 | +We start by reviewing the interaction Hamiltonians between the EM environment and atom/cavity |
| 150 | + |
| 151 | +- Atom: $$ |
| 152 | + \hat{H}_{\text{a}}^\text{int} = \sum_l \alpha_l \left( \hat{b}_l + \hat{b}_l^\dagger \right) \left( \hat{\sigma} + \hat{\sigma}^\dagger \right) |
| 153 | + $$ |
| 154 | +- Cavity: $$ |
| 155 | + \hat{H}_{\text{c}}^\text{int} = \sum_l \beta_l \left( \hat{b}_l + \hat{b}_l^\dagger \right) \left( \hat{a} + \hat{a}^\dagger \right) |
| 156 | + $$ |
| 157 | + |
| 158 | +where for the $l$-th mode |
| 159 | + |
| 160 | +- $\alpha_l$ is the coupling strength with the atom |
| 161 | +- $\beta_l$ is the coupling strength with the cavity |
| 162 | +- $\hat{b}_l$ is the annihilation operator |
| 163 | + |
| 164 | + |
| 165 | +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. |
| 166 | + |
| 167 | +Therefore, the time evolution of the dissipative JC model can be described by the master equation |
| 168 | + |
| 169 | +$$ |
| 170 | +\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) |
| 171 | +$$ |
| 172 | + |
| 173 | +where the $(\Gamma_i, \hat{S}_i)$ pairs are |
| 174 | + |
| 175 | +1. $\kappa \cdot n(\omega_c, T), \hat{a}^\dagger$ |
| 176 | +2. $\kappa \cdot (1 + n(\omega_c, T)), \hat{a}$ |
| 177 | +3. $\gamma \cdot n(\omega_a, T), \hat{\sigma}^\dagger$ |
| 178 | +4. $\gamma \cdot (1 + n(\omega_a, T)), \hat{\sigma}$ |
| 179 | + |
| 180 | +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. |
| 181 | + |
| 182 | +### Solve for evolutions in dissipative case |
| 183 | + |
| 184 | +We can now define variables in `julia` and solve the evolution of dissipative JC model |
| 185 | +```{julia} |
| 186 | +# Bose-Einstein distribution |
| 187 | +function nb(_ω, _KT) |
| 188 | + if _KT < 0 |
| 189 | + println("Temperature should not go below absolute zero, returned 0") |
| 190 | + _n = 0.0 |
| 191 | + else |
| 192 | + _n = (_KT == 0) ? 0 : (1 / (exp(_ω / _KT) - 1)) |
| 193 | + end |
| 194 | + return _n |
| 195 | +end |
| 196 | +
|
| 197 | +# Collapse operators for interaction with the environment with variable |
| 198 | +# dissipation rates, cavity frequency, and thermal energy of the environment |
| 199 | +cop_ls(_γ, _κ, _ωc, _KT) = ( |
| 200 | + √(_κ * nb(_ωc, _KT)) * a', |
| 201 | + √(_κ * (1 + nb(_ωc, _KT))) * a, |
| 202 | + √(_γ * nb(ωa, _KT)) * σ', |
| 203 | + √(_γ * (1 + nb(ωa, _KT))) * σ, |
| 204 | +) |
| 205 | +``` |
| 206 | + |
| 207 | +```{julia} |
| 208 | +# use the same ψ0, tlist, and eop_ls from isolated case |
| 209 | +γ = 4e-3 |
| 210 | +κ = 0.005 |
| 211 | +KT = 0 # for theoretical vacuum EM field |
| 212 | +
|
| 213 | +# `mesolve()` only has one additional keyword argument `c_ops` from `sesolve()` |
| 214 | +sol_me = mesolve(Htot, ψ0, tlist, cop_ls(γ, κ, ωc, KT), e_ops = eop_ls) |
| 215 | +sol_me_ = mesolve(Htot_, ψ0, tlist, cop_ls(γ, κ, ωc_, KT), e_ops = eop_ls) |
| 216 | +
|
| 217 | +print(sol_me) |
| 218 | +println("\n####################\n") |
| 219 | +print(sol_me_) |
| 220 | +``` |
| 221 | + |
| 222 | +```{julia} |
| 223 | +# Create the first figure for the expectation value of photon number |
| 224 | +n_me = real.(sol_me.expect[1, :]) |
| 225 | +n_me_ = real.(sol_me_.expect[1, :]) |
| 226 | +fign_me = Figure(size = (600, 350)) |
| 227 | +axn_me = Axis( |
| 228 | + fign_me[1, 1], |
| 229 | + title = "⟨n⟩", |
| 230 | + xlabel = "time [1/ωa]", |
| 231 | + ylabel = "expectation value", |
| 232 | + xlabelsize = 15, |
| 233 | + ylabelsize = 15, |
| 234 | + width = 400, |
| 235 | + height = 220 |
| 236 | + ) |
| 237 | +lines!(axn_me, tlist, n_me, label = "resonant") |
| 238 | +lines!(axn_me, tlist, n_me_, label = "off resonant") |
| 239 | +axislegend(axn_me; position = :rt, labelsize = 12) |
| 240 | +display(fign_me); |
| 241 | +``` |
| 242 | + |
| 243 | + |
| 244 | +```{julia} |
| 245 | +# Create the second figure for atom excited state population |
| 246 | +e_me = real.(sol_me.expect[2, :]) |
| 247 | +e_me_ = real.(sol_me_.expect[2, :]) |
| 248 | +fige_me = Figure(size = (600, 350)) |
| 249 | +axe_me = Axis( |
| 250 | + fige_me[1, 1], |
| 251 | + title = "⟨e⟩", |
| 252 | + xlabel = "time [1/ωa]", |
| 253 | + ylabel = "probability", |
| 254 | + xlabelsize = 15, |
| 255 | + ylabelsize = 15, |
| 256 | + width = 400, |
| 257 | + height = 220 |
| 258 | + ) |
| 259 | +ylims!(axe_me, 0, 1) |
| 260 | +lines!(axe_me, tlist, e_me, label = "resonant") |
| 261 | +lines!(axe_me, tlist, e_me_, label = "off resonant") |
| 262 | +axislegend(axe_me; position = :rb, labelsize = 12) |
| 263 | +display(fige_me); |
| 264 | +``` |
| 265 | + |
| 266 | +### Note for RWA |
| 267 | + |
| 268 | +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()`](). |
| 269 | + |
| 270 | + |
| 271 | + |
| 272 | +## Version Information |
| 273 | +```{julia} |
| 274 | +QuantumToolbox.versioninfo() |
| 275 | +``` |
| 276 | + |
0 commit comments