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]   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 $\hat{U}$ for the Schrodinger equation with a small time step $\mathrm{\Delta }t$ can be used as

$$|\psi (t+\mathrm{\Delta }t)⟩=\hat{U}(\mathrm{\Delta }t)|\psi (t)⟩$$

For a Hamiltonian H, the corresponding propagator $\hat{U}(\mathrm{\Delta }t)$ 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]   size=(2, 2)   ishermitian=false
2×2 Matrix{ComplexF64}:
   -0.817499-0.449725im  -5.64934e-17-0.35978im
 5.64934e-17-0.35978im      -0.817499+0.449725im

Thus, the state at time $t=n\mathrm{\Delta }t$ can be computed by applying the propagator $\hat{U}(\mathrm{\Delta }t)$ repeatedly $n$ times on the initial state $|\psi (0)⟩$:

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

ψt = U_se^n * ψ0

Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
    0.9893929007493156 - 0.113432308543756im
 5.551115123125783e-17 - 0.09074584683500489im

Similarly, the propagator $\mathcal{U}$ for Lindblad master equation with a small time step $\mathrm{\Delta }t$ can be used as

$$|\hat{\rho }(t+\mathrm{\Delta }t)⟩⟩=\mathcal{U}(\mathrm{\Delta }t)|\hat{\rho }(t)⟩⟩,$$

where $|\hat{\rho }(t)⟩⟩$ is the vectorized density operator.

For a Liouvillian superoperator L, the corresponding propagator $\mathcal{U}(\mathrm{\Delta }t)$ 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]   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-5.11789e-17im     -0.161802+0.29412im
 0.129442+0.0im      -0.161802-0.29412im          0.870558+0.0im

Thus, the state at time $t=n\mathrm{\Delta }t$ can be computed by applying the propagator $\mathcal{U}(\mathrm{\Delta }t)$ repeatedly $n$ times on the initial vectorized density operator $|\hat{\rho }(0)⟩⟩$:

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

ρt_vec = U_me^n * ρ0_vec

Quantum Object:   type=OperatorKet()   dims=[2]   size=(4,)
4-element Vector{ComplexF64}:
   0.9917633321003388 + 1.3313997427657217e-17im
  0.01029583487457792 - 0.0897836419924306im
 0.010295834874577928 + 0.08978364199243057im
  0.00823666789966114 + 1.1986845980022588e-17im