最近研究qutip库中的一些函数，对我这个量子力学基础并不牢靠的人来说，用程序中的矩阵辅助理解是非常有效果的。下面是我自己在玩儿这个程序的时候得到的一些结果以及对物理的理解。
一、基本函数
首先先给出一些qutip库中物理量的基本矩阵形式，
basis(2,0)
上面的结果是[1 0]T，因为只是自己做个记录，没想着有人看，所以后面的列矢量我就都不加“T”了，默认是列矢量就好。
create(2)
上面的结果是[[0 0] [1 0]]
二、激发单粒子
首先在qutip中有一个函数可以生成fock态，basis(N,m)，N的意思是fock态的数量，m的意思是第几个态。
那么我就先尝试了一下：一个很简单的哈密顿量作用在一个很简单的单粒子基态上。哈密顿量是产生算符，他的作用就是让基态粒子激发到激发态。让我们看看时间演化的结果怎么样吧，代码如下，结果随后
N=1 # 粒子数是1，后面会用到
psi = basis(2,0) # 用basis产生一个fock态
ad_list = create(2) # 定义产生算符
a_list = destroy(2) # 定义湮灭算符
hal = ad_list # 哈密顿量
tlist = np.linspace(0, 10, 10) # 演化时间
res = sesolve(hal,psi,tlist) # 解薛定谔方程
states = [s * s.dag() for s in res.states] # 计算密度矩阵
np.set_printoptions(linewidth=10000) # 设置jupyter显示宽度
states = np.real(np.array(states)) # 提取密度矩阵的实部
rounded_states = np.round(states, decimals=4) # 只显示四位小数
np.set_printoptions(suppress=True) # 不使用科学计数法
print(rounded_states) # 输出每一时刻的密度矩阵
上面可以看到，初始的密度矩阵为[[1 0] [0 0]]，随着时间的演化，最后的密度矩阵为[[0.0099 0] [0 0.9901]]，也就是说态从初始的|0>态变成了最后的99%|1>态，相信随着时间的演化最终结果将全变为100%的|1>态。
这一猜想也可以通过直接将哈密顿量作用在初态上得到：
s = ad_list * basis(2,0)
这个矩阵结果手算也可以得到，[[0 0] [1 0]]*[1 0]=[0 1]。
三、激发双粒子
三粒子及以后的规律和之前的一样，可以自己试试。
主要代码和上面的还是差不多的，只是修改一下总粒子数，产生湮灭算符，初态。代码如下，结果随后。
N=2 # 粒子数是2
state = []
for i in range(1, N+1): # 写一个小循环，将每个粒子的fock态写在一起state.append(basis(2,0))
psi = tensor(state) # 用tensor将所有粒子初态直积在一起a_list, ad_list = [], []
for i in range(N): # 用循环定义每个单格点上的产生湮灭算符op_list = [qeye(2)] * Nop_list[i] = destroy(2)a_list.append(tensor(op_list))op_list[i] = create(2)ad_list.append(tensor(op_list))hal = 0
for i in range(N): # 哈密顿量改为每个格点上都有产生算符hal += ad_list[i]tlist = np.linspace(0, 10, 10) # 演化时间
res = sesolve(hal,psi,tlist) # 解薛定谔方程
states = [s * s.dag() for s in res.states] # 计算密度矩阵
np.set_printoptions(linewidth=10000) # 设置jupyter显示宽度
states = np.real(np.array(states)) # 提取密度矩阵的实部
rounded_states = np.round(states, decimals=4) # 只显示四位小数
np.set_printoptions(suppress=True) # 不使用科学计数法
print(rounded_states) # 输出每一时刻的密度矩阵
这就是每个时刻的密度矩阵。我怕我以后看不懂密度矩阵所以多解释一下。就用第一个密度矩阵举例子，能看到除了左上角是1外其他全是0，这是因为我们的初态是|0>×|0>，直积符号不想打，懂得都懂。那么我们密度矩阵的基矢就从上到下依次为|00>，|01>，|10>，|11>，那么左上角的1自然而然就代表初态只有一个|000>态，也就是说初态三个粒子全是基态|0>。那么最后一个矩阵的意思是什么呢，显然，代表1%|01>+1%|10>+98%|11>。
那既然我研究的是动力学，肯定要关注稳态之前的内容，那么用上面的规律我们可以看到态的演化大概经过了一个从100%|00>到25%|00>+25%|01>+25%|10>+25%|11>再到7.5%|01>+7.5%|10>+85%|11>最后到达1%|01>+1%|10>+98%|11>的状态，那么我们可以相信时间无穷长之后一定会变成完美的100%|11>的。这是一个全激发。
四、三粒子相互作用哈密顿量
前面的哈密顿量都是简单的激发每个格点上的粒子，没有考虑格点之间的相互作用，那么接下来我们就改一下哈密顿量，使得其具有相互作用。
这里需要修改一下初态和哈密顿量，哈密顿量里加入相互作用，具体形式是能量向右传播，既然这样那我们必须给予体系一个初始能量，这一点可以通过修改初态来实现。代码如下，结果随后
N=3 # 粒子数是3psi = tensor(basis(2,1),basis(2,0),basis(2,0)) # 制备一个第一个粒子是激发态其他都是基态的初态a_list, ad_list = [], []
for i in range(N): # 用循环定义每个单格点上的产生湮灭算符op_list = [qeye(2)] * Nop_list[i] = destroy(2)a_list.append(tensor(op_list))op_list[i] = create(2)ad_list.append(tensor(op_list))hal = 0
for i in range(N): # 哈密顿量改为每个格点上都有产生算符hal += ad_list[i+1] * a_list[i]tlist = np.linspace(0, 10, 10) # 演化时间
res = sesolve(hal,psi,tlist) # 解薛定谔方程
states = [s * s.dag() for s in res.states] # 计算密度矩阵
np.set_printoptions(linewidth=10000) # 设置jupyter显示宽度
states = np.real(np.array(states)) # 提取密度矩阵的实部
rounded_states = np.round(states, decimals=4) # 只显示四位小数
np.set_printoptions(suppress=True) # 不使用科学计数法
print(rounded_states) # 输出每一时刻的密度矩阵
通过以上结果我们可以看到，初态是|100>，位于第一个密度矩阵的第五行第五列，这是因为密度矩阵的基矢是fock态，排列顺序是|000>，|001>，|010>，|011>，|100>，|101>，|110>，|111>。随着密度矩阵的演化我们可以看到最后到达了96%|001>+3.8%|010>。这非常符合我们的物理期望：第一个格点上是激发态，随着哈密顿量的作用能量依次往右传递，最终到达第三个粒子。这是一个单向能流。
五、三粒子空间对称哈密顿量
我们进一步修改我们的哈密顿量，前一节讲了哈密顿量的作用效果是使得能量自左向右流动，那么如果我们给哈密顿量加一个空间对称性，也就是复共轭，结果会怎么样的。代码如下，结果随后
N=3 # 粒子数是3psi = tensor(basis(2,1),basis(2,0),basis(2,0)) # 制备一个第一个粒子是激发态其他都是基态的初态a_list, ad_list = [], []
for i in range(N): # 用循环定义每个单格点上的产生湮灭算符op_list = [qeye(2)] * Nop_list[i] = destroy(2)a_list.append(tensor(op_list))op_list[i] = create(2)ad_list.append(tensor(op_list))hal = 0
for i in range(N): # 哈密顿量改为每个格点上都有产生算符hal += ad_list[i+1] * a_list[i] + a_list[i+1] * ad_list[i]tlist = np.linspace(0, 10, 10) # 演化时间
res = sesolve(hal,psi,tlist) # 解薛定谔方程
states = [s * s.dag() for s in res.states] # 计算密度矩阵
np.set_printoptions(linewidth=10000) # 设置jupyter显示宽度
states = np.real(np.array(states)) # 提取密度矩阵的实部
rounded_states = np.round(states, decimals=4) # 只显示四位小数
np.set_printoptions(suppress=True) # 不使用科学计数法
print(rounded_states) # 输出每一时刻的密度矩阵
发现了没有，发现了没有。这两张图片貌似没有区别，其实不是的，你会发现激发态密度在周期性振荡。为了方便观看，我把fock态的排序复制下来：|000>，|001>，|010>，|011>，|100>，|101>，|110>，|111>。那么密度矩阵从100%|100>到25%|001>+50%|010>+25%|100>到100%|001>再到25%|001>+50%|010>+25%|100>再到100%|100>。发现没有，形成了一个周期振荡。也就是说初始第一个粒子激发态的能量通过第二个传给了第三个，然后第三个又把能量按照空间反演对称性传了回来，完成了一个周期。貌似状态没有发生改变。
六、展望
先写这么多，后续可能会加入其他更加复杂的物理过程，到时候再慢慢写吧。