为了理解开放量子系统，我们需要能够描述状态不完全已知系统的方法。这属于统计力学范畴，而密度矩阵则是解决这一问题的核心工具。

主方程

若已知所有粒子的哈密顿量，理论上可以通过薛定谔方程或海森堡方程得到所有粒子的动力学行为，但对于一个自由度很高的开放系统，薛定谔方程或海森堡方程一般并不实用，因此我们需要一个热库模型
对于腔内光子的一个算符，有广义量子朗之万方程，对应于封闭系统的海森堡方程

$\dot{c}=i[H_{\mathbf{sys}},c]-(\frac{\gamma}{2}a^\dagger+\sqrt{\gamma}A_{\mathbf{in}}^\dagger)[a,c]+[a^\dagger,c]\left(\frac{\gamma}{2}a+\sqrt{\gamma}A_{\mathbf{in}}\right)$

对于原子，我们通常使用他们的态而不是对应的算符来描述系统，因为原子态更加直观，因此我们无法使用量子朗之万方程
密度矩阵常用于描述态，而主方程是密度矩阵的方程，这就是采用主方程的原因，其在热库模型中经常使用。

经典刘维尔方程

为了建立直觉上的理解，我们首先回顾其经典对应。考虑一个 $N$ 粒子经典系统，系统在任意时刻的全部状态可由 $6N$ 维相空间中的一个点描述，这个相空间由其正则坐标和正则动量给出

$(q_1,q_2,\ldots q_N,p_1,p_2,\ldots p_N)$

在经典统计力学中，我们不再使用相空间的单点 $(q_i,p_i)$ 来描述具有多个可能状态的系统，而是使用概率分布函数

$\rho(q_1,\ldots,q_N,p_1,\ldots,p_N,t)$

其中 $\rho(q,p,t)dqdp$ 表示在 $t$ 时刻，系统状态落在点 $(q,p)$ 附近极微小 $dqdp$ 体积内的概率。若将这个分布解释为系综中总系统数量的分布，则 $\rho(q,p,t)$ 的归一化条件为

$\int\rho(q,p,t)dqdp=1$

有时候，我们也会将 $\rho(q,p,t)$ 解释为相空间中某点附近的粒子数密度，此时归一化条件变为

$\int\rho(q,p,t)dqdp=N$

如果已知系统的哈密顿量 $H$ ，那么相空间分布函数的时间演化由刘维尔方程给出

$\frac{\partial}{\partial t}\rho=-\{\rho,H\}=-\sum_{i=1}^N(\frac{\partial H}{\partial p_i}\frac{\partial\rho}{\partial q_i}-\frac{\partial H}{\partial q_i}\frac{\partial\rho}{\partial p_i})$

我们可以定义刘维尔算符或刘维尔量为

$\hat{L}=\sum_{i=1}^N\frac{\partial H}{\partial p_i}\frac{\partial}{\partial q_i}-\frac{\partial H}{\partial q_i}\frac{\partial}{\partial p_i}$

这样刘维尔方程可以简洁地写成

$\frac{\partial}{\partial t}\rho=-\{\rho,H\}=-\hat{L}\rho$

这个体系是经典统计力学的基石，成功描述了从热力学到流体力学内一切系综的演化。但是我们在量子力学中面临着类似的挑战：如何描述一个与一个不可控的大热库相互作用的量子系统？我们不可能去追踪热库的无数个自由度，因此很自然地需要从统计上描述一个热库。而由于系统和热库之间的相互作用，我们也必须放弃用一个纯态波函数 $\left|\psi\right\rangle$ 来描述系统的想法，转而使用密度矩阵 $\rho$ 来描述系统的状态。

经典分布函数 $\rho_{cl}$ 的量子对应就是密度矩阵 $\rho$ 。由于经典刘维尔方程决定了经典系综的演化，那么量子主方程就应决定密度矩阵 $\rho$ 的演化，使得我们能够捕捉到系统的内部量子行为和其与热库的耗散相互作用。

接下来我们将定义量子密度矩阵并导出其运动方程，有利于我们理解开放量子系统的动力学。

密度矩阵

既然 $\ket{\psi}$ 已经作为量子力学的基础表示，那为什么还要引入 $\rho$ 呢？因为密度矩阵在实际应用中的范围更大并具有计算优势。

波函数 $\ket{\psi}$ 只能描述纯态孤立系统，而密度矩阵 $\rho$ 可以描述混态系统或者与环境纠缠的子系统。例如，考虑等概率处于量子态 $\ket{a}$ 和 $\ket{b}$ 的一个系统，这个统计上的混合态无法用单个的波函数表示：叠加态 $\ket{\psi}=(\ket{a}+\ket{b})/\sqrt{2}$ 描述的是一个确定的相干态，而不是一个非相干的概率混合态。

而密度矩阵为纯态和混态提供了统一的框架。密度矩阵得名字来源于，他在量子力学中扮演了概率分布的角色，同时包含了统计混合和量子相干性。虽然在纯态下密度矩阵等价于波函数，但 $\rho$ 使得我们能够严格处理退相干、弛豫和其他不可逆过程。因此，密度矩阵是高级时变问题中不可或缺的工具，尤其在弛豫动力学和非线性光谱信号分析中，系统统计特性至关重要。

考虑一个算符 $A$ 的期望值 $\langle A\rangle=\langle\psi|A|\psi\rangle$ ，我们可以将波函数 $\ket{\psi}$ 展开在某个完备正交归一基底 $\{|n\rangle\}$ 上

$\langle A\rangle=\langle\psi|A|\psi\rangle=\sum_m\sum_nu_n^*u_m\langle n|A|m\rangle$

这表明，期望就是所有矩阵元的加权求和，其中权重由 $u_n^*u_m$ 给出。我们来看 $u_n^*u_m$ 的物理意义，定义密度矩阵为

$\begin{aligned} \rho & =|\psi(t)\rangle\langle\psi(t)| \\[1em] & =\sum_{m}u_{m}(t)|x_{n}\rangle\sum_{n}u_{n}(t)^{*}\langle x_{n}| \\[1em] & =\sum_{m}\sum_{n}u_{m}(t)u_{n}^{*}(t)|x_{m}\rangle\langle x_{n}| \\[1em] & =\sum_m\sum_n\rho_{mn}(t)|x_m\rangle\langle x_n| \end{aligned} \qquad \Longrightarrow \qquad \rho= \begin{pmatrix} \rho_{11} & \rho_{12} & \cdots & \rho_{1N} \\[1em] \rho_{21} & \rho_{22} & \cdots & \rho_{2N} \\[1em] \vdots & \vdots & \ddots & \vdots \\[1em] \rho_{N1} & \rho_{N2} & \cdots & \rho_{NN} \end{pmatrix}$

由于量子态由一组 $\{u_n\}$ 唯一确定，那么密度矩阵的矩阵元（也就是权重 $u_n^*u_m$ ）包含了这个量子态的所有信息：对角元为概率幅，非对角元为量子相干。

通过定义矩阵的迹为

$\mathrm{Tr}(\rho)=\sum_n\rho_{nn}=\sum_nu_nu_n^*$

可以得到

