|
| 1 | +# Quick Start |
| 2 | + |
| 3 | +### Content |
| 4 | +- [Import HierarchicalEOM.jl](#Import-HierarchicalEOM.jl) |
| 5 | +- [System and Bath](#System-and-Bath) |
| 6 | +- [HEOM Liouvillian superoperator](#HEOM-Liouvillian-superoperator) |
| 7 | +- [Time Evolution](#Time-Evolution) |
| 8 | +- [Stationary State](#Stationary-State) |
| 9 | +- [Reduced Density Operator](#Reduced-Density-Operator) |
| 10 | +- [Expectation Value](#Expectation-Value) |
| 11 | +- [Multiple Baths](#Multiple-Baths) |
| 12 | + |
| 13 | +### Import HierarchicalEOM.jl |
| 14 | + |
| 15 | +Here are the functions in `HierarchicalEOM.jl` that we will use in this tutorial (Quick Start): |
| 16 | + |
| 17 | +```@example quick_start |
| 18 | +import HierarchicalEOM |
| 19 | +import HierarchicalEOM: Boson_DrudeLorentz_Pade, M_Boson, HEOMsolve, getRho, BosonBath |
| 20 | +``` |
| 21 | + |
| 22 | +Note that you can also type `using HierarchicalEOM` to import everything you need in `HierarchicalEOM.jl`. |
| 23 | +To check the versions of dependencies of `HierarchicalEOM.jl`, run the following function |
| 24 | + |
| 25 | +```@example quick_start |
| 26 | +HierarchicalEOM.versioninfo() |
| 27 | +``` |
| 28 | + |
| 29 | +### System and Bath |
| 30 | + |
| 31 | +Let us consider a simple two-level system coupled to a Drude-Lorentz bosonic bath. The system Hamiltonian, ``H_{sys}``, and the bath spectral density, ``J_D``, are |
| 32 | + |
| 33 | +```math |
| 34 | +H_{sys}=\frac{\epsilon \sigma_z}{2} + \frac{\Delta \sigma_x}{2} ~~\text{and} |
| 35 | +``` |
| 36 | + |
| 37 | +```math |
| 38 | +J_{D}(\omega)=\frac{2\lambda W\omega}{W^2+\omega^2}, |
| 39 | +``` |
| 40 | +#### System Hamiltonian and initial state |
| 41 | + |
| 42 | +You must construct system hamiltonian, initial state, and coupling operators by [`QuantumToolbox`](https://github.com/qutip/QuantumToolbox.jl) framework. It provides many useful functions to create arbitrary quantum states and operators which can be combined in all the expected ways. |
| 43 | + |
| 44 | +```@example quick_start |
| 45 | +import QuantumToolbox: Qobj, sigmaz, sigmax, basis, ket2dm, expect, steadystate |
| 46 | +``` |
| 47 | + |
| 48 | +```@example quick_start |
| 49 | +# The system Hamiltonian |
| 50 | +ϵ = 0.5 # energy of 2-level system |
| 51 | +Δ = 1.0 # tunneling term |
| 52 | +
|
| 53 | +Hsys = 0.5 * ϵ * sigmaz() + 0.5 * Δ * sigmax() |
| 54 | +
|
| 55 | +# System initial state |
| 56 | +ρ0 = ket2dm(basis(2, 0)); |
| 57 | +
|
| 58 | +# Define the operators that measure the populations of the two system states: |
| 59 | +P00 = ket2dm(basis(2, 0)) |
| 60 | +P11 = ket2dm(basis(2, 1)) |
| 61 | +
|
| 62 | +# Define the operator that measures the 0, 1 element of density matrix |
| 63 | +# (corresponding to coherence): |
| 64 | +P01 = basis(2, 0) * basis(2, 1)' |
| 65 | +``` |
| 66 | + |
| 67 | +#### Bath Properties |
| 68 | + |
| 69 | +Now, we demonstrate how to describe the bath using the built-in implementation of ``J_D(\omega)`` under Pade expansion by calling [`Boson_DrudeLorentz_Pade`](@ref) |
| 70 | + |
| 71 | +```@example quick_start |
| 72 | +λ = 0.1 # coupling strength |
| 73 | +W = 0.5 # band-width (cut-off frequency) |
| 74 | +kT = 0.5 # the product of the Boltzmann constant k and the absolute temperature T |
| 75 | +
|
| 76 | +Q = sigmaz() # system-bath coupling operator |
| 77 | +
|
| 78 | +N = 2 # Number of expansion terms to retain: |
| 79 | +
|
| 80 | +# Padé expansion: |
| 81 | +bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) |
| 82 | +``` |
| 83 | + |
| 84 | +For other different expansions of the different spectral density correlation functions, please refer to [Bosonic Bath](@ref doc-Bosonic-Bath) and [Fermionic Bath](@ref doc-Fermionic-Bath). |
| 85 | + |
| 86 | +### HEOM Liouvillian superoperator |
| 87 | + |
| 88 | +For bosonic bath, we can construct the HEOM Liouvillian superoperator matrix by calling [`M_Boson`](@ref) |
| 89 | + |
| 90 | +```@example quick_start |
| 91 | +tier = 5 # maximum tier of hierarchy |
| 92 | +L = M_Boson(Hsys, tier, bath) |
| 93 | +``` |
| 94 | + |
| 95 | +To learn more about the HEOM Liouvillian superoperator matrix (including other types: `M_Fermion`, `M_Boson_Fermion`), please refer to [HEOMLS Matrices](@ref doc-HEOMLS-Matrix). |
| 96 | + |
| 97 | +### Time Evolution |
| 98 | + |
| 99 | +Next, we can calculate the time evolution for the entire auxiliary density operators (ADOs) by calling [`HEOMsolve`](@ref) |
| 100 | + |
| 101 | +```@example quick_start |
| 102 | +tlist = 0:0.2:50 |
| 103 | +sol = HEOMsolve(L, ρ0, tlist; e_ops = [P00, P11, P01]) |
| 104 | +``` |
| 105 | + |
| 106 | +To learn more about `HEOMsolve`, please refer to [Time Evolution](@ref doc-Time-Evolution). |
| 107 | + |
| 108 | +### Stationary State |
| 109 | + |
| 110 | +We can also solve the stationary state of the auxiliary density operators (ADOs) by calling [`steadystate`](@ref). |
| 111 | + |
| 112 | +```@example quick_start |
| 113 | +ados_steady = steadystate(L) |
| 114 | +``` |
| 115 | + |
| 116 | +To learn more about `steadystate`, please refer to [Stationary State](@ref doc-Stationary-State). |
| 117 | + |
| 118 | +### Reduced Density Operator |
| 119 | + |
| 120 | +To obtain the reduced density operator, one can either access the first element of auxiliary density operator (`ADOs`) or call [`getRho`](@ref): |
| 121 | + |
| 122 | +```@example quick_start |
| 123 | +# reduce density operator in the final time (`end`) of the evolution |
| 124 | +ados_list = sol.ados |
| 125 | +ρ = ados_list[end][1] # index `1` represents the reduced density operator |
| 126 | +ρ = getRho(ados_list[end]) |
| 127 | +
|
| 128 | +# reduce density operator in stationary state |
| 129 | +ρ = ados_steady[1] |
| 130 | +ρ = getRho(ados_steady) |
| 131 | +``` |
| 132 | + |
| 133 | +One of the great features of `HierarchicalEOM.jl` is that we allow users to not only considering the density operator of the reduced state but also easily take high-order terms into account without struggling in finding the indices (see [Auxiliary Density Operators](@ref doc-ADOs) and [Hierarchy Dictionary](@ref doc-Hierarchy-Dictionary) for more details). |
| 134 | + |
| 135 | +### Expectation Value |
| 136 | + |
| 137 | +We can now compare the results obtained from `HEOMsolve` and `steadystate`: |
| 138 | + |
| 139 | +```@example quick_start |
| 140 | +# for time evolution |
| 141 | +p00_e = real(sol.expect[1, :]) # P00 is the 1st element in e_ops |
| 142 | +p01_e = real(sol.expect[3, :]); # P01 is the 3rd element in e_ops |
| 143 | +
|
| 144 | +# for steady state |
| 145 | +p00_s = expect(P00, ados_steady) |
| 146 | +p01_s = expect(P01, ados_steady); |
| 147 | +``` |
| 148 | + |
| 149 | +### Plot the results |
| 150 | + |
| 151 | +```@example quick_start |
| 152 | +using Plots, LaTeXStrings |
| 153 | +
|
| 154 | +plot(tlist, p00_e, label = L"\textrm{P}_{00}", linecolor = :blue, linestyle = :solid, linewidth = 3, grid = false) |
| 155 | +plot!(tlist, p01_e, label = L"\textrm{P}_{01}", linecolor = :red, linestyle = :solid, linewidth = 3) |
| 156 | +plot!( |
| 157 | + tlist, |
| 158 | + ones(length(tlist)) .* p00_s, |
| 159 | + label = L"\textrm{P}_{00} \textrm{(Steady State)}", |
| 160 | + linecolor = :blue, |
| 161 | + linestyle = :dash, |
| 162 | + linewidth = 3, |
| 163 | +) |
| 164 | +plot!( |
| 165 | + tlist, |
| 166 | + ones(length(tlist)) .* p01_s, |
| 167 | + label = L"\textrm{P}_{01} \textrm{(Steady State)}", |
| 168 | + linecolor = :red, |
| 169 | + linestyle = :dash, |
| 170 | + linewidth = 3, |
| 171 | +) |
| 172 | +
|
| 173 | +xlabel!("time") |
| 174 | +ylabel!("Population") |
| 175 | +``` |
| 176 | + |
| 177 | +### Multiple Baths |
| 178 | + |
| 179 | +`HierarchicalEOM.jl` also supports for system to interact with multiple baths. All you need to do is to provide a list of baths instead of a single bath |
| 180 | + |
| 181 | +```@example quick_start |
| 182 | +# The system Hamiltonian |
| 183 | +Hsys = Qobj([ |
| 184 | + 0.25 1.50 2.50 |
| 185 | + 1.50 0.75 3.50 |
| 186 | + 2.50 3.50 1.25 |
| 187 | +]) |
| 188 | +
|
| 189 | +# System initial state |
| 190 | +ρ0 = ket2dm(basis(3, 0)); |
| 191 | +
|
| 192 | +# Projector for each system state: |
| 193 | +P00 = ket2dm(basis(3, 0)) |
| 194 | +P11 = ket2dm(basis(3, 1)) |
| 195 | +P22 = ket2dm(basis(3, 2)); |
| 196 | +
|
| 197 | +# Construct one bath for each system state: |
| 198 | +# note that `BosonBath[]` make the list created in type: Vector{BosonBath} |
| 199 | +baths = BosonBath[] |
| 200 | +for i in 0:2 |
| 201 | + # system-bath coupling operator: |i><i| |
| 202 | + Q = ket2dm(basis(3, i)) |
| 203 | + push!(baths, Boson_DrudeLorentz_Pade(Q, λ, W, kT, N)) |
| 204 | +end |
| 205 | +
|
| 206 | +L = M_Boson(Hsys, tier, baths) |
| 207 | +
|
| 208 | +tlist = 0:0.025:5 |
| 209 | +sol = HEOMsolve(L, ρ0, tlist; e_ops = [P00, P11, P22]) |
| 210 | +
|
| 211 | +# calculate population for each system state: |
| 212 | +p0 = real(sol.expect[1, :]) |
| 213 | +p1 = real(sol.expect[2, :]) |
| 214 | +p2 = real(sol.expect[3, :]) |
| 215 | +
|
| 216 | +plot(tlist, p0, linewidth = 3, linecolor = "blue", label = L"P_0", grid = false) |
| 217 | +plot!(tlist, p1, linewidth = 3, linecolor = "orange", label = L"P_1") |
| 218 | +plot!(tlist, p2, linewidth = 3, linecolor = :green, label = L"P_2") |
| 219 | +xlabel!("time") |
| 220 | +ylabel!("Population") |
| 221 | +``` |
| 222 | + |
| 223 | +Note that this example can also be found in [qutip documentation](https://qutip.org/docs/latest/guide/heom/bosonic.html). |
0 commit comments