add QobjEvo documentation page

Author: ytdHuang

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

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)
 

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.
+```@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
+ge = tensor(qeye(N), g * e')  # |g><e|
+ue = tensor(qeye(N), u * e')  # |u><e|
+uu = tensor(qeye(N), u * u')  # |u><u|
+gg = tensor(qeye(N), g * g')  # |g><g|
+a = tensor(destroy(N), qeye(3))
+
+# Hamiltonian
+g = 5  # coupling strength
+H0 = -g * (ge' * a + a' * ge)
+H1 = ue' + ue
+f1(p, t) = 9 * exp(-(t / 5)^2)
+H_t = QobjEvo(
+    (
+        H0,
+        (H1, f1)
+    )
+)
+
+# Build collapse operators
+κ = 1.5 # Cavity decay rate
+γ0 = 6   # Atomic decay rate
+c_ops = [
+    sqrt(κ) * a,
+    sqrt(4 * γ0 / 9) * ge, # Use Rb branching ratio of 4/9 e -> g
+    sqrt(5 * γ0 / 9) * ue  # 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,
+    uu,
+    gg
+]
+
+# 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
+```