Solving Problems with Time-dependent Hamiltonians
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, mesolve, etc., are all capable of handling time-dependent Hamiltonians and collapse operators.
QuantumToolbox.jl uses QuantumObjectEvolution (or the abbreviated synonym: QobjEvo) to represent time-dependent quantum operators, and its data field (attribute) takes advantage of SciMLOperators.jl.
To generate a time-dependent operator QobjEvo, one needs to specify a time-independent Qobj together with a time-dependent coefficient function with the form coef(p, t), namely
coef(p, t) = sin(t)
H_t = QobjEvo(sigmax(), coef)
Quantum Object Evo.: type=Operator() dims=[2] size=(2, 2) ishermitian=true isconstant=false
ScalarOperator(0.0 + 0.0im) * MatrixOperator(2 × 2)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 front of t) in the declaration of coef. We will describe how to use the parameter p in the section Using parameters.
The QobjEvo can also be generated by specifying many pairs of time-independent Qobj and time-dependent coefficient function. For instance, we will look at a case with the total Hamiltonian
The following code sets up this problem by using a Tuple which contains many time-dependent pairs (H_j, f_j), namely
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)
Quantum Object Evo.: type=Operator() dims=[2] size=(2, 2) ishermitian=true isconstant=false
(DiagonalOperator(2 × 2) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(2 × 2) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(2 × 2) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(2 × 2))This is equivalent by generating the QobjEvo separately and + (or sum) them up:
H_t == H0 + QobjEvo(H1, f1) + QobjEvo(H2, f2) + QobjEvo(H3, f3)trueMost solvers will accept any format that could be made into a QobjEvo for the Hamiltonian. All of the followings are equivalent:
sol = mesolve(H_t, ...)
sol = mesolve(H_tuple, ...)Collapse operators c_ops only accept a list where each element is either a Qobj or a QobjEvo. For example, in the following call:
γ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 will see 2 collapse operators:
As an example, we will look at a case with a time-dependent Hamiltonian of the form
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><e|
u_e = tensor(qeye(N), u * e') # |u><e|
u_u = tensor(qeye(N), u * u') # |u><u|
g_g = tensor(qeye(N), g * g') # |g><g|
a = tensor(destroy(N), qeye(3))
# Hamiltonian
g = 5 # coupling strength
H0 = -g * (g_e' * a + a' * g_e)
H1 = u_e' + u_e
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) * g_e, # Use Rb branching ratio of 4/9 e -> 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, ntraj = 100, 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 is identical to that shown in the examples, the mcsolve however will be noticeably off, suggesting we should increase the number of trajectories ntraj = 100 for this example.
In addition, we can also consider the decay of a simple Harmonic oscillator with time-varying decay rate
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)
QuantumObjectEvolution as a time-dependent quantum system, as its main functionality creates a QuantumObject at a specified time
coef(p, t) = sin(π * t)
Ht = QobjEvo(sigmaz(), coef)
Ht(0.25) # t = 0.25
Quantum Object: type=Operator() dims=[2] size=(2, 2) ishermitian=true
2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
0.707107+0.0im ⋅
⋅ -0.707107+0.0imQuantumObjectEvolution shares a lot of properties with the QuantumObject:
Ht.dataScalarOperator(0.7071067811865475 + 0.0im) * MatrixOperator(2 × 2)Ht.typeOperator()Ht.dims1-element StaticArraysCore.SVector{1, Int64} with indices SOneTo(1):
2The properties of a QuantumObjectEvolution can also be retrieved using several functions with inputting QuantumObjectEvolution:
size(Ht)(2, 2)shape(Ht) # synonym of size(Ht)(2, 2)length(Ht)4eltype(Ht) # element typeComplexF64 (alias for Complex{Float64})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(Ht)
println(issymmetric(Ht)) # symmetric
println(isposdef(Ht)) # positive definite (and Hermitian)false
false
true
false
false
false
false
true
true
true
falseQobjEvo follow the same mathematical operations rules as Qobj. They can be added, subtracted and multiplied with scalar, Qobj and QobjEvo. They also support adjoint, transpose, and conj methods, and can be used for SuperOperator transformation:
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)
Quantum Object Evo.: type=SuperOperator() dims=[2] size=(4, 4) ishermitian=true isconstant=false
(ScalarOperator(0.0 + 0.0im) * MatrixOperator(4 × 4) + (ScalarOperator(-1) * ScalarOperator(0.0 - 0.0im)) * MatrixOperator(4 × 4) + (ScalarOperator(0.0 - 0.0im) * ScalarOperator(0.0 + 0.0im)) * MatrixOperator(4 × 4))Or equivalently:
L2 = liouvillian(Ht, [c_op])
Quantum Object Evo.: type=SuperOperator() dims=[2] size=(4, 4) ishermitian=true isconstant=false
(ScalarOperator(0.0 + 0.0im) * MatrixOperator(4 × 4) + (ScalarOperator(0.0 - 0.0im) * ScalarOperator(0.0 + 0.0im)) * MatrixOperator(4 × 4))t = rand()
L1(t) == L2(t)trueOptimization 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 or lindblad_dissipator to generate time-dependent SuperOperator because these functions have been optimized.
Using parameters
Until now, the coefficients were only functions of time t. In the definition of QuantumToolbox.jl allows you to accomplish this by adding extra parameters p to coefficient functions that make the QobjEvo. For instance, instead of explicitly writing 9 for 5 for p as a Vector or NamedTuple:
# 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)
Quantum Object: type=Operator() dims=[2] size=(2, 2) ishermitian=true
2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
8.6471+0.0im ⋅
⋅ -8.6471+0.0im# 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)
Quantum Object: type=Operator() dims=[2] size=(2, 2) ishermitian=true
2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
8.6471+0.0im ⋅
⋅ -8.6471+0.0imt = rand()
H_t_v(p_v, t) == H_t_n(p_n, t)trueCustom 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:
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))Solution of time evolution
(return code: Success)
--------------------------
num_states = 1
num_expect = 1
ODE alg.: OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false))
abstol = 1.0e-8
reltol = 1.0e-6Similarly, with p::NamedTuple:
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))Solution of time evolution
(return code: Success)
--------------------------
num_states = 1
num_expect = 1
ODE alg.: OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false))
abstol = 1.0e-8
reltol = 1.0e-6sol1.expect == sol2.expecttrue