|
| 1 | +--- |
| 2 | +title: "Cavity QED system" |
| 3 | +author: Shen-Liang Yang, Yi-Te Huang |
| 4 | +date: last-modified |
| 5 | +date-format: iso |
| 6 | + |
| 7 | +engine: julia |
| 8 | +--- |
| 9 | + |
| 10 | +## Package Imports |
| 11 | +```{julia} |
| 12 | +using HierarchicalEOM |
| 13 | +using LaTeXStrings |
| 14 | +import Plots |
| 15 | +``` |
| 16 | + |
| 17 | +## Introduction |
| 18 | + |
| 19 | +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. |
| 20 | + |
| 21 | +## Hamiltonian |
| 22 | + |
| 23 | +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}$). |
| 24 | + |
| 25 | +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. |
| 26 | + |
| 27 | +$$ |
| 28 | +\begin{aligned} |
| 29 | +H_{\textrm{s}}&=H_{\textrm{A}}+H_{\textrm{c}}+H_{\textrm{int}},\\ |
| 30 | +H_{\textrm{A}}&=\frac{\omega_A}{2}\sigma_z,\\ |
| 31 | +H_{\textrm{c}}&=\omega_{\textrm{c}} a^\dagger a,\\ |
| 32 | +H_{\textrm{int}}&=g (a^\dagger\sigma^-+a\sigma^+), |
| 33 | +\end{aligned} |
| 34 | +$$ |
| 35 | + |
| 36 | +where $\sigma^-$ ($\sigma^+$) is the annihilation (creation) operator of the atom, and $a$ ($a^\dagger$) is the annihilation (creation) operator of the cavity. |
| 37 | + |
| 38 | +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 |
| 39 | + |
| 40 | +$$ |
| 41 | +\begin{aligned} |
| 42 | +H_{\textrm{b}} &=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k},\\ |
| 43 | +H_{\textrm{sb}} &=(a+a^\dagger)\sum_{k}g_{k}(b_k + b_k^{\dagger}). |
| 44 | +\end{aligned} |
| 45 | +$$ |
| 46 | + |
| 47 | +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. |
| 48 | + |
| 49 | +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. |
| 50 | + |
| 51 | +```{julia} |
| 52 | +N = 3 ## system cavity Hilbert space cutoff |
| 53 | +ωA = 2 |
| 54 | +ωc = 2 |
| 55 | +g = 0.1 |
| 56 | +
|
| 57 | +## operators |
| 58 | +a_c = destroy(N) |
| 59 | +I_c = qeye(N) |
| 60 | +σz_A = sigmaz() |
| 61 | +σm_A = sigmam() |
| 62 | +I_A = qeye(2) |
| 63 | +
|
| 64 | +## operators in tensor-space |
| 65 | +a = tensor(a_c, I_A) |
| 66 | +σz = tensor(I_c, σz_A) |
| 67 | +σm = tensor(I_c, σm_A) |
| 68 | +
|
| 69 | +## Hamiltonian |
| 70 | +H_A = 0.5 * ωA * σz |
| 71 | +H_c = ωc * a' * a |
| 72 | +H_int = g * (a' * σm + a * σm') |
| 73 | +
|
| 74 | +H_s = H_A + H_c + H_int |
| 75 | +
|
| 76 | +## initial state |
| 77 | +ψ0 = tensor(basis(N, 0), basis(2, 0)); |
| 78 | +``` |
| 79 | + |
| 80 | +## Construct bath objects |
| 81 | +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: |
| 82 | + |
| 83 | +- the coupling strength $\Gamma$ between system and reservoir |
| 84 | +- the band-width $W$ |
| 85 | +- the product of the Boltzmann constant $k$ and the absolute temperature $T$ : $kT$ |
| 86 | +- the total number of exponentials for the reservoir $(N + 1)$ |
| 87 | + |
| 88 | +```{julia} |
| 89 | +Γ = 0.01 |
| 90 | +W = 1 |
| 91 | +kT = 0.025 |
| 92 | +N = 20 |
| 93 | +Bath = Boson_DrudeLorentz_Pade(a + a', Γ, W, kT, N) |
| 94 | +``` |
| 95 | + |
| 96 | +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. |
| 97 | + |
| 98 | +```{julia} |
| 99 | +tlist_test = 0:0.1:10; |
| 100 | +
|
| 101 | +Bath_test = Boson_DrudeLorentz_Pade(a + a', Γ, W, kT, 1000); |
| 102 | +Ct = correlation_function(Bath, tlist_test); |
| 103 | +Ct2 = correlation_function(Bath_test, tlist_test); |
| 104 | +
|
| 105 | +Plots.plot(tlist_test, real(Ct), label = "N=20 (real part )", linestyle = :dash, linewidth = 3) |
| 106 | +Plots.plot!(tlist_test, real(Ct2), label = "N=1000 (real part)", linestyle = :solid, linewidth = 3) |
| 107 | +Plots.plot!(tlist_test, imag(Ct), label = "N=20 (imag part)", linestyle = :dash, linewidth = 3) |
| 108 | +Plots.plot!(tlist_test, imag(Ct2), label = "N=1000 (imag part)", linestyle = :solid, linewidth = 3) |
| 109 | +
|
| 110 | +Plots.xaxis!("t") |
| 111 | +Plots.yaxis!("C(t)") |
| 112 | +``` |
| 113 | + |
| 114 | +## Construct HEOMLS matrix |
| 115 | + |
| 116 | +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). |
| 117 | + |
| 118 | +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`. |
| 119 | + |
| 120 | +```{julia} |
| 121 | +pump = 0.01 |
| 122 | +J_pump = sqrt(pump) * σm' |
| 123 | +
|
| 124 | +tier = 2 |
| 125 | +M_Heom = M_Boson(H_s, tier, threshold = 1e-6, Bath) |
| 126 | +M_Heom = addBosonDissipator(M_Heom, J_pump) |
| 127 | +``` |
| 128 | + |
| 129 | +## Solve time evolution of ADOs |
| 130 | + |
| 131 | +```{julia} |
| 132 | +t_list = 0:1:500 |
| 133 | +sol_H = HEOMsolve(M_Heom, ψ0, t_list; e_ops = [σz, a' * a]) |
| 134 | +``` |
| 135 | + |
| 136 | +## Solve stationary state of ADOs |
| 137 | + |
| 138 | +```{julia} |
| 139 | +steady_H = steadystate(M_Heom); |
| 140 | +``` |
| 141 | + |
| 142 | +## Expectation values |
| 143 | + |
| 144 | +observable of atom: $\sigma_z$ |
| 145 | + |
| 146 | +```{julia} |
| 147 | +σz_evo_H = real(sol_H.expect[1, :]) |
| 148 | +σz_steady_H = expect(σz, steady_H) |
| 149 | +``` |
| 150 | + |
| 151 | +observable of cavity: $a^\dagger a$ (average photon number) |
| 152 | + |
| 153 | +```{julia} |
| 154 | +np_evo_H = real(sol_H.expect[2, :]) |
| 155 | +np_steady_H = expect(a' * a, steady_H) |
| 156 | +``` |
| 157 | + |
| 158 | +plot results |
| 159 | + |
| 160 | +```{julia} |
| 161 | +p1 = Plots.plot( |
| 162 | + t_list, |
| 163 | + [σz_evo_H, ones(length(t_list)) .* σz_steady_H], |
| 164 | + label = [L"\langle \sigma_z \rangle" L"\langle \sigma_z \rangle ~~(\textrm{steady})"], |
| 165 | + linewidth = 3, |
| 166 | + linestyle = [:solid :dash], |
| 167 | +) |
| 168 | +p2 = Plots.plot( |
| 169 | + t_list, |
| 170 | + [np_evo_H, ones(length(t_list)) .* np_steady_H], |
| 171 | + label = [L"\langle a^\dagger a \rangle" L"\langle a^\dagger a \rangle ~~(\textrm{steady})"], |
| 172 | + linewidth = 3, |
| 173 | + linestyle = [:solid :dash], |
| 174 | +) |
| 175 | +Plots.plot(p1, p2, layout = [1, 1]) |
| 176 | +Plots.xaxis!("t") |
| 177 | +``` |
| 178 | + |
| 179 | +## Power spectrum |
| 180 | + |
| 181 | +```{julia} |
| 182 | +ω_list = 1:0.01:3 |
| 183 | +psd_H = PowerSpectrum(M_Heom, steady_H, a, ω_list) |
| 184 | +
|
| 185 | +Plots.plot(ω_list, psd_H, linewidth = 3) |
| 186 | +Plots.xaxis!(L"\omega") |
| 187 | +``` |
| 188 | + |
| 189 | +## Compare with Master Eq. approach |
| 190 | + |
| 191 | +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 |
| 192 | + |
| 193 | +```{julia} |
| 194 | +## Drude_Lorentzian spectral density |
| 195 | +Drude_Lorentz(ω, Γ, W) = 4 * Γ * W * ω / ((ω)^2 + (W)^2) |
| 196 | +
|
| 197 | +## Bose-Einstein distribution |
| 198 | +n_b(ω, kT) = 1 / (exp(ω / kT) - 1) |
| 199 | +
|
| 200 | +## build the jump operators |
| 201 | +jump_op = |
| 202 | + [sqrt(Drude_Lorentz(ωc, Γ, W) * (n_b(ωc, kT) + 1)) * a, sqrt(Drude_Lorentz(ωc, Γ, W) * (n_b(ωc, kT))) * a', J_pump]; |
| 203 | +
|
| 204 | +## construct the HEOMLS matrix for master equation |
| 205 | +M_master = M_S(H_s) |
| 206 | +M_master = addBosonDissipator(M_master, jump_op) |
| 207 | +
|
| 208 | +## time evolution |
| 209 | +sol_M = HEOMsolve(M_master, ψ0, t_list; e_ops = [σz, a' * a]); |
| 210 | +
|
| 211 | +## steady state |
| 212 | +steady_M = steadystate(M_master); |
| 213 | +
|
| 214 | +## expectation value of σz |
| 215 | +σz_evo_M = real(sol_M.expect[1, :]) |
| 216 | +σz_steady_M = expect(σz, steady_M) |
| 217 | +
|
| 218 | +## average photon number |
| 219 | +np_evo_M = real(sol_M.expect[2, :]) |
| 220 | +np_steady_M = expect(a' * a, steady_M) |
| 221 | +
|
| 222 | +p1 = Plots.plot( |
| 223 | + t_list, |
| 224 | + [σz_evo_M, ones(length(t_list)) .* σz_steady_M], |
| 225 | + label = [L"\langle \sigma_z \rangle" L"\langle \sigma_z \rangle ~~(\textrm{steady})"], |
| 226 | + linewidth = 3, |
| 227 | + linestyle = [:solid :dash], |
| 228 | +) |
| 229 | +p2 = Plots.plot( |
| 230 | + t_list, |
| 231 | + [np_evo_M, ones(length(t_list)) .* np_steady_M], |
| 232 | + label = [L"\langle a^\dagger a \rangle" L"\langle a^\dagger a \rangle ~~(\textrm{steady})"], |
| 233 | + linewidth = 3, |
| 234 | + linestyle = [:solid :dash], |
| 235 | +) |
| 236 | +Plots.plot(p1, p2, layout = [1, 1]) |
| 237 | +Plots.xaxis!("t") |
| 238 | +``` |
| 239 | + |
| 240 | +We can also calculate the power spectrum |
| 241 | + |
| 242 | +```{julia} |
| 243 | +ω_list = 1:0.01:3 |
| 244 | +psd_M = PowerSpectrum(M_master, steady_M, a, ω_list) |
| 245 | +
|
| 246 | +Plots.plot(ω_list, psd_M, linewidth = 3) |
| 247 | +Plots.xaxis!(L"\omega") |
| 248 | +``` |
| 249 | + |
| 250 | +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. |
| 251 | + |
| 252 | +## Version Information |
| 253 | +```{julia} |
| 254 | +HierarchicalEOM.versioninfo() |
| 255 | +``` |
0 commit comments