Schrödinger Equation Solver

Unitary evolution

The dynamics of a closed (pure) quantum system is governed by the Schrödinger equation

$$i\hslash \frac{\partial }{\partial t}\mathrm{\Psi }(\overrightarrow{x},t)=\hat{H}\mathrm{\Psi }(\overrightarrow{x},t),$$

where $\mathrm{\Psi }(\overrightarrow{x},t)$ is the wave function, $\hat{H}$ is the Hamiltonian, and $\hslash$ is reduced Planck constant. In general, the Schrödinger equation is a partial differential equation (PDE) where both $\mathrm{\Psi }$ and $\hat{H}$ are functions of space $\overrightarrow{x}$ and time $t$. For computational purposes it is useful to expand the PDE in a set of basis functions that span the Hilbert space of the Hamiltonian, and to write the equation in matrix and vector form, namely

$$i\hslash \frac{d}{dt}|\psi (t)⟩=\hat{H}|\psi (t)⟩,$$

where $|\psi (t)⟩$ is the state vector, and the Hamiltonian $\hat{H}$ is now under matrix representation. This matrix equation can, in principle, be solved by diagonalizing the Hamiltonian matrix $\hat{H}$. In practice, however, it is difficult to perform this diagonalization unless the size of the Hilbert space (dimension of the matrix $\hat{H}$) is small. Analytically, it is a formidable task to calculate the dynamics for systems with more than two states. If, in addition, we consider dissipation due to the inevitable interaction with a surrounding environment, the computational complexity grows even larger, and we have to resort to numerical calculations in all realistic situations. This illustrates the importance of numerical calculations in describing the dynamics of open quantum systems, and the need for efficient and accessible tools for this task.

The Schrödinger equation, which governs the time-evolution of closed quantum systems, is defined by its Hamiltonian and state vector. In the previous sections, Manipulating States and Operators and Tensor Products and Partial Traces, we showed how Hamiltonians and state vectors are constructed in QuantumToolbox.jl. Given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state vector $|\psi (0)⟩$ using the QuantumToolbox time evolution problem sesolveProblem or directly call the function sesolve. It evolves the state vector $|\psi (t)⟩$ and evaluates the expectation values for a set of operators e_ops at each given time points, using an ordinary differential equation solver provided by the powerful julia package DifferentialEquation.jl.

Example: Spin dynamics

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

H = 2 * π * 0.1 * sigmax()
ψ0 = basis(2, 0) # spin-up
tlist = LinRange(0.0, 10.0, 20)

prob = sesolveProblem(H, ψ0, tlist, e_ops = [sigmaz()], progress_bar = Val(false))
sol = sesolve(prob)

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-6

Note

Here, we generate the time evolution problem by sesolveProblem first and then put it into the function sesolve. One can also directly call sesolve, which we also demonstrates in below.

The function returns TimeEvolutionSol, as described in the previous section Time Evolution Solutions. The attribute expect in solution is a list of expectation values for the operator(s) that are passed to the e_ops keyword argument.

sol.expect

1×20 Matrix{ComplexF64}:
 1.0+0.0im  0.789141+0.0im  0.245485+0.0im  …  0.789141+0.0im  1.0+0.0im

Passing multiple operators to e_ops as a Vector results in the expectation values for each operators at each time points.

tlist = LinRange(0.0, 10.0, 100)
sol = sesolve(H, ψ0, tlist, e_ops = [sigmaz(), sigmay()])

Solution of time evolution
(return code: Success)
--------------------------
num_states = 1
num_expect = 2
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-6

Note

Here, we call sesolve directly instead of pre-defining sesolveProblem first (as shown previously).

println(size(sol.times)) # time points corresponds to stored expectation values
println(size(sol.times_states)) # time points corresponds to stored states

(100,)
(1,)

expt = sol.expect
print(size(expt))

(2, 100)

We can therefore plot the expectation values:

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, tlist, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
lines!(ax, tlist, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)

axislegend(ax, position = :rb)

fig

If the keyword argument e_ops is not specified (or given as an empty Vector), the sesolve and functions return a TimeEvolutionSol that contains a list of state vectors which corresponds to the time points specified in tlist:

tlist = [0, 10]
sol = sesolve(H, ψ0, tlist) # or specify: e_ops = []

println(size(sol.times))
println(size(sol.times_states))

[sesolve] 100%|███████████████████████████| Time: 0:00:00 (15.97 μs/it)
(2,)
(2,)

sol.states

2-element Vector{QuantumObject{Ket, Dimensions{1, Tuple{Space}}, Vector{ComplexF64}}}:
 
Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 0.0 + 0.0im
 
Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 0.9999999543220187 + 0.0im
                0.0 - 5.6377384923229734e-8im

This is because the states will be saved depend on the keyword argument saveat in kwargs. If e_ops is empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state).

One can also specify e_ops and saveat separately:

tlist = [0, 5, 10]
sol = sesolve(H, ψ0, tlist, e_ops = [sigmay()], saveat = tlist)

Solution of time evolution
(return code: Success)
--------------------------
num_states = 3
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-6

print(size(sol.times))

(3,)

sol.expect

1×3 Matrix{ComplexF64}:
 0.0+0.0im  -7.18972e-8+0.0im  -1.12755e-7+0.0im

print(size(sol.times_states))

(3,)

sol.states

3-element Vector{QuantumObject{Ket, Dimensions{1, Tuple{Space}}, Vector{ComplexF64}}}:
 
Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 0.0 + 0.0im
 
Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 -0.9999999984043708 + 0.0im
                 0.0 + 3.594862168420766e-8im
 
Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 0.9999999543220187 + 0.0im
                0.0 - 5.6377384923229734e-8im