Wick’s theorem is a powerful mathematical tool in quantum field theory, providing a systematic method for evaluating time-ordered products of field operators [1, 2]. By expressing such products in terms of normal-ordered components and contractions, the theorem enables the derivation of Feynman rules and facilitates the computation of scattering amplitudes [3–7]. In finite-temperature field theory, Wick’s theorem plays a central role in computing thermal Green’s functions, which are essential for analyzing equilibrium and nonequilibrium thermodynamic properties [8–11]. Furthermore, in statistical mechanics and condensed matter systems, it is routinely employed to decompose multi-point correlation functions into sums of products of two-point correlation functions, thereby enabling analytical progress in the study of interacting many-body systems [12–17]. Further applications include the evaluation of quantum coherence [18], the characterization of order parameters [19], and the analysis of spin correlation functions [20].

Beyond the standard formulation, Wick’s theorem for superoperators has been proposed in Liouville space, avoiding the need for backward propagations and analytic continuations using artificial times required in Hilbert space [21]. Wick’s theorem for superoperators allows the factorization of higher-order response functions into products of two fundamental Green’s functions in Liouville space [21]. This method motivates us to propose a derivation of Wick’s theorem for superoperators in the study of reduced dynamics which involves higher-order superoperators based on operator algebra.

In this work, we prove Wick’s theorem for superoperators that are linear in bosonic operators rigorously. This theorem enables the exact solution of the reduced dynamics for systems interacting with a bosonic bath. Our approach employs superoperators labeled by upper indices to explicitly record whether system and bath operators multiply the density matrix from the left or the right, thus making their relative ordering with respect to the density matrix manifest. We apply the stable partition theory, which rearranges a sequence while preserving the relative ordering of elements with identical labels [22]. Using this framework, we show that any higher-order time-ordered superoperator term can be decomposed into sums of products of second-order terms. We derive explicit expressions for all even-order terms of the reduced density matrix, while all odd-order terms vanish. Higher-order terms of reduced density matrix can be expressed in terms of second-order terms. Consequently, the reduced dynamics can be constructed by summing these even-order contributions, yielding a formally exact expression.

In section 2, we study the reduced dynamics of the system coupled to a bosonic bath and highlight the challenge of evaluating time-ordered higher-order superoperator terms. In section 3, we prove Wick’s theorem for superoperators rigorously by applying a stable partition theory and demonstrate how higher-order time-ordered superoperators can be systematically decomposed into second-order terms. In section 4, we apply Wick’s theorem for superoperators, leading to an exact expression for the reduced density matrix. In section 5, we make our conclusions.

We consider the system interacting with a bosonic bath, described by the total Hamiltonian where $H_{S}$ and $H_{B}$ are the free Hamiltonians for the system and bath, respectively. The system-bath interaction Hamiltonian $H_{SB}$ in the interaction picture becomes and the total density matrix in the interaction picture is where the superoperator is defined by $A^{\times}B\equiv\left[A,\, B\right]$ and ${\rho}_{\text{tot}}$ is the total density matrix in the Schrödinger picture. The time evolution of $\tilde{\rho}_{\text{tot}}\left(t\right)$ follows the equation Therefore, $\tilde{\rho}_{\text{tot}}\left(t\right)$ can be written as where $\mathcal{T}$ is the time-ordering operator.

In general, we assume a factorized initial state of the system and bath, $\rho_{ \text{tot}}\left(0\right) = \rho_{S}\left(0\right)\rho_{B}\left(0\right)$. Substituting this into equation (5) and tracing out the bath degrees of freedom, $ \tilde{\rho}_{S}\left(t\right) = \text{Tr}_{B}\left[\tilde{\rho}_{\text{tot}}\left(t\right)\right] $, we obtain the reduced density matrix as follows: where the $n$th order contribution $\varrho_{n}$ is given by The system-bath interaction is bilinear, where $V\left(t\right)$ and $B\left(t\right)$ are the system and bath operators, respectively. The system operator $V\left(t\right)$ could be generalized coordinates or spin operators. Below, we take the bath operator $B\left(t\right)$ to have the form where $a_{k}$ and $a_{k}^{\dagger}$ are the annihilation and creation bosonic operators in mode $k$. $\xi_{k}(t)$ is the complex coefficient that determines the contribution of each mode $k$ to $B(t)$. To perform the partial trace of the bath, we assume the bath is in thermal equilibrium. The density matrix of the bath is $\rho_{B} = \frac{\exp\left(-\beta H_{B}\right)}{\text{Tr} \exp\left(-\beta H_{B}\right)},$ where $ H_B = \sum_{k}\omega _{k}a_{k}^{\dagger }a_{k}$. Since $\text{Tr}\left[a\rho_{B}\right] = \text{Tr}\left[a^\dagger\rho_{B}\right] = 0$ and more generally $\text{Tr}\left[(a+a^\dagger)^{2n+1}\rho_{B}\right] = 0$, we conclude that all odd-order terms $\varrho_{2n+1}(t) = 0$. Thus, we are only required to compute the even-order terms $\varrho_{2n}(t)$.

