13

很长一段时间以来，我一直致力于使用Lindblad 方程对开放量子系统进行建模。哈密​​顿量如下：

哈密​​顿量

然而，另外两个矩阵被添加到哈密顿量中。其中一个的所有对角项都等于 -33.3333i，其他所有项都为零。另一个是第三个对角项等于 -0.033333i 的矩阵。

Lindblad 方程是这样的：

林布拉德方程

其中 L_i 是矩阵（在列表中：[L1,L2,L3,L4,L5,L6,L7]）。L_i 的矩阵只是一个 7x7 矩阵，除了 L_(ii)=1 之外全为零。H 是总哈密顿量，$\rho$是密度矩阵，并且$\伽马$是一个常数，$2\pi kT/\hbar*E_{R}/(\hbar\omega_{c})$其中 T 是温度，k 是玻尔兹曼常数，并且$\hbar$ = $h/2\pi$，其中 h 是普朗克常数。（请注意，伽玛是自然单位）

以下代码求解 Lindblad 方程，从而计算密度矩阵。然后它计算并绘制它与时间的关系：

人口

这被称为站点 3 人口。胸罩被称为胸罩，酮被称为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'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()

此代码使用称为qutip的 Python 模块。它有一个使用 scipy.integrate.odeint 的内置 Lindblad 方程求​​解器。

目前，该程序显示以下内容：

结果

但是，站点 3 人口的限制应该是 0。因此，它应该慢慢减少到零。特别是到 t=75 时，应该开始下降。

该代码运行，但没有产生我解释的正确结果。那么现在，为什么它不能产生正确的结果呢？我的代码有问题吗？

我查看了我的代码，每一行看它是否与我正在使用的模型匹配。他们完美匹配。问题必须出在代码中，而不是物理上。

我做了一些调试提示，所有的矩阵和伽玛都是正确的。然而，我仍然怀疑trap矩阵中的某些东西。我这么认为的原因是因为情节看起来像没有矩阵的系统的动力学trap，我没有注意到的陷阱矩阵的定义是否有问题？

请注意，代码需要几分钟才能运行。运行代码时请耐心等待！

4

1 回答 1

6

（注意：这是我希望在编程意义上的答案，但不是物理答案。）

我独立运行了你的模拟，而不是使用qutip，我得到了基本相同的结果。所以好消息（也许？:)）这不是您的编程问题，而是物理问题或至少是“选择参数”问题。这是我的结果： 在此处输入图像描述

和我工作的笔记本，参数都和你的一样，除了不同的时间尺度（解释如下）。我使用的是与 qutip 本身相同的集成方法，qutip但不是 qutip 本身：Notebook Link。

几点注意事项：

  1. 当您执行from math import *导入函数gamma，然后命名变量gamma时，这给我带来了问题，您将来可能要小心。

  2. 当您将 linblad 运算符乘以\gamma而不是总和时，它们将在主方程中出现两次，因此您实际上是\gamma^2在此处指定。这会影响时间尺度。

  3. <3|rho(t)|3>只是第三个对角矩阵元素，这里不需要内积。

在物理/参数方面要检查一些事情。

  1. 从您链接的论文中，
    • 肯定是\Gamma100/3？
    • 绝对是\kappa_30.1/3 而其他都是 0？
    • 初始状态肯定是所有人口都处于 0 状态吗？
  2. 我不是最新的能量转移模型，但这里的哈密顿是非厄米特和非平凡的虚数部分（尽管仍然很小）是在密度矩阵对角线上生成的。确保您完全了解这些人如何以及为什么使用此模型，因为这对我来说似乎很奇怪！
于 2016-04-15T12:51:54.147 回答