模拟一个五量子比特系统，其中四个量子比特（编号为1, 2, 3, 4）分别与第五个量子比特（编号为5）耦合，形成一个星型结构。分析量子比特1和2的纠缠熵随时间的变化。

系统的哈密顿量H描述了量子比特间的相互作用，这里使用sigmax()表示泡利X矩阵，用于描述量子比特间的耦合。qeye(2)是2x2的单位矩阵，用于扩展哈密顿量的维度。

import numpy as np
import matplotlib.pyplot as plt
from qutip import (basis, mesolve, tensor, sigmax, qeye, destroy,entropy_vn, concurrence, ptrace, Options, sigmam)# 配置全局绘图参数
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['axes.unicode_minus'] = Falsedef compute_Gamma(t_points, gamma, omega, kBT):"""根据文档中的公式(39)，计算耗散系数Γ(t)。"""Gamma_values = []dt = t_points[1] - t_points[0]for t in t_points:if t == 0:Gamma_values.append(0)else:# 计算积分 ∫_0^t ⟨{f̂(t), f̂(t')}⟩ dt'integral = 0for t_prime in np.linspace(0, t, int(t / dt)):correlation = 2 * gamma * omega * kBT * np.exp(-gamma * abs(t - t_prime)) * np.cos(omega * (t - t_prime))integral += correlation * dtGamma_values.append(np.max([integral, 0]))  # 确保Gamma值非负return Gamma_values# 第一部分：五量子比特星型耦合系统
def simulate_five_qubits_star():J = 0.3  # 相邻量子比特耦合强度gamma = 0.1  # 耗散率omega = 1.0  # 环境频率kBT = 1.0  # 温度乘以玻尔兹曼常数total_time = 20num_steps = 1000  # 增加时间步数以提高精度t_points = np.linspace(0, total_time, num_steps + 1)# 构建星型哈密顿量 (1-5, 2-5, 3-5, 4-5)H = J * (tensor(sigmax(), qeye(2), qeye(2), qeye(2), sigmax()) +tensor(qeye(2), sigmax(), qeye(2), qeye(2), sigmax()) +tensor(qeye(2), qeye(2), sigmax(), qeye(2), sigmax()) +tensor(qeye(2), qeye(2), qeye(2), sigmax(), sigmax()))psi0 = tensor(basis(2, 0), basis(2, 0), basis(2, 0), basis(2, 0), basis(2, 0))# 计算每个时间点上的耗散系数Γ(t)Gamma_values = compute_Gamma(t_points, gamma, omega, kBT)states = [psi0]for i in range(num_steps):t_start = t_points[i]t_end = t_points[i + 1]Gamma_t = Gamma_values[i]if Gamma_t < 0:Gamma_t = 0  # 确保Gamma值非负c_ops = [np.sqrt(Gamma_t) * tensor(sigmam(), qeye(2), qeye(2), qeye(2), qeye(2))]result = mesolve(H, states[-1], [t_start, t_end], c_ops, [],options={'store_states': True, 'atol': 1e-8, 'rtol': 1e-6})states.append(result.states[-1])entropies_1 = [entropy_vn(ptrace(s, [0])) for s in states]  # Qubit 1entropies_2 = [entropy_vn(ptrace(s, [1])) for s in states]  # Qubit 2return t_points, entropies_1, entropies_2# 第二部分：增强型量子关联分析（含纠缠熵和纠缠度）
def enhanced_quantum_correlation_analysis():J = 0.5gamma = 0.1# 两量子比特系统H = J * tensor(sigmax(), sigmax())psi0 = tensor(basis(2, 0), basis(2, 0))c_ops = [np.sqrt(gamma) * tensor(sigmam(), qeye(2))]t_points = np.linspace(0, 10, 50)opts = {'store_states': True}result = mesolve(H, psi0, t_points, c_ops, [], options=opts)entropies = [entropy_vn(ptrace(s, [0])) for s in result.states]concurrences = [concurrence(s) for s in result.states]return t_points, entropies, concurrences# =============================================================================
# 执行主程序
# =============================================================================
if __name__ == "__main__":print("五量子比特星型系统模拟中...")t_points_star, entropies_1, entropies_2 = simulate_five_qubits_star()print("\n量子关联综合分析中...")t_points_corr, entropies_corr, concurrences = enhanced_quantum_correlation_analysis()print("所有模拟完成！")# 绘制所有图表plt.figure(figsize=(16, 12))# 五量子比特星型系统图表plt.subplot(2, 2, 1)plt.plot(t_points_star, entropies_1, 'b--', label='Qubit 1 Entropy')plt.plot(t_points_star, entropies_2, 'g-.', label='Qubit 2 Entropy')plt.xlabel('Time')plt.ylabel('Entropy')plt.title('Five-Qubit Star: Entanglement Propagation')plt.legend()plt.grid(alpha=0.3)# 增强型量子关联分析图表plt.subplot(2, 2, 2)plt.plot(t_points_corr, entropies_corr, 'b-o', label='Entanglement Entropy')plt.plot(t_points_corr, concurrences, 'r-s', label='Concurrence')plt.ylabel('Quantum Correlation')plt.title('Comparative Analysis of Quantum Correlations')plt.legend()plt.grid(alpha=0.3)plt.tight_layout()plt.show()