To express the commutation structure resulting from the form of the interaction Hamiltonian, we introduce the following superoperators: Then $H_{SB}^{\times}$ acting on the density matrix can be written as Substituting these definitions into equation (7), the second-order contribution $\varrho_2(t)$ becomes where $\mathcal{T}_{S}$ and $\mathcal{T}_{B}$ are the time-ordering operators for the system and bath, respectively. Instead of integrating over the full domain $0 \unicode{x2A7D} \tau_1, \tau_2 \unicode{x2A7D} t$ with the time-ordering operator, which requires the use of Heaviside step functions to enforce proper ordering, we simplify the calculation by restricting the integration to the ordered domain $0 \unicode{x2A7D} \tau_1 \unicode{x2A7D} \tau_2 \unicode{x2A7D} t$ after removing the time-ordering operator. Thus, we have For convenience, we introduce the symmetric correlation function $G\left( \tau \right) $, and the response function $ \chi \left( \tau \right) $, given by Indeed, $G$ and $\chi $ are related to the real and imaginary parts of the bath time correlation function $C\left( \tau _{2}-\tau _{1}\right) = \left\langle B_{2}B_{1}\right\rangle _{B}. $ Then where $V^{\circ }\rho = \left\{V,\,\rho \right\} $. Compared to second-order term $\varrho _{2}$, higher-order terms are extremely difficult to express using symmetric correlation functions and response functions directly. In the next section, we propose a method to simplify higher-order terms.

Since all odd-order contributions vanish, we focus exclusively on even-order terms. This naturally suggests applying Wick’s theorem, which is commonly employed to evaluate expectation values of time-ordered even numbers of operators [1]. However, the conventional form of Wick’s theorem, which applies to operator products, is not applicable to the superoperators involved in our formulation.

To overcome this limitation, we propose Wick’s theorem for superoperators that simplifies the analysis of higher-order contributions. This theorem states that the bath trace of any time-ordered product of $2n$ superoperators can be expressed as a sum of products over all pairings of the $2n$ superoperators into $n$ pairs, with each pair taken in internal time order, where $m = \frac{(2n)!}{2^nn!}$ is the total number of such pairings. Following the left-hand side of equation (16), we have the base index sequence Here $Q_{l,j}$ (for $l = 1,2,3,4$ and $j = 1,\dots,m$) denotes the $j$th way of pairing up all $2n$ elements of the $l$th sequence $\mathrm{SEQ}_l$ into $n$ pairs $(u_l, v_l)$. In the theorem we only consider the case $l = 1$. The indices $l = 2,3,4$ are used only to describe the sequences $\mathrm{SEQ}_l$ that appear in the proof. For simplicity, we write $u,v$ in place of $u_1, v_1$. The case $2n = 4$ is provided as an example, There are $m = 3$ ways of picking pairs, We now consider the case of $2n$ elements. Without loss of generality, one way to pair up all elements in this sequence is After performing the time ordering, the left side of equation (16) can be expressed as where $\sigma$ is a permutation that sorts the time labels in nonincreasing order: Equation (21) involves $2n$ superoperators, among which $p$ carries upper index 2 and the remaining $d = 2n-p$ carry upper index 1. The time-ordering operator arranges these superoperators chronologically, resulting in a sequence where $L^{(1)}$ and $L^{(2)}$ terms are mixed. So from equation (21), after the action of the superoperators on the density matrix, whether each $B_\sigma$ is on the right side of $\rho$ or on the left is obscured.

To address this, we gather superoperators with the same upper indices together while keeping their relative order unchanged. This rearrangement is justified by the superoperator’s property: superoperators with different upper indices commute, whereas those with the same upper indices do not.

Formally, this process corresponds to a stable partition of the superoperator sequence based on the upper indices. A stable partition algorithm is presented in appendix A. Without a stable partition, the proof of Wick’s theorem for superoperators can hardly be presented clearly. Stable partition groups terms with the same index into contiguous blocks while preserving the original relative order within each group. In this way, we introduce a canonical reordering of a general sequence under a partial commutativity constraint: commuting factors can be permuted to form contiguous index blocks, while the relative order of non-commuting factors is preserved.

