Lindblad Master Equation Solver

Von Neumann equation

While the evolution of the state vector in a closed quantum system is deterministic (as we discussed in the previous section: Schrödinger Equation Solver), open quantum systems are stochastic in nature. The effect of an environment on the system of interest is to induce stochastic transitions between energy levels, and to introduce uncertainty in the phase difference between states of the system. The state of an open quantum system is therefore described in terms of ensemble averaged states using the density matrix formalism. A density matrix $\hat{\rho }$ describes a probability distribution of quantum states $|{\psi }_n⟩$ in a matrix representation, namely

$$\hat{\rho }=\sum _n{p}_n|{\psi }_n⟩⟨{\psi }_n|,$$

where $p_n$ is the classical probability that the system is in the quantum state $|{\psi }_n⟩$. The time evolution of a density matrix $\hat{\rho }$ is the topic of the remaining portions of this section.

The time evolution of the density matrix $\hat{\rho }(t)$ under closed system dynamics is governed by the von Neumann equation:

$$\frac{d}{dt}\hat{\rho }(t)=-\frac{i}{\hslash }[\hat{H},\hat{\rho }(t)],$$

where $[\cdot ,\cdot ]$ represents the commutator. The above equation is equivalent to the Schrödinger equation described in the previous section under the density matrix formalism.