$\langle A\rangle=\mathrm{Tr}(\rho A)$

密度矩阵的性质

$\rho_{nn} \geq 0$ ，因为 $\rho_{nn}$ 代表处于态 $\ket{n}$ 的概率
定义： $\text{Tr}(\rho)=\sum_n \rho_{nn}$ ，那么有 $\text{Tr}(\rho)=1$ ，因为波函数满足归一化 $\left|\psi\left(t\right)\right|^2=1$ ，可得 $\mathrm{Tr}(\rho)=\sum_n\rho_{nn}=\sum_n|u_n(t)|^2=1$

纯态密度矩阵的运动方程

在定义纯态密度矩阵后，我们想要得到其运动方程，这就是经典刘维尔方程的量子对应，先从时间依赖的薛定谔方程出发：

$i\frac{\partial}{\partial t}|\psi(t)\rangle=H|\psi(t)\rangle$

对这个方程取厄米共轭可得到右矢的演化

$-i\frac{\partial}{\partial t}\langle\psi(t)|=H\langle\psi(t)|$

对上两式两边相乘得到密度矩阵的时间导数

$\begin{aligned} i\frac{\partial}{\partial t}\rho & =i(|\psi(t)\rangle\langle\psi(t)|+|\psi(t)\rangle\langle\dot{\psi(t)}|) \\[1em] & =i(\frac{1}{i}H|\psi(t)\rangle\langle\psi(t)|-\frac{1}{i}|\psi(t)\rangle\langle\psi(t)|H) \\[1em] & = \begin{bmatrix} H,\rho \end{bmatrix} \end{aligned}$

因此得到密度矩阵的运动方程为

$i\frac{\partial}{\partial t}\rho=[H,\rho]$

对于一个纯态来说，这个方程与薛定谔方程完全等价。

混态密度矩阵及其推导

纯态下的密度矩阵运动方程似乎和薛定谔方程没什么区别，然而，在一些量子系综中，系统并非都被制备在相同的量子态上，还有在一些我们只知道其统计性质的量子系统（例如热库）中，密度矩阵才真正彰显其神威。

考虑一个具有 $N$ 个量子系统的巨大系综，假设其中有 $N_i$ 个系统被制备在纯态 $\ket{\psi_i}$ 上，从该系综中随机抽取一个系统，其量子态为 $\ket{\psi_i}$ 的概率为

$p_i=\frac{N_i}{N},\quad\sum p_i=1$

这是一个混态，其表示给定系统处于某个纯态经典统计信息。我们定义混态的密度矩阵为各纯态投影算符的加权求和

$\rho^\mathrm{mix}=\sum p_i\rho_i^\mathrm{pure}=\sum p_i|\psi_i\rangle\langle\psi_i|$

这个定义概括了纯态情形，并同时包含了系综的经典信息和量子信息。我们接下来将看到，混态密度矩阵的运动方程与纯态情形是相同的。

$\rho$ 的时间依赖性体现在纯态 $\ket{\psi_i(t)}$ 上，因此我们通过乘积法则和每个纯态分量的已知演化来计算混态密度矩阵的时间导数

$\begin{aligned} i\hbar\frac{\partial\rho}{\partial t} & =i\hbar\frac{\partial}{\partial t}\left(\sum_{i}p_{i}|\psi_{i}\rangle\langle\psi_{i}|\right)=\sum_{i}p_{i}i\hbar\left(\frac{\partial|\psi_{i}\rangle}{\partial t}\langle\psi_{i}|+|\psi_{i}\rangle\frac{\partial\langle\psi_{i}|}{\partial t}\right) \\[1em] &=\sum_ip_i(H|\psi_i\rangle\langle\psi_i|-|\psi_i\rangle\langle\psi_i|H)=\sum_ip_i\left[H,\rho_i\right] =\left[H,\sum_ip_i\rho_i\right] \\[1em] & =[H,\rho] \end{aligned}$

因此，我们得到了混态密度矩阵的冯诺依曼方程

$i\hbar\frac{\partial\rho}{\partial t}=[H,\rho]$

所以，纯态和混态下算符 $A$ 的期望值都可以由下式给出

$\langle A\rangle=\mathrm{Tr}(A\rho)=\mathrm{Tr}(\rho A)$

密度矩阵体系的核心是冯·诺依曼方程，它既能够描述孤立量子系统的相干动力学，也能够描述纯态和混合态系综的统计特性。

旋转框架

在研究光与物质的相互作用（例如激光驱动双能级原子）时，我们经常面临着一个困难：哈密顿量是时变的。旋转框架可以将在光学频率处的快速、平凡的振荡与动力学过程分离，从而简化这一问题。

消除时变哈密顿量

考虑一个标准的双能级原子，基态为 $\ket{1}$ ，激发态为 $\ket{2}$ ，跃迁频率为 $\omega_2$ 。当系统被频率为 $\omega$ 的单色激光场驱动时，哈密顿量为

$\begin{aligned} H(t)=\frac{\omega_2}{2}\sigma_z+er_{12}\cdot Ee^{-i\omega t}\cdot\sigma_{21}+er_{12}\cdot Ee^{i\omega t}\cdot\sigma_{12} \end{aligned}$

其中 $\sigma_{12}=\ket{1}\bra{2}$ 和 $\sigma_{21}=\ket{2}\bra{1}$ 分别是原子的降算符和升算符。

由于振荡项 $e^{\pm i\omega t}$ 的存在，哈密顿量是时间依赖的，直接求解薛定谔方程是十分困难的。然而，直觉和计算都告诉我们：系统的响应（例如诱导的原子偶极矩）将主要在驱动频率 $\omega$ 下振荡。这表明频率为 $\omega$ 的快速振荡往往代表动力学中的平凡部分，从而掩盖了一些潜在的物理现象。

我们的核心思想是找到一种能够消除这种显式时间依赖的方法，从而简化分析。一种简单的假设是通过替换： $\sigma\to\sigma e^{-i\omega t}，\sigma^{\dagger}\to\sigma^{\dagger}e^{i\omega t}$ 将时间依赖性吸收到算符中来简单处理，那么就会得到与时间无关的哈密顿量

$H=\omega_2\sigma_{22}+er_{12}\cdot E\sigma_{21}+er_{12}\cdot E\cdot\sigma_{12}\quad\mathrm{X}$

看起来是更简单了，但在数学上是不成立的，因为在其变换算符的同时，没有对态矢量进行变换，我们应该对系统的整体量子描述进行幺正变换：从静止参考系到一个旋转频率为 $\omega$ 的旋转参考系。

我们先写出哈密顿量的矩阵形式

$H= \begin{pmatrix} 0 & \Omega e^{i\omega t} \\[1em] \Omega^*e^{-i\omega t} & \omega_2 \end{pmatrix}$

再定义一个幺正算符为

$R= \begin{pmatrix} e^{if(t)} & 0 \\[1em] 0 & e^{-if(t)} \end{pmatrix}$

计算 $R\cdot H\cdot R^{\dagger}$ 得到

$R\cdot H\cdot R^\dagger= \begin{pmatrix} 0 & \Omega e^{i\omega t-2if(t)} \\[1em] \Omega^*e^{-i\omega t+2if(t)} & \omega_2 \end{pmatrix}$

