十七、多粒子激发+振荡+相互作用
这节将是截至目前位置最复杂的一个哈密顿量，它综合了多粒子的激发、振荡和相互作用。废话不多说，代码如下，结果随后（说明：参数比较多，为了控制我们一律将t1设置为1）
N=3
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))
sigmap_list, sigmam_list = [], []
for i in range(N):op_list = [qeye(2)] * Nop_list[i] = sigmap()sigmap_list.append(tensor(op_list))op_list[i] = sigmam()sigmam_list.append(tensor(op_list))
ttotal = 100
tsteps = 1000
tlist = np.linspace(0, ttotal, tsteps)
psi = tensor(basis(2,0),basis(2,0),basis(2,0))
hal = 0
t_1 = 1
t_2 = 1
omega = 1 # 原子失谐
for i in range(N-1):if i%2 == 0: hal += -t_1 * (ad_list[i+1] * a_list[i] + ad_list[i] * a_list[i+1])else:hal += -t_2 * (ad_list[i+1] * a_list[i] + ad_list[i] * a_list[i+1])hal += ad_list[i]hal += omega * (ad_list[i] * a_list[i])
hal += ad_list[-1]
hal += omega * (ad_list[-1] * a_list[-1])
res = sesolve(hal,psi,tlist)
np.set_printoptions(threshold=np.inf) # 遇到显示不全的用
states = [s * s.dag() for s in res.states]
# np.set_printoptions(linewidth=1000000)
states = np.real(np.array(states))
rounded_states = np.round(states, decimals=4)
np.set_printoptions(suppress=True)
# print(rounded_states)
psi000 = []
psi001 = []
psi010 = []
psi011 = []
psi100 = []
psi101 = []
psi110 = []
psi111 = []
for i in range(tsteps):psi000.append(rounded_states[i][0][0])
for i in range(tsteps):psi001.append(rounded_states[i][1][1])
for i in range(tsteps):psi010.append(rounded_states[i][2][2])
for i in range(tsteps):psi011.append(rounded_states[i][3][3])
for i in range(tsteps):psi100.append(rounded_states[i][4][4])
for i in range(tsteps):psi101.append(rounded_states[i][5][5])
for i in range(tsteps):psi110.append(rounded_states[i][6][6])
for i in range(tsteps):psi111.append(rounded_states[i][7][7])
plt.plot(tlist,psi000, label='000')
plt.plot(tlist,psi001, label='001')
plt.plot(tlist,psi010, label='010')
plt.plot(tlist,psi100, label='100')
plt.plot(tlist,psi011, label='011') # 顺序换了一下，为了按激发数归组
plt.plot(tlist,psi101, label='101')
plt.plot(tlist,psi110, label='110')
plt.plot(tlist,psi111, label='111')
plt.xlabel('Time')
plt.ylabel('Occuption')
plt.ylim(-0.1,1.1)
# plt.xlim(37,40)
text = "Evolution of states"
plt.title(text)
plt.legend()
只能说，结果乱中有序吧。不出意外，基态|000>占据还是很快衰减为0，其他单、双和全激发态各领风骚。但是我把时间调至无穷大的时候，出现了一起基态的死灰复燃，应该是计算问题或者是不稳定的态。
接下来我们逐渐改变参数，看看会不会发现一些有趣的事情。
我们将t2=0后，就是一些简单的振荡。没什么特别的东西。但是当我们把t2调到0.5-0.55范围内的时候出现了一些相对稳定的态分布，如图
在一段时间内态分布是相对稳定的，或许这个参数区间可以达到相对的稳态。比起它，后面这个参数区间更加稳定
这个稳态持续的时间更长，而且态的波动更小。t2大于2之后基本就没什么稳态了，但是也有可能是我没试到
接下来开始调整omega的值：
我们在omega=0.7附近也找到了一个类似的稳态，甚至他的稳定时间更长。这个参数区间的态组合为15%|011>+30%|101>+15%|110>+40%|111>，这样组合出来的激发数密度为|0.85 0.7 0.85>。
当omega再增大的时候我没有发现什么特别的区间
在一些其他的参数区间偶尔发现了一个稳态
看来当t1=1，t2的取值在2附近的时候更容易有稳态吧（只是个非常不成熟的猜测）
打脸来的如此之快
这里也有一组类似稳态的参数区间。试参数是一项非常耗时耗力的工作，今天就先这样吧。
十八、展望
后续还会继续改变哈密顿量，希望发现更有趣的物理现象。