模拟一个五量子比特系统，其中四个量子比特（编号为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'] = False


def 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 = 0
            for 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 * dt
            Gamma_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 = 20
    num_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 1
    entropies_2 = [entropy_vn(ptrace(s, [1])) for s in states]  # Qubit 2

    return t_points, entropies_1, entropies_2


# 第二部分：增强型量子关联分析（含纠缠熵和纠缠度）
def enhanced_quantum_correlation_analysis():
    J = 0.5
    gamma = 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()