diff --git a/docs/src/users_guide/QuantumObject/QuantumObject.md b/docs/src/users_guide/QuantumObject/QuantumObject.md index 01a86914f..58377dcd7 100644 --- a/docs/src/users_guide/QuantumObject/QuantumObject.md +++ b/docs/src/users_guide/QuantumObject/QuantumObject.md @@ -174,6 +174,7 @@ println(isoper(a)) # operator println(isoperket(a)) # operator-ket println(isoperbra(a)) # operator-bra println(issuper(a)) # super operator +println(isconstant(a)) # time-independent or not println(ishermitian(a)) # Hermitian println(isherm(a)) # synonym of ishermitian(a) println(issymmetric(a)) # symmetric diff --git a/docs/src/users_guide/steadystate.md b/docs/src/users_guide/steadystate.md index 2c53efb03..52e64419c 100644 --- a/docs/src/users_guide/steadystate.md +++ b/docs/src/users_guide/steadystate.md @@ -104,7 +104,7 @@ exp_mc = real(sol_mc.expect[1, :]) sol_me = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops, progress_bar=false) exp_me = real(sol_me.expect[1, :]) -# plot the results +# plot by CairoMakie.jl fig = Figure(size = (500, 350)) ax = Axis(fig[1, 1], title = L"Decay of Fock state $|10\rangle$ in a thermal environment with $\langle n\rangle=2$", diff --git a/docs/src/users_guide/time_evolution/intro.md b/docs/src/users_guide/time_evolution/intro.md index 2aaf3eeb0..3b9ccf985 100644 --- a/docs/src/users_guide/time_evolution/intro.md +++ b/docs/src/users_guide/time_evolution/intro.md @@ -19,6 +19,9 @@ - [Monte-Carlo Solver](@ref doc-TE:Monte-Carlo-Solver) - [Stochastic Solver](@ref doc-TE:Stochastic-Solver) - [Solving Problems with Time-dependent Hamiltonians](@ref doc-TE:Solving-Problems-with-Time-dependent-Hamiltonians) + - [Generate QobjEvo](@ref doc-TE:Generate-QobjEvo) + - [QobjEvo fields (attributes)](@ref doc-TE:QobjEvo-fields-(attributes)) + - [Using parameters](@ref doc-TE:Using-parameters) # [Introduction](@id doc-TE:Introduction) diff --git a/docs/src/users_guide/time_evolution/mesolve.md b/docs/src/users_guide/time_evolution/mesolve.md index 6c06c1568..9ad924b96 100644 --- a/docs/src/users_guide/time_evolution/mesolve.md +++ b/docs/src/users_guide/time_evolution/mesolve.md @@ -2,6 +2,9 @@ ```@setup mesolve using QuantumToolbox + +using CairoMakie +CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) ``` ## [Von Neumann equation](@id doc-TE:Von-Neumann-equation) @@ -123,13 +126,11 @@ sol = mesolve(H, ψ0, tlist, c_ops, e_ops = [sigmaz(), sigmay()]) We can therefore plot the expectation values: ```@example mesolve -using CairoMakie -CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) - times = sol.times expt_z = real(sol.expect[1,:]) expt_y = real(sol.expect[2,:]) +# plot by CairoMakie.jl fig = Figure(size = (500, 350)) ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values") lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid) @@ -166,7 +167,7 @@ e_ops = [a' * a] sol = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops) Num = real(sol.expect[1, :]) -# plot the results +# plot by CairoMakie.jl fig = Figure(size = (500, 350)) ax = Axis(fig[1, 1], title = L"Decay of Fock state $|10\rangle$ in a thermal environment with $\langle n\rangle=2$", @@ -213,6 +214,7 @@ times = sol.times N_atom = real(sol.expect[1,:]) N_cavity = real(sol.expect[2,:]) +# plot by CairoMakie.jl fig = Figure(size = (500, 350)) ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values") lines!(ax, times, N_atom, label = "atom excitation probability", linestyle = :solid) diff --git a/docs/src/users_guide/time_evolution/sesolve.md b/docs/src/users_guide/time_evolution/sesolve.md index f41df647d..d3ba8a21c 100644 --- a/docs/src/users_guide/time_evolution/sesolve.md +++ b/docs/src/users_guide/time_evolution/sesolve.md @@ -23,6 +23,9 @@ The Schrödinger equation, which governs the time-evolution of closed quantum sy ```@setup sesolve using QuantumToolbox + +using CairoMakie +CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) ``` For example, the time evolution of a quantum spin-``\frac{1}{2}`` system (initialized in spin-``\uparrow``) with tunneling rate ``0.1`` is calculated, and the expectation values of the Pauli-Z operator ``\hat{\sigma}_z`` is also evaluated, with the following code @@ -68,12 +71,10 @@ print(size(expt)) We can therefore plot the expectation values: ```@example sesolve -using CairoMakie -CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) - expt_z = real(expt[1,:]) expt_y = real(expt[2,:]) +# plot by CairoMakie.jl fig = Figure(size = (500, 350)) ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values") lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid) diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md index 455c94991..5294ef214 100644 --- a/docs/src/users_guide/time_evolution/solution.md +++ b/docs/src/users_guide/time_evolution/solution.md @@ -2,6 +2,9 @@ ```@setup TE-solution using QuantumToolbox + +using CairoMakie +CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) ``` ## [Solution](@id doc-TE:Solution) @@ -61,9 +64,7 @@ nothing # hide we can plot the resulting expectation values: ```@example TE-solution -using CairoMakie -CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) - +# plot by CairoMakie.jl fig = Figure(size = (500, 350)) ax = Axis(fig[1, 1], xlabel = L"t") lines!(ax, times, expt1, label = L"\langle 0 | \rho(t) | 0 \rangle") diff --git a/docs/src/users_guide/time_evolution/time_dependent.md b/docs/src/users_guide/time_evolution/time_dependent.md index f9346e884..ba9c0042c 100644 --- a/docs/src/users_guide/time_evolution/time_dependent.md +++ b/docs/src/users_guide/time_evolution/time_dependent.md @@ -1,3 +1,353 @@ # [Solving Problems with Time-dependent Hamiltonians](@id doc-TE:Solving-Problems-with-Time-dependent-Hamiltonians) -This page is still under construction, please visit [API](@ref doc-API) first. \ No newline at end of file +```@setup QobjEvo +using QuantumToolbox + +using CairoMakie +CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) +``` + +## [Generate QobjEvo](@id doc-TE:Generate-QobjEvo) + +In the previous examples of solving time evolution, we assumed that the systems under consideration were described by time-independent Hamiltonians. However, many systems have explicit time dependence in either the Hamiltonian, or the collapse operators (`c_ops`) describing coupling to the environment, and sometimes both components might depend on time. The time evolution solvers such as [`sesolve`](@ref), [`mesolve`](@ref), etc., are all capable of handling time-dependent Hamiltonians and collapse operators. + +`QuantumToolbox.jl` uses [`QuantumObjectEvolution`](@ref) (or the abbreviated synonym: [`QobjEvo`](@ref)) to represent time-dependent quantum operators, and its `data` field (attribute) takes advantage of [`SciMLOperators.jl`](https://docs.sciml.ai/SciMLOperators/stable/). + +To generate a time-dependent operator [`QobjEvo`](@ref), one needs to specify a time-independent [`Qobj`](@ref) together with a time-dependent coefficient function with the form `coef(p, t)`, namely + +```@example QobjEvo +coef(p, t) = sin(t) + +H_t = QobjEvo(sigmax(), coef) +``` + +!!! warning "The inputs of coefficient function" + Please note that although we didn't use the argument `p` in the definition of `coef`, we still need to put a dummy input `p` in the declaration of `coef`. We will describe how to use the parameter `p` in the section [Using parameters](@ref doc-TE:Using-parameters). + +The [`QobjEvo`](@ref) can also be generated by specifying many pairs of time-independent [`Qobj`](@ref) and time-dependent coefficient function. For instance, we will look at a case with the total Hamiltonian ``\hat{H}(t)`` can be separate into time-independent part (``\hat{H}_0``) and a summation of many time-dependent operators, which takes the form: + +```math +\hat{H}(t) = \hat{H}_0 + \sum_j f_j(t) \hat{H}_j. +``` + +The following code sets up this problem by using a `Tuple` which contains many time-dependent pairs `(H_j, f_j)`, namely + +```@example QobjEvo +H0 = qeye(2) + +H1 = sigmax() +f1(p, t) = sin(t) + +H2 = sigmay() +f2(p, t) = cos(t) + +H3 = sigmaz() +f3(p, t) = 9 * exp(-(t / 5)^2) + +H_tuple = ( + H0, + (H1, f1), + (H2, f2), + (H3, f3), +) + +H_t = QobjEvo(H_tuple) +``` + +This is equivalent by generating the [`QobjEvo`](@ref) separately and `+` (or `sum`) them up: + +```@example QobjEvo +H_t == H0 + QobjEvo(H1, f1) + QobjEvo(H2, f2) + QobjEvo(H3, f3) +``` + +Most solvers will accept any format that could be made into a [`QobjEvo`](@ref) for the Hamiltonian. All of the followings are equivalent: + +```julia +sol = mesolve(H_t, ...) +sol = mesolve(H_tuple, ...) +``` + +Collapse operators `c_ops` only accept a list where each element is either a [`Qobj`](@ref) or a [`QobjEvo`](@ref). For example, in the following call: + +```julia +γ1 = sqrt(0.1) +γ2(p, t) = sqrt(sin(t)) +c_ops = ( + γ1 * sigmaz(), + QobjEvo(sigmax() + sigmay(), γ2) +) + +sol = mesolve(H_t, ..., c_ops, ...) +``` + +[`mesolve`](@ref) will see `2` collapse operators: ``\sqrt{0.1} \hat{\sigma}_z`` and ``\sqrt{\sin(t)} (\hat{\sigma}_x + \hat{\sigma}_y)``. + +As an example, we will look at a case with a time-dependent Hamiltonian of the form ``\hat{H}(t) = \hat{H}_0 + f_1(t) \hat{H}_1``, where ``f_1(t) = A \exp \left[ -(t / \sigma)^2 \right]`` is the time-dependent driving strength. The following code sets up the problem + +```@example QobjEvo +N = 2 # Set where to truncate Fock state for cavity + +# basis states for Atom +u = basis(3, 0) # u state +e = basis(3, 1) # excited state +g = basis(3, 2) # ground state + +# operators +g_e = tensor(qeye(N), g * e') # |g> g + sqrt(5 * γ0 / 9) * u_e # 5/9 e -> u +] + +tlist = LinRange(0, 4, 200) # Define time vector +ψ0 = tensor(basis(N, 0), u) # Define initial state + +# Build observables +e_ops = [ + a' * a, + u_u, + g_g +] + +# solve dynamics +exp_me = mesolve(H_t, ψ0, tlist, c_ops, e_ops = e_ops; progress_bar = Val(false)).expect +exp_mc = mcsolve(H_t, ψ0, tlist, c_ops, e_ops = e_ops; progress_bar = Val(false)).expect + +# plot by CairoMakie.jl +fig = Figure(size = (500, 350)) +ax = Axis(fig[1, 1], xlabel = L"Time $t$", ylabel = "expectation value") +lines!(ax, tlist, real(exp_me[1,:]), label = L"$\hat{a}^\dagger \hat{a}$ (mesolve)", color = :red, linestyle = :solid) +lines!(ax, tlist, real(exp_mc[1,:]), label = L"$\hat{a}^\dagger \hat{a}$ (mcsolve)", color = :red, linestyle = :dash) +lines!(ax, tlist, real(exp_me[2,:]), label = L"$|u \rangle\langle u|$ (mesolve)", color = :green, linestyle = :solid) +lines!(ax, tlist, real(exp_mc[2,:]), label = L"$|u \rangle\langle u|$ (mcsolve)", color = :green, linestyle = :dash) +lines!(ax, tlist, real(exp_me[3,:]), label = L"$|g \rangle\langle g|$ (mesolve)", color = :blue, linestyle = :solid) +lines!(ax, tlist, real(exp_mc[3,:]), label = L"$|g \rangle\langle g|$ (mcsolve)", color = :blue, linestyle = :dash) + +axislegend(ax, position = :rc) + +fig +``` + +The result from [`mesolve`](@ref) is identical to that shown in the examples, the [`mcsolve`](@ref) however will be noticeably off, suggesting we should increase the number of trajectories `ntraj` for this example. + +In addition, we can also consider the decay of a simple Harmonic oscillator with time-varying decay rate ``\gamma_1(t)`` + +```@example QobjEvo +N = 10 # number of basis states +a = destroy(N) +H = a' * a + +ψ0 = basis(N, 9) # initial state + +γ0 = 0.5 +γ1(p, t) = sqrt(γ0 * exp(-t)) +c_ops_ti = [sqrt(γ0) * a] # time-independent collapse term +c_ops_td = [QobjEvo(a, γ1)] # time-dependent collapse term + +e_ops = [a' * a] + +tlist = LinRange(0, 10, 100) +exp_ti = mesolve(H, ψ0, tlist, c_ops_ti, e_ops = e_ops; progress_bar = Val(false)).expect[1,:] +exp_td = mesolve(H, ψ0, tlist, c_ops_td, e_ops = e_ops; progress_bar = Val(false)).expect[1,:] + +# plot by CairoMakie.jl +fig = Figure(size = (500, 350)) +ax = Axis(fig[1, 1], xlabel = L"Time $t$", ylabel = L"\langle \hat{a}^\dagger \hat{a} \rangle") +lines!(ax, tlist, real(exp_ti), label = L"\gamma_0", linestyle = :solid) +lines!(ax, tlist, real(exp_td), label = L"\gamma_1(t) = \gamma_0 e^{-t}", linestyle = :dash) + +axislegend(ax, position = :rt) + +fig +``` + +## [QobjEvo fields (attributes)](@id doc-TE:QobjEvo-fields-(attributes)) + +[`QuantumObjectEvolution`](@ref) as a time-dependent quantum system, as its main functionality creates a [`QuantumObject`](@ref) at a specified time ``t``: + +```@example QobjEvo +coef(p, t) = sin(π * t) +Ht = QobjEvo(sigmaz(), coef) + +Ht(0.25) +``` + +[`QuantumObjectEvolution`](@ref) shares a lot of properties with the [`QuantumObject`](@ref): + +```@example QobjEvo +Ht.data +``` + +```@example QobjEvo +Ht.type +``` + +```@example QobjEvo +Ht.dims +``` + +The properties of a [`QuantumObjectEvolution`](@ref) can also be retrieved using several functions with inputting [`QuantumObjectEvolution`](@ref): + +```@example QobjEvo +size(Ht) +``` + +```@example QobjEvo +shape(Ht) # synonym of size(H) +``` + +```@example QobjEvo +length(Ht) +``` + +```@example QobjEvo +eltype(Ht) # element type +``` + +```@example QobjEvo +println(isket(Ht)) # ket +println(isbra(Ht)) # bra +println(isoper(Ht)) # operator +println(isoperket(Ht)) # operator-ket +println(isoperbra(Ht)) # operator-bra +println(issuper(Ht)) # super operator +println(isconstant(Ht)) # time-independent or not +println(ishermitian(Ht)) # Hermitian +println(isherm(Ht)) # synonym of ishermitian(a) +println(issymmetric(Ht)) # symmetric +println(isposdef(Ht)) # positive definite (and Hermitian) +``` + +[`QobjEvo`](@ref) follow the same mathematical operations rules as [`Qobj`](@ref). They can be added, subtracted and multiplied with scalar, [`Qobj`](@ref) and [`QobjEvo`](@ref). They also support [`adjoint`](@ref), [`transpose`](@ref), and [`conj`](@ref) methods, and can be used for [`SuperOperator`](@ref) transformation: + +```@example QobjEvo +coef(p, t) = sin(π * t) +Ht = QobjEvo(sigmaz(), coef) + +coef2(p, t) = exp(-t) +c_op = QobjEvo(destroy(2), coef2) + +L1 = -1im * (spre(Ht) - spost(Ht')) +L1 += lindblad_dissipator(c_op) +``` + +Or equivalently: +```@example QobjEvo +L2 = liouvillian(Ht, [c_op]) +``` + +```@example QobjEvo +t = rand() +L1(t) == L2(t) +``` + +!!! note "Optimization for superoperator transformation" + Although the value of `L1` and `L2` here is equivalent, one can observe that the structure of `L2.data` is much simpler than `L1.data`, which means that it will be more efficient to solve the time evolution using `L2` instead of `L1`. Therefore, we recommend to use [`liouvillian`](@ref) or [`lindblad_dissipator`](@ref) to generate time-dependent [`SuperOperator`](@ref) because these functions have been optimized. + +## [Using parameters](@id doc-TE:Using-parameters) + +Until now, the coefficients were only functions of time `t`. In the definition of ``f_1(t) = A \exp \left[ -(t / \sigma)^2 \right]`` in the previous example, the driving amplitude ``A`` and width of the gaussian driving term ``\sigma`` were hardcoded with their numerical values. This is fine for problems that are specialized, or that we only want to run once. However, in many cases, we would like to study the same problem with a range of parameters and don't have to worry about manually changing the values on each run. `QuantumToolbox.jl` allows you to accomplish this by adding extra parameters `p` to coefficient functions that make the [`QobjEvo`](@ref). For instance, instead of explicitly writing `9` for ``A`` and `5` for ``\sigma``, we can either specify them with `p` as a `Vector` or `NamedTuple`: + +```@example QobjEvo +# specify p as a Vector +f1_v(p, t) = p[1] * exp(-(t / p[2])^2) + +H_t_v = QobjEvo(sigmaz(), f1_v) + +p_v = [9, 5] # 1st element represents A; 2nd element represents σ +t = 1 +H_t_v(p_v, t) +``` + +```@example QobjEvo +# specify p as a NamedTuple +f1_n(p, t) = p.A * exp(-(t / p.σ)^2) + +H_t_n = QobjEvo(sigmaz(), f1_n) + +p_n = (A = 9, σ = 5) +t = 1 +H_t_n(p_n, t) +``` + +```@example QobjEvo +t = rand() +H_t_v(p_v, t) == H_t_n(p_n, t) +``` + +!!! note "Custom structure of parameter" + For more advanced usage, any custom `struct` of parameter `p` can be used. + +When solving time evolutions, the solvers take an single keyword argument `params`, which will be directly passed to input `p` of the coefficient functions to build the time dependent Hamiltonian and collapse operators. For example, with `p::Vector`: + +```@example QobjEvo +f1(p, t) = p[1] * cos(p[2] * t) +f2(p, t) = p[3] * sin(p[4] * t) +γ(p, t) = sqrt(p[5] * exp(-p[6] * t)) +p_total = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6] + +H_t = sigmaz() + QobjEvo(sigmax(), f1) + QobjEvo(sigmay(), f2) + +c_ops = [ + QobjEvo(destroy(2), γ) +] + +e_ops = [sigmaz()] + +ψ0 = basis(2, 0) +tlist = LinRange(0, 10, 20) +sol1 = mesolve(H_t, ψ0, tlist, c_ops, params = p_total, e_ops = e_ops; progress_bar = Val(false)) +``` + +Similarly, with `p::NamedTuple`: + +```@example QobjEvo +f1(p, t) = p.A1 * cos(p.ω1 * t) +f2(p, t) = p.A2 * sin(p.ω2 * t) +γ(p, t) = sqrt(p.A3 * exp(-p.γ0 * t)) +p_total = ( + A1 = 0.1, + ω1 = 0.2, + A2 = 0.3, + ω2 = 0.4, + A3 = 0.5, + γ0 = 0.6 +) + +H_t = sigmaz() + QobjEvo(sigmax(), f1) + QobjEvo(sigmay(), f2) + +c_ops = [ + QobjEvo(destroy(2), γ) +] + +e_ops = [sigmaz()] + +ψ0 = basis(2, 0) +tlist = LinRange(0, 10, 20) +sol2 = mesolve(H_t, ψ0, tlist, c_ops, params = p_total, e_ops = e_ops; progress_bar = Val(false)) +``` + +```@example QobjEvo +sol1.expect == sol2.expect +```