使用Qutip计算JC哈密顿量本征值在特定耦合常数处异常

问题根源 #

你遇到的锯齿状曲线，本质是本征值索引跳变导致的：

• Qutip的eigenstates()方法会返回按升序排列的本征值，但当耦合常数g变化时，相邻g点的能级顺序可能发生交换（尤其在delta=0时，能级会出现避免交叉现象）。
• 如果你固定取第i个索引的本征值绘制曲线，当能级顺序交换时，该索引对应的实际能级会突然切换，导致曲线出现断裂/锯齿。

解决方案：本征态跟踪 #

要让曲线连续，需要确保每条曲线对应同一个量子态随g的演化，而非固定索引的本征值。这里提供两种可行方法：

方法1：使用Qutip内置的连续本征态计算函数 #

Qutip的continuous_eigenstates()专门用于处理随参数变化的哈密顿量，会自动通过本征态重叠匹配相邻参数点的同一态，保持本征值连续。

修改后的完整代码：

import qutip as qt
import numpy as np
import matplotlib.pyplot as plt

def JC_hamiltonian(N, omega, delta):
    # 定义以g为参数的哈密顿量函数
    def H(g):
        sz = qt.tensor(qt.qeye(N), qt.sigmaz())
        a = qt.tensor(qt.destroy(N), qt.qeye(2))
        sp = qt.tensor(qt.qeye(N), qt.sigmap())
        sm = qt.tensor(qt.qeye(N), qt.sigmam())
        return omega*a.dag()*a + 0.5*(delta + omega)*sz + g*(a*sp + a.dag()*sm)
    return H

# 参数设置
N = 10  # 腔模截断维度，建议取足够大避免截断误差
omega = 1.0
delta = 0.0
g_list = np.linspace(0, 2, 200)  # 扫描g的范围和点数

# 计算随g变化的连续本征值
H_func = JC_hamiltonian(N, omega, delta)
eigenvals, _ = qt.continuous_eigenstates(H_func, g_list, sparse=False)

# 绘制前7个能级（i=0到6）
plt.figure(figsize=(10,6))
for i in range(7):
    plt.plot(g_list, eigenvals[:, i], label=f"能级 {i}")
plt.xlabel("耦合常数 g")
plt.ylabel("本征能量")
plt.legend()
plt.title("Jaynes-Cummings模型本征能量随g变化曲线（delta=0）")
plt.show()

方法2：手动跟踪本征态（原理展示） #

如果想理解底层逻辑，可以通过计算相邻g点本征态的重叠，手动匹配同一态并重新排序本征值：

import qutip as qt
import numpy as np
import matplotlib.pyplot as plt

def JC(N, g, omega, delta):
    sz = qt.tensor(qt.qeye(N), qt.sigmaz())
    a = qt.tensor(qt.destroy(N), qt.qeye(2))
    sp = qt.tensor(qt.qeye(N), qt.sigmap())
    sm = qt.tensor(qt.qeye(N), qt.sigmam())
    H = omega*a.dag()*a + 0.5*(delta + omega)*sz + g*(a*sp + a.dag()*sm)
    return H.eigenstates()

# 参数设置
N = 10
omega = 1.0
delta = 0.0
g_list = np.linspace(0, 2, 200)

# 初始化跟踪后的本征值列表
tracked_eigenvals = []

# 计算第一个g点的本征值和本征态
eig_vals_prev, eig_vecs_prev = JC(N, g_list[0], omega, delta)
tracked_eigenvals.append(eig_vals_prev)

# 遍历后续g点，手动匹配本征态
for g in g_list[1:]:
    eig_vals_curr, eig_vecs_curr = JC(N, g, omega, delta)
    # 计算前一g点所有态与当前g点所有态的重叠绝对值
    overlap_matrix = np.abs(np.array([
        [prev_vec.overlap(curr_vec) for curr_vec in eig_vecs_curr]
        for prev_vec in eig_vecs_prev
    ]))
    # 找到每个前态对应的当前态（重叠最大的索引）
    match_indices = np.argmax(overlap_matrix, axis=1)
    # 按匹配结果重新排序当前本征值
    sorted_eig_vals = eig_vals_curr[match_indices]
    tracked_eigenvals.append(sorted_eig_vals)
    # 更新前态为当前匹配后的态
    eig_vecs_prev = [eig_vecs_curr[idx] for idx in match_indices]

# 转换为数组并绘图
tracked_eigenvals = np.array(tracked_eigenvals)
plt.figure(figsize=(10,6))
for i in range(7):
    plt.plot(g_list, tracked_eigenvals[:, i], label=f"能级 {i}")
plt.xlabel("耦合常数 g")
plt.ylabel("本征能量")
plt.legend()
plt.title("Jaynes-Cummings模型本征能量随g变化曲线（delta=0）")
plt.show()

额外说明 #

• 你的原始哈密顿量写法是正确的，当delta=0时符合JC模型的标准形式。
• 腔模截断数N建议取足够大（比如N=10），避免截断误差影响能级形状。

内容的提问来源于stack exchange，提问作者Mhd Mahmoud Yehya