With these explicit upper indices, we can then determine the ordering of the operators resulting from the action of the superoperators in the subsequent calculation. Specifically, the sequence $\mathrm{SEQ}_1$ can be partitioned into two subsequences, $\mathrm{SUBSEQ}_1$ and $\mathrm{SUBSEQ}_2$, where Here $k_{j}$ (respectively $k_{j}^{^{\prime}}$) denotes the index in the permutation $\sigma$ such that the superoperator $L_{B_{\sigma(k_{j})}}^{i_{\sigma(k_{j})}}$ $\left({\text{respectively}}L_{B_{\sigma(k^{^{\prime}}_{j})}}^{i_{\sigma(k^{^{\prime}}_{j})}}\right)$ carries upper index 1 (respectively 2). Under this reordering, equation (21) can be expressed as This allows us to express the corresponding expression in terms of operators rather than superoperators. When these superoperators act on the density matrix, we obtain Applying Wick’s theorem, we have where However, our target is to pair superoperators. By applying the cyclic property of trace again, equation (27) can be rewritten into the superoperator representation,

where Next, we consider the right-hand side of equation (16). By applying the time-ordering operator, we obtain

where After applying the stable partition operation, equation (29) is consistent with equation (31). This completes the proof of equation (16).

Having established the theorem, we now apply Wick’s theorem for superoperators to evaluate the even-order contributions $\varrho_{2n}(t)$. First, we consider the fourth-order term $\varrho_{4}(t)$. From equation (7) we have The integrand can be written as where we have identified three terms $A_1, A_2, A_3$ corresponding to the different ways of pairing the bath superoperators. These terms are given by Up to this point, we have factorized the bath superoperator contributions using Wick’s theorem for superoperators. We now take into account the system superoperators. Importantly, system superoperators can be reordered freely without altering the result because under the time-ordering operator, system superoperators commute. Therefore, the system parts of $A_1$, $A_2$, and $A_3$ can be rearranged to align with the corresponding bath superoperator orderings. Implementing this reordering, $A_1$ can be written as and $A_2,A_3$ can be written analogously. The order of integration can be exchanged to align with the superoperator order without altering the result. Therefore, the integrals for $A_1$, $A_2$, and $A_3$ are identical, since they differ only in dummy indices, By removing the time-ordering operator $\mathcal{T}_{B}$ acting on bath superoperators in $A_1$, we obtain An interesting feature of equation (40) is that the double integrals over $(\tau_{2}, \tau_{1})$ and $(\tau_{4}, \tau_{3})$ play identical roles. Equation (40) can be factorized as the square of a single double-integral expression. In fact, by using the second-order result equation (15), equation (40) can be rewritten in a compact form, where $G_{kl} = \text{Re}\left\langle B_{k}B_{l}\right\rangle ,\;\chi_{kl} = -\frac{2}{ \hbar}\text{Im}\left\langle B_{k}B_{l}\right\rangle $.

Using a similar procedure, we obtain the $2n$-th term in appendix B.For the $2n$-th term, we can express it with a power of $n$, establishing a correspondence between the term order $2n$ and the power exponent $n$. Finally, we derive the exact expression for the reduced dynamics, Our result, derived using Wick’s theorem for superoperators, is consistent with that obtained from the Feynman and Vernon influence functional method, which yields a continued fraction expression for the reduced density matrix [23–25].

We have rigorously proved Wick’s theorem for superoperators that are linear in bosonic operators that decompose higher-order time-ordered superoperators into a sum of products of second-order terms. This derivation employs stable partition theory to reorder superoperators according to their upper indices, thereby enabling a systematic factorization of higher-order contributions.

Applying Wick’s theorem for superoperators, explicit expressions for all even-order terms of the reduced density matrix are derived. We have derived a formally exact expression for the reduced density matrix of a system coupled to a bosonic bath initially in thermal equilibrium. The Wick’s theorem for superoperators offers new insights into the reduced dynamics of open quantum systems. We believe there will be other applications for the Wick’s theorem for superoperators.

A stable partition algorithm rearranges the elements of a sequence such that elements satisfying a given condition precede those that do not, while preserving the original relative order of the elements within each subsequence.For example, consider the sequence: Applying the condition that an element is odd, we identify two subsequences. The subsequence of the elements satisfying the condition is: and the subsequence of elements not satisfying the condition is: A stable partition of the original sequence $S$ with respect to this condition yields the sequence:

The expression for the $2n$-th order terms is Employing Wick’s theorem for superoperators and equation (11), the integrand can be written as where Due to the presence of the time-ordering operator, we can always rearrange the order of the superoperators $L_{S}^{\left( i\right) }$ corresponding to the pairing of superoperators $ L_{B}^{\left( i\right)} $, we have Since the $A_{j} ( 1\unicode{x2A7D} j \unicode{x2A7D} m)$ differ only by dummy indices, they become identical after integration, Therefore, $\varrho _{2n}$ can be derived as Removing the time-ordering operator $\mathcal{T}_B$ acting on the bath superoperator, we have