@@ -10,8 +10,7 @@ engine: julia
1010## Package Imports
1111``` {julia}
1212using HierarchicalEOM
13- using LaTeXStrings
14- import Plots
13+ using CairoMakie
1514```
1615
1716## Introduction
@@ -53,22 +52,26 @@ N = 3 ## system cavity Hilbert space cutoff
5352ωA = 2
5453ωc = 2
5554g = 0.1
56- ## operators
55+
56+ # operators
5757a_c = destroy(N)
5858I_c = qeye(N)
5959σz_A = sigmaz()
6060σm_A = sigmam()
6161I_A = qeye(2)
62- ## operators in tensor-space
62+
63+ # operators in tensor-space
6364a = tensor(a_c, I_A)
6465σz = tensor(I_c, σz_A)
6566σm = tensor(I_c, σm_A)
66- ## Hamiltonian
67+
68+ # Hamiltonian
6769H_A = 0.5 * ωA * σz
6870H_c = ωc * a' * a
6971H_int = g * (a' * σm + a * σm')
7072H_s = H_A + H_c + H_int
71- ## initial state
73+
74+ # initial state
7275ψ0 = tensor(basis(N, 0), basis(2, 0));
7376```
7477
@@ -94,13 +97,19 @@ Before incorporating the correlation function into the HEOMLS matrix, it is esse
9497tlist_test = 0:0.1:10;
9598Bath_test = Boson_DrudeLorentz_Pade(a + a', Γ, W, kT, 1000);
9699Ct = correlation_function(Bath, tlist_test);
97- Ct2 = correlation_function(Bath_test, tlist_test);
98- Plots.plot(tlist_test, real(Ct), label = "N=20 (real part )", linestyle = :dash, linewidth = 3)
99- Plots.plot!(tlist_test, real(Ct2), label = "N=1000 (real part)", linestyle = :solid, linewidth = 3)
100- Plots.plot!(tlist_test, imag(Ct), label = "N=20 (imag part)", linestyle = :dash, linewidth = 3)
101- Plots.plot!(tlist_test, imag(Ct2), label = "N=1000 (imag part)", linestyle = :solid, linewidth = 3)
102- Plots.xaxis!("t")
103- Plots.yaxis!("C(t)")
100+ Ct2 = correlation_function(Bath_test, tlist_test)
101+
102+ # plot
103+ fig = Figure(size = (500, 350))
104+ ax = Axis(fig[1, 1], xlabel = L"t", ylabel = L"C(t)")
105+ lines!(ax, tlist_test, real(Ct2), label = L"$N=1000$ (real part)", linestyle = :solid)
106+ lines!(ax, tlist_test, real(Ct), label = L"$N=20$ (real part)", linestyle = :dash)
107+ lines!(ax, tlist_test, imag(Ct2), label = L"$N=1000$ (imag part)", linestyle = :solid)
108+ lines!(ax, tlist_test, imag(Ct), label = L"$N=20$ (imag part)", linestyle = :dash)
109+
110+ axislegend(ax, position = :rt)
111+
112+ fig
104113```
105114
106115## Construct HEOMLS matrix
@@ -149,83 +158,102 @@ np_steady_H = expect(a' * a, steady_H)
149158plot results
150159
151160``` {julia}
152- p1 = Plots.plot(
153- t_list,
154- [σz_evo_H, ones(length(t_list)) .* σz_steady_H],
155- label = [L"\langle \sigma_z \rangle" L"\langle \sigma_z \rangle ~~(\textrm{steady})"],
156- linewidth = 3,
157- linestyle = [:solid :dash],
158- )
159- p2 = Plots.plot(
160- t_list,
161- [np_evo_H, ones(length(t_list)) .* np_steady_H],
162- label = [L"\langle a^\dagger a \rangle" L"\langle a^\dagger a \rangle ~~(\textrm{steady})"],
163- linewidth = 3,
164- linestyle = [:solid :dash],
165- )
166- Plots.plot(p1, p2, layout = [1, 1])
167- Plots.xaxis!("t")
161+ fig = Figure(size = (600, 350))
162+
163+ ax1 = Axis(fig[1, 1], xlabel = L"t")
164+ lines!(ax1, t_list, σz_evo_H, label = L"\langle \sigma_z \rangle", linestyle = :solid)
165+ lines!(ax1, t_list, ones(length(t_list)) .* σz_steady_H, label = L"\langle \sigma_z \rangle ~~(\textrm{steady})", linestyle = :dash)
166+ axislegend(ax1, position = :rt)
167+
168+ ax2 = Axis(fig[2, 1], xlabel = L"t")
169+ lines!(ax2, t_list, np_evo_H, label = L"\langle a^\dagger a \rangle", linestyle = :solid)
170+ lines!(ax2, t_list, ones(length(t_list)) .* np_steady_H, label = L"\langle a^\dagger a \rangle ~~(\textrm{steady})", linestyle = :dash)
171+ axislegend(ax2, position = :rt)
172+
173+ fig
168174```
169175
170176## Power spectrum
171177
172178``` {julia}
173179ω_list = 1:0.01:3
174180psd_H = PowerSpectrum(M_Heom, steady_H, a, ω_list)
175- Plots.plot(ω_list, psd_H, linewidth = 3)
176- Plots.xaxis!(L"\omega")
181+
182+ # plot
183+ fig = Figure(size = (500, 350))
184+ ax = Axis(fig[1, 1], xlabel = L"\omega")
185+ lines!(ax, ω_list, psd_H)
186+
187+ fig
177188```
178189
179190## Compare with Master Eq. approach
180191
181192The 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
182193
183194``` {julia}
184- ## Drude_Lorentzian spectral density
195+ # Drude_Lorentzian spectral density
185196Drude_Lorentz(ω, Γ, W) = 4 * Γ * W * ω / ((ω)^2 + (W)^2)
186- ## Bose-Einstein distribution
197+
198+ # Bose-Einstein distribution
187199n_b(ω, kT) = 1 / (exp(ω / kT) - 1)
188- ## build the jump operators
189- jump_op =
190- [sqrt(Drude_Lorentz(ωc, Γ, W) * (n_b(ωc, kT) + 1)) * a, sqrt(Drude_Lorentz(ωc, Γ, W) * (n_b(ωc, kT))) * a', J_pump];
191- ## construct the HEOMLS matrix for master equation
200+
201+ # build the jump operators
202+ jump_op = [
203+ sqrt(Drude_Lorentz(ωc, Γ, W) * (n_b(ωc, kT) + 1)) * a,
204+ sqrt(Drude_Lorentz(ωc, Γ, W) * (n_b(ωc, kT))) * a',
205+ J_pump
206+ ];
207+
208+ # construct the HEOMLS matrix for master equation
192209M_master = M_S(H_s)
193210M_master = addBosonDissipator(M_master, jump_op)
194- ## time evolution
211+
212+ # time evolution
195213sol_M = HEOMsolve(M_master, ψ0, t_list; e_ops = [σz, a' * a]);
196- ## steady state
214+
215+ # steady state
197216steady_M = steadystate(M_master);
198- ## expectation value of σz
217+
218+ # expectation value of σz
199219σz_evo_M = real(sol_M.expect[1, :])
200220σz_steady_M = expect(σz, steady_M)
201- ## average photon number
221+
222+ # average photon number
202223np_evo_M = real(sol_M.expect[2, :])
203- np_steady_M = expect(a' * a, steady_M)
204- p1 = Plots.plot(
205- t_list,
206- [σz_evo_M, ones(length(t_list)) .* σz_steady_M],
207- label = [L"\langle \sigma_z \rangle" L"\langle \sigma_z \rangle ~~(\textrm{steady})"],
208- linewidth = 3,
209- linestyle = [:solid :dash],
210- )
211- p2 = Plots.plot(
212- t_list,
213- [np_evo_M, ones(length(t_list)) .* np_steady_M],
214- label = [L"\langle a^\dagger a \rangle" L"\langle a^\dagger a \rangle ~~(\textrm{steady})"],
215- linewidth = 3,
216- linestyle = [:solid :dash],
217- )
218- Plots.plot(p1, p2, layout = [1, 1])
219- Plots.xaxis!("t")
224+ np_steady_M = expect(a' * a, steady_M);
225+ ```
226+
227+ plot results
228+
229+ ``` {julia}
230+ fig = Figure(size = (600, 350))
231+
232+ ax1 = Axis(fig[1, 1], xlabel = L"t")
233+ lines!(ax1, t_list, σz_evo_M, label = L"\langle \sigma_z \rangle", linestyle = :solid)
234+ lines!(ax1, t_list, ones(length(t_list)) .* σz_steady_M, label = L"\langle \sigma_z \rangle ~~(\textrm{steady})", linestyle = :dash)
235+ axislegend(ax1, position = :rt)
236+
237+ ax2 = Axis(fig[2, 1], xlabel = L"t")
238+ lines!(ax2, t_list, np_evo_M, label = L"\langle a^\dagger a \rangle", linestyle = :solid)
239+ lines!(ax2, t_list, ones(length(t_list)) .* np_steady_M, label = L"\langle a^\dagger a \rangle ~~(\textrm{steady})", linestyle = :dash)
240+ axislegend(ax2, position = :rt)
241+
242+ fig
220243```
221244
222245We can also calculate the power spectrum
223246
224247``` {julia}
225248ω_list = 1:0.01:3
226249psd_M = PowerSpectrum(M_master, steady_M, a, ω_list)
227- Plots.plot(ω_list, psd_M, linewidth = 3)
228- Plots.xaxis!(L"\omega")
250+
251+ # plot
252+ fig = Figure(size = (500, 350))
253+ ax = Axis(fig[1, 1], xlabel = L"\omega")
254+ lines!(ax, ω_list, psd_M)
255+
256+ fig
229257```
230258
231259Due 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.
0 commit comments