如果我们令 $f(t)=\omega t/2$ ，那么 $R\cdot H\cdot R^{\dagger}$ 就与时间无关，这就是旋转框架的核心：一个幺正变换在驱动频率下通过旋转基矢来消除哈密顿量的显式时间依赖性。旋转框架无非就是找到合适的 $R(t)$ ，使得在新的参考系下，哈密顿量变为时间无关的形式。

二能级系统的旋转框架

旋转框架来源于核磁共振，用于简化在振荡磁场下的自旋分析。旋转框架的目的：进行一个幺正变换，将系统转换到一个以驱动场频率旋转的参考系中，从而消除薛定谔方程的时间依赖性。

我们先将实验室参考系中的态矢量哈密顿量分为自由项和相互作用项

$H_\mathrm{lab}=H_0+V(t)$

其中

$\begin{aligned} H_0 & =\omega|e\rangle\langle e|=\frac{\omega}{2}\sigma_z+\mathrm{const.}\\[1em] V(t) & =\left(\Omega e^{-i\omega t}\sigma_{21}+\Omega^*e^{i\omega t}\sigma_{12}\right) \end{aligned}$

其中 $\Omega=E\cdot er_{12}$ 是拉比频率。系统的时间演化由含时薛定谔方程控制

$i\frac{\partial}{\partial t}|\psi(t)\rangle=(H_0+V(t))|\psi(t)\rangle$

现在开始进行旋转框架变换，定义一个旋转算符并施加到态矢量上，将其变换到新的框架中

$|\phi(t)\rangle=R(t)|\psi(t)\rangle$

那么逆变换为

$|\psi(t)\rangle=R^\dagger(t)|\phi(t)\rangle$

我们的目的就是找到一个有效的哈密顿量，其能够通过类似于薛定谔方程的形式来描述 $|\phi(t)\rangle$ 的演化。为了做到这一点，我们对 $|\phi(t)\rangle$ 求时间导数

$\begin{aligned} i\frac{\partial}{\partial t}|\phi\rangle & =i\frac{\partial}{\partial t}R(t)|\psi\rangle \\[1em] & =[i\frac{\partial}{\partial t}R(t)]|\psi\rangle+iR(t)\frac{\partial}{\partial t}|\psi\rangle \\[1em] & =i\dot{R(t)}|\psi(t)\rangle+R(t)(H_{0}+V)|\psi(t)\rangle \\[1em] & =[i\dot{R(t)}+R(t)(H_0+V)]|\psi(t)\rangle \end{aligned}$

将 $|\psi(t)\rangle=R^\dagger(t)|\phi(t)\rangle$ 代入上式，得到

$i\frac{\partial}{\partial t}|\phi(t)\rangle=\left[i\dot{R}(t)R^\dagger(t)+R(t)(H_0+V(t))R^\dagger(t)\right]|\phi(t)\rangle$

将这个方程与薛定谔方程进行比较，我们可以提取出旋转框架下的有效哈密顿量

$H_{\mathrm{rot}}=[i\dot{R(t)}R^\dagger(t)+R(t)(H_0+V)R^\dagger(t)]$

因此，在旋转框架下，得到一个类似于薛定谔方程的形式

$i\frac{\partial}{\partial t}|\phi\rangle=H_{\mathrm{rot}}|\phi\rangle$

这个结果是具有一般性的，关键在于有策略地选择 $R(t)$ 来简化 $H$ 。我们选择 $R(t) = e^{iH_0t/\hbar}$ ，是自由演化算符的逆变换，这样就能够精确地“逆转”由 $H_0$ 产生的自由演化，从而分离出自由演化的影响，使得态矢量的演化仅由相互作用项 $V(t)$ 决定。

接下来我们来看这个幺正算符的用处。对于双能级系统与单色驱动场相互作用的情形，我们选择

$R(t)=e^{i\omega\sigma_{22}t}$

这表示以驱动激光场 $\omega$ 的频率进行“旋转”。这个旋转算符具有以下性质

$R^\dagger(t)=e^{-i\omega\sigma_{22}t} \qquad \qquad \dot{R}(t)=i\omega\sigma_{22}e^{i\omega\sigma_{22}t}$

将这些结果代入有效哈密顿量 $H_{\mathrm{rot}}$ 的表达式中，得到

