Computing propagators

Using a solver to compute a propagator

Many solvers accept an identity matrix as the initial state. When such QuantumObject is passed as the initial state, the propagator is computed. This can be used to compute a propagator for sesolve, mesolve, brmesolve, etc. For example,

ϵ0 = 1.0 * 2π
Ω = 0.8 * 2π
Δt = π/5
H = (ϵ0/2) * sigmaz() + (Ω/2) * sigmax()

Quantum Object:   type=Operator()   dims=([2], [2])   size=(2, 2)   ishermitian=true
2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 4 stored entries:
 3.14159+0.0im   2.51327+0.0im
 2.51327+0.0im  -3.14159+0.0im

The propagator for the Schrodinger equation with a small time step can be used as

For a Hamiltonian H, the corresponding propagator can be generated by using sesolve with the initial state given as qeye_like(H), which is an identity Operator (same as H.type), namely

U_se = sesolve(H, qeye_like(H), [0, Δt]; progress_bar = Val(false)).states[end]

Quantum Object:   type=Operator()   dims=([2], [2])   size=(2, 2)   ishermitian=false
2×2 Matrix{ComplexF64}:
    -0.817499-0.449725im  5.77598e-17-0.35978im
 -5.77598e-17-0.35978im     -0.817499+0.449725im

Thus, the state at time   can be computed by applying the propagator repeatedly times on the initial state :

n = 10
ψ0 = basis(2, 0)

ψt = U_se^n * ψ0

Quantum Object:   type=Ket()   dims=([2], [1])   size=(2,)
2-element Vector{ComplexF64}:
 0.9893929007501819 - 0.113432308543722im
                0.0 - 0.09074584683497754im

Similarly, the propagator for Lindblad master equation with a small time step can be used as

where is the vectorized density operator.

For a Liouvillian superoperator L, the corresponding propagator can be generated by using mesolve with the initial state given as qeye_like(L), which is an identity SuperOperator (same as L.type), namely

L = liouvillian(H)
U_me = mesolve(L, qeye_like(L), [0, Δt]; progress_bar = Val(false)).states[end]

Quantum Object:   type=SuperOperator()   dims=([2], [2])   size=(4, 4)
4×4 Matrix{ComplexF64}:
 0.870558+0.0im       0.161802+0.29412im      …   0.129442+0.0im
 0.161802+0.29412im   0.466053-0.7353im          -0.161802-0.29412im
 0.161802-0.29412im   0.129442-2.62112e-17im     -0.161802+0.29412im
 0.129442+0.0im      -0.161802-0.29412im          0.870558+0.0im

Thus, the state at time   can be computed by applying the propagator repeatedly times on the initial vectorized density operator :

n = 10
ρ0_vec = mat2vec(ket2dm(basis(2, 0))) # type=OperatorKet()

ρt_vec = U_me^n * ρ0_vec

Quantum Object:   type=OperatorKet()   dims=([2], [1])   size=(4,)
4-element Vector{ComplexF64}:
   0.9917633321003378 - 9.379593720852414e-18im
 0.010295834874578129 - 0.08978364199243058im
 0.010295834874578115 + 0.08978364199243065im
 0.008236667899661237 + 1.2297230239048356e-17im