Open quantum system modelling

Question

I have been working for a long time now on modelling an open quantum system using the Lindblad Equation. The Hamiltonian is the following:

However, two other matrices are added to the Hamiltonian. One of them has all the diagonal terms equal to -33.3333i and everything else zero. Another is a matrix with the third diagonal term equaling -0.033333i.

The Lindblad Equation is this:

where L_i are matrices (in the list: [L1,L2,L3,L4,L5,L6,L7]). The matrix for L_i is simply a 7x7 matrix with all zeros except L_(ii)=1. H is the total Hamiltonian, is the density matrix, and is a constant equal to where T is the temperature, k is the Boltzmann constant, and , where h is Planck's constant. (Note that gamma is in natural units)

The following codes solves the Lindblad Equation, therefore calculating the density matrix. It then calculates and plots this versus time:

This is known as the site 3 population. is called a bra and is called a ket. Both are vectors. See the code for their definition in this case.

Here is the code:

from qutip import Qobj, Options, mesolve
import numpy as np
import scipy
from math import *
import matplotlib.pyplot as plt

hamiltonian = np.array([
    [215, -104.1, 5.1, -4.3, 4.7, -15.1, -7.8],
    [-104.1, 220.0, 32.6, 7.1, 5.4, 8.3, 0.8],
    [5.1, 32.6, 0.0, -46.8, 1.0, -8.1, 5.1],
    [-4.3, 7.1, -46.8, 125.0, -70.7, -14.7, -61.5],
    [4.7, 5.4, 1.0, -70.7, 450.0, 89.7, -2.5],
    [-15.1, 8.3, -8.1, -14.7, 89.7, 330.0, 32.7],
    [-7.8, 0.8, 5.1, -61.5, -2.5, 32.7, 280.0]
])

recomb = np.zeros((7, 7), dtype=complex)
np.fill_diagonal(recomb, 33.33333333)
recomb = recomb * -1j
trap = np.zeros((7, 7), complex)
trap[2][2] = -0.033333333333j
hamiltonian = recomb + trap + hamiltonian
H = Qobj(hamiltonian)

# Note the extra .0 on the end to convert to float
gamma = (2 * pi) * (296 * 0.695) * (35.0 / 150)

L1 = np.array([
    [1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]
])

L2 = np.array([
    [0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]
])

L3 = np.array([
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]
])      

L4 = np.array([
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]
])

L5 = np.array([
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]
])

L6 = np.array([
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0],
    [0, 0, 0, 0, 0, 0, 0]
])

L7 = np.array([
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 1]
])

# Since our gamma variable cannot be directly applied onto
# the Lindblad operator, we must multiply it with
# the collapse operators:  

rho0=Qobj(L1)

L1 = Qobj(gamma * L1)
L2 = Qobj(gamma * L2)
L3 = Qobj(gamma * L3)
L4 = Qobj(gamma * L4)
L5 = Qobj(gamma * L5)
L6 = Qobj(gamma * L6)
L7 = Qobj(gamma * L7)

options = Options(nsteps=1000000, atol=1e-5)

bra3 = [[0, 0, 1, 0, 0, 0, 0]]
bra3q = Qobj(bra3)

ket3 = [[0], [0], [1], [0], [0], [0], [0]]
ket3q = Qobj(ket3)

starttime = 0
# this is effectively just a label - `mesolve` alwasys starts from `rho0` -
# it's just saying what we're going to call the time at t0
endtime = 100
# Arbitrary - this solves with the options above
# (max 1 million iterations to converge - tolerance 1e-10)
num_intermediate_state = 100

state_evaluation_times = np.linspace(
    starttime,
    endtime,
    num_intermediate_state
)

result = mesolve(
    H,
    rho0,
    state_evaluation_times,
    [L1, L2, L3, L4, L5, L6, L7],
    [],
    options=options
)

number_of_interest = bra3q * (result.states * ket3q)

points_to_plot = []
for number in number_of_interest:
    if number == number_of_interest[0]:
        points_to_plot.append(0)
    else:
        points_to_plot.append(number.data.data.real[0])

plt.plot(state_evaluation_times, points_to_plot)
plt.show()
exit()

This code uses a Python module known as qutip. It has a builtin Lindblad Equation solver using scipy.integrate.odeint.

Currently, this program displays this:

However, the limit of the site 3 population is supposed to be 0. Therefore, it should decrease slowly down to zero. Especially by t=75, the decrease should start.

This code runs, but does not produce the correct result as I explained. So now, why doesn't it produce the right result? Is something wrong with my code?

I have looked at my code, each line to see if it matches up with the model I am using. They match up perfectly. The issue must be in the code, not the physics.

I have done some debug prompts, and all the matrices and the gamma is correct. I still, however, suspect something in the trap matrix. The reason I think so is because the plot looks like the dynamics of the system without the trap matrix, Could there be something wrong with the definition of the trap matrix that I am not noticing?

---

Note, the code takes a few minutes to run. Be patient while running the code!

Answer 1

(NB: this is an answer I hope in the programming sense, but not a physics answer.)

I ran your simulations independently, not using qutip, and I get essentially the same results. So the good news (maybe? :)) is that it's not an issue with your programming but a physics problem or at least a "choosing parameters" problem. Here are my results:

And a notebook with my working, the parameters are all the same as yours, apart from a different time-scale (explained below). I am using the same integration method asqutip but not qutip itself: Notebook Link.

A few notes:

1. When you perform from math import * you import a function gamma, and then you name a variable gamma, this caused problems for me and you might want to be careful about that in future.

2. When you multiply your linblad operators by \gamma rather than the sum they will appear twice in the master equation, so you are really specifying \gamma^2 here. This affects the time-scale.

3. <3|rho(t)|3> is just the third diagonal matrix element, you don't really need an inner product here.

A few things to check on the physics / parameter side.

1. From the paper you linked,
   • is \Gamma definitely 100/3?
   • is \kappa_3 definitely 0.1/3 and all others 0?
   • is the initial state definitely all population in 0 state?
2. I am not up to date with energy transfer models but the hamiltonian here is non-hermitian and non-trivial imaginary parts (though still small) are generated on the density matrices diagonal. Make sure you understand exactly how and why these guys are using this model as it seems odd to me!