开放量子系统建模

我一直在使用Lindblad方程对开放量子系统建模工作很长时间.汉密尔顿主义者如下:

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

Lindblad方程是这样的:

其中L_i是矩阵(在列表中:[L1,L2,L3,L4,L5,L6,L7]).L_i的矩阵只是一个7x7矩阵,除了L_(ii)= 1之外全零.H是总哈密顿量, 是密度矩阵,和 是一个等于的常数 其中T是温度,k是玻尔兹曼常数,和 ,其中h是普朗克常数.(请注意,gamma是自然单位)

以下代码解决了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模块.它有一个内置的Lindblad方程求解器,使用scipy.integrate.odeint.

目前,该程序显示如下:

但是,站点3总体的限制应该是0.因此,它应该缓慢降低到零.特别是在t = 75时,应该开始减少.

此代码运行,但没有像我解释的那样产生正确的结果.所以现在,为什么它不能产生正确的结果呢？我的代码有问题吗？

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

我做了一些调试提示,所有的矩阵和伽玛都是正确的.然而,我仍然怀疑trap矩阵中的某些东西.我想是这样的原因是因为情节貌似系统的动态无的trap矩阵,难道还有什么问题,我没有注意到的陷阱矩阵的定义是什么？

注意,代码需要几分钟才能运行.运行代码时要耐心等待!

(注意:这是我希望在编程意义上的答案,但不是物理答案.)

我独立运行你的模拟,而不是使用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状态？

我没有与能量转移模型约会,但这里的哈密顿函数是非埃尔米特和非平凡的虚部(尽管仍然很小)是在密度矩阵对角线上生成的.确保你完全理解这些人使用这个模型的方式和原因,因为这对我来说似乎很奇怪!