$\begin{aligned} H_{\mathrm{rot}} & =[i\dot{R(t)}R^{\dagger}(t)+R(t)(H_{0}+V)R^{\dagger}(t)] \\[1em] & =[i\cdot i\omega\sigma_{22}e^{i\omega\sigma_{22}t}e^{-i\omega\sigma_{22}t}+e^{i\omega\sigma_{22}t}H_{0}e^{-i\omega\sigma_{22}t}+e^{i\omega\sigma_{22}t}Ve^{-\omega\sigma_{22}t} \\[1em] & =[-i\omega\sigma_{22}+e^{i\omega\sigma_{22}t}\omega_2\sigma_{22}e^{-i\omega\sigma_{22}t}+e^{i\omega\sigma_{22}t}Ve^{-i\omega\sigma_{22}t}] \\[1em] & =-\omega\sigma_{22}+\omega_{2}\sigma_{22}+RVR^{\dagger} \\[1em] & =\delta\sigma_{22}+RVR^{\dagger} \end{aligned}$

其中 $\delta=\omega_2-\omega$ 是原子跃迁频率与驱动频率之间的失谐。利用 $\text{Baker-Campbell-Hausdorff}$ 公式，计算第二项

$\begin{aligned} R(t)V(t)R^{\dagger}(t) & =e^{i\omega\sigma_{22}t}\left(\Omega e^{-i\omega t}\sigma_{21}+\Omega^{*}e^{i\omega t}\sigma_{12}\right)e^{-i\omega\sigma_{22}t} \\[1em] & =\Omega e^{-i\omega t}\left(e^{i\omega\sigma_{22}t}\sigma_{21}e^{-i\omega\sigma_{22}t}\right)+\Omega^{*}e^{i\omega t}\left(e^{i\omega\sigma_{22}t}\sigma_{12}e^{-i\omega\sigma_{22}t}\right) \\[1em] & =\Omega e^{-i\omega t}\left(\sigma_{21}e^{+i\omega t}\right)+\Omega^{*}e^{i\omega t}\left(\sigma_{12}e^{-i\omega t}\right) \\[1em] & =\Omega\sigma_{21}+\Omega^{*}\sigma_{12} \end{aligned}$

此时复杂的时间依赖项被消除，留下了两个原子能级之间的静态耦合，这就是旋转框架的威力所在。那么有效哈密顿量最终为

$\begin{aligned} H & =-\omega\sigma_{22}+\omega_2\sigma_{22}+RVR^\dagger=\delta\sigma_{22}+\Omega\sigma_{21}+\Omega^*\sigma_{12} \\[1em] & =\delta\sigma_{22}+\Omega\sigma_{21}+\Omega^*\sigma_{12} \end{aligned}$

旋转框架的作用

简化哈密顿量：使用 $\sigma_{12}e^{-i\omega t}$ 替换 $\sigma_{12}$ 抵消了哈密顿量中的显式时间依赖，使得求解薛定谔方程更加容易，从而能够找到动力学、稳态和特征值的精确解。
旋转框架哈密顿量 $H_{\mathrm{rot}}$ 的物理解释：失谐 $\delta$ 表示原子和激光的能量不匹配，而拉比频率 $\Omega$ 则像一个代表了激光场与原子之间耦合强度的静态有效场。

旋转框架的核心是一种变量代换的数学方法，用于简便计算，例如用 $t=x-a$ 代换 $x$ 求解 $f(x)=x^2-2ax+a^2$ 时，可以简化为 $f(t)=t^2$ 。同理，旋转框架对量子态施加了一个“替换”（幺正变换 $R(t)$ ），能够消除动力学中的复杂部分（频率 $\omega$ 的振荡），从而得到一个更简单的方程。旋转框架中的动力学现象和随驱动场旋转时所观测到的结果一致。

相互作用绘景下的主方程

密度矩阵主方程

我们现在将目光从腔体转向与热库耦合的原子（或其他量子系统）动力学，我们希望得到仅针对系统自身的封闭运动方程，从而得到主方程。

系统加上环境的哈密顿量可分为三个部分

$H=H_s+H_b+V$

$H_s$ 是系统自身的特征频率

$H_s=\omega_s\sigma_{22}+\cdots$

$H_b$ 是一个具有大量环境模式的热库的本征哈密顿量，例如自由空间的电磁模式或固体中的声子模式

$H_b=\sum_{\omega_i}b^\dagger(\omega_i)b(\omega_i)=\int_{-\infty}^{+\infty}d\omega D(\omega)\omega b^\dagger(\omega)b(\omega)$

$V$ 是系统与热库之间的相互作用

$V=\sum_{\omega_i}\kappa(\omega_i)[b^\dagger(\omega_i)\sigma_{12}+b(\omega_i)\sigma_{21}]=\int_{-\infty}^{+\infty}d\omega D(\omega)\kappa(\omega)\left(b^\dagger(\omega)\sigma_{12}+b(\omega)\sigma_{21}\right)$

整个宇宙（系统+热库）的总密度矩阵 $\rho$ 的演化由主方程控制

$\frac{d}{dt}\rho=-i[H_s+H_b+V,\rho]$

但是我们只关心系统的状态，这就需要用到约化密度矩阵，将热库的自由度平均掉（迹掉）便可以得到系统的约化密度矩阵

$\rho_s=\mathrm{Tr}_b(\rho)$

其中 $\rho_s$ 是约化密度矩阵，与热库无关。我们希望导出 $\rho_s$ 的一个自洽、封闭运动方程，其不显式依赖于热库。接下来，我们将通过旋转框架和一系列合理的近似来推导 $\mathcal{L}[\rho_{s}]$ 的具体形式，最终得到著名的林德布拉德形式。这使得我们能够仅根据系统的自身性质及其与热库的耦合来预测系统的耗散动力学，例如自发辐射和退相干。

弱耦合近似

相互作用绘景就是一个旋转框架，我们选择将旋转算符 $R$ 定义为

$R=e^{i(H_s+H_b)t}\quad\qquad\quad R^\dagger=e^{-i(H_s+H_b)t}$

通过该旋转算符，可以转换到相互作用绘景下，之前得到哈密顿量的通用变换为

$\begin{array} {c}|\phi\rangle=R|\psi\rangle \\[1em] H_{\mathrm{rot}}=i\dot{R(t)}R^\dagger(t)+R(t)(H_0+V)R^\dagger(t) \end{array}$

因此薛定谔方程，再按照上述方法得到冯诺依曼方程

$i\frac{\partial}{\partial t}|\phi\rangle=H_{\mathsf{rot}}|\phi\rangle\quad\to\quad\frac{d}{dt}\rho_{\mathsf{rot}}=-i[H_{\mathsf{rot}},\rho_{\mathsf{rot}}]$

计算主方程的哈密顿量

$\begin{aligned} H_{\mathrm{rot}} & =i[R(t)R^{\dagger}(t)+R(t)(H_{s}+H_{b}+V)R^{\dagger}(t)] \\ & =-\omega_sa^\dagger a-\int_{-\infty}^{+\infty}\omega b^\dagger(\omega)b(\omega)d\omega+R(H_s+H_b)R^\dagger+RVR^\dagger \\ & =-\omega_sa^\dagger a-\int_{-\infty}^{+\infty}\omega b^\dagger(\omega)b(\omega)d\omega+H_s+H_b+RVR^\dagger \\ & =-\omega_sa^\dagger a-\int_{-\infty}^{+\infty}\omega b^\dagger(\omega)b(\omega)d\omega+H_s+H_b+RVR^\dagger \\ & =RVR^{\dagger} \\ & =V_{\mathrm{rot}} \end{aligned}$

因此在旋转框架（在这个特例中就是相互作用绘景）下有

$\frac{d}{dt}\rho_{\mathrm{rot}}=-i[V_{\mathrm{rot}},\rho_{\mathrm{rot}}]$

其中 $\rho_{\mathrm{rot}}$ 是相互作用绘景下的总密度矩阵。对上式从 $0$ 积分到 $t$ ，得到

$\rho_{\mathsf{rot}}(t)=\rho_{\mathsf{rot}}(0)-i\int_0^tdt^{\prime}[V_{\mathsf{rot}}(t^{\prime}),\rho_{\mathsf{rot}}(t^{\prime})]$

这说明 $t$ 时刻的密度矩阵 $\rho_{\mathrm{rot}}(t)$ 由从 $0$ 到 $t$ 的全部积分决定，这虽然自洽但是很难求解，我们可以使用迭代法解决这个问题：将上式右侧的表达式代代回到其自身的对易子中

$\rho_{\mathbf{rot}}(t)=\rho_{\mathbf{rot}}(0)-i\int_0^tdt^{\prime}\left[V_{\mathbf{rot}}(t^{\prime}),\quad\rho_{\mathbf{rot}}(0)-i\int_0^{t^{\prime}}dt^{\prime\prime}[V_{\mathbf{rot}}(t^{\prime\prime}),\rho_{\mathbf{rot}}(t^{\prime\prime})]\right]$

假设 $1$ ：弱耦合近似

我们引入一个重要的近似：假设系统与热库的耦合很弱。这就允许我们将相互作用 $V$ 当作微扰处理，并将 $\tilde {\rho}\left(t\right)$ 的迭代展开截断在二阶项。

取迭代一次后的表达式

$\rho_{\mathrm{rot}}(t)-\rho_{\mathrm{rot}}(0)=-i\int_0^tdt^{\prime}[V_{\mathrm{rot}}(t^{\prime}),\rho_{\mathrm{rot}}(0)]-\int_0^tdt^{\prime}\int_0^{t^{\prime}}dt^{\prime\prime}\left[V_{\mathrm{rot}}(t^{\prime}),\quad[V_{\mathrm{rot}}(t^{\prime\prime}),\rho_{\mathrm{rot}}(t^{\prime\prime})]\right]$

那么密度矩阵的运动方程可以重写为

$\frac{d}{dt}\rho_{\mathbf{rot}}(t)=-i[V_{\mathbf{rot}}(t),\rho_{\mathbf{rot}}(0)]-\int_0^tdt^{\prime}\left[V_{\mathbf{rot}}(t),\quad[V_{\mathbf{rot}}(t^{\prime}),\rho_{\mathbf{rot}}(t^{\prime})]\right]$

接下来我们推导系统约化密度矩阵的运动方程，将热库的自由度从总密度矩阵中迹掉

$\rho_{\mathrm{s,rot}}(t)=\mathrm{Tr}_b(\rho_{\mathrm{rot}}(t))$

式中下标 $s$ 表示系统， $b$ 表示热库， $rot$ 表示旋转框架。将这个偏迹应用至旋转框架的密度矩阵运动方程中，得到

$\frac{d}{dt}\rho_{\mathrm{s,rot}}(t)=\frac{d}{dt}\mathrm{Tr}_b(\rho_{\mathrm{rot}}(t))=\mathrm{Tr}_b(\frac{d}{dt}\rho_{\mathrm{rot}}(t))$

为了符号上的简便，接下里将省略 $\text{rot}$ 下标，除非另有说明，否则所有算符和密度矩阵都处于相互作用表象下。所以约化密度矩阵的运动方程为

$\frac{d}{dt}\rho_s(t)=-i\mathrm{Tr}_b[V(t),\rho(0)]-\mathrm{Tr}_b\int_0^tdt^{\prime}\left[V(t),\left[V(t),\rho(t^{\prime})\right] \right]$

在薛定谔绘景中，系统（例如一个原子）和热库（例如真空电磁场模式）的相互作用由下式给出

$V=\int d\omega\kappa(\omega)(ab^\dagger(\omega)+a^\dagger b(\omega))$

在旋转框架下，从薛定谔表象到相互作用表象的变换由系统和热库的自由演化 $R(t)=e^{i(H_{s}+H_{b})t}$ 完成，变换之后得到了两个代换

$a\to ae^{-i\omega_st},\quad b(\omega)\to b(\omega,0)e^{-i\omega t}$

$V=\int d\omega \quad \kappa(\omega) \left( a b^{\dagger}(\omega)e^{-i\omega_{s}t}e^{i\omega t}+ a^{\dagger} b(\omega)e^{i\omega_{s}t}e^{-i\omega t} \right)$

主方程的推导

接下来我们将计算下式的两个部分

$\frac{d}{dt}\rho_s(t)=\textcolor{blue}{-i\mathrm{Tr}_b[V(t),\rho(0)]}-\textcolor{red}{\mathrm{Tr}_b\int_0^tdt^{\prime}\left[V(t),\left[V(t),\rho(t^{\prime})\right] \right]}$

大均匀热库近似

为了计算第一项（蓝色部分），我们需要做一些近似

假设 $2$ ：可分离态

$\rho(t)=\rho_s(t)\otimes\rho_b(0)$

我们假设热库过于庞大，以至于系统与热库的耦合不会显著改变热库的状态。

$\begin{aligned} & \mathsf{Tr}_b\big([V(t),\rho(0)]\big) \\ & =\kappa\cdot\text{Tr}_b\big([\int d\omega a^\dagger b(\omega)e^{i\omega_st}e^{-i\omega t}+h.c.,\rho_s\otimes\rho_b\big]\big) \\ & =\kappa\cdot\mathsf{Tr}_{b}\big(\int d\omega a^{\dagger}b(\omega)e^{i\omega_{s}t}e^{-i\omega t}\rho_{s}\otimes\rho_{b}-\int d\omega \rho_{s}\otimes\rho_{b}a^{\dagger}b(\omega)e^{i\omega_{s}t}e^{-i\omega t}+h.c.\big) \\ & =\kappa\cdot\mathsf{Tr}_b\big(\int d\omega e^{i\omega_st}e^{-i\omega t}a^\dagger\rho_s\otimes b(\omega)\rho_b-\int d\omega e^{i\omega_st}e^{-i\omega t}\rho_sa^\dagger\otimes\rho_bb(\omega)+h.c.\big) \\ & =\kappa\cdot\int d\omega e^{i\omega_{s}t}e^{-i\omega t}a^{\dagger}\rho_{s}\cdot\mathsf{Tr}_{b}\big(b(\omega)\rho_{b}\big)-\kappa\cdot\int d\omega e^{i\omega_{s}t}e^{-i\omega t}\rho_{s}a^{\dagger}\cdot\mathsf{Tr}_{b}\big(\rho_{b}b(\omega)\big)+h.c \end{aligned}$

我们注意到，对热库的偏迹引入了期望值，例如 $\mathrm{Tr}_b\left(b^\dagger(\omega)\rho_b\right)=\langle b^\dagger(\omega)\rangle$ 和 $\mathrm{Tr}_{b}\left([b(\omega),\rho_{b}]\right)=\langle b(\omega)\rangle$ ，那么我们正好可以引入热库的统计性质

假设 $3$ ：热库的均匀性

对于一个处在静态的热库（例如真空态或热平衡态），其模式的期望值为零

$\langle b(\omega)\rangle=\mathrm{Tr}_b\left(\rho_bb(\omega)\right)=\mathrm{Tr}_b\left(b(\omega)\rho_b\right)=0$

这意味着热库在某种程度上是空的，例如，能使平均电磁场强度为零。

因此第一项（蓝色项）的每一部分都为零

$\mathrm{Tr}_b\left([V(t),\rho(0)]\right)=0$

第二项计算

接下来我们计算第二项（红色部分）

这一项代表系统-热库相互作用随时间的累积效应，包含了系统的历史记忆，是我们接下来推导主方程的重点近似对象。在开始计算前，我们需要一些必要的数学关系

$\begin{aligned} & \int_{-\infty}^{+\infty}e^{i\omega t^{\prime}}dt^{\prime}=2\pi\delta(\omega) \\[1em] & \int_{-\infty}^{+\infty}e^{-i\omega(t-t^{\prime})}dt^{\prime}=2\pi\delta(\omega) \\[1em] & \int_{-\infty}^te^{-i\omega(t-t^{\prime})}dt^{\prime}=\frac{1}{2}2\pi\delta(\omega)+\frac{P.V.}{\omega} \end{aligned}$

此外我们还引由热库的密度矩阵取表示热库的一些本征性质，对于热库模的产生和湮灭算符，分别有以下关系

$\begin{aligned} &\mathrm{Tr}_b\left(\rho_bb^\dagger(\omega)b(\omega)\right)=\mathrm{Tr}_b\left(b(\omega)\rho_bb^\dagger(\omega)\right)=\mathrm{Tr}_b\left(b^\dagger(\omega)b(\omega)\rho_b\right)=n(\omega) \\[1em] &\mathsf{Tr}_b\left(b(\omega)b^\dagger(\omega)\rho_b\right)=\mathsf{Tr}_b\left((\delta(\omega)+b^\dagger(\omega)b(\omega))\rho_b\right)=n(\omega)+\delta(\omega) \end{aligned}$

其中 $n(\omega)$ 是频率为 $\omega$ 的热库模的平均光子数， $\delta(\omega)$ 函数表示与算符序或恒等相关的函数。在热库模型背景下，对 $\delta$ 函数项的解释需要谨慎考量。我们正式开始第二项的计算

$\begin{aligned} & \mathsf{Tr}_b\int_0^tdt^{\prime}\left[V(t),[V(t^{\prime}),\rho(t^{\prime})]\right] \\ & =\mathsf{Tr}_b\int_0^tdt^{\prime}\quad[V(t),V(t^{\prime})\rho(t^{\prime})-\rho(t^{\prime})V(t^{\prime})] \\ & =\mathsf{T}r_b\int_0^tdt^{\prime}\left(\quad V(t)V(t^{\prime})\rho(t^{\prime})-V(t)\rho(t^{\prime})V(t^{\prime})-V(t^{\prime})\rho(t^{\prime})V(t)+\rho(t^{\prime})V(t^{\prime})V(t)\right) \\ & =\mathrm{Tr}_b\int_0^tdt^{\prime}\left(\quad\underbrace{V(t)V(t^{\prime})\rho(t^{\prime})}_{\large1}-\underbrace{V(t)\rho(t^{\prime})V(t^{\prime})}_{\large2}-\underbrace{V(t^{\prime})\rho(t^{\prime})V(t)}_{\large3}+\underbrace{\rho(t^{\prime})V(t^{\prime})V(t)}_{\large4}\right) \end{aligned}$

其中

$\begin{aligned} V(t)&=\kappa\int d\omega b^{\dagger}e^{i\omega t}ae^{-i\omega_{s}t}+be^{-i\omega t}a^{\dagger}e^{i\omega_{s}t} \\[1em] &=\kappa\int d\omega b^{\dagger}ae^{i(\omega-\omega_{s})t}+ba^{\dagger}e^{-i(\omega-\omega_{s})t} \end{aligned}$

第一项 $①V(t)V(t^{\prime})\rho(t^{\prime})$

我们将逐项计算上式中的四个部分，首先计算第一项

$\begin{aligned} & \mathrm{Tr}_{b}\int_{0}^{t}dt^{\prime}V(t)V(t^{\prime})\rho(t^{\prime}) \\[1em] & =\kappa^{2}\cdot\mathrm{Tr}_{b}\int_{0}^{t}dt^{\prime}\sum_{\omega}[b^{\dagger}ae^{i(\omega-\omega_{s})t}+ba^{\dagger}e^{-i(\omega-\omega_{s})t}]\sum_{\omega}^{\prime}[b^{\dagger}ae^{i(\omega^{\prime}-\omega_{s})t^{\prime}}+ba^{\dagger}e^{-i(\omega^{\prime}-\omega_{s})t^{\prime}}]\rho(t^{\prime}) \\[1em] & =\kappa^{2}\mathrm{Tr}_{b}\int_{0}^{t}dt^{\prime}\sum_{\omega}\sum_{\omega}^{\prime}[b^{\dagger}ab^{\dagger}ae^{i(\omega t+\omega^{\prime}t^{\prime})-i\omega_{s}(t+t^{\prime})}+b^{\dagger}aba^{\dagger}e^{i(\omega t-\omega^{\prime}t^{\prime})-i\omega_{s}(t-t^{\prime})}\\[1em] &\qquad\qquad +ba^\dagger b^\dagger ae^{i(-\omega t+\omega^{\prime}t^{\prime})+i\omega_s(t-t^{\prime})}+ba^\dagger ba^\dagger e^{i(-\omega t-\omega^{\prime}t^{\prime})+i\omega_s(t+t^{\prime})}]\rho_s(t^{\prime})\otimes\rho_b \end{aligned}$

上式包含了系统和热库算符的四种乘积。为了计算第 ① 项，我们必须对热库求迹，这将会涉及到热库算符对的期望，例如 $\mathrm{Tr}_b\left(b^\dagger(\omega)b^\dagger(\omega^\prime)\rho_b\right),\mathrm{Tr}_b\left(b^\dagger(\omega)b(\omega^\prime)\rho_b\right)$ ，等等。

第一项就这么复杂？别急，我们有近似大法。我们应用一个标准热库来简化这个表达式，核心猜想就是：热库处于静态（例如真空态或热平衡态），没有压缩态这样的关联。

假设 $4$ ：热库性质

对于一个热库

$\mathrm{Tr}_b(b^\dagger b^\dagger\rho_b)=0\quad\mathrm{and}\quad\mathrm{Tr}_b(bb\rho_b)=0$

这个条件适用于典型的热库，例如真空电磁场或热声子浴，但是对于压缩态或其他非经典态的热库，但对于压缩态热库，这些条件则不成立，因为其反常关联为非零值。光子数由下式给出

$\mathrm{Tr}_b\left(b^\dagger(\omega)b(\omega)\rho_b\right)=N(\omega)$

或者以更一般的形式表示

$\mathrm{Tr}_b\left(b^\dagger(\omega)b(\omega^{\prime})\rho_b\right)=N(\omega)\delta_{\omega\omega^{\prime}}$

应用上述假设至第 ① 项，大多数项都为零，只保留了涉及 $b^{\dagger}b$ 和 $bb^{\dagger}$ 的交叉项

$\begin{aligned} & \mathrm{Tr}_b\int_0^tdt^{\prime}V(t)V(t^{\prime})\rho(t^{\prime}) \\[1em] & =\kappa^2\mathrm{Tr}_b\int_0^tdt^{\prime}\sum_\omega\sum_{\omega^{\prime}}[b^\dagger aba^\dagger e^{i(\omega t-\omega^{\prime}t^{\prime})-i\omega_s(t-t^{\prime})}+ba^\dagger b^\dagger ae^{i(-\omega t+\omega^{\prime}t^{\prime})+i\omega_s(t-t^{\prime})}]\rho(t^{\prime})_s\otimes\rho_b \\[1em] & =\kappa^{2}\cdot\sum_{\omega}\sum_{\omega^{\prime}}\mathrm{Tr}_{b}\left(b^{\dagger}b\rho_{b}\right)\int_{0}^{t}dt^{\prime}e^{i(\omega t-\omega^{\prime}t^{\prime})-i\omega_{s}(t-t^{\prime})}aa^{\dagger}\rho_{s} \\[1em] &\quad +\kappa^{2}\cdot\sum_{\omega}\sum_{\omega^{\prime}}\mathrm{Tr}_{b}\left(bb^{\dagger}\rho_{b}\right)\int_{0}^{t}dt^{\prime}e^{i(-\omega t+\omega^{\prime}t^{\prime})+i\omega_{s}(t-t^{\prime})}a^{\dagger}a\rho_{s}\\[1em] & =\kappa^{2}\cdot\int_{0}^{t}dt^{\prime}\sum_{\omega}N(\omega)e^{i(\omega-\omega_{s})(t-t^{\prime})}aa^{\dagger}\rho_{s}(t^{\prime}) \\[1em] &\quad +\kappa^{2}\cdot\int_{0}^{t}dt^{\prime}\sum_{\omega}(N(\omega)+1)e^{-i(\omega-\omega_{s})(t-t^{\prime})}a^{\dagger}a\rho_{s}(t^{\prime}) \end{aligned}$

当热库模和系统频率失谐时，时间依赖只出现在振荡的相位因子 $e^{\pm i(\omega-\omega_s)(t-t^{\prime})}$ 中。

在实际物理系统中，热库通常具有连续模，也就是说，一个给定的 $\omega$ 有许多个空间模式，求和 $\sum_\omega$ 实际上是 $\sum_{\omega(\mathbf{k})}$ ，其中 $\mathbf{k}$ 是不同空间模式的三维指标。为了建立模型，我们模式的离散和转变为连续频率分布的积分，这就涉及到两个关键变化，首先需要用积分替换求和

$\sum_\omega\cdots\to\int_{-\infty}^\infty d\omega D(\omega)\ldots$

其中态密度 $D(\omega)$ 频率 $\omega$ 处单位频率范围内的热库模式数，将这个替换应用到上式中，得到

$\begin{aligned} & \mathrm{Tr}_b\int_0^tdt^{\prime}V(t)V(t^{\prime})\rho(t^{\prime}) \\[1em] & =\kappa^{2}\cdot\int_{0}^{t}dt^{\prime}\int D(\omega)d\omega N(\omega)e^{i(\omega-\omega_{s})(t-t^{\prime})}aa^{\dagger}\rho_{s}(t^{\prime}) \\[1em] &\quad +\kappa^{2}\cdot\int_{0}^{t}dt^{\prime}\int D(\omega)d\omega(N(\omega)+1)e^{-i(\omega-\omega_{s})(t-t^{\prime})}a^{\dagger}a\rho_{s}(t^{\prime}) \end{aligned}$

对于大部分实际物理系统， $D(\omega)N(\omega)$ 在系统频率 $\omega_s$ 附近近似于常数，我们现在引入一个基本近似，这是推导马尔可夫主方程的核心

假设 $5$ ：缓变（马尔可夫）近似

在系统频率 $\omega_s$ 附近，态密度 $D(\omega)$ 和热占据数 $N(\omega)$ 随频率变化缓慢。

驻相法是一种非常有效的积分近似方法，用于渐进近似以下形式积分

$I(\omega)=\int_a^bg(t)e^{i\omega f(t)}dt$

其中 $\omega$ 是一个大的实数参数， $f(t)$ 是一个实值相位函数， $g(t)$ 是一个缓慢变化的振幅函数。其核心思想：当 $\omega\gg1$ 时，被积函数 $e^{i\omega f(t)}$ 会迅速振荡，导致正负贡献相互抵消，这使得在大部分区间上积分的净贡献非常小，积分的主要贡献来自于相位平稳（即振荡变慢）的点： $f'(t_s) = 0$ ，这个点称为稳相点，点 $t_s$ 被称为驻点，该方法通过在 $t_s$ 附近展开函数来近似计算积分。最后得到的近似结果为

$I(\omega)\approx g(t_s)\sqrt{\frac{2\pi i}{\omega f^{\prime\prime}(t_s)}}e^{i\omega f(t_s)}$

缓变近似（或马尔可夫近似）非常适用于像原子与电磁场相互作用这样的系统，因为在这些系统中，相关带宽（约等于原子线宽）远窄于光学态密度或真空拉比频率发生显著变化的尺度。

在此近似下，常数 $D(\omega_s)$ 和 $N(\omega_s)$ 被提取到积分号之外；剩下的任务是计算频率和时间的积分，这将决定主方程中最终的耗散率。

$\begin{aligned} & \mathrm{Tr}_b\int_0^tdt^{\prime}V(t)V(t^{\prime})\rho(t^{\prime}) \\ & =\kappa^{2}\cdot\int_{0}^{t}dt^{\prime}D(\omega_{s})N(\omega_{s})\int d\omega e^{i(\omega-\omega_{s})(t-t^{\prime})}aa^{\dagger}\rho_{s}(t^{\prime}) \\ &\quad +\kappa^{2}\cdot\int_{0}^{t}dt^{\prime}D(\omega_{s})(N(\omega_{s})+1)\int d\omega e^{-i(\omega-\omega_{s})(t-t^{\prime})}a^{\dagger}a\rho_{s}(t^{\prime}) \\ & =\kappa^{2}\cdot\int_{0}^{t}dt^{\prime}D(\omega_{s})N(\omega_{s})2\pi\delta(t-t^{\prime})aa^{\dagger}\rho_{s}(t^{\prime}) \\ &\quad +\kappa^{2}\cdot\int_{0}^{t}dt^{\prime}D(\omega_{s})(N(\omega_{s})+1)2\pi\delta(t-t^{\prime})a^{\dagger}a\rho_{s}(t^{\prime}) \end{aligned}$

如果定义

$\begin{aligned} \frac{\gamma}{2\pi}=\kappa^2D(\omega_s) \end{aligned}$

那么我们总算得到第 ① 项的最终结果

$\begin{aligned} \mathrm{Tr}_b\int_0^tdt^{\prime}V(t)V(t^{\prime})\rho(t^{\prime}) = \frac{\gamma}{2}\cdot N\cdot aa^\dagger\rho_s+\frac{\gamma}{2}\cdot(N+1)\cdot a^\dagger a\rho_s \end{aligned}$

其中 $N$ 为系统频率 $\omega_s$ 处热库模的光子数密度。这一结果代表了热库相互作用对系统密度矩阵演化的二阶贡献，其中包含耗散项 $N a a^\dagger \rho_s$ 和涨落项 $(N+1) a^\dagger a \rho_s$ ，这些项是由处于腔频率 $\omega_s$ 的热库驱动的。

第二项 $②V(t)\rho(t^{\prime})V(t^{\prime})$

还在还只是计算完四项中的第一项，接下来我们计算第 ② 项

我们现在利用 $V(t)$ 的显式形式，并对热库的自由度进行求迹运算，来计算（标记为 ② 的）第二项，利用相同的近似与方法，同理我们得到

$\begin{aligned} & \mathsf{Tr}_b\int_0^tdt^{\prime}V(t)\rho(t^{\prime})V(t^{\prime}) \\ & =\kappa^2\mathrm{Tr}_b\int_0^tdt^{\prime}\int d\omega\left[b^\dagger ae^{i(\omega-\omega_s)t}+ba^\dagger e^{-i(\omega-\omega_s)t}\right]\rho_s(t^\prime)\otimes\rho_b\int d\omega^\prime[b^\dagger ae^{i(\omega^\prime-\omega_s)t^\prime}+ba^\dagger e^{-i(\omega^\prime-\omega_s)t} \\[1em] & =\kappa^{2}\mathrm{Tr}_{b}\int_{0}^{t}dt^{\prime}\sum_{\omega}\sum_{\omega}[b^{\dagger}\rho_{b}b^{\dagger}\cdot a\rho_{s}a\cdot e^{i(\omega t+\omega^{\prime}t^{\prime})-i\omega_{s}(t+t^{\prime})}+b^{\dagger}\rho_{b}b\cdot a\rho_{s}a^{\dagger}\cdot e^{i(\omega t-\omega^{\prime}t^{\prime})-i\omega_{s}(t-t^{\prime})} \\[1em] &\qquad\qquad\qquad\qquad\quad +b\rho_{b}b^{\dagger}\cdot a^{\dagger}\rho_{s}a\cdot e^{-i(\omega t-\omega^{\prime}t^{\prime})+i\omega_{s}(t-t^{\prime})}+b\rho_{b}b\cdot a^{\dagger}\rho_{s}a^{\dagger}\cdot e^{-i(\omega t+\omega^{\prime}t^{\prime})+i\omega_{s}(t+t^{\prime})}]\\[1em] &=\kappa^{2}\int_{0}^{t}dt^{\prime}\sum_{\omega}\sum_{\omega^{\prime}}[0+\mathrm{Tr}_{b}(b^{\dagger}\rho_{b}b)a\rho_{s}a^{\dagger}\cdot e^{i(\omega-\omega_{s})(t-t^{\prime})}+\mathrm{Tr}_{b}(b\rho_{b}b^{\dagger})a^{\dagger}\rho_{s}a\cdot e^{-i(\omega-\omega_{s})(t-t^{\prime})}+0] \\[1em] &=\kappa^{2}\int_{0}^{t}dt^{\prime}\sum_{\omega}\sum_{\omega^{\prime}}[0+\mathrm{Tr}_{b}(bb^{\dagger}\rho_{b})a\rho_{s}a^{\dagger}\cdot e^{i(\omega-\omega_{s})(t-t^{\prime})}+\mathrm{Tr}_{b}(b^{\dagger}b\rho_{b})a^{\dagger}\rho_{s}a\cdot e^{-i(\omega-\omega_{s})(t-t^{\prime})}+0]\\[1em] &=\frac{\gamma}{2}\cdot(N+1)\cdot a\rho_sa^\dagger+\frac{\gamma}{2}\cdot N\cdot a^\dagger\rho_s a \end{aligned}$

第三项 $③V(t^{\prime})\rho(t^{\prime})V(t)$

第三项与第二项非常相似，只是 $V(t)$ 和 $V(t^{\prime})$ 的顺序交换了，类似地

$\begin{aligned} & \mathrm{Tr}_b\int_0^tdt^{\prime}\quad V(t^{\prime})\rho(t^{\prime})V(t) \\[1em] & =\cdots\quad[b^{\dagger}a\cdots+ba^{\dagger}\cdots]\rho_{s}(t^{\prime})\otimes\rho_{b}[b^{\dagger}a\cdots+ba^{\dagger}\cdots] \\[1em] & =\cdots[b^{\dagger}\rho_{b}b^{\dagger}\cdot a\rho_{s}a\cdots+b^{\dagger}\rho_{b}b\cdot a\rho_{s}a^{\dagger}\cdots+b\rho_{b}b^{\dagger}\cdot a^{\dagger}\rho_{s}a\cdots+b \\[1em] & =0+\frac{\gamma}{2}\cdot(N+1)\cdot a\rho_{s}a^{\dagger}+\frac{\gamma}{2}\cdot N\cdot a^{\dagger}\rho_{s}a+0 \end{aligned}$

第四项 $④\rho(t^{\prime})V(t^{\prime})V(t)$

同理可得

$\begin{aligned} & \mathrm{Tr}_b\int_0^tdt^{\prime}\quad\rho(t^{\prime})V(t^{\prime})V(t) \\[1em] & =\cdots\quad\rho_{s}(t^{\prime})\otimes\rho_{b}\quad[b^{\dagger}a\cdots+ba^{\dagger}\cdots]\quad[b^{\dagger}a\cdots+ba^{\dagger}\cdots] \\[1em] & =\cdots[\rho_{b}b^{\dagger}b^{\dagger}\cdot\rho_{s}aa\cdots+\rho_{b}b^{\dagger}b\cdot\rho_{s}aa^{\dagger}\cdots+\rho_{b}bb^{\dagger}\cdot\rho_{s}a^{\dagger}a\cdots+\rho_{b}bb\cdot\rho_{s}a^{\dagger}a^{\dagger}\cdots]\\[1em] &=0+\frac{\gamma}{2}\cdot N\cdot\rho_saa^\dagger+\frac{\gamma}{2}\cdot(N+1)\cdot\rho_sa^\dagger a+0 \end{aligned}$

终于，经过计算并对热库进行求迹后，总算得到四个单独的项结果如下

$\begin{gathered} V(t)V(t^{\prime})\rho(t^{\prime})\to\frac{\gamma}{2}\cdot N\cdot aa^{\dagger}\rho_{s}+\frac{\gamma}{2}\cdot(N+1)\cdot a^{\dagger}a\rho_{s} \\[1em] V(t)\rho(t^{\prime})V(t^{\prime})\to\frac{\gamma}{2}\cdot(N+1)\cdot a\rho_{s}a^{\dagger}+\frac{\gamma}{2}\cdot N\cdot a^{\dagger}\rho_{s}a \\[1em] V(t^{\prime})\rho(t^{\prime})V(t)\to\frac{\gamma}{2}\cdot(N+1)\cdot a\rho_{s}a^{\dagger}+\frac{\gamma}{2}\cdot N\cdot a^{\dagger}\rho_{s}a \\[1em] \rho(t^{\prime})V(t^{\prime})V(t)\to\frac{\gamma}{2}\cdot N\cdot\rho_{s}aa^{\dagger}+\frac{\gamma}{2}\cdot(N+1)\cdot\rho_{s}a^{\dagger}a \end{gathered}$

其中 $N$ 表示 $N(\omega_s)$ ，即位于腔频率 $\omega_s$ 处的热库模式的平均光子数。我们四项的结果代入到密度矩阵运动方程的第二项（红色部分）中，得到

$\begin{aligned} & \mathrm{Tr}_{b}\int_{0}^{t}dt^{\prime}\left[V(t),[V(t^{\prime}),\rho(t^{\prime})]\right] \\[1em] & =\frac{\gamma}{2}N\cdot(aa^{\dagger}\rho_{s}-2a^{\dagger}\rho_{s}a+\rho_s aa^{\dagger})+\frac{\gamma}{2}(N+1)\cdot(a^{\dagger}a\rho_{s}-2a^{\dagger}\rho_{s}a+\rho_s a^{\dagger}a) \end{aligned}$

旋转框架下的系统约化密度矩阵 $\rho_s(t)$ 的演化受主方程控制

$\begin{aligned} \frac{d}{dt}\rho_{s}(t) & =-i\mathrm{Tr}_b[V(t),\rho(0)]-\mathrm{Tr}_b\int_0^tdt^{\prime}\left[V(t),\quad[V(t^{\prime}),\rho(t^{\prime})]\right] \\[1em] & =\frac{\gamma}{2}N\cdot(2a^{\dagger}\rho_{s}a-aa^{\dagger}\rho_{s}-\rho_s aa^{\dagger})+\frac{\gamma}{2}(N+1)\cdot(2a\rho_{s}a^{\dagger}-a^{\dagger}a\rho_{s}-\rho_s a^{\dagger}a) \end{aligned}$