Открытое квантовое моделирование системы

Я долгое время работал над моделированием открытой квантовой системы с помощью Уравнение Линдблада. Гамильтониан следующий:

Однако к гамильтониану добавляются две другие матрицы. Один из них имеет все диагональные члены, равные -33.3333i, и все остальное ноль. Другая - матрица с третьим диагональным членом, равным -0.033333i.

Уравнение Линдблада таково:

где L_i - матрицы (в списке: [L1, L2, L3, L4, L5, L6, L7]). Матрица для L_i является просто матрицей 7x7 со всеми нулями, кроме L_ (ii) = 1. H - полный гамильтониан, $\rho$ - матрица плотности, а $\gamma$ - постоянная, равная $2\pi kT/\hbar*E_{R}/(\hbar\omega_{c})$ , где T - температура, k - постоянная Больцмана и $\hbar$ = $h/2\pi$ , где h - постоянная Планка. (Обратите внимание, что гамма находится в естественных единицах)

Следующие коды решают уравнение Линдблада, поэтому вычисляют матрицу плотности. Затем он вычисляет и вычисляет это в зависимости от времени:

Это называется популяцией сайта 3. называется bra и называется ket. Оба являются векторами. См. Код для их определения в этом случае.

Вот код:

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 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()

В этом коде используется модуль Python, известный как qutip. Он имеет встроенный решатель уравнения Линдбада с использованием scipy.integrate.odeint.

В настоящее время эта программа отображает следующее:

Однако предел популяции сайта 3 должен быть равен 0. Следовательно, он должен медленно уменьшаться до нуля. Особенно при t = 75, должно начаться снижение.

Этот код работает, но не дает правильного результата, как я объяснил. Итак, почему же он не дает правильного результата? Что-то не так с моим кодом?

Я просмотрел свой код, каждая строка, чтобы увидеть, совпадает ли он с моделью, которую я использую. Они отлично сочетаются. Проблема должна быть в коде, а не в физике.

Я сделал некоторые отладочные подсказки, и все матрицы и гамма правильны. Тем не менее, я все же подозреваю что-то в матрице trap. Причина, по которой я так думаю, состоит в том, что сюжет выглядит как динамика системы без матрицы trap, может быть что-то не так с определением матрицы ловушки, которую я не замечаю?

Обратите внимание, что для выполнения кода требуется несколько минут. Будьте терпеливы при запуске кода!

Ответ 1

(NB: это ответ, который я надеюсь в смысле программирования, но не физический ответ.)

Я запускал ваши симуляции самостоятельно, не используя qutip, и получаю практически те же результаты. Значит, хорошая новость (может быть?:)) заключается в том, что это не проблема с вашим программированием, а проблема физики или, по крайней мере, проблема "выбора параметров". Вот мои результаты:

И ноутбук с моей работой, параметры все те же, что и у вас, кроме другого временного масштаба (поясняется ниже). Я использую тот же метод интеграции, что и qutip, но не сам qutip: Notebook Link.

Несколько примечаний:

Когда вы выполняете from math import *, вы импортируете функцию gamma, а затем вы указываете переменную gamma, это вызывает проблемы для меня, и вы можете быть осторожны в этом будущем.
Когда вы умножаете свои операторы linblad на \gamma, а не на сумму, они будут отображаться дважды в главном уравнении, поэтому вы действительно указываете \gamma^2 здесь. Это влияет на шкалу времени.
<3|rho(t)|3> - это только третий элемент диагональной матрицы, вам действительно не нужен внутренний продукт.

Несколько вещей, которые нужно проверить на стороне физики/параметра.

Из документа, который вы связали,
- \gamma определенно 100/3?
- \kappa_3 определенно 0,1/3 и все остальные 0?
- является исходным состоянием определенно всей совокупности в состоянии 0?
Я не обновляюсь с моделями переноса энергии, но гамильтониан здесь является неэрмитовым, а нетривиальные мнимые части (хотя и небольшие) порождаются на диагональной матрице плотности. Убедитесь, что вы точно понимаете, как и почему эти ребята используют эту модель, как мне кажется странным!