In QuantumToolbox, given a Hamiltonian, we can calculate the unitary (non-dissipative) time-evolution of an arbitrary initial state using the QuantumToolbox time evolution problem mesolveProblem or directly call the function mesolve. It evolves the density matrix $\hat{\rho }(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.

H = 0.5 * sigmax()
state0 = basis(2, 0)         # state vector
state0 = ket2dm(basis(2, 0)) # density matrix
tlist = LinRange(0.0, 10.0, 20)

sol = mesolve(H, state0, tlist, e_ops = [sigmaz()])

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

Type of initial state

The initial state state0 here can be given as a state vector $|\psi (0)⟩$ (in the type of Ket) or a density matrix $\hat{\rho }(0)$ (in the type of Operator). If it is given as a Ket, it will be transformed to density matrix $\hat{\rho }(0)=|\psi (0)⟩⟨\psi (0)|$ internally in mesolve.

The function returns TimeEvolutionSol, as described in the previous section Time Evolution Solutions. The stored states will always be in the type of Operator (density matrix).

sol.states

1-element Vector{QuantumObject{OperatorQuantumObject, Dimensions{1, Tuple{Space}}, Matrix{ComplexF64}}}:
 
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 Matrix{ComplexF64}:
 0.0804643+0.0im            0.0-0.272011im
       0.0+0.272011im  0.919536+0.0im

Here, only the final state is stored 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 = mesolve(H, state0, 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

sol.expect

1×3 Matrix{ComplexF64}:
 0.0+0.0im  0.958924+0.0im  0.544021+0.0im

sol.states

3-element Vector{QuantumObject{OperatorQuantumObject, Dimensions{1, Tuple{Space}}, Matrix{ComplexF64}}}:
 
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
 
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 Matrix{ComplexF64}:
 0.641831+0.0im            0.0-0.479462im
      0.0+0.479462im  0.358169+0.0im
 
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 Matrix{ComplexF64}:
 0.0804643+0.0im            0.0-0.272011im
       0.0+0.272011im  0.919536+0.0im

The Lindblad master equation

The standard approach for deriving the equations of motion for a system interacting with its environment is to expand the scope of the system to include the environment. The combined quantum system is then closed, and its evolution is also governed by the von Neumann equation

$$\frac{d}{dt}{\hat{\rho }}_{\text{tot}}(t)=-\frac{i}{\hslash }[{\hat{H}}_{\text{tot}},{\hat{\rho }}_{\text{tot}}(t)].$$

Here, the total Hamiltonian

$${\hat{H}}_{\text{tot}}={\hat{H}}_{\text{sys}}+{\hat{H}}_{\text{env}}+{\hat{H}}_{\text{int}},$$

includes the original system Hamiltonian ${\hat{H}}_{\text{sys}}$, the Hamiltonian for the environment ${\hat{H}}_{\text{env}}$, and a term representing the interaction between the system and its environment ${\hat{H}}_{\text{int}}$. Since we are only interested in the dynamics of the system, we can, perform a partial trace over the environmental degrees of freedom, and thereby obtain a master equation for the motion of the original system density matrix ${\hat{\rho }}_{\text{sys}}(t)={\text{Tr}}_{\text{env}}[{\hat{\rho }}_{\text{tot}}(t)]$. The most general trace-preserving and completely positive form of this evolution is the Lindblad master equation for the reduced density matrix, namely

$$\frac{d}{dt}{\hat{\rho }}_{\text{sys}}(t)=-\frac{i}{\hslash }[{\hat{H}}_{\text{sys}},{\hat{\rho }}_{\text{sys}}(t)]+\sum _n{\hat{C}}_n{\hat{\rho }}_{\text{sys}}(t){\hat{C}}_n^{†}-\frac{1}{2}{\hat{C}}_n^{†}{\hat{C}}_n{\hat{\rho }}_{\text{sys}}(t)-\frac{1}{2}{\hat{\rho }}_{\text{sys}}(t){\hat{C}}_n^{†}{\hat{C}}_n$$

where ${\hat{C}}_n\equiv \sqrt{{\gamma }_n}{\hat{A}}_n$ are the collapse operators, ${\hat{A}}_n$ are the operators acting on the system in ${\hat{H}}_{\text{int}}$ which describes the system-environment interaction, and ${\gamma }_n$ are the corresponding rates. The derivation of Lindblad master equation may be found in several sources, and will not be reproduced here. Instead, we emphasize the approximations that are required to arrive at the above Lindblad master equation from physical arguments, and hence perform a calculation in QuantumToolbox:

• Separability: At $t=0$, there are no correlations between the system and environment, such that the total density matrix can be written as a tensor product, namely ${\hat{\rho }}_{\text{tot}}(0)={\hat{\rho }}_{\text{sys}}(0)\otimes {\hat{\rho }}_{\text{env}}(0)$.

• Born approximation: Requires: (i) the state of the environment does not significantly change as a result of the interaction with the system; (ii) the system and the environment remain separable throughout the evolution. These assumptions are justified if the interaction is weak, and if the environment is much larger than the system. In summary, ${\hat{\rho }}_{\text{tot}}(t)\approx {\hat{\rho }}_{\text{sys}}(t)\otimes {\hat{\rho }}_{\text{env}}(0)$.

• Markov approximation: The time-scale of decay for the environment ${\tau }_{\text{env}}$ is much shorter than the smallest time-scale of the system dynamics, i.e., ${\tau }_{\text{sys}}\gg {\tau }_{\text{env}}$. This approximation is often deemed a “short-memory environment” as it requires the environmental correlation functions decay in a fast time-scale compared to those of the system.

• Secular approximation: Stipulates that elements in the master equation corresponding to transition frequencies satisfy $|{\omega }_{ab}-{\omega }_{cd}|\ll 1/{\tau }_{\text{sys}}$, i.e., all fast rotating terms in the interaction picture can be neglected. It also ignores terms that lead to a small renormalization of the system energy levels. This approximation is not strictly necessary for all master-equation formalisms (e.g., the Block-Redfield master equation), but it is required for arriving at the Lindblad form in the above equation which is used in mesolve.

For systems with environments satisfying the conditions outlined above, the Lindblad master equation governs the time-evolution of the system density matrix, giving an ensemble average of the system dynamics. In order to ensure that these approximations are not violated, it is important that the decay rates ${\gamma }_n$ be smaller than the minimum energy splitting in the system Hamiltonian. Situations that demand special attention therefore include, for example, systems strongly coupled to their environment, and systems with degenerate or nearly degenerate energy levels.

What is new in the master equation compared to the Schrödinger equation (or von Neumann equation) are processes that describe dissipation in the quantum system due to its interaction with an environment. For example, evolution that includes incoherent processes such as relaxation and dephasing. These environmental interactions are defined by the operators ${\hat{A}}_n$ through which the system couples to the environment, and rates ${\gamma }_n$ that describe the strength of the processes.

In QuantumToolbox, the function mesolve can also be used for solving the master equation. This is done by passing a list of collapse operators (c_ops) as the fourth argument of the mesolve function in order to define the dissipation processes of the Lindblad master equation. As we mentioned above, each collapse operator ${\hat{C}}_n$ is the product of $\sqrt{{\gamma }_n}$ (the square root of the rate) and ${\hat{A}}_n$ (operator which describes the dissipation process).

Furthermore, QuantumToolbox solves the master equation in the SuperOperator formalism. That is, a Liouvillian SuperOperator will be generated internally in mesolve by the input system Hamiltonian ${\hat{H}}_{\text{sys}}$ and the collapse operators ${\hat{C}}_n$. One can also generate the Liouvillian SuperOperator manually for special purposes, and pass it as the first argument of the mesolve function. To do so, it is useful to read the section Superoperators and Vectorized Operators, and also the docstrings of the following functions:

• spre

• spost

• sprepost

• liouvillian

• lindblad_dissipator

Example: Dissipative Spin dynamics

Using the example with the dynamics of spin-$\frac{1}{2}$ from the previous section (Schrödinger Equation Solver), we can easily add a relaxation process (describing the dissipation of energy from the spin to the environment), by adding [sqrt(γ) * sigmax()] in the fourth parameter of the mesolve function.

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

γ = 0.05
c_ops = [sqrt(γ) * sigmax()]

sol = mesolve(H, ψ0, tlist, c_ops, 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

We can therefore plot the expectation values:

times = sol.times
expt_z = real(sol.expect[1,:])
expt_y = real(sol.expect[2,:])

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)

axislegend(ax, position = :rt)

fig

Example: Harmonic oscillator in thermal bath

Consider a harmonic oscillator (single-mode cavity) couples to a thermal bath. If the single-mode cavity initially is in a 10-photon fock state, the dynamics is calculated with the following code:

# Define parameters
N = 20  # number of basis states to consider
a = destroy(N)
H = a' * a
ψ0 = fock(N, 10)  # initial state
κ = 0.1  # coupling to oscillator
n_th = 2 # temperature with average of 2 excitations
tlist = LinRange(0, 50, 100)

# collapse operators
c_ops = [
    sqrt(κ * (n_th + 1)) * a, # emission
    sqrt(κ *  n_th     ) * a' # absorption
]

# find expectation value for particle number
e_ops = [a' * a]

sol = mesolve(H, ψ0, tlist, c_ops, e_ops=e_ops)
Num = real(sol.expect[1, :])

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1],
    title = L"Decay of Fock state $|10\rangle$ in a thermal environment with $\langle n\rangle=2$",
    xlabel = "Time",
    ylabel = "Number of excitations",
)
lines!(ax, tlist, Num)

fig

Example: Two-level atom coupled to dissipative single-mode cavity

Consider a two-level atom coupled to a dissipative single-mode cavity through a dipole-type interaction, which supports a coherent exchange of quanta between the two systems. If the atom initially is in its ground state and the cavity in a 5-photon fock state, the dynamics is calculated with the following code:

Generate Liouviilian superoperator manually

In this example, we demonstrate how to generate the Liouvillian SuperOperator manually and pass it as the first argument of the mesolve function.

# two-level atom
σm = tensor(destroy(2), qeye(10))
H_a = 2 * π * σm' * σm

# dissipative single-mode cavity
γ = 0.1 # dissipation rate
a  = tensor(qeye(2), destroy(10))
H_c = 2 * π * a' * a
c_ops = [sqrt(γ) * a]

# coupling between two-level atom and single-mode cavity
g = 0.25 # coupling strength
H_I = 2 * π * g * (σm * a' + σm' * a)

ψ0 = tensor(basis(2,0), fock(10, 5)) # initial state
tlist = LinRange(0.0, 10.0, 200)

# generate Liouvillian superoperator manually
L = liouvillian(H_a + H_c + H_I, c_ops)
sol = mesolve(L, ψ0, tlist, e_ops=[σm' * σm, a' * a])

times = sol.times

# expectation value of Number operator
N_atom   = real(sol.expect[1,:])
N_cavity = real(sol.expect[2,:])

# plot by CairoMakie.jl
fig = Figure(size = (500, 350))
ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
lines!(ax, times, N_atom, label = "atom excitation probability", linestyle = :solid)
lines!(ax, times, N_cavity, label = "cavity photon number", linestyle = :dash)

axislegend(ax, position = :rt)

fig