Here, $\Omega(F)$ is the set of all unique eigenvalues, and $\hat P(f)$ are the projector operators [i.e., $\hat P(f)\hat P(f') = \delta_{f,f'}\hat P(f)$] onto the degenerate subspaces of $\hat F$, i.e., the subspaces spanned by the eigenkets of $\hat F$ that correspond to the same eigenvalue $f$. Of course, partitioning the basis of eigenkets into degenerate subspaces does not interfere with the decomposition of identity,

\begin{align}
\label{eq:subspace_identity}
\sum_{f\in\Omega(F)}\hat P(f) = \sum_n\left(\sum_{f\in\Omega(F)}\delta_{f,f(n)}\right)|n\rangle\langle n| = \sum_n |n\rangle\langle n| = \hat{1}.
\end{align}

The dynamical map can then be written as a quasi-moment series,

\begin{align}
\nonumber
\langle\mathcal{U}\rangle_I(t,0) &= \bullet + \sum_{k=1}^\infty (-i\lambda)^k \int_0^tds_1\cdots\int_0^{s_{k-1}}ds_k\\
\nonumber
\!\!\sum_{f_1,\bar{f}_1\in\Omega(F)}\!\!\cdots\!\!\sum_{f_k,\bar{f}_k\in\Omega(F)}\!\! Q_F^{(k)}(f_1,\bar{f}_1,s_1;\ldots;f_k,\bar{f}_k,s_k)\,
\mathcal{W}(f_1,\bar{f}_1,s_1)\cdots\mathcal{W}(f_k,\bar{f}_k,s_k)\\
\nonumber
&\equiv \sum_{k=0}^\infty(-i\lambda)^k \left\langle
\int\limits_0^t\!\!ds_1\mathcal{W}(F(s_1),\bar{F}(s_1),s_1)\cdots\!\!\!\int\limits_0^{s_{k-1}}\!\!\!\!ds_k\mathcal{W}(F(s_k),\bar{F}(s_k),s_k)
\right\rangle\\
\label{eq:Q_rep}
&=\left\langle T_\mathcal{W}e^{-i\lambda\int_0^t\mathcal{W}(F(s),\bar{F}(s),s)ds}\right\rangle,
\end{align}

where the super-operator in $S$-operator space is defined as

\begin{align}
\label{eq:super-W}
\mathcal{W}(f,\bar{f},s) &= f\hat V_I(s)\bullet - \bullet\hat V_I(s)\bar{f},
\end{align}

and the joint quasiprobability distribution of order $k$ is given by

\begin{align}
\label{eq:Q_def}
Q_F^{(k)}(f_1,\bar{f}_1,s_1;\ldots;f_k,\bar{f}_k,s_k)&= \operatorname{tr}_E\left[
\hat P_I(f_1,s_1)\cdots \hat P_I(f_k,s_k)\ \hat\rho_E\
\hat P_I(\bar{f}_k,s_k)\cdots \hat P_I(\bar{f}_1,s_1)\right].
\end{align}

The subscript $F$ indicates that the quasiprobability is associated with the coupling operator $\hat F_I(s)$.

[[sol:quasi-moment_series]Solution] Prove the joint quasiprobability representation of the dynamical map (<ref>).

To hammer home the analogy with stochastic dynamics, compare the quasiprobability representation with the moment series for $\overline{\mathcal{U}}_I$,

\begin{align}
\nonumber
\overline{\mathcal{U}}_I(t,0) &=\bullet +\sum_{k=1}^\infty(-i\lambda)^k\int_0^t\!\!ds_1\cdots\int_0^{s_{k-1}}\!\!\!\!ds_k\sum_{f_1}\cdots\sum_{f_k}\\
\nonumber
&\phantom{=}\times P_F^{(k)}(f_1,s_1;\ldots;f_k,s_k)[f_1\hat V_I(s_1),\bullet]\cdots[f_k\hat V_I(s_k),\bullet]\\
\nonumber
&= \sum_{k=0}^\infty (-i\lambda)^k \overline{
\int\limits_0^t\!\!ds_1 [F(s_1)\hat V_I(s_1),\bullet]\cdots\!\!\!\int\limits_0^{s_{k-1}}\!\!\!\!ds_k [F(s_k)\hat V_I(s_k),\bullet]
\nonumber
&=\sum_{k=0}^\infty (-i\lambda)^k \overline{
\int\limits_0^t\!\!ds_1 \mathcal{W}(F(s_1),F(s_1),s_1)\cdots\!\!\!\int\limits_0^{s_{k-1}}\!\!\!\!ds_k \mathcal{W}(F(s_k),F(s_k),s_k)
\end{align}

where we assumed that process $F(t)$ in the stochastic Hamiltonian is discrete (sums instead of integrals over $f_i$) and we have used the fact that

\begin{align*}
[f\hat V_I(s),\bullet] &= \mathcal{W}(f,f,s).
\end{align*}

The structure of the moment and quasi-moment series are essentially identical; formally, they differ by the method of averaging the exp: we have the quasi-average $\langle \ldots \rangle$ for dynamical map and the stochastic average $\overline{(\ldots)}$ for stochastic map. In fact, if we were to substitute an ordinary probability distributions for joint quasiprobability distributions,

\begin{align*}
Q_F^{(k)}(f_1,\bar{f}_1,s_1;\ldots) \to \delta_{f_1,\bar{f}_1}\cdots\delta_{f_k,\bar{f}_k} P_F^{(k)}(f_1,s_1;\ldots),
\end{align*}

then the dynamical map would simply turn into the stochastic map.

The members of the family of joint quasiprobability distributions $\{Q_F^{(k)}\}_{k=1}^\infty$, analogously to joint probability distributions, are all related to each other via their version of Chapman-Kolmogorov equation (the consistency relation),

\begin{align}\label{eq:Q_consistency}
\sum_{f_i,\bar{f}_i}Q_F^{(k)}(f_1,\bar{f}_1,s_1;\ldots;f_k,\bar{f}_k,s_k) = Q_F^{(k-1)}(\ldots;f_{i-1},\bar{f}_{i-1},s_{i-1};f_{i+1},\bar{f}_{i+1},s_{i+1};\ldots),
\end{align}

for $k>1$ and $\sum_{f,\bar{f}}Q_F^{(1)}(f,\bar{f},s) = 1$ for $k=1$.

[[sol:Q_consistency]Solution] Prove consistency relation (<ref>).

Find the explicit form of quasi-moments (assume $t>s$):

\begin{align*}
\langle F(t) \rangle &= \sum_{f_1,\bar{f}_1}Q_F^{(1)}(f_1,\bar{f}_1,t)f_1;\\
\langle \bar{F}(t) \rangle &= \sum_{f_1,\bar{f}_1}Q_F^{(1)}(f_1,\bar{f}_1,t)\bar{f}_1;\\
\langle F(t)F(s) \rangle&= \sum_{f_1,\bar{f}_1}\sum_{f_2,\bar{f}_2}Q_F^{(2)}(f_1,\bar{f}_1,t;f_2,\bar{f}_2,s)f_1 f_2;\\
\langle F(t)\bar{F}(s) \rangle &=\sum_{f_1,\bar{f}_1}\sum_{f_2,\bar{f}_2}Q_F^{(2)}(f_1,\bar{f}_1,t;f_2,\bar{f}_2,s)f_1 \bar{f}_2;\\
\langle \bar{F}(t)F(s) \rangle &=\sum_{f_1,\bar{f}_1}\sum_{f_2,\bar{f}_2}Q_F^{(2)}(f_1,\bar{f}_1,t;f_2,\bar{f}_2,s)\bar{f}_1 f_2;\\
\langle \bar{F}(t)\bar{F}(s)\rangle &= \sum_{f_1,\bar{f}_1}\sum_{f_2,\bar{f}_2}Q_F^{(2)}(f_1,\bar{f}_1,t;f_2,\bar{f}_2,s)\bar{f}_1 \bar{f}_2.
\end{align*}

Find the explicit form of super-quasi-moments,

\begin{align*}
\langle \mathcal{W}(F(t),\bar{F}(t),t)\rangle,\quad\langle \mathcal{W}(F(t),\bar{F}(t),t)\mathcal{W}(F(s),\bar{F}(s),s) \rangle\ (t>s).
\end{align*}

§.§ Super-cumulant expansion for dynamical maps

Previously, we have shown that stochastic maps and dynamical maps are structurally analogous; we have reached this conclusion by comparing the moment and quasi-moment series representation of the respective maps. Essentially, the two maps differ by the method of averaging: the stochastic average $\overline{(\ldots)}$ over joint probability distributions $P_F^{(k)}$ that produces moments, versus the quasi-average $\langle\ldots\rangle$ over joint quasiprobability distributions $Q_F^{(k)}$ that produces quasi-moments. At this point it should not be surprising that there is also a “quasi” analogue of cumulants—quasi-cumulants, or qumulants, if you will.

First, note that the stochastic cumulants can be written as a kind of average, e.g.,

\begin{align}
\nonumber
C_F^{(2)}(s_1,s_2) &=\overline{F(s_1)F(s_2)}-\overline{F(s_1)}\ \overline{F(s_2)}\\
\nonumber
&\equiv\sum_{f_1,f_2}\tilde{P}_F^{(2)}(f_1,s_1;f_2,s_2)f_1 f_2.
\end{align}

In general, we have

\begin{align}
C_F^{(k)}(s_1,\ldots,s_k) &=\sum_{f_1,\ldots,f_k}\tilde{P}_F^{(k)}(f_1,s_1;\ldots;f_k,s_k)f_1\cdots f_k \equiv \overline{\overline{F(s_1)\cdots F(s_k)}},
\end{align}

where the symbol $\overline{\overline{(\ldots)}}$ stands for the cumulant average, and the cumulant distributions $\tilde{P}_F^{(k)}$ the cumulant average is carried out with, can be found by inspecting the the corresponding cumulants $C_F^{(k)}$ [see Eq. (<ref>)]; e.g., the first three distributions are defined as,

\begin{align}
\tilde{P}_F^{(1)}(f,s) &= P_F^{(1)}(f,s),\\
\tilde{P}_F^{(2)}(f_1,s_1;f_2,s_2) &= P_F^{(2)}(f_1,s_1;f_2,s_2) - P_F^{(1)}(f_1,s_1)P_F^{(1)}(f_2,s_2),\\
\nonumber
\tilde{P}_F^{(3)}(f_1,s_1;f_2,s_2;f_3,s_3) &= P_F^{(3)}(f_1,s_1;f_2,s_2;f_3,s_3)+ P_F^{(1)}(f_1,s_1)P_F^{(1)}(f_2,s_2)P_F^{(1)}(f_3,s_3)\\
\nonumber
\end{align}

Using the concept of cumulant average we can rewrite the cumulant expansion of the stochastic map,

\begin{align}
\nonumber
\overline{\mathcal{U}}_I(t,0) &=T_{\mathcal{V}_I}e^{\sum_{k=1}^\infty(-i\lambda)^k\int_0^tds_1\cdots\int_0^{s_{k-1}}ds_k C_F^{(k)}(s_1,\ldots,s_k)\mathcal{V}_I(s_1)\cdots\mathcal{V}_I(s_k)}\\
\nonumber
&= T_{\mathcal{V}_I}e^{\sum_{k=1}^\infty(-i\lambda)^k\int_0^tds_1\cdots\int_0^{s_{k-1}}ds_k \overline{\overline{F(s_1)\cdots F(s_k)}}\mathcal{V}_I(s_1)\cdots\mathcal{V}_I(s_k)}\\
\label{eq:cum_avg_form}
&=T_{\mathcal{W}}e^{\sum_{k=1}^\infty(-i\lambda)^k \overline{\overline{\int_0^t ds_1 \mathcal{W}(F(s_1),F(s_1),s_1)\cdots\int_0^{s_{k-1}}ds_k\mathcal{W}(F(s_k),F(s_k),s_k)}}}.
\end{align}

For clarity, this is how one should interpret the cumulant average of super-operators,

\begin{align}
\nonumber
\nonumber
\sum_{f_1,\ldots,f_k}\tilde{P}_F^{(k)}(f_1,s_1;\ldots;f_k,s_k)\mathcal{W}(f_1,f_1,s_1)\cdots\mathcal{W}(f_k,f_k,s_k)\\
\sum_{f_1,\ldots,f_k}\tilde{P}_F^{(k)}(f_1,s_1;\ldots;f_k,s_k)f_1\cdots f_k\mathcal{V}_I(s_1)\cdots\mathcal{V}_I(s_k).
\end{align}

If we now take the quasi-moment representation of a dynamical map (<ref>), and we replicate on it all the steps that previously lead from the moment series to the cumulant expansion (<ref>) (see assignment <ref> of Sec. <ref>) while replacing all stochastic averages with quasi-averages, we will inevitable end up with the qumulant expansion,

\begin{align}
&\langle \mathcal{U}\rangle_I(t,0)= T_{\mathcal{W}}e^{\sum_{k=1}^\infty(-i\lambda)^k\left\langle\!\left\langle\int_0^tds_1\mathcal{W}(F(s_1),\bar{F}(s_1),s_1)\cdots\int_0^{s_{k-1}}ds_k\mathcal{W}(F(s_k),\bar{F}(s_k),s_k)\right\rangle\!\right\rangle}
\end{align}

where the quasi-cumulant average (or qumulant average) $\langle\!\langle\ldots\rangle\!\rangle$ operates as follows

\begin{align}
\nonumber
\sum_{f_1,\bar{f}_1}\cdots\sum_{f_k,\bar{f}_k}\tilde{Q}_F^{(2)}(f_1,\bar{f}_1,s_1;\ldots;f_k,\bar{f}_k,s_k)\mathcal{W}(f_1,\bar{f}_1,s_1)\cdots\mathcal{W}(f_k,\bar{f}_k,s_k),
\end{align}

and the qumulant distributions $\tilde{Q}_F^{(k)}$ are constructed according the same pattern as cumulant distributions $\tilde{P}_F^{(k)}$, e.g.,

\begin{align}
\tilde{Q}_F^{(1)}(f,\bar{f},s) &= Q_F^{(1)}(f,\bar{f},s),\\
\tilde{Q}_F^{(2)}(f_1,\bar{f}_1,s_1;f_2,\bar{f}_2,s_2) &= Q_F^{(2)}(f_1,\bar{f}_1,s_1;f_2,\bar{f}_2,s_2) - Q_F^{(1)}(f_1,\bar{f}_1,s_1)Q_F^{(1)}(f_2,\bar{f}_2,s_2),
\end{align}

etc.. Essentially, the qumulant expansion is obtained form the cumulant expansion by replacing the cumulant average $\overline{\overline{(\ldots)}}$ with the qumulant average $\langle\!\langle\ldots\rangle\!\rangle$; this is just another consequence of the structural parallelism between stochastic and dynamical maps.

Find the first two qumulant distributions $\tilde{Q}_F^{(1)}$ and $\tilde{Q}_F^{(2)}$ without invoking the structural parallelism between stochastic and dynamical maps.

The main purpose behind introducing the concept of qumulants is to use them as a jump-off point for the super-qumulant expansion—the analogue of super-cumulant expansion for stochastic map. Consider now an implicit definition of super-qumulants $\{\mathcal{L}^{(k)}\}_{k=1}^\infty$,

\begin{align}
\nonumber
\langle\mathcal{U}\rangle_I(t,0) &= \left\langle T_\mathcal{W}e^{-i\lambda\int_0^t\mathcal{W}(F(s),\bar{F}(s),s)ds}\right\rangle
\equiv T_{\mathcal{L}^{(k)}}e^{\sum_{k=1}^\infty (-i\lambda)^k\int_0^t\mathcal{L}^{(k)}(s)ds}.
\end{align}

Of course, the explicit form of $\mathcal{L}^{(k)}$'s can be found in a standard way: expand the RHS into powers of $\lambda$ and compare them with the terms of corresponding order in the quasi-moment series on the LHS. However, we do not have to do all that labor! Due to structural parallelism, we can obtain $\mathcal{L}^{(k)}$'s by simply taking the super-cumulants $\mathcal{C}^{(k)}$ and replacing the cumulant averages with the qumulant versions. In assignment <ref> of Sec. <ref> you have already found first four super-cumlants,

\begin{align}
\mathcal{C}^{(1)}(s) &= C_F^{(1)}(s)\mathcal{V}_I(s) = \overline{\overline{\mathcal{W}(F(s),F(s),s)}},\\
\mathcal{C}^{(2)}(s) &= \mathcal{V}_I(s)\int_0^sC_F^{(2)}(s,u)\mathcal{V}_I(u)du = \int_0^s\overline{\overline{\mathcal{W}(F(s),F(s),s)\mathcal{W}(F(u),F(u),u)}}du,\\
\nonumber
\mathcal{C}^{(3)}(s) &= \int\limits_0^s\!\!du_1\!\!\int\limits_0^{u_1}\!\!du_2\,\overline{\overline{\mathcal{W}(F(s),F(s),s)\mathcal{W}(F(u_1),F(u_1),u_1)\mathcal{W}(F(u_2),F(u_2),u_2)}}\\
\overline{\overline{\mathcal{W}(F(s),F(s),s)\left[\overline{\mathcal{W}(F(u_1),F(u_1),u_1)},\mathcal{W}(F(u_2),F(u_2),u_2)\right]}}\\
\nonumber
\end{align}

This gives us immediately the corresponding super-qumulants,

\begin{align}
\mathcal{L}^{(1)}(s) &= \left\langle\!\left\langle\mathcal{W}(F(s),\bar{F}(s),s)\right\rangle\!\right\rangle,\\
\mathcal{L}^{(2)}(s) &= \int_0^s\left\langle\!\left\langle\mathcal{W}(F(s),\bar{F}(s),s)\mathcal{W}(F(u),\bar{F}(u),u)\right\rangle\!\right\rangle du,\\
\nonumber
\mathcal{L}^{(3)}(s) &= \int\limits_0^s\!\!du_1\!\!\int\limits_0^{u_1}\!\!du_2\,
\left\langle\!\left\langle\mathcal{W}(F(s),\bar{F}(s),s)\mathcal{W}(F(u_1),\bar{F}(u_1),u_1)\mathcal{W}(F(u_2),\bar{F}(u_2),u_2)\right\rangle\!\right\rangle\\
\label{eq:L3}
\left\langle\!\left\langle\mathcal{W}(F(s),\bar{F}(s),s)\left[\langle\mathcal{W}(F(u_1),\bar{F}(u_1),u_1)\rangle,\mathcal{W}(F(u_2),\bar{F}(u_2),u_2)\right]\right\rangle\!\right\rangle\\
\nonumber
\end{align}

[[sol:explicit_super_qum]Solution] Calculate the explicit forms of super-qumulants $\mathcal{L}^{(1)}(s)$ and $\mathcal{L}^{(2)}(s)$.

[[sol:Q_for_general_coupling]Solution] Find the quasiprobability representation of the dynamical map for general $SE$ coupling

\begin{align*}
\hat V_{SE} &= \lambda\sum_\alpha \hat V_\alpha\otimes\hat F_\alpha,
\end{align*}

where $\hat V_\alpha$ and $\hat F_\alpha$ are hermitian operators.

Note that this also includes couplings of form

\begin{align*}
\hat V_{SE} &= \lambda(\hat S\otimes\hat E + \hat S^\dagger\otimes\hat E^\dagger),
\end{align*}

where $\hat S$ and $\hat E$ are non-hermitian. Operator like this can always be transformed into a combination of hermitian products:

\begin{align*}
&\hat S\otimes\hat E+\hat S^\dagger\otimes\hat E^\dagger
= \frac{\hat S+\hat S^\dagger}{\sqrt 2}\otimes\frac{\hat E+\hat E^\dagger}{\sqrt 2} + \frac{\hat S-\hat S^\dagger}{i\sqrt 2}\otimes\frac{\hat E^\dagger-\hat E}{i\sqrt 2}.
\end{align*}

[[sol:general_super_qum_explicit]Solution] Calculate the explicit forms of super-qumulants $\mathcal{L}^{(1)}(s)$ and $\mathcal{L}^{(2)}(s)$ for the general $SE$ coupling,

\begin{align*}
\hat V_{SE} &= \lambda\sum_\alpha \hat V_\alpha\otimes\hat F_\alpha,
\end{align*}

where $\hat V_\alpha$ and $\hat F_\alpha$ are hermitian operators.

§.§ An interpretation of joint quasi-probability distributions

So far we have been single-mindedly exploiting the structural parallelism between dynamical and stochastic maps to arrive as efficiently as possible at the super qumulant expansion. On our way there, we have identified the joint quasiprobability distributions $\{Q_F^{(k)}\}_{k=1}^\infty$; although, surely, they are a new and interesting object, up to this point we have looked at them only from purely utilitarian perspective. Now is the moment to pause for a minute and consider possible physical interpretations.

The choice of the prefix “quasi” in the name have been very much deliberate: joint quasiprobability distributions are similar to joint probability distributions in certain ways (both function as distributions to integrate with, and both satisfy analogous consistency relations) but—and this is vital—there are also key differences that completed change their respective natures. Let us illustrate what we have in mind here with an example. First and foremost, joint probability distributions $P_F^{(k)}(f_1,s_1;\ldots;f_k,s_k)$ are probability distributions for chronological sequences of form $(f_k,\ldots,f_1)$. The most natural physical interpretation is that such a sequence is a record of trajectory of process $F(t)$ captured by repeated observations at the consecutive moments in time $0 < s_k < s_{k-1} < \cdots < s_1$; hence, the sequence can be though of as a sample trajectory of $F(t)$. This observation has important practical implications. On the one hand, if one has access to functions $P_F^{(k)}$, then one can generate trajectories of $F(t)$ by drawing sample sequences from the probability distributions (see assignment <ref> of Sec. <ref>). On the other hand, when one has access to sufficiently large ensemble of sample trajectories, then one can approximate the stochastic average over distribution $P_F^{(k)}$ with the sample average (see Sec. <ref>); the equivalence of sample average and average over distribution follows from the very definition of probability. Since functions $Q_F^{(k)}$ are complex-valued, they are quasi-probability distributions—they appear in formulas as if they were probability distributions, but cannot be actual probabilities because they are not nonnegative. Consequently, it is impossible to use $Q_F^{(k)}$'s to generate something analogous to trajectories, nor to approximate the quasi-average with some sort of sample average. This is the key difference between the two types of distributions and the reason why, despite all the formal similarities, they are of a completely different natures.

What is the nature of joint quasiprobability distributions, then? Obviously, it is a tricky question to answer, but we may try to find an inspiration by looking for other contexts where quasi-probabilities appear naturally. One example of such context is the classic double-slit-type of experiments that have been often used to demonstrate, among other thing, the super-position principle and the effects of quantum interference. The most famous type of double-slit experiment is of the Young's experiment; it typically involves an emitter of particle beam shot towards a screen that has the ability to detect the position of particles hitting it. However, the direct path to the detector is blocked by an impassible barrier with a couple of slits punctured through. As a result, there are only a number of alliterative paths—each passing through different slit—a particle can take to reach the screen. For quantum particles, those different paths are in super-position and they can interfere with each other. To avoid technical complications, inevitable when attempting to describe such a setup, we will consider instead a much simpler variant of double-slit experiment—the Stern-Gerlach (SG) experiment.

The SG experiment involves particles with spin degree of freedom, described with operator

\begin{align}\label{eq:spin_op}
\hat S = \sum_{m=-s}^s m |m\rangle\langle m|
\end{align}

representing, say, the $z$-component of the spin (e.g., for $s=1/2$ we would have $\hat S=\hat\sigma_z/2$). The beam of those particles is emitted towards the SG apparatus that employs the magnetic field gradient to deflect moving particles in the directions dependent on the value of $m$—the higher the value the steeper the deflection angle. Consequently, when the initial beam passes through the apparatus, it leaves split into $2s+1$ distinct beams, each characterized by a definite value of $m$. The detector screen is placed right behind the output ports of the SG apparatus; since the beam has been split, one is now able to correlate the position of the detected hit with the value of $m$, thus realizing the measurement of the spin. According to the rules of quantum mechanics, the probability of measuring a particular value $m$ is given by the modulus square of the probability amplitude for the corresponding process leading to the detection event,

\begin{align}
\operatorname{Prob}(m) &= |\operatorname{Amp}(m)|^2.
\end{align}

Assuming that the travel time from the emitter to the screen is $t$, and the initial state of the spin is given by $|\Psi\rangle =\sum_{m'} c_{m'}|m'\rangle$, and the Hamiltonian describing the evolution of the beam is $\hat H_\mathrm{beam}$ (it includes the beam-splitting SG apparatus), we can write the amplitude as

\begin{align}
\operatorname{Amp}(m) &= \langle m|e^{-it\hat H_\mathrm{beam}}|\Psi\rangle = \sum_{m'=-s}^sc_{m'}\langle m|e^{-it\hat H_\mathrm{beam}}|m'\rangle \equiv \sum_{m'=-s}^s c_{m'}\phi_t(m|m').
\end{align}

We are not indicating here the spatial degrees of freedom because they do not have any impact on our discussion.

To turn this basic setup into proper double-slit experiment, we insert an additional SG apparatus somewhere in between the emitter and the detection screen (say, it takes time $u<t$ for particles to reach it). Right after the apparatus we place an impassible barrier with two slits punctured through in such a way that only $m=\pm s$ beams are not being blocked. Then, we install one more SG apparatus with the gradient orientation in the opposite direction so that the split beams recombine into a single beam. In this setup, the amplitude for the spin to be found in state $m$, provided that it has passed trough one of the slits, is calculated according to the product rule,

\begin{align}
\operatorname{Amp}(m|{\pm s}) &= \langle m|e^{-i(t-u)\hat H_\mathrm{beam}}|{\pm s}\rangle\langle {\pm s}|e^{-iu\hat H_\mathrm{beam}}|\Psi\rangle
= \sum_{m'=-s}^sc_{m'}\phi_{t-u}(m|{\pm s})\phi_u({\pm s}|m').
\end{align}

The total amplitude is a sum of amplitudes for passing through each slit,

\begin{align}
\operatorname{Amp}(m) &= \operatorname{Amp}(m|{+s})+\operatorname{Amp}(m|{-s})
=\sum_{m'=-s}^sc_{m'}\sum_{m''=\pm s}\phi_{t-u}(m|m'')\phi_u(m''|m'),
\end{align}

and thus, the probability of measuring $m$ reads

\begin{align}
\nonumber
\operatorname{Prob}(m) &=|\operatorname{Amp}(m)|^2 = \big|\sum_{m'=\pm s}\operatorname{Amp}(m|m')\big|^2\\
\nonumber
&=\sum_{m'=\pm s}|\operatorname{Amp}(m|m')|^2
+ \operatorname{Amp}(m|{+s})\operatorname{Amp}(m|{-s})^*
+ \operatorname{Amp}(m|{-s})\operatorname{Amp}(m|{+s})^*\\
&\equiv \sum_{m'=\pm s}\operatorname{Prob}(m|m') + \Phi(m|s,{-s}).
\end{align}

We got here the standard result for double-slit experiment: $\operatorname{Prob}(m|{+s})+\operatorname{Prob}(m|{-s})$ is the classical probability for reaching the screen through two alternative paths, and the interference term $\Phi(m|s,{-s})$ is responsible for all the quantum effects visible in the measured results.

Now we establish the connection with quasi-probabilities. Given the spectral decomposition of the spin operator $\hat S$, the Hamiltonian of the spin particle beam $\hat H_\mathrm{beam}$, and the initial pure state, $\hat\rho = |\Psi\rangle\langle\Psi|$, we can define a joint quasiprobability distribution associated with $\hat S$,

\begin{align}
\nonumber
\nonumber
&\phantom{==} = \operatorname{tr}\left(
\hat P_I(m_1,s_1)\cdots\hat P_I(m_k,s_k)\hat\rho\hat P_I(\bar{m}_k,s_k)\cdots\hat P_I(\bar{m}_1,s_1)
\right)\\
\nonumber
&\phantom{==} = \operatorname{tr}\left(
e^{is_1\hat H_\mathrm{beam}}|m_1\rangle\langle m_1|e^{-is_1\hat H_\mathrm{beam}}\cdots e^{is_k\hat H_\mathrm{beam}}|m_k\rangle\langle m_k|e^{-is_k\hat H_\mathrm{beam}}|\Psi\rangle
\right.\\
\nonumber
\langle\Psi|e^{is_k\hat H_\mathrm{beam}}|\bar{m}_k\rangle\langle \bar{m}_k|e^{-is_k\hat H_\mathrm{beam}}\cdots e^{is_1\hat H_\mathrm{beam}}|\bar{m}_1\rangle\langle \bar{m}_1|e^{-is_k\hat H_\mathrm{beam}}
\right)\\
\nonumber
&\phantom{==} = \delta_{m_1,\bar{m}_1}\langle m_1|e^{-i(s_1-s_2)\hat H_\mathrm{beam}}|m_2\rangle\cdots\langle m_k|e^{-is_k\hat H_\mathrm{beam}}|\Psi\rangle\\
\nonumber
\langle m_1|e^{-i(s_1-s_2)\hat H_\mathrm{beam}}|\bar{m}_2\rangle\cdots\langle \bar{m}_k|e^{-is_k\hat H_\mathrm{beam}}|\Psi\rangle
\right)^*\\
&\phantom{==} = \delta_{m_1,\bar{m}_1}\sum_{m_{k+1}}\sum_{\bar{m}_{k+1}}c_{m_{k+1}}c_{\bar{m}_{k+1}}^*
\prod_{i=1}^k\phi_{s_i-s_{i+1}}(m_i|m_{i+1})\phi_{s_i-s_{i+1}}(\bar{m}_i|\bar{m}_{i+1})^*,
\end{align}

and we have set here $s_{k+1} = 0$. Observe that $Q_S^{(k)}$ has a form of a product of the elementary amplitudes $\phi_{s_i-s_{i+1}}(m_i|m_{i+1})$. To clarify this point, let us focus on the second-order quasiprobability,

\begin{align}
\nonumber
Q_S^{(2)}(m_1,m_1,t;m_2,\bar{m}_2,u) &= \sum_{m_0}c_{m_0}\phi_{t-u}(m_1|m_2)\phi_u(m_2|m_0)
\sum_{\bar{m}_0}c_{\bar{m}_0}^*\phi_{t-u}(m_1|\bar{m}_2)^*\phi_u(\bar{m}_2|\bar{m}_0)^*\\ &=\operatorname{Amp}(m_1|m_2)\operatorname{Amp}(m_1|\bar{m}_2)^*.
\end{align}

We see here how exactly a quasiprobability relates to probability amplitudes $\operatorname{Amp}$ we have been discussing so far. Consequently, both parts of the measurement probability $\operatorname{Prob}(m)$ in our double-slit experiment can be rewritten in terms of quasiprobability distribution,

\begin{align}
\operatorname{Prob}(m|{\pm s}) &= Q_S^{(2)}(m,m,t;\pm s,\pm s, u);\\
\Phi(m|{s},{-s}) &= Q_S^{(2)}(m,m,t;{+}s,{-}s,u)+Q_S^{(2)}(m,m,t;{-}s,{+}s,u).
\end{align}

Of course, this result can be generalized: Each additional SG apparatus and barrier-with-slits combo increases the order $k$ of used quasi-probabilities $Q^{(k)}_S$; each additional slit punctured through barrier adds more possible values of argument pair $m_i,\bar{m}_i$. In any case, the overall pattern always remains unchanged. First, the quasiprobability with matching sets of arguments, $\bar{m}_i = m_i$ for all $i$, equals the classical conditional probability for particles to get to the detector through alternative paths defined by the sequence of traversed slits $(m_1,m_2,\ldots,m_k)$,

\begin{align}
Q_S^{(k)}(m_1,m_1,t;m_2,m_2,u_2;\ldots;m_k,m_k,u_k) &= \operatorname{Prob}(m_1|m_2,\ldots,m_k).
\end{align}

Second, all the quasi-probabilities with one path $(m_1,m_2,\ldots,m_k)$ not matching the other path $(m_1,\bar{m}_2,\ldots,\bar{m}_k)$ (even if the only difference is one pair of arguments), contribute collectively to the interference term. One could even say that each individual mismatched quasiprobability, $Q_S^{(k)}(m_1,m_1,t;m_2,\bar{m}_2,u_2;\ldots)$, quantifies the quantum interference between the two alternative paths.

Finally, note how in our example the joint quasiprobability distributions appeared even though the SG experiment is a closed quantum system. This underlines an important feature of $\{Q_F^{(k)}\}_{k=1}^\infty$: joint quasi-probabilities are defined exclusively by the dynamical properties of the system the observable $\hat F$ belongs to. Therefore, in the context of dynamical map, the quasi-probabilities that holistically describe the influence $E$ has on the dynamics of open system $S$, are completely independent of $S$ itself.

§.§ Central super-qumulants

While studying qumulants and super-qumulants (or cumulants and super-cumulants), it is not difficult to notice that the vast majority of formulas simplify tremendously when the average value $\langle F(t)\rangle = \operatorname{tr}_E(\hat F_I(t)\hat\rho_E)$ (or $\overline{F(t)}$) is set to zero. Here, we will show how to eliminate the contribution from the average even when $\langle F(t)\rangle\neq 0$ by switching the dynamical map to central picture.

We start by rewriting the total Hamiltonian of the $SE$ system,

\begin{align}
\nonumber
\hat H_{SE} &= \hat H_S\otimes\hat 1 + \hat 1\otimes\hat H_E + \lambda \hat V\otimes\hat F +\left(\lambda\langle F(t)\rangle\hat V\otimes\hat{1}-\lambda\langle F(t)\rangle\hat V\otimes\hat{1}\right)\\
\nonumber
&= \left(\hat H_S+\lambda\langle F(t)\rangle\hat V\right)\otimes\hat 1 + \hat 1\otimes\hat H_E + \lambda\hat V\otimes(\hat F - \langle F(t)\rangle\hat 1)\\
&= \hat H_S'(t)\otimes\hat 1 + \hat 1\otimes\hat H_E + \lambda\hat V\otimes\hat X(t),
\end{align}

\begin{align}
\hat H_S'(t) &= \hat H_S + \lambda\langle F(t)\rangle\hat V,\\
\hat X(t) &= \hat F - \langle F(t)\rangle \hat 1 = \sum_{f}(f-\langle F(t)\rangle)\hat P(f).
\end{align}

The operator $\hat H_S'$ can be though of as a modified Hamiltonian of the open system $S$, and the “centered” coupling $\hat X$ has zero average,

\begin{align}
\langle X(t) \rangle &=\operatorname{tr}_E\left(e^{it\hat H_E}\hat X(t)e^{-it\hat H_E}\hat\rho_E\right) = \operatorname{tr}_E\left(\hat F_I(t)\hat\rho_E\right)-\langle F(t)\rangle\operatorname{tr}_E(\hat\rho_E)
= 0.
\end{align}

The next step is to disentangle from the dynamical map the $S$-only part of the total Hamiltonian, $\mathcal{H}_S'(t)=[\hat H_S'(t),\bullet]$,

\begin{align}
\langle\mathcal{U}\rangle(t,0) &= \mathcal{U}_S'(t,0)\langle\mathcal{U}\rangle_c(t,0),
\end{align}

where $\mathcal{U}_S'(t,0) = T_{\mathcal{H}_S'}\exp[-i\int_0^t\mathcal{H}_S'(s)ds]$ and $\langle\mathcal{U}\rangle_c$ is the map in central picture, that we got instead of the interaction picture which we would have obtained if we had disentangle only $[\hat H_S,\bullet]$. Since the eigenkets of $\hat X$ and $\hat F$ are identical, we can express the central picture map as a quasi-moment series with only slight modifications (compare with the assignment <ref> of Sec. <ref>),

\begin{align}
\nonumber
\langle\mathcal{U}\rangle_c(t,0) &= \sum_{k=0}^\infty (-i\lambda)^k\int_0^t\!ds_1\cdots\int_0^{s_{k-1}}\!\!\!ds_k
\sum_{f_1,\bar{f}_1}\cdots\sum_{f_k,\bar{f}_k}Q_F^{(k)}(f_1,\bar{f}_1,s_1;\ldots;f_k,\bar{f}_k,s_k)\\
\nonumber
\mathcal{W}_c(f_1-\langle F(s_1)\rangle,\bar{f}_1-\langle F(s_1)\rangle,s_1)\cdots\mathcal{W}_c(f_k-\langle F(s_k)\rangle,\bar{f}_k-\langle F(s_k)\rangle,s_k)\\
\nonumber
\int\limits_0^t\!\!ds_1\mathcal{W}_c(F(s_1)-\langle F(s_1)\rangle,\bar{F}(s_1)-\langle F(s_1)\rangle,s_1)\cdots\right.\\
\nonumber
\cdots\!\!\!\int\limits_0^{s_{k-1}}\!\!\!\!ds_k\mathcal{W}_c(F(s_k)-\langle F(s_k)\rangle,\bar{F}(s_k)-\langle F(s_k)\rangle,s_k)
\right\rangle\\
\nonumber
\int\limits_0^t\!\!ds_1\mathcal{W}_c(X(s_1),\bar{X}(s_1),s_1)\cdots\!\!\!\!\int\limits_0^{s_{k-1}}\!\!\!\!ds_k\mathcal{W}_c(X(s_k),\bar{X}(s_k),s_k)
\right\rangle\\
&=\left\langle T_\mathcal{W}e^{-i\lambda\int_0^t\mathcal{W}_c(X(s),\bar{X}(s),s)ds}\right\rangle,
\end{align}

where the central picture of super-operators $\mathcal{W}$ are given by

\begin{align}
&\hat V_c(s) = \left(T_{H_S'}e^{-i\int_0^s\hat H_S'(u)du}\right)^\dagger \hat V\,T_{H_S'}e^{-i\int_0^s\hat H_S'(u)du}= \hat U_S'(0,s)\hat V\hat U_S'(s,0),\\[.2cm]
&\mathcal{W}_c(x,\bar{x},s) = x\hat V_c(s)\bullet-\bullet\hat V_c(s)\bar{x}.
\end{align}

Note that the central and interaction pictures are related via the following transformation,

\begin{align}
\hat V_c(s) &= \hat U_S'(0,s) e^{-is\hat H_S}\hat V_I(s)\,e^{is\hat H_S}\hat U_S'(s,0),\\
\mathcal{W}_c(x,\bar{x},s) &=\mathcal{U}_S'(0,s)e^{-is[\hat H_S,\bullet]}\mathcal{W}(x,\bar{x},s)e^{is[\hat H_S,\bullet]}\mathcal{U}_S'(s,0).
\end{align}

The final step is to introduce the central super-qumulants $\mathcal{L}_c^{(k)}(s)$ in the usual fashion,

\begin{align}
\left\langle T_{\mathcal{W}_c}e^{-i\lambda\int_0^t\mathcal{W}_c(X(s),\bar{X}(s),s)ds}\right\rangle &\equiv T_{\mathcal{L}_c^{(k)}}e^{\sum_{k=1}^\infty (-i\lambda)^k\int_0^t\mathcal{L}_c^{(k)}(s)ds}.
\end{align}

We, then, find the consecutive $\mathcal{L}_c^{(k)}$ by comparing the terms of the corresponding orders in $\lambda$ on LHS and RHS.

As advertised, the first-order central super-qumulant vanishes,

\begin{align}
\nonumber
\mathcal{L}_c^{(1)}(s) &= \langle \mathcal{W}_c(X(s),\bar{X}(s),s)\rangle = \langle F(s)-\langle F(s)\rangle\rangle \hat V_c(s)\bullet-\bullet\hat V_c(s)\langle \bar{F}(s)-\langle F(s)\rangle\rangle\\
&=\left(\langle F(s)\rangle-\langle F(s)\rangle\right)\hat V_c(s)\bullet-\bullet\hat V_c(s)\left(\langle \bar{F}(s)\rangle -\langle F(s)\rangle\right) = 0.
\end{align}

This, of course, has ripple effects on all subsequent orders. The second central super-qumulant reduces to

\begin{align}
\nonumber
\mathcal{L}_c^{(2)}(s) &= \int_0^s\langle\mathcal{W}_c(X(s),\bar{X}(s),s)\mathcal{W}_c(X(u),\bar{X}(u),u)\rangle du - \int_0^s\mathcal{L}_c^{(1)}(s)\mathcal{L}^{(1)}_c(u)du\\
\nonumber
&=\int_0^s\Big\langle\!\big(F(s)-\langle F(s)\rangle\big)\big(F(u)-\langle F(u)\rangle\big)\!\Big\rangle\left(\hat V_c(s)\hat V_c(u){\bullet}-\hat V_c(u){\bullet}\hat V_c(s)\right)du\\
\nonumber
&\phantom{=}{+}\int_0^s\Big\langle\!\big(\bar{F}(s)-\langle F(s)\rangle\big)\big(\bar{F}(u)-\langle F(u)\rangle\big)\!\Big\rangle\left({\bullet}\hat V_c(u)\hat V_c(s)-\hat V_c(s){\bullet}\hat V_c(u)\right)du\\
\nonumber
&=\int_0^s\langle\!\langle F(s)F(u)\rangle\!\rangle\left(\hat V_c(s)\hat V_c(u){\bullet}-\hat V_c(u){\bullet}\hat V_c(s)\right)du\\
\nonumber
&\phantom{=}{+}\int_0^s\langle\!\langle \bar{F}(s)\bar{F}(u)\rangle\!\rangle\left({\bullet}\hat V_c(u)\hat V_c(s)-\hat V_c(s){\bullet}\hat V_c(u)\right)du\\
&=\int_0^s \langle\!\langle \mathcal{W}_c(F(s),\bar{F}(s),s)\mathcal{W}_c(F(u),\bar{F}(u),u)\rangle\!\rangle du,
\end{align}

or, using the shorthand notation,

\begin{align}
&\tilde{Q}_{i_1\cdots i_k}\mathcal{W}_{i_1}\cdots\mathcal{W}_{i_k}
= \sum_{f_{i_1},\bar{f}_{i_1}}\!\!\cdots\!\!\sum_{f_{i_k},\bar{f}_{i_k}}
\tilde{Q}_F^{(k)}(f_{i_1},\bar{f}_{i_1},s_{i_1};\ldots;f_{i_k},\bar{f}_{i_k},s_k)
\mathcal{W}_c(f_{i_1},\bar{f}_{i_1},s_{i_1})\cdots\mathcal{W}_c(f_{i_k},\bar{f}_{i_k},s_{i_k}),
\end{align}

we can write it as

\begin{align}
\mathcal{L}_c^{(2)}(s_1) &= \int_0^{s_1}\!\!\!\!ds_2\, \tilde{Q}_{12}\mathcal{W}_1\mathcal{W}_2.
\end{align}

Then, normally very complex higher order super-qumulants, simplify significantly in the central picture,

\begin{align}
\mathcal{L}_c^{(3)}(s_1) &= \int_0^{s_1}\!\!\!\!ds_2\int_0^{s_2}\!\!\!\!ds_3\,\tilde{Q}_{123}\mathcal{W}_1\mathcal{W}_2\mathcal{W}_3,\\
\nonumber
\mathcal{L}_c^{(4)}(s_1) &= \int_0^{s_1}\!\!\!\!ds_2\cdots\int_0^{s_4}\!\!\!\!ds_4\,\tilde{Q}_{1234} \mathcal{W}_1\mathcal{W}_2\mathcal{W}_3\mathcal{W}_4\\
\label{eq:Lc4}
\tilde{Q}_{13}\tilde{Q}_{24}\mathcal{W}_1[\mathcal{W}_2,\mathcal{W}_3]\mathcal{W}_4
\end{align}

The cost of switching to the central picture is rather minuscule: computing the transformation operator $\hat U_S'$ is a trivial problem (it satisfies a standard dynamical equation). The readily apparent benefits we have just showcased above are definitely worth this additional bit of work.

§.§ Computability of dynamical maps

Due to complexities of the environment $E$, the computation of dynamical maps is a notoriously difficult problem (even more so than it is for stochastic maps). Since, in general, dynamical maps do not satisfy any form of dynamical equation—the structural parallelism with stochastic maps can be invoked here to explain why it is the case—a systematic way forward is to approximate the central super-qumulant (cSQ) series with a closed-form super-operator, and by doing so, force a dynamical equation upon the map. Once we have the equation, we can integrate it numerically, and thus, compute a map approximating the real dynamical map. This is not the only way to achieve computability, mind you; there are myriad specialized methods to compute an approximated map, but, typically, they work only in a narrowly defined contexts. Here, we want to discuss a general purpose approach, hence the focus on the cSQ series.

A good starting point for most approximation schemes is to assume that $E$ is Gaussian, i.e., to neglect the contribution from qumulant distributions $\tilde{Q}_F^{(k)}$ with $k>2$—an approximation analogous to cutting the cumulant series of stochastic process at the second order. The Gaussian environment means that the totality of the information about $E$ is contained in the trivial expectation value of the coupling $\langle F(t)\rangle$ and the second-order qumulants,

\begin{align}
\langle\!\langle F(t)F(s) \rangle\!\rangle &= \langle\!\langle \bar{F}(t)F(s) \rangle\!\rangle = \langle F(t)F(s)\rangle - \langle F(t)\rangle\langle F(s)\rangle = C_F(t,s)-i K_F(t,s),\\
\langle\!\langle \bar{F}(t)\bar{F}(s)\rangle\!\rangle &= \langle\!\langle F(t)\bar{F}(s)\rangle\!\rangle = \langle \bar{F}(t)\bar{F}(s)\rangle - \langle F(t)\rangle\langle F(s)\rangle = C_F(t,s)+i K_F(t,s);
\end{align}

that can be expressed through two physically significant quantities,

\begin{align}
\nonumber
C_F(s,u) &= \frac{1}{2}\left(\langle F(t)F(s)\rangle+\langle \bar{F}(t)\bar{F}(s)\rangle \right) - \langle F(t)\rangle\langle F(s)\rangle\\
&=\frac{1}{2}\operatorname{tr}_E\left(\Big\{ \hat F_I(t)-\langle F(t)\rangle\ ,\ \hat F_I(s)-\langle F(s)\rangle\Big\}\hat\rho_E\right);\\[.2cm]
\nonumber
K_F(s,u) &= \frac{1}{2}i\left(\langle F(t)F(s)\rangle - \langle \bar{F}(t)\bar{F}(s)\rangle\right)\\
&=\frac{1}{2}i\operatorname{tr}_E\left(\Big[\hat F_I(t)-\langle F(t)\rangle\ ,\ \hat F_I(s)-\langle F(s)\rangle\Big]\hat\rho_E\right).
\end{align}

These two objects, the susceptibility $K_F$ and the correlation function $C_F$ (and the spectral density $S(\omega)$, the Fourier transform of $C_F$), can often be measured—with the linear response technique to obtain the former, and with noise spectroscopy to get the latter. Moreover, when $E$ is in thermal equilibrium (i.e., $\hat\rho_E \propto e^{-\beta \hat H_E}$), then the spectral density and the susceptibility are related through the Fluctuation-Dissipation theorem, so that one can be obtained from the other, given the temperature of the environment.

[[sol:FDT]Solution] Prove the Fluctuation-Dissipation theorem,

\begin{align*}
S(\omega) = i\left( \frac{1+e^{-\beta\omega}}{1-e^{-\beta\omega}}\right)\kappa(\omega),
\end{align*}

where $\beta = 1/k_\mathrm{B}T$, and

\begin{align*}
S(\omega) &= \int_{-\infty}^\infty C_F(\tau,0)e^{-i\omega \tau}d\tau;\
\kappa(\omega) =\int_{-\infty}^\infty K_F(\tau,0)e^{-i\omega \tau}d\tau.
\end{align*}

The Gaussianity is a natural approximation—or even an exact property, in some cases—for environments composed of indistinguishable bosons; some prominent examples of such systems include phononic excitations in solids and electromagnetic quantum fields. Gaussianity can also be expected in $E$ where the central limit theorem-style of argument is applicable, i.e., systems that are composed of large number of roughly independent and similar subsystems,

\begin{align}
\hat  F \approx \frac{1}{\sqrt N}\sum_{i=1}^{N}\hat 1^{\otimes (i-1)}\otimes\hat f\otimes\hat 1^{\otimes(N-i)};
\ \hat H_E \approx \sum_{i=1}^N\hat 1^{\otimes(i-1)}\otimes\hat h\otimes\hat 1^{\otimes(N-i)};
\ \hat\rho_E\approx \hat \rho^{\otimes N}.
\end{align}

An example of such $E$ would be a bath of nuclear spins in crystal lattices. Finally, note that the Gaussianity, or any other approximation to quasi-moments and qumulants, can be argued for independently of any properties of the open system $S$ because $Q_F^{(k)}$'s are defined exclusively by the coupling operator $\hat F$ and the dynamical properties of $E$.

Assuming that $E$ is Gaussian, the simplest approximation to cSQ series on the $S$ side is to truncate it at some finite order $n$ (note that for Gaussian $E$ odd-order cSQs are zero), so that

\begin{align}
\langle\mathcal{U}\rangle_c(t,0) &\approx T_{\mathcal{L}_c^{(k)}}e^{\sum_{k=1}^n(-i\lambda)^{2k} \int_0^t\mathcal{L}_c^{(2k)}(s)ds} \equiv \langle\mathcal{U}\rangle^{(2n)}(t,0),
\end{align}

where, of course, the approximating map $\langle\mathcal{U}\rangle^{(2n)}$ is computable because it satisfies the dynamical equation,

\begin{align}
\frac{d}{dt}\langle\mathcal{U}\rangle^{(2n)}(t,s) &=\left(\sum_{k=1}^n(-i\lambda)^{2k} \mathcal{L}_c^{(2k)}(t)\right)\langle\mathcal{U}\rangle^{(2n)}(t,s);\quad\langle\mathcal{U}\rangle^{(2n)}(s,s)=\bullet.
\end{align}

The justification for this strategy is based on the estimate of the order of magnitude for cSQs,

\begin{align}\label{eq:L_magnitude}
\lambda^{2k}\int_0^t\mathcal{L}_c^{(2k)}(s)ds &\sim \lambda^{2k}t\tau_c^{2k-1}\left(\frac{\tau_c}{\tau_S}\right)^{k-1};
\end{align}

let us now explain where this estimation comes from. Here, it is assumed the coupling operators has been defined in such a way that $||\mathcal{W}_c||\sim 1$ so that the overall magnitude of the coupling is quantified by $\lambda$. Next, the correlation time $\tau_c$ is defined as a range of the correlation function,

\begin{align}
C_F(t,s) \xrightarrow{|t-s|\gg\tau_c} 0.
\end{align}

More generally, $\tau_c$ is a time-scale on which the quasi-process $(F(t),\bar{F}(t))$ decorrelates,

\begin{align}
\left.\begin{array}{c}
\langle\!\langle \cdots F(t)\cdots F(s)\cdots\rangle\!\rangle \\[.2cm]
\langle\!\langle \cdots \bar{F}(t)\cdots F(s)\cdots\rangle\!\rangle \\[.2cm]
\langle\!\langle \cdots F(t)\cdots \bar{F}(s)\cdots\rangle\!\rangle \\[.2cm]
\langle\!\langle \cdots \bar{F}(t)\cdots \bar{F}(s)\cdots\rangle\!\rangle \\
\end{array}\right\}
&\xrightarrow{|t-s| \gg \tau_c} 0.
\end{align}

This follows from the analogy with the statistical irreducibility of cumulants (see Sec. <ref>), that can also be applied to the the process evaluated at different points in time. Typically, the correlation between the values of stochastic process at different time points, say $F(t)$ and $F(s)$, weaken as the distance $|t-s|$ increases; when the correlations drop to zero, it means that $F(t)$ and $F(s)$ become statistically independent, and thus, the cumulant $\overline{\overline{\cdots F(t)\cdots F(s)\cdots}}$ must vanish because it is statistically irreducible. Since the irreducibility is a consequence of the cumulant structure (see, e.g., Sec. <ref>), due to parallelism, it also applies to quasi-cumulants. This explains the factor $\lambda^{2k}t\tau_c^{2k-1}$ which comes for the $2k$-fold time integrations that are being restricted by the range of correlations withing the central qumulants constituting the cSQ. Note that this estimation would not work for non-central SQs because the contributions form the average $\langle F(s)\rangle$ interwoven into (non-central) qumulants would not restrict the time-integrals.

The second factor in Eq (<ref>), $(\tau_c/\tau_S)^{k-1}$, estimates the contribution from commutators between $\mathcal{W}_c$'s at different times that appear in $\mathcal{L}_c^{(2k)}$'s [see, e.g., Eq. (<ref>)]; these commutators have the following general form,

\begin{align}
\sum_{f_2,\bar{f}_2}\sum_{f_3,\bar{f}_3}\tilde{Q}_F^{(2)}(f_1,\bar{f}_1,s_1;f_3,\bar{f}_3,s_3)\tilde{Q}_F^{(2)}(f_2,\bar{f}_2,s_2;f_4,\bar{f}_4,s_4)
[\mathcal{W}_c(f_2,\bar{f}_2,s_2),\mathcal{W}_c(f_3,\bar{f}_3,s_3) ].
\end{align}

Of course, the distance between time arguments of $\mathcal{W}_c$ in question ($|s_2-s_3|$ in the above example) is capped by $\tau_c$ due to the range of correlations in $E$ that we have discussed previously. Naturally, if the central picture does not change appreciably on the time-scale $\tau_c$, i.e.,

\begin{align}
\mathcal{W}_c(f,\bar{f},s+\tau_c) &= \mathcal{U}_S'(s,s+\tau_c)\mathcal{W}_c(f,\bar{f},s)\mathcal{U}_S'(s+\tau_c,s)\approx \mathcal{W}_c(f,\bar{f},s),
\end{align}

then the commutator will tend to zero because $\mathcal{W}_c$ at $s_2$ and $s_3$ will, basically, be the same super-operator. Here, we want to estimate the magnitude of the commutator when the change is small but still finite. To this end, we assume $E$ is stationary, i.e., $[\hat\rho_E,\hat H_E] = 0$, so that $\langle F(t)\rangle$ is time-independent, and thus, so is $\hat H_S'$; then, the generator of the central picture has spectral decomposition

\begin{align}
\mathcal{H}_S' &= \sum_{a,b} (\epsilon_a-\epsilon_b)|a\rangle\langle b|\operatorname{tr}_S\left(|b\rangle\langle a|\bullet\right).
\end{align}

defined by the modified energy levels of $S$, $\hat H_S'|a\rangle =\epsilon_a|a\rangle$. Now, let $\Omega$ be the maximal spacing between the levels, i.e.,

\begin{align}\Omega = \max_{a,b}|\epsilon_a-\epsilon_b| \equiv \frac{2\pi}{\tau_S},
\end{align}

which makes $\tau_S$ the fastest time scale associated with the central picture. The assumption that the central picture changes very little can be quantified as $\Omega\tau_c\ll 1$ or $\tau_c\ll \tau_S$, and it implies that the evolution of super-operators $\mathcal{W}_c$ on time-scale $\tau_c$ can be approximated with the lowest-order correction,

\begin{align}
\nonumber
[\mathcal{W}_c(s+\tau_c),\mathcal{W}_c(s)] &\approx [\mathcal{W}_c(s)+\tau_c\frac{d\mathcal{W}_c(s)}{ds},\mathcal{W}_c(s)] = [ -i\tau_c[\mathcal{H}_S',\mathcal{W}_c(s)] , \mathcal{W}_c(s)]\\
\end{align}

where we used the fact that $\mathcal{U}_S'$ is a time-ordered exp generated by $\mathcal{H}_S'$. The matrix elements, $\mathcal{A}_{ab,a'b'} = \operatorname{tr}(|b\rangle\langle a|\mathcal{A}|a'\rangle\langle b'|)$, of the commutator then read

\begin{align}
|([\mathcal{W}_c(s+\tau_c),\mathcal{W}_c(s)])_{ab,a'b'}|\approx \tau_c|(\epsilon_a-\epsilon_b+\epsilon_{a'}-\epsilon_{b'})(\mathcal{W}_c(s)^2)_{ab,a'b'}|\sim \Omega\tau_c
\sim \frac{\tau_c}{\tau_S},
\end{align}

which gives us our estimate.

A more advanced approximations to the cSQ series can be obtained with diagramatic techniques. The idea is to represent the series, and cSQs themselves, as a sum of irreducible diagrams. These are highly sophisticated techniques because inventing the set of rules for translating mathematical formulas into diagrams is more of an art form than an engineering problem; how one deals with it depends on the particular context and one's skill, intuition and overall technical aptitude. Historical examples show that a well formulated diagram rules can potentially expose some normally hidden patterns that could be then exploited to, e.g., analytically sum up a whole sub-series of diagrams. Such finds are the Holy Grail of diagramatic methods, as they allow to incorporate corrections of all orders of magnitudes.

§ MASTER EQUATION

§.§ Completely positive maps

As we have defined it here, the objective of the open system theory is to find an efficient and accurate description of the dynamics of system $S$ induced by its coupling with the environment $E$. When the theory is viewed from this perspective, the central problem to solve is to derive from “first principles” the dynamical map that encapsulates all the relevant influence $E$ exerts onto $S$ and the effects the dynamics has on the state are simply the consequence of the map's action. Then, the main practical concern is to achieve the computability of the map while maintaining a high degree of accuracy and, simultaneously, minimizing the necessary input of information about $E$.

A different way to interpret the goals of open system theory is to focus on the effects dynamical maps have on the state of $S$ rather than the physical origins of the map. In such a framework, instead of deriving the dynamical map from first principles, one's objective is to determine what kind of map, treated as an abstract transformation, produces a desired effect on a state of system $S$. There are two main practical concerns regarding this approach. The first is to define the constraints on the form of an abstract map so that it could, at least in principle, correspond to some real dynamical map induced by interactions with actual physical environment. The second is to find a convenient parametrization of the valid abstract maps that would allow one to clearly define the expected effects it has on the states.

The industry standard answer to the first point is the so-called Kraus representation; consider the following relation,

\begin{align}
\nonumber
\langle \mathcal{U}\rangle(t,0) &= \operatorname{tr}_E\left(e^{-it[\hat H_{SE},\bullet]}\bullet\otimes\hat\rho_E\right)
= \operatorname{tr}_E\left(e^{-it\hat H_{SE}}\bullet\otimes\left(\sum_{i} \rho_j |j\rangle\langle j|\right) e^{it\hat H_{SE}}\right)\\
&=\sum_{j,j'} \sqrt{\rho_j}\langle j'|e^{-it\hat H_{SE}}|j\rangle \bullet \langle j|e^{it\hat H_{SE}}|j'\rangle \sqrt{\rho_j}
\equiv \sum_{j,j'} \hat K_{jj'} {\bullet}\hat K_{jj'}^\dagger,
\end{align}

where $B_E=\{ |j\rangle\}_j$ is a basis in $E$ ket-space made out of eigenkets of $\hat\rho_E$ and the Kraus operators are defined as

\begin{align}
\nonumber
\hat K_{jj'} &= \sqrt{\rho_j}\langle j'|e^{-it\hat H_{SE}}|j\rangle
= \sum_{n,n'}\sqrt{\rho_j}\left(\langle n'|\otimes\langle j'|e^{-it\hat H_{SE}}|n\rangle\otimes|j\rangle\right)|n'\rangle\langle n|\\
\label{eq:derived_Kraus}
&=\sum_{n,n'}|n'\rangle\sqrt{\rho_j}\left(e^{-it\hat H_{SE}}\right)_{n'j',nj}\langle n|,
\end{align}

($\{|n\rangle\}_n$ is a basis in $S$). The key property of the set $\{\hat K_{jj'}\}_{j,j'}$ is the relation,

\begin{align}
\nonumber
\sum_{j,j'}\hat K_{jj'}^\dagger\hat K_{jj'} &=\sum_{j,j'}\rho_j
\sum_{n,n'}|n'\rangle\left(e^{-it\hat H_{SE}}\right)_{nj',n'j}^*\langle n|
\sum_{n'',n'''}|n'''\rangle\left(e^{-it\hat H_{SE}}\right)_{n'''j',n''j}\langle n''|\\
\nonumber
&=\sum_{j,j'}\rho_j \sum_{n,n',n''}|n'\rangle\left(e^{it\hat H_{SE}}\right)_{n'j,nj'}\left(e^{-it\hat H_{SE}}\right)_{nj',n''j}\langle n''|\\
\nonumber
&= \sum_{j}\rho_j\sum_{n',n''}|n'\rangle\left(e^{it\hat H_{SE}}e^{-it\hat H_{SE}}\right)_{n'j,n''j}\langle n''|
=\sum_j \rho_j \sum_{n',n''}|n'\rangle\delta_{n',n''}\delta_{j,j}\langle n''|\\
&=\Big(\sum_j \rho_j\Big)\sum_{n'} |n'\rangle\langle n'| = \hat 1.
\end{align}

The representation of the dynamical map with the sum over Kraus operators is not unique, however, e.g., we could have picked a different basis $B_E' = \{ |\mu\rangle\}_\mu$ to represent the partial trace,

\begin{align*}
\operatorname{tr}_E(\hat A_{SE}) = \sum_\mu\langle \mu|\hat A_{SE}|\mu\rangle,
\end{align*}

or a different decomposition of the density matrix,

\begin{align*}
\hat\rho_E = \sum_\alpha p_\alpha |\Psi_\alpha\rangle\langle \Psi_\alpha|
\end{align*}

($p_\alpha\geqslant 0$ and $\sum_\alpha p_\alpha = 1$, but the kets $\{|\Psi_\alpha\rangle\}_\alpha$ are not necessarily mutually orthogonal, see Sec. <ref>), to obtain a different set of Kraus operators,

\begin{align}
\hat K_{\alpha\mu} = \sqrt{p_\alpha}\langle \mu|e^{-it\hat H_{SE}}|\Psi_\alpha\rangle,\quad \sum_{\alpha,\mu}\hat K_{\alpha\mu}^\dagger\hat K_{\alpha\mu}= \hat 1,
\end{align}

and a different representation of the same dynamical map. In summary, any dynamical map evaluated at a specific time $t$ has Kraus representation(s),

\begin{align}
\langle\mathcal{U}\rangle(t,0) = \sum_{j,j'}\hat K_{jj'}{\bullet}\hat K_{jj'}^\dagger = \sum_{\alpha,\mu}\hat K_{\alpha\mu}\bullet\hat K^\dagger_{\alpha\mu}=\ldots\,.
\end{align}

Conversely, according to the Naimark's dilation theorem, given a set of operators $\{\hat K_m\}_m$ in $S$ that satisfy $\sum_m\hat K_m^\dagger\hat K_m = \hat 1$, it is always possible to find an unitary operator $\hat U_{SS'}$ that acts in a composite space $SS'$, and a ket in $S'$ space $|\Psi_{S'}\rangle$, such that

\begin{align}
\sum_m \hat K_m{\bullet}\hat K_m^\dagger = \operatorname{tr}_{S'}\left(\hat U_{SS'}{\bullet}\otimes|\Psi_{S'}\rangle\langle \Psi_{S'}|\hat U_{SS'}^\dagger\right).
\end{align}

Of course, in such a setting, $S'$ can be formally treated as the environment to $S$ and $\hat U_{SS'}$ as an evolution operator generated by a corresponding $SS'$ Hamiltonian over some evolution time. Even though the system $S'$ and the evolution operator $\hat U_{SS'}$ are only postulated here, they could exist in principle, and that is enough. Indeed, for as long as one only observes system $S$ after the map $\sum_m\hat K_m{\bullet}\hat K_m^\dagger$ produces its effects, it is impossible to verify if $S'$ is real or not.

Kraus representation is also an answer to the question of parametrization for the effects-oriented view of dynamical maps. In particular, for qubit system $S$, any operator—including Kraus operators—can be decomposed into linear combination of Pauli matrices. It is easy to categorize what kind of effects can be produced by Kraus representations made out of simple combinations, e.g., $\{\hat K_1 =\sqrt{(1+p)/2}\,\hat 1, \hat K_2 =\sqrt{(1-p)/2}  \hat\sigma_z\}$ results in

\begin{align}
\sum_{m=1}^2\hat K_m\hat\rho_S\hat K_m^\dagger &= \left(\frac{1+p}{2}\right)\hat\rho_S+\left(\frac{1-p}2\right) \hat\sigma_z\hat\rho_S\hat\sigma_z
\langle{\uparrow}|\hat\rho_S|{\uparrow}\rangle & p\langle{\uparrow}|\hat\rho_S|{\downarrow}\rangle \\
p\langle{\downarrow}|\hat\rho_S|{\uparrow}\rangle & \langle{\downarrow}|\hat\rho_S|{\downarrow}\rangle \\
\end{array}\right].
\end{align}

Kraus representation made out of more complex combinations would, in a sense, combine those simple effects. Thus, there are only relatively few options for Kraus representations that can produce distinct and unique effects. Unfortunately, for non-qubit $S$, the categorization of possible Kraus representations by their effects is not so clear cut. Historically, this was not a big problem because, originally, the development of the open system theory as a distinct field was, almost wholly, motivated by the context of quantum computers and quantum information processing, and thus, non-qubit systems were not often considered.

Traditionally, valid abstract maps, that is, maps with Kraus representations, are referred to as completely positive (CP) maps; therefore, “CP map” can be used as a shorthand for “dynamical map in $S$ obtained as a partial trace of the total unitary map in $SE$, acting on an uncorrelated initial state with fixed $\hat\rho_E$ and variable $\hat\rho_S$.” However, as one might have guessed, there is some history behind this very specific nomenclature. The concept of complete positivity arose as quintessentially effects-oriented way to solve the problem of constraints for abstract maps. Here, we will review the basics of reasoning that leads naturally to the definition of CP.

Consider a physical system $S$ that has its dynamics defined by the map $\mathcal{M}_t$, i.e., if $\hat\rho_S$ is an arbitrary initial state of $S$, then $\mathcal{M}_t\hat\rho_S$ is the state at later time $t>0$, $\hat\rho_S(t)$. In what follows we fully embrace the effects-oriented reading of the theory, and thus, we are willfully ignoring the first-principles origin of the map as a reduced unitary evolution of the system coupled to its environment. In that case, we have to somehow deduce which properties we should demand of $\mathcal{M}_t$ so that it can function as a valid dynamical map. The most obvious first step is to make sure that the output operator $\mathcal{M}_t\hat\rho_S$ is a proper density matrix. Therefore, at the very least, the map $\mathcal{M}_t$ should preserve the trace and positivity of the input state,

* $\mathcal{M}_t$ is trace-preserving, i.e., $\forall_{t>0}\ \operatorname{tr}_S(\mathcal{M}_t\hat\rho_S) = \operatorname{tr}_S(\hat\rho_S)$;

* $\mathcal{M}_t$ is positive, i.e., if $\hat\rho_S> 0$, then $\forall_{t>0}\ \mathcal{M}_t\hat\rho_S > 0$.

However, positivity and trace-preservation (PTP) alone turns out to be insufficient. To see what is missing, we have to take into consideration other physical systems governed by their respective dynamical laws. In particular, it is feasible to have two systems—our initial $S$ and the other $S'$—such that, starting from $t=0$, their dynamical laws are independent of each other, i.e., if $\mathcal{M}_t$ is the PTP map of $S$ and $\mathcal{M}_t'$ is the PTP map of $S'$, then the total dynamics for $t>0$ is described by $\mathcal{M}_t\otimes\mathcal{M}_t'$. We can leave unspecified what happened in the past, t<0; for now, all we need to know is the initial (as in, $t=0$) state of the total $SS'$ system, $\hat\rho_{SS'}$. Even though we are assuming the dynamical laws are independent, the total map should, at least, be PTP; after all, the composition of physical systems is a physical system, and thus, $\mathcal{M}_t\otimes\mathcal{M}_t'\hat\rho_{SS'}$ should result in a proper density matrix of the composite system $SS'$. Nothing seems amiss if we assume that $S$ and $S'$ also start as uncorrelated systems, $\hat\rho_{SS'} = \hat\rho_S\otimes\hat\rho_{S'}$; indeed, the trace is preserved,

\begin{align*}
&\operatorname{tr}_{SS'}\left(\mathcal{M}_t\otimes\mathcal{M}_t'\hat\rho_S\otimes\hat\rho_{S'}\right) = \operatorname{tr}_S\left(\mathcal{M}_t\hat\rho_S\right)\operatorname{tr}_{S'}\left(\mathcal{M}_t'\hat\rho_{S'}\right) = \operatorname{tr}_S(\hat\rho_S)\operatorname{tr}_{S'}(\hat\rho_{S'}) = \operatorname{tr}_{SS'}\left(\hat\rho_S\otimes\hat\rho_{S'}\right);
\end{align*}

and so is positivity,

\begin{align*}
&\mathcal{M}_t\otimes\mathcal{M}_t'\hat\rho_S\otimes\hat\rho_{S'} = \left(\mathcal{M}_t\hat\rho_S\right)\otimes\left(\mathcal{M}_t'\hat\rho_{S'}\right) > 0,
\end{align*}

because the outer product $\otimes$ of positive operators is positive. We can even allow for classical correlations between $S$ and $S'$, $\hat\rho_{SS'} = \sum_i p_i \hat\rho_S^{(i)}\otimes\hat\rho_{S'}^{(i)}$ ($p_i \geqslant 0$, $\sum_i p_i = 1$), and still everything looks fine. The inadequacy of PTP is exposed only once we introduce non-classical correlations—the quantum entanglement—between $S$ and $S'$, i.e., we set $\hat\rho_{SS'}$ that is not of the form $\sum_i p_i \hat\rho_S^{(i)}\otimes\hat\rho_{S'}^{(i)}$. We can demonstrate the breakdown of PTP as the definite constraint with the famous Peres-Horodecki criterion for non-entangled states. The test is defined as follows: First, one constructs the following transformation in $SS'$,

\begin{align}
\mathcal{M}_{SS'} = \bullet\otimes(\bullet^T),
\end{align}

where $(\bullet^T)\hat A = \hat A^T$ is the matrix transposition—note here that the identity $\bullet$ and the transposition $\bullet^T$ are both PTP transformations. Second, one applies $\mathcal{M}_{SS'}$ to the subject of the test, the $SS'$ state $\hat\rho_{SS'}$. In the last step, one verifies if the resultant output is a positive operator. As we have discussed above, if the tested state $\hat\rho_{SS'}$ was only classically correlated, than the product of PTP transformations would produce a positive output. Therefore, a non-positive output $\bullet\otimes(\bullet^T)\hat\rho_{SS'}$ implies that $\hat\rho_{SS'} \neq \sum_i p_i\hat\rho_S^{(i)}\otimes\hat\rho_{S'}^{(i)}$ (although, the converse is true only for qubit $S$, $S'$ or for qubit $S$ and qutrit $S'$). So, the PTP requirement alone is so transparently inadequate to define the constraints on abstract maps that its failure is used to witness non-classical correlations! Clearly, there is more than PTP that one should require of physical dynamical maps $\mathcal{M}_t$; one way to "fix" the failure point with entangled composite systems is to introduce the third condition:

* $\mathcal{M}_t$ is such that $\mathcal{M}_{SS_d} = \mathcal{M}_t\otimes\bullet$, a transformation constructed to act in composite system $SS_d$ with $d$-dimensional $S_d$, is positive for any $d$, $t>0$.

The PTP transformation $\mathcal{M}_t$ that also satisfies the above is called a CP map—this is the standard mathematical definition of complete positivity. Of course, a total map $\mathcal{M}_t\otimes\mathcal{M}_t'$ made out of CP maps will produce a proper output density matrix even if the initial state is entangled. Also, a product of CP maps is CP itself, so the condition of complete positivity is consistent in this sense (unlike maps that are only PTP).

The definition of complete positivity presented here (i.e., PTP plus the condition 3) is mathematically equivalent to a plethora of many different sets of conditions; a good portion of those alternative definitions make an explicit reference to entangled states (see [7] for detailed discussion on this subject). However, from the point of view of our current discussion, there is one alternative definition that is of particular interest: the Choi's theorem asserts that all CP maps have Kraus representation and all maps with Kraus representation are CP. This is the origin of the nomenclature; since CP is mathematically equivalent to Kraus representation, it is okay to call maps with Kraus representation “CP maps”.

§.§ Time-dependent Kraus representations

The set of Kraus operators constructed according to the “first principles” blueprint (<ref>) represents a “snapshot” of the dynamical map taken at certain point in time $t$. When one consistently follows the blueprint for all points in time,

\begin{align}
\forall_{t>0}\ \hat K_{jj'}(t) \equiv \sqrt{\rho_j}\langle j'|e^{-it\hat H_{SE}}|j\rangle,
\end{align}

then these $t$-dependent Kraus representations can be logically conjoined together under one family parametrized by the passage of time, $\{\hat K_{jj'}(t) \}_{j,j'; t\geqslant 0}$, such that $ \forall_{j,j'}\,\hat K_{jj'}(s)\to\hat K_{jj'}(t)$ when $s\to t$. As a result, one naturally upgrades the collection of disjointed temporal snapshots of the dynamical map into the Kraus representation of the whole dynamics of the open system.

But what about the effects-oriented point of view? How to introduce time-dependence when Kraus representations are abstracted from any underlying physical $SE$ dynamics? Of course, without the underpinning of unitary evolution of the total system, there is no reason for any form of temporal relation between two abstract representations; such a relation can be established only by introducing additional constraints on the form of Kraus operators. One such constraint, that is known to lead to useful parametrization, is to postulate that the map represented by the $t$-dependent family of Kraus operators $\{\hat K_m(t)\}_{m,t\geqslant 0}$ obeys the simplest form of the composition rule,

\begin{align}
\left(\sum_m \hat K_m(t)\bullet\hat K_m^\dagger(t)\right)\left(\sum_m\hat K_m(s)\bullet\hat K_m^\dagger(s)\right) = \sum_m \hat K_m(t+s)\bullet\hat K_m^\dagger(t+s).
\end{align}

As we have learned previously (e.g., Sec. <ref>), this composition rule implies time-independent dynamical equation (in fact, the two are equivalent),

\begin{align}
\frac{d}{dt}\sum_m\hat K_m(t)\bullet\hat K_m^\dagger(t) = \mathcal{L}_\text{GKLS}\sum_m\hat K_m(t)\bullet\hat K_m^\dagger(t);\quad \sum_m\hat K_m(0)\bullet\hat K_m^\dagger(0) = \bullet,
\end{align}

where the super-operator $\mathcal{L}_\text{GKLS}$ is constant and it acts as a generator of the map, i.e.,

\begin{align}
\sum_{m}\hat K_m(t)\bullet\hat K_m^\dagger(t) &= e^{t\mathcal{L}_\text{GKLS}}.
\end{align}

The enforced composition rule by itself does not tell us anything else about $\mathcal{L}_\text{GKLS}$. However, $\hat K_m(t)$'s also satisfy the standard Kraus condition, $\forall_{t\geqslant 0}\,\sum_m\hat K_m^\dagger(t)\hat K_m(t)= \hat {1}$, which, of course, constraints the allowed forms of the generator; in particular, it was shown by the trio of V. Gorini, A. Kossakowski, and G. Sudarshan [8] and independently by G. Lindblad [9] (GKLS in alphabetic order) that $\mathcal{L}_\text{GKLS}$ must be of the GKLS form:

\begin{align}
\mathcal{L}_\text{GKLS} = -i[\hat H,\bullet] +\sum_{\alpha,\beta}\gamma_{\alpha,\beta}\left(\hat L_\alpha{\bullet}\hat L_\beta^\dagger - \frac{1}{2}\{\hat L_\beta^\dagger\hat L_\alpha,\bullet\}\right),
\end{align}

where $\hat H$ is an arbitrary hermitian operator, $\{\hat L_\alpha\}_\alpha$ is an arbitrary set of operators, and—this is key—$\gamma_{\alpha,\beta}$ is a positive matrix. The resultant dynamical equation for the density matrix of $S$,

\begin{align}
\frac{d}{dt}\hat\rho_S(t) = -i[\hat H,\hat\rho_S(t)] + \sum_{\alpha,\beta}\gamma_{\alpha,\beta}\left(\hat L_\alpha\hat\rho_S(t)\hat L_\beta^\dagger-\frac{1}{2}\{\hat L_\beta^\dagger\hat L_\alpha,\hat\rho_S(t)\}\right)
\end{align}

is called the GKLS equation or simply the master equation (see, e.g., [10] for a simplified proof).

The theorem that super-operator in the GKLS form generates CP map can be considered as the foundational result for the theory of open systems. It would not be an exaggeration to say that the discovery of GKLS form was the moment when the open system theory was born as a distinct field in physics. Indeed, see [11, 12] as an example of works on open systems where GKLS maps are treated as the central concept of the theory.

§.§ Master equation as a physical model

At the time, the result obtained by GKLS captured imagination of many; a mathematical characterization of valid dynamical maps (at least, of a class of maps with constant generators) is a supremely strong result indeed. Of course, people started to wonder how this result would apply to physics: can it be argued that dynamical maps in GKLS form are sufficiently accurate models of real open systems? And if they are, then what is the physical meaning behind operators $\hat L_\alpha$ and the matrix $\gamma_{\alpha,\beta}$? Essentially, this is a question of link between the effect-oriented and first principles approaches.

The first person to answer these questions was E.B. Davies, who demonstrated that the dynamics induced by the coupling with a generic environment tends to the GKLS form when one can justify the application of a sequence of intuitive approximations; here, we recreate their reasoning using the modern tools we have developed previously. To begin with, we assume that $E$ is initialized in the state of thermal equilibrium, $\hat\rho_E=e^{-\beta\hat H_E}/Z$—the immediate consequence is that $\hat H_S' = \hat H_S + \lambda\langle F\rangle\hat V$ is constant. Then, we apply the following approximations,

* Born approximation: relaying on the estimate (<ref>), we assume the weak coupling $\lambda\tau_c\ll 1$ and the separation of time scales $\tau_c\ll \tau_S$, so that we can truncate the cSQ series at the lowest order,

\begin{align}
\langle\mathcal{U}\rangle_c(t,0) &\approx \langle\mathcal{U}\rangle^{(2)}(t,0) = T_{\mathcal{L}_c^{(k)}}e^{-\lambda^2\int_0^t \mathcal{L}_c^{(2)}(s)ds}.
\end{align}

* Markov approximation: First, we exploit the stationarity of $E$, brought about by the initial thermal state, to shift the integration variable in the remaining cSQ,

\begin{align}
\nonumber
\int_0^t\mathcal{L}_c^{(2)}(s)ds &=\int_0^t\!\!ds[\hat V_c(s) ,\bullet] \int_0^s\!\!du
\left(\langle\!\langle F(s)F(u)\rangle\!\rangle \hat V_c(u){\bullet}-{\bullet}\hat V_c(u)\langle\!\langle \bar{F}(s)\bar{F}(u)\rangle\!\rangle\right)\\
\nonumber
&=\int_0^t\!\!ds\,\mathcal{U}_S'(-s)[\hat V,\bullet]\\
\nonumber
\int_0^s\!\!du\left(\langle\!\langle F(s-u)F\rangle\!\rangle \hat V_c(u-s){\bullet}-{\bullet}\hat  V_c(u-s)\langle\!\langle F(s-u)F\rangle\!\rangle^*\right)
\mathcal{U}_S'(s)\\
&=\int_0^t\!\!\!ds\,\mathcal{U}_S'(-s)[\hat V,\bullet]
\int_0^s\!\!\!d\tau\left(\langle\!\langle F(\tau)F\rangle\!\rangle \hat V_c(-\tau){\bullet}-{\bullet}\hat  V_c(-\tau)\langle\!\langle F(\tau)F\rangle\!\rangle^*\right)
\mathcal{U}_S'(s),
\end{align}

where we used the relations

\begin{align*}
&\langle\!\langle \bar{F}(s)\bar{F}(u)\rangle\!\rangle = \langle\!\langle F(s)F(u)\rangle\!\rangle^*,\\
&\mathcal{U}_S'(s) = e^{-is\mathcal{H}_S'} = e^{-is[\hat H_S',\bullet]} = \hat U_S'(s){\bullet}\hat U_S'(-s),\\
&\hat U_S'(s) = e^{-is\hat H_S'} = e^{-is(\hat H_S+\lambda\langle F\rangle\hat V)},\\
&[\hat V_c(s),\bullet] = \mathcal{U}_S'(-s)[\hat V,\bullet]\mathcal{U}_S'(s),\\
&\hat V_c(s-u)\bullet = \mathcal{U}_S'(s)(\hat V_c(u)\bullet)\mathcal{U}_S'(-s),\\
&{\bullet}\hat V_c(s-u) = \mathcal{U}_S'(s)(\bullet\hat V_c(u))\mathcal{U}_S'(-s),\\
&{\bullet} =\mathcal{U}_S'(-s)\mathcal{U}_S'(s).
\end{align*}

Then, we apply the approximation where we extend the upper limit of the integral over $\tau$ to infinity,

\begin{align}
\nonumber
\int_0^t\mathcal{L}_c^{(2)}(s)ds & \approx \int_0^t\!\!ds\, \mathcal{U}_S'(-s)[ \hat V , \bullet]\\
\nonumber
\left(\int_0^\infty\!\!\!\!\langle\!\langle F(\tau)F\rangle\!\rangle\hat V_c(-\tau)d\tau\right){\bullet} - {\bullet}\left(\int_0^\infty\!\!\!\! \langle\!\langle F(\tau)F\rangle\!\rangle\hat V_c(-\tau)d\tau\right)^\dagger
\right)
\mathcal{U}_S'(s)\\
&\equiv \int_0^t\!\!ds\,\mathcal{U}_S'(-s)[\hat V, \hat\Lambda{\bullet} - {\bullet}\hat\Lambda^\dagger ]\mathcal{U}_S'(s) \equiv \int_0^t \mathcal{L}_R(s)ds,
\end{align}

where we have defined, now time-independent, operator

\begin{align}
\hat\Lambda = \int_0^\infty  \langle\!\langle F(\tau)F \rangle\!\rangle \hat V_c(-\tau)d\tau.
\end{align}

Given that we are already assuming the weak coupling, $\lambda\tau_c\ll 1$, the approximation is justified under the condition that $t$ is much larger than the correlation time $\tau_c$. When $t\gg\tau_c$, the effective upper limit of the integral over $\tau$ is not $t$ but rather $\tau_c$—the range of qumulants. Then, we can use a crude estimate $|\langle\!\langle F(\tau)F\rangle\!\rangle|\sim \theta(\tau_c-\tau)$, and we write

\begin{align}
\nonumber
&\lambda^2\left|\,\Big|\Big|\int_0^t\mathcal{L}_c^{(2)}(s)ds\Big|\Big|-\Big|\Big|\int_0^t\mathcal{L}_R(s)ds\Big|\Big|	\,\right|\\
\nonumber
&\phantom{=}	\sim\lambda^2\left|\int_0^t\!\!ds\left(\int_0^s\theta(\tau_c-\tau)d\tau -\int_0^\infty \theta(\tau_c-\tau)d\tau\right)\right|\\
\nonumber
=\lambda^2\left|\left(t\tau_c - \frac{1}{2}\tau_c^2\right) - t\tau_c\right|
= \frac{1}{2}\lambda^2\tau_c^2\ll 1.
\end{align}

Hence, when $t\gg\tau_c$ and $\lambda\tau_c\ll 1$, the difference between the original cSQ and the approximating super-operator is negligible.

If we were to stop here, we would obtain the so-called Redfield equation for density matrix of $S$,

\begin{align}
\nonumber
\frac{d}{dt}\hat\rho_S(t)&=\frac{d}{dt}\langle\mathcal{U}\rangle(t,0)\hat\rho_S =\frac{d}{dt}\mathcal{U}_S'(t)\langle\mathcal{U}\rangle_c(t,0)\hat\rho_S\\
\nonumber
&=	-i\mathcal{H}_S'\mathcal{U}_S'(t)\langle\mathcal{U}\rangle_c(t,0)\hat\rho_S+\mathcal{U}_S'(t)\frac{d}{dt}\langle \mathcal{U}\rangle_c(t,0)\hat\rho_S\\
\nonumber
&\approx -i\mathcal{H}_S'\langle\mathcal{U}\rangle(t,0)\hat\rho_S
-\lambda^2\mathcal{U}_S'(t)\mathcal{U}_S'(-t)[\hat V,\hat\Lambda{\bullet}-{\bullet}\hat\Lambda^\dagger]\mathcal{U}_S'(t)\langle\mathcal{U}_c\rangle(t,0)\hat\rho_S\\
&= -i[\hat H_S+\lambda\langle F\rangle\hat V,\hat\rho_S(t)] -\lambda^2 [\hat V,\hat\Lambda\hat\rho_S(t)-\hat\rho_S(t)\hat\Lambda^\dagger],
\end{align}

and its formal solution reads

\begin{align}
\hat\rho_S(t) = \exp\left(-it[\hat H_S+\lambda\langle F\rangle\hat V,\bullet]-\lambda^2t[\hat V,\hat\Lambda{\bullet}-{\bullet}\hat\Lambda^\dagger]\right)\hat\rho_S.
\end{align}

* Secular approximation (also known as rotating wave approximation): we introduce the frequency decomposition of the coupling operator obtained using the eigenkets of $\hat H_S' = \sum_a \epsilon_a |a\rangle\langle a|$,

\begin{align}
\nonumber
\hat V_c(s) &= \hat U_S'(-s)\hat V\hat U_S'(s) = \sum_{a,b}\hat U_S'(-s)|a\rangle\langle a|\hat V|b\rangle\langle b|\hat U_S'(s)\\
\nonumber
& = \sum_{a,b}e^{is(\epsilon_a-\epsilon_b)}|a\rangle V_{ab}\langle b|	= \sum_{\omega}e^{is\omega}\sum_{a,b}\delta(\omega-\epsilon_a+\epsilon_b)|a\rangle V_{ab}\langle b|\\
\label{eq:freq_decomp}
&\equiv \sum_\omega e^{is\omega}\hat V_\omega = \sum_{\omega}e^{-is\omega}\hat V_\omega^\dagger.
\end{align}

We use the frequency decomposition to transform the previously obtained operators $\hat\Lambda$,

\begin{align}
\hat\Lambda &= \int_0^\infty d\tau\langle\!\langle F(\tau)F\rangle\!\rangle \hat V_c(-\tau) = \sum_\omega \left(\int_0^\infty \langle\!\langle F(\tau)F\rangle\!\rangle e^{-i\omega\tau}d\tau\right)\hat V_\omega
\equiv \sum_\omega \gamma_\omega\hat V_\omega,
\end{align}

as well as the super-qumulant,

\begin{align}
\nonumber
\int_0^t\!\!\mathcal{L}_R(s)ds &= \int_0^t\!\!ds\,\mathcal{U}_S'(-s)[\hat V,\hat\Lambda{\bullet}-{\bullet}\hat\Lambda^\dagger]\mathcal{U}_S'(s)\\
\nonumber
&=\int_0^t\!\!ds\left(\hat U_S'(-s)\hat V\hat U_S'(s)\hat U_S'(-s)\hat\Lambda\hat U_S'(s){\bullet}{-}\hat U_S'(-s)\hat\Lambda\hat U_S'(s){\bullet}\hat U_S'(-s)\hat V\hat U_S'(s)\right.\\
\nonumber
{\bullet}\hat U_S'(-s)\hat\Lambda^\dagger\hat U_S'(s)\hat U_S'(-s)\hat V\hat U_S'(s){-}\hat U_S'({-}s)\hat V\hat U_S'(s){\bullet}\hat U_S'(-s)\hat \Lambda^\dagger\hat U_S'(s)\right)\\
\gamma_\omega\hat V^\dagger_{\omega'}\hat V_\omega{\bullet}+{\bullet}\hat V^\dagger_{\omega'}\hat V_\omega\gamma_{\omega'}^*
-\left(\gamma_\omega+\gamma_{\omega'}^*\right)\hat V_\omega{\bullet}\hat V^\dagger_{\omega'}
\right)\int_0^te^{is(\omega-\omega')}ds.
\end{align}

As a result, we have confined the totality of $t$-dependence to phase factors.

The approximation we are about to implement is based on the observation that when the frequencies match up, $\omega = \omega'$, the phase factor simplifies to $t$,

\begin{align}
\int_0^t e^{is(\omega-\omega')}ds\bigg|_{\omega=\omega'} = \int_0^t ds = t,
\end{align}

which is a very different behavior than the mismatched case, $\omega\neq\omega'$,

\begin{align}
\int_0^t e^{i(\omega-\omega')s}ds &= e^{\frac{1}{2}i(\omega-\omega')t}t\operatorname{sinc}\left[\frac{(\omega-\omega')t}{2}\right].
\end{align}

Therefore, the phase factors with matching frequencies—on-resonance terms—overshadow the off-resonance terms on a sufficiently long time scale,

\begin{align}
\frac{\left|\int_0^t e^{is(\omega-\omega')}ds\right|}{\left|\int_0^t ds\right|} = \left|\operatorname{sinc}\left[\frac{(\omega-\omega')t}{2}\right]\right|
\xrightarrow{|\omega-\omega'|t\gg 2\pi} 0.
\end{align}

To obtain this effect for all off-resonance terms, $t$ has to be long in comparison to the slowest time-scale in $S$, or

\begin{align}
\min_{\omega\neq\omega'}|\omega-\omega'| t  \equiv \frac{2\pi}{T_S}t\gg 2\pi \ \Rightarrow\  t\gg T_S.
\end{align}

When this condition is in effect, then all “off-diagonal” terms where $\omega\neq\omega'$ can be neglected in comparison with the “diagonal” $\omega=\omega'$ terms,

\begin{align}
\nonumber
\int_0^t\!\!\mathcal{L}_R(s)ds&\approx t\sum_{\omega}\left(
\gamma_\omega\hat V^\dagger_\omega\hat V_\omega{\bullet}+{\bullet}\hat V^\dagger_\omega\hat V_\omega\gamma_\omega^*
-\left(\gamma_\omega+\gamma_\omega^*\right)\hat V_\omega{\bullet}\hat V^\dagger_\omega
\right)
\\
\nonumber
&=it\left[ \sum_\omega \operatorname{Im}\{\gamma_\omega\}\hat V^\dagger_\omega\hat V_\omega, \bullet\right]
-t\sum_\omega 2\operatorname{Re}\{\gamma_\omega\}\left(\hat V_\omega{\bullet}\hat V^\dagger_\omega-\frac{1}{2}\{\hat V^\dagger_\omega\hat V_\omega,\bullet\}\right)\\
&\equiv -t \mathcal{L}_D.
\end{align}

Now note the following: (i) the operator $\sum_\omega \operatorname{Im}\{\gamma_\omega\}\hat V^\dagger_\omega\hat V_\omega$ is hermitian; (ii) the coefficients $\operatorname{Re}\{\gamma_\omega\}$ are positive because

\begin{align}
\nonumber
2\operatorname{Re}\{\gamma_\omega\} &= \int_0^\infty \left(\langle\!\langle F(\tau)F\rangle\!\rangle e^{-i\omega\tau}+\langle\!\langle F(\tau)F\rangle\!\rangle^* e^{i\omega\tau}\right)d\tau\\
\nonumber
&=\int_0^\infty \langle\!\langle F(\tau)F\rangle\!\rangle e^{-i\omega\tau}d\tau + \int_{-\infty}^0 \langle\!\langle F(-\tau)F\rangle\!\rangle^* e^{-i\omega\tau}d\tau\\
\nonumber
&= \int_0^\infty\langle\!\langle F(\tau)F\rangle\!\rangle e^{-i\omega\tau}d\tau+\int_{-\infty}^0\langle\!\langle F(\tau)F\rangle\!\rangle e^{-i\omega\tau}d\tau\\
\nonumber
&=\int_{-\infty}^\infty\langle\!\langle F(\tau)F\rangle\!\rangle e^{-i\omega\tau}d\tau = \int_{-\infty}^\infty\left[C_F(\tau,0)-iK_F(\tau,0)\right]e^{-i\omega\tau}d\tau\\
\label{eq:gamma}
&= S(\omega)-i\kappa(\omega) = S(\omega)\left(1+\frac{e^{-\beta\omega}-1}{e^{-\beta\omega}+1}\right) =\frac{2S(\omega)}{e^{\beta\omega}+1},
\end{align}

and the spectral density $S(\omega)$ cannot be negative or, otherwise, it would violate the second law of thermodynamics (we used the fluctuation-dissipation theorem to eliminate susceptibility $\kappa(\omega)$, see assignment <ref> of Sec. <ref>). Note that it is possible to prove the positivity of $\operatorname{Re}\{\gamma_\omega\}$ without resorting to physical arguments [13].

Thus, the approximated map we get at the end,

\begin{align}
\langle\mathcal{U}\rangle_c(t,0) \approx e^{\lambda^2 t\mathcal{L}_D},
\end{align}

has the GKLS form; indeed, we see that $\lambda^2\mathcal{L}_D = \mathcal{L}_\mathrm{GKLS}$ when we identify

\begin{align*}
\hat H &\leftrightarrow -\lambda^2\sum_\omega\operatorname{Im}\{\gamma_\omega\}\hat V^\dagger_\omega\hat V_\omega,\\
\{\hat L_\alpha\}_\alpha &\leftrightarrow \{\hat V_\omega\}_\omega,\\
\gamma_{\alpha,\beta} &\leftrightarrow \gamma_{\omega,\omega'} = 2\lambda^2\operatorname{Re}\{\gamma_\omega\}\delta_{\omega,\omega'}.
\end{align*}

The matrix $\gamma_{\omega,\omega'}$ is, of course, positive because its eigenvalues, $2\lambda^2\operatorname{Re}\{\gamma_\omega\}$, are nonnegative numbers.

The corresponding master equation for the open system density matrix,

\begin{align}
\nonumber
\frac{d}{dt}\hat\rho_S(t) &= \frac{d}{dt}\mathcal{U}_S'(t)e^{\lambda^2t\mathcal{L}_D}\hat\rho_S = -i\mathcal{H}_S'\hat\rho_S(t)
\nonumber
&=\left(-i\mathcal{H}_S' + \lambda^2\mathcal{L}_D\right)\hat\rho_S(t)\\
\nonumber
&=-i\left[\hat H_S +\lambda\langle F\rangle\hat V -\lambda^2\sum_\omega\operatorname{Im}\{\gamma_\omega\}\hat V_\omega^\dagger\hat V_\omega\ ,\ \hat\rho_S(t)\right]\\
\label{eq:Davies_eqn}
{+}\lambda^2\sum_\omega\operatorname{Re}\{\gamma_\omega\}\left(2\hat V_\omega\hat\rho_S(t)\hat V_\omega^\dagger-\{\hat V_\omega^\dagger\hat V_\omega,\hat\rho_S(t)\}\right),
\end{align}

is called the Davies master equation.

[[sol:RWA]Solution] Show that

\begin{align*}
\mathcal{U}_S'(t)\mathcal{L}_D\mathcal{U}_S'(-t) = \mathcal{L}_D.
\end{align*}

[[sol:pauli_ME]Solution] Assuming that the Hamiltonian $\hat H_S' = \sum_a\epsilon_a|a\rangle\langle a|$ is nondegenerate,

\begin{align*}
\forall_{a\neq b}\epsilon_a\neq \epsilon_b,
\end{align*}

derive the Pauli master equation: the rate equation obtained from Davies equation by restricting it to the diagonal elements of the density matrix,

\begin{align*}
\frac{d}{dt}\langle a|\hat\rho_S(t)|a\rangle &=\sum_b M_{a,b}\langle b|\hat\rho_S(t)|b\rangle\\
&= \sum_{b}\lambda^2\left(\gamma_{a,b} - \delta_{a,b}\sum_c \gamma_{c,a}\right)\langle b|\hat\rho_S(t)|b\rangle,
\end{align*}

with $\gamma_{a,b} = 2|V_{ab}|^2 S(\omega_{ab})/(e^{\beta\omega_{ab}}+1)$.

Note that the thermal equilibrium state with respect to $\hat H_S'$,

\begin{align*}
\hat p &= \frac{e^{-\beta\hat H_S'}}{\operatorname{tr}_S\left(e^{-\beta\hat H_S'}\right)} = \frac{1}{Z}\sum_a e^{-\beta\epsilon_a}|a\rangle\langle a|,
\end{align*}

is a steady state solution to Pauli equation. This can be easily showed by using the detailed balance property of the transition matrix,

\begin{align*}
\gamma_{b,a} &= |V_{ba}|^2\frac{2S(\omega_{ba})}{e^{\beta\omega_{ba}}+1}
= e^{-\beta\omega_{ba}}|V_{ab}|^2 \frac{2S(-\omega_{ab})}{1+e^{-\beta\omega_{ba}}}
= e^{\beta\omega_{ab}}|V_{ab}|^2\frac{2S(\omega_{ab})}{e^{\beta\omega_{ab}}+1}
= e^{\beta\omega_{ab}}\gamma_{a,b},
\end{align*}

to transform the Pauli equation to

\begin{align*}
\frac{d}{dt}\langle a|\hat\rho_S(t)|a\rangle &= \lambda^2\langle a|\mathcal{L}_D\hat\rho_S(t)|a\rangle \\
&= \lambda^2\sum_b\left(\gamma_{a,b}\langle b|\hat\rho_S(t)|b\rangle-\gamma_{b,a}\langle a|\hat\rho_S(t)|a\rangle \right)\\
&= \lambda^2\sum_b\gamma_{a,b}\left(\langle b|\hat\rho_S(t)|b\rangle - e^{\beta\omega_{ab}}\langle a|\hat\rho_S(t)|a\rangle\right).
\end{align*}

Then, when we substitute $\hat p$ for $\hat\rho_S(t)$ on the RHS we get

\begin{align*}
\lambda^2\langle a|\mathcal{L}_D\hat p|a\rangle = \frac{\lambda^2}{Z}\sum_b \gamma_{a,b}\left(e^{-\beta\epsilon_b}-e^{\beta(\epsilon_a-\epsilon_b)}e^{-\beta\epsilon_a}\right) = 0.
\end{align*}

In addition, it can be showed that any solution of Pauli master equation satisfies

\begin{align*}
\lim_{t\to\infty}\langle a|\hat\rho_S(t)|a\rangle = \langle a|\hat p|a\rangle.
\end{align*}

Moreover, with some additional assumptions about operator $\hat V$, one can prove that all solutions of the Davies equation tend to the equilibrium state as the time goes to infinity, which gives some insight into the process of thermalization. We will talk about this topic in more detail later in Chapter <ref>.

§.§ Nakajima–Zwanzig equation

In previous section we have derived the Davies master equation starting from the super-qumulant expansion; this is not a standard textbook route, however. The typical textbook derivation (e.g., see [14]) starts with the von Neumann equation for the density matrix of the total system,

\begin{align}
\frac{d}{dt}\hat R_{SE}(t) &= -i\lambda[ \hat V_c(t)\otimes\hat X_I(t),\hat R_{SE}(t)],
\end{align}

where $\hat R_{SE}(t) = \mathcal{U}_S'(-t)\otimes e^{it[\hat H_E,\bullet]}\hat\rho_{SE}(t)$ is the total density matrix in the central picture. The equation can be formally integrated with the assumed initial condition $\hat R_{SE}(0) = \hat\rho_S\otimes\hat\rho_E$,

\begin{align}
\hat R_{SE}(t) = \hat\rho_S\otimes\hat\rho_E -i \lambda\int_0^t [ \hat V_c(s)\otimes\hat X_I(s),\hat R_{SE}(s) ]ds.
\end{align}

This formal solution is then substituted back into RHS of the equation and the partial trace over $E$ is taken,

\begin{align}
\label{eq:iterated_vonNeumann}
\frac{d}{dt}\hat\rho_c(t) = -\lambda^2\int_0^t \operatorname{tr}_E\left([\hat V_c(t)\otimes\hat X_I(t),[\hat V_c(s)\otimes\hat X_I(s),\hat R_{SE}(s)]]\right)ds,
\end{align}

where, of course, $\operatorname{tr}_E(\hat R_{SE}(t)) = \hat\rho_c(t)$. The issue is that the above formula is not a dynamical equation because there is no $\hat\rho_c$ on the RHS; before we can proceed, something has to be done about it.

Here comes the pivotal—and the most suspect—moment of the derivation: to “close” the equation it is assumed that the total density matrix preserves the product form at all times, plus, that the state on the $E$-side remains static,

\begin{align}\label{eq:BM_ansatz}
\hat R_{SE}(s) \approx \hat\rho_c(s)\otimes\hat\rho_E.
\end{align}

This step is typically referred to as the Born-Markov approximation—we will discuss later how one might try to justify this ansatz and if “approximation” is a correct descriptor for it. The ansatz itself, of course, closes the equation,

\begin{align}
\nonumber
\frac{d}{dt}\hat\rho_c(t) &\approx -\lambda^2\int_0^t \operatorname{tr}_E\left([\hat V_c(t)\otimes\hat X_I(t),[\hat V_c(s)\otimes\hat X_I(s),\hat\rho_c(s)\otimes\hat\rho_E]]\right)ds\\
\label{eq:closed_eqn}
&= -\lambda^2\int_0^t \langle\!\langle \mathcal{W}_c(X(t),\bar{X}(t),t)\mathcal{W}_c(X(s),\bar{X}(s),s)\rangle\!\rangle\hat\rho_c(s) ds.
\end{align}

It is not over yet, as the equation is non-local in time, i.e., $\hat\rho_c(s)$ on RHS depends on the integrated variable. To fix this problem, the density matrix has to be brought outside the integral,

\begin{align}
\frac{d}{dt}\hat\rho_c(t) &\approx -\lambda^2\int_0^t \langle\!\langle\mathcal{W}_c(X(t),\bar{X}(t),t)\mathcal{W}_c(X(s),\bar{X}(s),s)\rangle\!\rangle ds \hat\rho_c(t)
= -\lambda^2\mathcal{L}_c^{(2)}(t)\hat\rho_c(t).
\end{align}

A rough estimation shows that the resultant error is $\sim\lambda\tau_c$ times smaller than the leading term, hence, it can be neglected when the coupling is weak, $\lambda\tau_c\ll 1$. The equation—now local in time—can be solved without difficulty,

\begin{align}
\hat\rho_c(t) =T_{\mathcal{L}_c^{(k)}}e^{-\lambda^2\int_0^t \mathcal{L}_c^{(2)}(s)ds}\hat\rho_S = \langle\mathcal{U}\rangle^{(2)}(t,0)\hat\rho_S.
\end{align}

At last, we have arrived at the first step of the SQ-based derivation. From this point on, we can get the Davies equation by simply replicating all the steps described in the previous section.

Now let us get back to the Born-Markov approximation. Typically, textbooks try to argue that the form (<ref>) is a genuine approximation, and some of the arguments do sound semi-convincingly. The issue is that it has been demonstrated in numerical experiments that the assumption about product form is simply incorrect [15]—the Born-Markov ansatz does not describe the real form of the total density matrix. So, if it is not an adequate approximation, then why it is still invoked as such an integral part of the textbook derivation? The honest answer is that it is the simplest way to force the issue of closing the Eq. (<ref>). Remember, the SQ expansion—which does not even have this kind of problems because it produces dynamical equation by design—is too recent to make it to traditional textbooks.

A less problematic, but orders of magnitude more complicated, derivation still makes use of the Born-Markov ansatz but it abandons completely any notions that it is some form of approximation. Instead, the ansatz is implemented as a formal super-operator acting in the total $SE$ operator space,

\begin{align}
\mathcal{P}\hat A_{SE} &= \left(\operatorname{tr}_E({\bullet})\otimes\hat\rho_E\right)\hat A_{SE} = \operatorname{tr}_E(\hat A_{SE})\otimes\hat\rho_E.
\end{align}

For example, when $\mathcal{P}$ acts on the total density matrix $\hat R_{SE}(s)$, it effectively enforces the ansatz. The projection super-operator is then used on the von Neumann equation, and after some algebra and the most rudimentary calculus [16], one is lead to the Nakajima-Zwanzig equation (NZE),

\begin{align}
\frac{d}{dt}\mathcal{P}\hat R_{SE}(t) &= -\lambda^2\int_0^t  \mathcal{P}[\hat V_c(t)\otimes\hat X_I(t),\bullet]
T_{\mathcal{G}}e^{-i\lambda\int_s^t \mathcal{G}(u)du}\mathcal{G}(s)\mathcal{P}\hat R_{SE}(s)ds,
\end{align}

where $\mathcal{G}(s) = \left( {\bullet} - \mathcal{P} \right)[\hat V_c(s)\otimes\hat X_I(s),\bullet]$.

Taking the partial trace over $E$ transforms NZE into a time-non-local equation for the system's density matrix,

\begin{align}
\frac{d}{dt}\hat\rho_c(t) &= -\lambda^2 \int_0^t \operatorname{tr}_E\left([\hat V_c(t)\otimes\hat X_I(t), T_{\mathcal{G}}e^{-i\lambda\int_s^t\mathcal{G}(u)du}\mathcal{G}(s)\bullet\otimes\hat\rho_E]\right)
\hat\rho_c(s)ds.
\end{align}

Now, when we take the above equation as an entry point instead of Eq. (<ref>), then we can get to the closed equation (<ref>) by expanding the exp into the power series,

\begin{align}
T_\mathcal{G}e^{-i\lambda\int_s^t\mathcal{G}(u)du} &= \bullet -i\lambda\int_s^t\mathcal{G}(u)du+\ldots,
\end{align}

and then cutting off everything above the lowest $\lambda^2$-order,

\begin{align}
\nonumber
\frac{d}{dt}\hat\rho_c(t) &\approx -\lambda^2\int_0^t \operatorname{tr}_E\left([\hat V_c(t)\otimes\hat X_I(t),
\left(\bullet-\mathcal{P}\right)[\hat V_c(s)\otimes\hat X_I(s),\bullet\otimes\hat\rho_E]
\nonumber
&= -\lambda^2\int_0^t \langle\!\langle\mathcal{W}_c(X(t),\bar{X}(t),t)\mathcal{W}_c(X(s),\bar{X}(s),s)\rangle\!\rangle \hat\rho_c(s)ds\\
\nonumber
&\phantom{=} +\lambda^2\int_0^t \langle X(t)\rangle\langle X(s)\rangle [\hat V_c(t),[\hat V_c(s),\hat\rho_c(s)]ds\\
&=-\lambda^2\int_0^t \langle\!\langle\mathcal{W}_c(X(t),\bar{X}(t),t)\mathcal{W}_c(X(s),\bar{X}(s),s)\rangle\!\rangle \hat\rho_c(s)ds.
\end{align}

And thus, we get the same result as before but without invoking the Born-Markov “approximation”. Overall, it can be shown that the solution to the Davies master equation (<ref>)—the end point of all those deliberations—tends to the solution of Nakajima–Zwanzig equation asymptotically in the so-called weak coupling limit, where $\lambda\to 0$ while $t\to\infty$ at the rate $\sim\lambda^{-2}$ [17].

§ DYNAMICS OF THERMALISATION

§.§ Thermalisation

Thermalisation is a process where the state of system $S$ is brought over time to the thermal equilibrium due to interactions with the other system $E$ that already is in the state of equilibrium. An intuitive expectation is that the thermalisation occurs for most systems open to thermal baths (i.e., “large” $E$ initialized in thermal equilibrium), especially when the coupling between $S$ and $E$ is weak. It should be kept in mind, however, that thermalization is not guaranteed—it is quite possible that some peculiarities of the structure of energy levels of $S$ or the symmetries of the coupling will close off the path to the equilibrium. Here we will discuss the necessary conditions for the thermalisation and we will describe the dynamics of the process itself using the formalism of the open system theory. We will restrict our considerations to the weak coupling regime, and thus, the results and conclusions presented below do not exhaust the topic by any means, shape, or form. Thermalisation is still an open problem and a very much active field of research.

Formally, we say that the system thermalizes when

\begin{align}
\forall_{\hat\rho_S}\ \langle\mathcal{U}\rangle(t,0)\hat\rho_S \xrightarrow{t\to\infty} \hat p = \frac{1}{\sum_a e^{-\beta\epsilon_a}}\sum_a e^{-\beta \epsilon_a}|a\rangle\langle a|
= \frac{1}{Z}\sum_a e^{-\beta\epsilon_a}|a\rangle\langle a|,
\end{align}

where $\{\epsilon_a\}_a$ and $\{|a\rangle\}_a$ are the eigenvalues and eigenkets of the operator representing the energy in system $S$ (this operator does not have be equal to the free Hamiltonian $\hat H_S$ of $S$) and $\beta = (k_B T)^{-1}$ is the inverse temperature determined by the initial thermal state of $E$,

\begin{align}
\hat\rho_E =\frac{e^{-\beta\hat H_E}}{\operatorname{tr}(e^{-\beta\hat H_E})}.
\end{align}

In the case considered here, the dynamical map $\langle\mathcal{U}\rangle$ is induced by the total Hamiltonian

\begin{align}
\hat H_{SE} = \hat H_S\otimes\hat 1 + \hat 1\otimes\hat H_E + \lambda \hat V\otimes\hat F,
\end{align}

and we assume that the coupling is weak, so that the Born approximation (see Sec. <ref>) is sufficiently accurate,

\begin{align}
\langle\mathcal{U}\rangle(t,0) &\approx e^{-it[\hat H_S+\lambda\langle F\rangle \hat V,\bullet]}T_{\mathcal{L}_c^{(k)}}e^{-\lambda^2\int_0^t\mathcal{L}_c^{(2)}(s)ds}.
\end{align}

Since the process of thermalisation occurs over long period of time we also naturally satisfy the conditions for the Markov and secular approximations from Sec. <ref> that lead us to the Davies generator,

\begin{align}
\mathcal{L}_D &= -i[\hat H_D,\bullet]
+\sum_{\omega\in \Omega}\frac{2S(\omega)}{e^{\beta\omega}+1}\left(\hat V_\omega{\bullet}\hat V_\omega^\dagger-\frac{1}{2}\{\hat V_\omega^\dagger\hat V_\omega,\bullet\}\right)
\equiv i\mathcal{J}+\mathcal{D},
\end{align}

Here, $\Omega$ is the set of unique frequencies of form $\omega_{ab} = \epsilon_a-\epsilon_b$ where $\epsilon_a$'s are the eigenvalues of the shifted Hamiltonian $\hat H_S' = \hat H_S+\lambda\langle F\rangle\hat V$, $\hat H_S'|a\rangle = \epsilon_a|a\rangle$. The hermitian part of the generator reads

\begin{align}
\hat H_D &= - \sum_{\omega\in\Omega}\operatorname{Im}\left\{\int_0^\infty\langle\!\langle F(\tau)F\rangle\!\rangle e^{-i\omega\tau}d\tau\right\}\hat V_\omega^\dagger\hat V_\omega
= - \sum_{\omega\in\Omega}h(\omega)\hat V_\omega^\dagger\hat V_\omega,
\end{align}

and the frequency decomposition of the coupling operator is defined as

\begin{align}
\hat V_\omega &= \sum_{a,b}\delta(\omega-\omega_{ab})|a\rangle V_{ab} \langle b| = \sum_{(a,b)\in I(\omega)}|a\rangle\langle b|V_{ab} = \hat V^\dagger_{-\omega}.
\end{align}

Here, the index pair $(a,b)$ belonging to the set of pairs $I(\omega)$ (an energy shell), is such that $\omega_{ab} = \omega$; note, e.g., that if $(a,b)\in I(\omega)$, then $(b,a)\notin I(\omega)$, unless $\omega=0$.

§.§ Necessary conditions for thermalisation

There is a theorem cited in a number of textbooks [5, 18] asserting that the Davies map generated by $\mathcal{L}_D$ thermalizes the system with respect to the energy operator $\hat H_S' = \sum_a \epsilon_a |a\rangle\langle a|$,

\begin{align}
\forall_{\hat\rho_S}\ e^{\lambda^2t\mathcal{L}_D}\hat\rho_S \xrightarrow{t\to\infty}\hat p =\frac{e^{-\beta\hat H_S'}}{\operatorname{tr}(e^{-\beta\hat H_S'})} = \frac{1}{Z}\sum_a e^{-\beta\epsilon_a}|a\rangle\langle a|,
\end{align}

when the following conditions are met:

* the system $S$ is ergodic:

\begin{align}
\forall_{\omega\in\Omega}\ [\hat A , \hat V_\omega ]=0\quad\Rightarrow\quad \hat A \propto \hat 1,
\end{align}

* the dissipative part of the Davies generator, $\mathcal{D}$, satisfies the quantum detailed balance,

\begin{align}
\forall_{\hat A,\hat B}\ \operatorname{tr}\left(\hat p \hat A^\dagger \mathcal{D}^\dagger\hat B\right) &=\operatorname{tr}\left(\hat p (\mathcal{D}^\dagger\hat A)^\dagger \hat B\right).
\end{align}

These somewhat cryptic conditions can be clarified with a more down-to-earth parameterization. We will show here that the thermalisation occurs when

*  Hamiltonian $\hat H_S'$ is not degenerate,

\begin{align}\label{eq:condition_a}
\forall_{a\neq b}\ \epsilon_a\neq\epsilon_b,
\end{align}

*  the coupling operator $\hat V$ does not have zero off-diagonal matrix elements,

\begin{align}\label{eq:condition_b}
\forall_{a\neq b}\ V_{ab} = \langle a|\hat V|b\rangle \neq 0.
\end{align}

As we have shown in assignment <ref> of Sec <ref>, when $\hat H_S'$ is not degenerate, the matrix of $\mathcal{L}_D$ in basis $\mathcal{B}_S'=\{|a\rangle\langle b|\}_{a,b}$ has the block-diagonal form with respect to the split between subspaces spanned by the projectors $\mathcal{B}_P = \{ |a\rangle\langle a|\}_{a}$ and coherences $\mathcal{B}_C = \{|a\rangle\langle b|\}_{a\neq b}$,

\begin{align}
\mathcal{L}_D &=\left[\begin{array}{cc} [\mathcal{M}]_{\mathcal{B}_P} & 0 \\ 0 & [i\mathcal{J}+\mathcal{G}]_{\mathcal{B}_C} \\
\end{array}\right]_{\mathcal{B}_S'},
\end{align}

where the generator of the Pauli master equation (see assignment <ref>) is given by,

\begin{align}
M_{a,b} = \operatorname{tr}\left(|a\rangle\langle a|\mathcal{L}_D|b\rangle\langle b|\right) = \operatorname{tr}\left(|a\rangle\langle a|\mathcal{M}|b\rangle\langle b|\right) = \gamma_{a,b} - \delta_{a,b}\sum_c\gamma_{c,a},
\end{align}

and the two-part generator in the coherence subspace,

\begin{align}
i J_{ab,cd} &= \operatorname{tr}\left(|b\rangle\langle a|i\mathcal{J}|c\rangle\langle d|\right)
= i\delta_{a,c}\delta_{b,d}\sum_e\left(h(\omega_{ea})|V_{ea}|^2-h(\omega_{eb})|V_{eb}|^2\right);\\
\nonumber
G_{ab,cd} &= \operatorname{tr}\left(|b\rangle\langle a|\mathcal{D}|c\rangle\langle d|\right) = \operatorname{tr}\left(|b\rangle\langle a|\mathcal{G}|c\rangle\langle d|\right)\\
&= 2V_{ac}^*V_{bd}\delta(\omega_{ab}-\omega_{cd})\frac{S(\omega_{ac})}{e^{\beta\omega_{ac}}+1}-\delta_{a,b}\delta_{c,d}\sum_{e}\frac{\gamma_{e,a}+\gamma_{e,b}}{2} ,
\end{align}

with the transition matrix given by

\begin{align}
\gamma_{a,b} &= |V_{ab}|^2\frac{2S(\omega_{ab})}{e^{\beta\omega_{ab}}+1}.
\end{align}

Our set of assumptions is consistent with the assumptions of the general theorem. First, we can show that <ref> and <ref> lead to ergodic system:

\begin{align}
\nonumber
[ \hat A , \hat V_\omega ] &= \sum_{(a,b)\in I(\omega)}\sum_{c,d}A_{cd}V_{ab} [ |c\rangle\langle d| , |a\rangle\langle b| ]\\
&= \sum_{(a,b)\in I(\omega)}V_{ab}\left(|a\rangle\langle b|(A_{aa} - A_{bb})+\sum_{c\neq a}A_{ca}|c\rangle\langle b| - \sum_{d\neq b}A_{bd}|a\rangle\langle d|\right)
= 0,
\end{align}

which then implies

\begin{align}
\forall_{(a,b)\in I(\omega)}\ \left( A_{aa}=A_{bb}\ \mathrm{and}\ \forall_{c\neq a}\forall_{d\neq b}\ A_{ca} = A_{bd}=0 \right),
\end{align}

because $V_{ab}\neq 0$ (assumption <ref>) and, if $(a,b)\in I(\omega)$, then $(c,b),(a,d)\notin I(\omega)$ for $c\neq a, d\neq b$ (assumption <ref>). When we check the condition for commutation for every $\omega\in\Omega$, then we get

\begin{align}
\nonumber
&\forall_{\omega\in\Omega}\forall_{(a,b)\in I(\omega)}\ \left( A_{aa}=A_{bb}\ \mathrm{and}\ \forall_{c\neq a}\forall_{d\neq b}\ A_{ca} = A_{bd}=0 \right)\\
\nonumber
&\Leftrightarrow \forall_{a,b}\ \left( A_{aa}=A_{bb}\ \mathrm{and}\ \forall_{c\neq a}\forall_{d\neq b}\ A_{ca} = A_{bd}=0 \right)\\
&\Rightarrow \hat A \propto \hat 1.
\end{align}

That is, the condition of ergodicity is satisfied.

To show that the the quantum detailed balance is also satisfied, first we need to find the symmetries of generators,

\begin{align}
M_{b,a} &= \gamma_{b,a} - \delta_{b,a}\sum_{e}\gamma_{e,b} = \frac{e^{-\beta\epsilon_{b}}}{e^{-\beta\epsilon_a}}\left(\gamma_{a,b} - \delta_{a,b}\sum_{e}\gamma_{e,a}\right)
\nonumber
G_{cd,ab} &= 2 V_{ca}^*V_{db}\delta(\omega_{cd}-\omega_{ab})\frac{S(\omega_{ca})}{e^{\beta\omega_{ca}}+1}-\delta_{c,a}\delta_{d,b}\sum_e\frac{\gamma_{e,c}+\gamma_{e,d}}{2} \\
\nonumber
&= 2 V_{ac}V_{bd}^* \delta(\omega_{ab}-\omega_{cd})\frac{S(-\omega_{ac})}{e^{-\beta\omega_{ca}}+1}\frac{1}{e^{\beta\omega_{ca}}}
\nonumber
&=\left(2 V_{ac}V_{bd}^* \delta(\omega_{ab}-\omega_{cd})\frac{S(\omega_{ac})}{e^{\beta\omega_{ac}}+1}
- \delta_{a,c}\delta_{b,d}\sum_e\frac{\gamma_{e,a}+\gamma_{e,b}}{2}\right)\frac{e^{-\beta\epsilon_c}}{e^{-\beta\epsilon_a}}\\
\nonumber
i J_{cd,ab} &= i\delta_{a,c}\delta_{b,d}\sum_{e}\left(h(\omega_{ec})|V_{ec}|^2-h(\omega_{ed})|V_{ed}|^2\right)\\
&= i\frac{e^{-\beta\epsilon_c}}{e^{-\beta\epsilon_a}}\delta_{a,c}\delta_{b,d}\sum_{e}\left(h(\omega_{ea})|V_{ea}|^2-h(\omega_{ea})|V_{eb}|^2\right) = -\frac{p_c}{p_a}(i J_{ab,cd})^*.
\end{align}

Now, the proof is straightforward,

\begin{align}
\nonumber
\operatorname{tr}\left(\hat p \hat A^\dagger \mathcal{D}^\dagger\hat B\right) &= \sum_{a}p_a A_{aa}^*\left(\sum_b M_{b,a}B_{bb}\right)
+ \sum_{a\neq b}p_b A_{ab}^* \left(\sum_{c\neq d}G_{cd,ab}^* B_{cd}\right)\\
\nonumber
&=\sum_{a,b}p_b\frac{p_a}{p_b}M_{b,a}A_{aa}^*B_{bb} + \sum_{a\neq b}\sum_{c\neq d} p_d \frac{p_bp_c}{p_dp_a}\frac{p_a}{p_c}G^*_{cd,ab}A_{ab}^*B_{cd}\\
\nonumber
&=\sum_b p_b \left(\sum_a M_{a,b}A_{aa}\right)^* + \sum_{c\neq d}p_d \left(\sum_{a\neq b} e^{\beta(\omega_{ab} - \omega_{cd})}G_{ab,cd}^* A_{ab}\right)^*B_{cd}\\
\nonumber
&=\sum_b p_b \left(\sum_a M_{a,b}A_{aa}\right)^* + \sum_{c\neq d}p_d \left(\sum_{a\neq b} G_{ab,cd}^* A_{ab}\right)^*B_{cd}\\
&=\operatorname{tr}\left(\hat p (\mathcal{D}^\dagger\hat A)^\dagger \hat B\right),
\end{align}

where we have used the fact that $G_{ab,cd}\propto \delta(\omega_{ab}-\omega_{cd})$.

§.§ H-theorem

The diagonal elements of the system's density matrix,

\begin{align}
\rho_a(t) = \langle a|\hat\rho_S(t)|a\rangle = \langle a|e^{\lambda^2t\mathcal{L}_D}\hat\rho_S|a\rangle = \operatorname{tr}\left(|a\rangle\langle a| e^{\lambda^2t\mathcal{M}}\sum_b \rho_b(0)|b\rangle\langle b|\right),
\end{align}

are positive and normalized, $\sum_a \rho_a(t) = \operatorname{tr}(\exp(\lambda^2t\mathcal{L}_D)\hat\rho_S) = \operatorname{tr}(\hat\rho_S) = 1$—therefore, $\rho_a(t)$ can be treated as a time-dependent probability distribution.

One can define the nonnegative quantity called the relative entropy between the equilibrium distribution $p_a = \langle a|\hat p|a\rangle$ and $\rho_a(t)$,

\begin{align}
H(t) &= \sum_a \rho_a(t) \ln\left(\frac{\rho_a(t)}{p_a}\right)\geqslant 0.
\end{align}

It can be shown that $H(t)=0$, if and only if, $\rho_a = p_a$ for all $a$—i.e., $H$ vanishes only when $\rho_a(t)$ is identical with the equilibrium distribution.

The so-called H-theorem states that $H(t)$ is a monotonically decreasing function of time,

\begin{align}
\nonumber
\frac{1}{\lambda^2}\frac{dH(t)}{dt} &= \frac{1}{\lambda^2}\sum_a \frac{d\rho_a(t)}{dt}\ln\left(\frac{\rho_a(t)}{p_a}\right)+\frac{1}{\lambda^2}\sum_a \rho_a(t)\frac{p_a}{\rho_a(t)}\frac{1}{p_a}\frac{d\rho_a(t)}{dt} \\
\nonumber
&= \sum_{a,b}M_{a,b}\rho_b(t)\ln\left(\frac{\rho_a(t)}{p_a}\right) + \frac{1}{\lambda^2}\frac{d}{dt}\sum_{a}\rho_a(t)\\
\nonumber
\nonumber
\nonumber
&= \frac{1}{2}\sum_{a,b}\left(\gamma_{a,b}-\delta_{a,b}\sum_{c}\gamma_{c,a}\right)p_b
\left[\frac{\rho_b(t)}{p_b}\ln\left(\frac{\rho_a(t)}{p_a}\right)+\frac{\rho_a(t)}{p_a}\ln\left(\frac{\rho_b(t)}{p_b}\right)\right]\\
\nonumber
\nonumber
\nonumber
\frac{\rho_b(t)}{p_b}\ln\left(\frac{\rho_a(t)}{p_a}\right)+\frac{\rho_a(t)}{p_a}\ln\left(\frac{\rho_b(t)}{p_b}\right)\right.\\
\nonumber
\left.-\frac{\rho_a(t)}{p_a}\ln\left(\frac{\rho_a(t)}{p_a}\right)-\frac{\rho_b(t)}{p_b}\ln\left(\frac{\rho_b(t)}{p_b}\right)
\right]\\
\left(\ln\left(\frac{\rho_b(t)}{p_b}\right)-\ln\left(\frac{\rho_a(t)}{p_a}\right)\right)\leqslant 0.
\end{align}

This is because $\gamma_{a,b} > 0$ and

\begin{align}
(y-x)(\ln x - \ln y ) = (y-x)\ln\left(\frac{x}{y}\right) \leqslant (y-x)\left(\frac{x}{y}-1\right) = \frac{(y-x)(x-y)}{y} \leqslant 0,
\end{align}

for any $x\neq y$—the consequence of the concavity of $\ln$,

\begin{align}
\ln a - \ln b \leqslant  (a-b)\ln'b = \frac{a-b}b\ \xrightarrow{b=1, a=x/y}\ \ln\left(\frac{x}{y}\right) \leqslant \frac{x}{y} - 1.
\end{align}

The immediate conclusion from H-theorem is that $H(t)$ keeps decreasing until it reaches its $t\to\infty$ limit value $H(\infty) \geqslant 0$. We can determine when the limit is reached by checking at which point the derivative vanishes. There is only one circumstance when $dH/dt=0$: all ratios $\rho_a(t)/p_a$ are mutually equal. This can be achieved only by setting $\rho_a(t) = p_a$ for each $a$, i.e., when $\rho_a(t)$ becomes the equilibrium distribution. In other words, the distribution $\rho_a(t)$ evolves under the action of Davies map until it thermalises, then the evolution stops because $p_a$ is the stationary solution of Pauli master equation, $\sum_{b}M_{a,b}p_b = 0$ (see assignment <ref> of Sec. <ref>).

§.§ Eigensystem analysis

The H-theorem had shown us that the diagonal of the density matrix always ends up at the thermal equilibrium distribution under the action of Davies map. However, we also want to describe the dynamics of this process. The most transparent analysis is based on the eigensystem of generator $\mathcal{M}$.

First, observe that the matrix $M_{a,b}$ is hermitian, but in a specific sense. Consider a scalar product that “embraces” the special symmetry of the Davies generator (<ref>)

\begin{align}
\langle\!\langle \hat A | \hat B \rangle\!\rangle &\equiv \sum_{a,b}e^{\beta\epsilon_a}A_{ab}^*B_{ab}.
\end{align}

For this scalar product we have

\begin{align}
\nonumber
\langle\!\langle \hat A|\mathcal{M}\hat B\rangle\!\rangle &= \sum_{a}e^{\beta\epsilon_a}A_{aa}^*\sum_b M_{a,b} B_{bb}
= \sum_b e^{\beta \epsilon_b}\left(\sum_a \frac{e^{\beta\epsilon_a}}{e^{\beta\epsilon_b}}M_{a,b}A_{aa}\right)^*B_{bb}\\
\nonumber
&= \sum_b e^{\beta \epsilon_b}\left(\sum_a \frac{p_b}{p_a}M_{a,b}A_{aa}\right)^*B_{bb}
=\sum_b e^{\beta\epsilon_b}\left(\sum_a M_{b,a}A_{aa}\right)^*B_{bb}\\
&=\langle\!\langle \mathcal{M}\hat A|\hat B\rangle\!\rangle,
\end{align}

which is the definition of hermicity. From basic algebra we know that hermitian matrices have real eigenvalues and the corresponding eigenvectors (or eigenoperators in the case of matrices of super-operators) form an orthogonal basis in the subspace the matrix operates in, i.e.,

\begin{align}
\mathcal{M}\hat q(\mu) = \mu \hat q(\mu),
\end{align}

where $\mu\in\mathbb{R}$, the eigenoperators live in the subspace $\mathcal{B}_P$, $\hat q(\mu) = \sum_a |a\rangle q_a(\mu) \langle a|$, and they are mutually orthogonal $\langle\!\langle \hat q(\mu)|\hat q(\mu') \rangle\!\rangle \propto \delta_{\mu,\mu'}$, and they form basis,

\begin{align}
\bullet_{\mathcal{B}_P} = \sum_\mu \frac{\hat q(\mu)\langle\!\langle \hat q(\mu)|\bullet\rangle\!\rangle}{\langle\!\langle \hat q(\mu)|\hat q(\mu)\rangle\!\rangle},
\end{align}

where $\bullet_{\mathcal{B}_P}$ is the identity in subspace $\mathcal{B}_P$.

We already know that one of the eigenvalues equals zero and that it corresponds to the equilibrium state, $\mathcal{M}\hat p = 0$ (see assignment <ref> of Sec. <ref>). When we consider the projection onto $\hat p = \hat q(0)$ we get an interesting result,

\begin{align}
\hat p \frac{\langle\!\langle \hat p|\bullet\rangle\!\rangle}{\langle\!\langle\hat p|\hat p\rangle\!\rangle} &=
\hat p \frac{Z^{-1}\sum_a e^{\beta\epsilon_a}e^{-\beta\epsilon_a}\langle a|\bullet|a\rangle}{Z^{-2}\sum_a e^{\beta\epsilon_a}e^{-2\beta\epsilon_a}}
= \hat p \sum_a\langle a|\bullet|a\rangle
=\hat p\operatorname{tr}(\bullet),
\end{align}

which means that all other eigenoperators $\hat q(\mu)$ are traceless because of orthogonality,

\begin{align*}
0 = \langle\!\langle \hat p |\hat q(\mu)\rangle\!\rangle  = \frac{1}{Z}\operatorname{tr}(\hat q(\mu)).
\end{align*}

All non-zero eigenvalues are negative because $\mathcal{M}$ is a non-positive matrix,

\begin{align}
\nonumber
\langle\!\langle \hat A|\mathcal{M}\hat A\rangle\!\rangle &= \sum_a e^{\beta\epsilon_a}A_{aa}^*\left(\sum_b M_{a,b}A_{bb}\right)
= \sum_{a,b}e^{\beta\epsilon_a}\gamma_{a,b}A_{aa}^*A_{bb} - \sum_{a}e^{\beta\epsilon_a}|A_{aa}|^2\sum_e\gamma_{e,a}\\
\nonumber
&=	\sum_{a > b}\left(e^{\beta\epsilon_{a}}\gamma_{a,b}A_{aa}^*A_{bb}+e^{\beta\epsilon_{b}}\gamma_{b,a}A_{bb}^*A_{aa}\right)
-\sum_a e^{\beta\epsilon_a}|A_{aa}|^2\sum_{e\neq a}\gamma_{e,a}\\
\nonumber
-\sum_{c\neq b}e^{\beta\epsilon_c}|A_{cc}|^2\sum_{e\neq a}\gamma_{e,c}
-\sum_{c\neq a}e^{\beta\epsilon_c}|A_{cc}|^2\sum_{e\neq b}\gamma_{e,c}\right]\\
\nonumber
\nonumber
-\sum_{c\neq a,b}e^{\beta\epsilon_{c}}|A_{cc}|^2\sum_{e\neq a,b}\gamma_{e,c}\right]\\
&= -\sum_{a>b}\left(e^{\beta\epsilon_{a}}\gamma_{a,b}|A_{aa}e^{-\beta\epsilon_b}-A_{bb}e^{\beta\epsilon_b}|^2
+\sum_{c\neq a,b}e^{\beta\epsilon_{c}}|A_{cc}|^2\sum_{e\neq a,b}\gamma_{e,c}
\right) \leqslant 0.
\end{align}

With all these facts in hand, we can write the explicit form of the Davies map in the projector subspace,

\begin{align}
\nonumber
e^{\lambda^2t\mathcal{M}} &= \exp\left[\lambda^2t\sum_\mu \mu \frac{\hat q(\mu)\langle\!\langle \hat q(\mu)|\bullet\rangle\!\rangle}{\langle\!\langle \hat q(\mu)|\hat q(\mu)\rangle\!\rangle}\right]
=\sum_{\mu}e^{\lambda^2 \mu t}\frac{\hat q(\mu)\langle\!\langle \hat q(\mu)|\bullet\rangle\!\rangle}{\langle\!\langle \hat q(\mu)|\hat q(\mu)\rangle\!\rangle}\\
\label{eq:etM}
&=\hat p\operatorname{tr}(\bullet) + \sum_{\mu < 0}e^{-\lambda^2|\mu|t}\frac{\hat q(\mu)\langle\!\langle \hat q(\mu)|\bullet\rangle\!\rangle}{\langle\!\langle \hat q(\mu)|\hat q(\mu)\rangle\!\rangle},
\end{align}

and thus

\begin{align}
e^{\lambda^2t\mathcal{M}}\sum_{a}|a\rangle\rho_a(0)\langle a|
= \hat p \operatorname{tr}(\hat\rho_S) + \sum_{\mu<0}e^{-\lambda^2|\mu|t}\hat q(\mu) \frac{\langle\!\langle \hat q(\mu)|\hat\rho_S\rangle\!\rangle}{\langle\!\langle \hat q(\mu)|\hat q(\mu)\rangle\!\rangle}
\xrightarrow{t\to\infty} \hat p.
\end{align}

We knew from the H-theorem that there is only one zero eigenvalue: if there was another $\mu=0$, then the corresponding eigenoperator would also be a stationary solution, and thus, not all initial distributions would tend to $p_a$, which would contradict the theorem.

§.§ Dephasing

For the thermalisation to be complete, bringing the diagonal of the density matrix to the equilibrium distribution is not enough—we also need to have the total dephasing of the off-diagonal elements,

\begin{align}
e^{\lambda^2 t(i\mathcal{J}+\mathcal{G})}\sum_{a\neq b}|a\rangle\rho_{ab}(0)\langle b| \xrightarrow{t\to\infty} 0.
\end{align}

Since the off-diagonal elements,

\begin{align}
\rho_{ab}(t) = \operatorname{tr}(|b\rangle\langle a|e^{\lambda^2t(i\mathcal{J}+\mathcal{G})}\hat\rho_S),
\end{align}

are complex numbers, they cannot be interpreted as probability distribution, and thus, there is no H-theorem for the coherence subspace. Still, the best way to demonstrate that the map in $\mathcal{B}_C$ eliminates everything over time is to investigate the eigenvalues of the generator. However, there are some complications. Even though $\mathcal{G}$ is hermitian,

\begin{align}
\nonumber
\langle\!\langle  \hat A | \mathcal{G} \hat B\rangle\!\rangle &= \sum_{a\neq b}e^{\beta\epsilon_a} A_{ab}^*\left(\sum_{c\neq d}G_{ab,cd}B_{cd}\right)
= \sum_{c\neq d}e^{\beta\epsilon_c}\left(\sum_{a\neq b}\frac{e^{\beta\epsilon_a}}{e^{\beta\epsilon_c}}G_{ab,cd}A_{ab}^*\right)B_{cd}\\
\nonumber
&= \sum_{c\neq d}e^{\beta\epsilon_c}\left(\sum_{a\neq b}\frac{p_c}{p_a}G_{ab,cd}^*A_{ab}\right)^*B_{cd}
=\sum_{c\neq d}e^{\beta\epsilon_c}\left(\sum_{a\neq b}G_{cd,ab}A_{ab}\right)^*B_{cd}\\
& = \langle\!\langle \mathcal{G}\hat A | \hat B\rangle\!\rangle,
\end{align}

and $i\mathcal{J}$ is anti-hermitian

\begin{align}
\langle\!\langle \hat A|i\mathcal{J}\hat B\rangle\!\rangle &= -\langle\!\langle i\mathcal{J}\hat A|\hat B\rangle\!\rangle,
\end{align}

so that both have full set of orthogonal eigenoperators that can be used as a basis, the two generators do not commute,

\begin{align}
[ i\mathcal{J} ,\mathcal{G} ] \neq 0.
\end{align}

Therefore, we need to solve the eigenproblem for the total generator, $i\mathcal{J}+\mathcal{G}$, and, of course, the sum of hermitian and anti-hermitian matrix is neither of those. This means that there is no guarantee that the matrix of the total generator has enough eigenoperators to form a basis, and so, we cannot rely on getting a map with form similar to (<ref>). All is not lost, however; even if matrix does not have enough eigenoperators to form a basis, one can still express the map in terms of eigenvalues $\mu$ and exponential factors $\exp(\lambda^2\mu t)$.

Mathematically, eigenvalues $\mu$ are the roots of the characteristic polynomial,

\begin{align}
\operatorname{det}\left(i\mathcal{J}+\mathcal{G}-\mu\bullet_{\mathcal{B}_C}\right) = 0.
\end{align}

According to the fundamental theorem of algebra, every polynomial of degree $n$ has as many roots; however, it is always possible that some roots (or eigenvalues) have the same value,

\begin{align}
\operatorname{det}\left(i\mathcal{J}+\mathcal{G}-\mu\bullet_{\mathcal{B}_C}\right) &= \prod_k (\mu_k-\mu) = (\mu_1-\mu)^{n_1}(\mu_2-\mu)^{n_2}\cdots .
\end{align}

When $n_k > 1$ we say that the eigenvalue $\mu_k$ is degenerate and $n_k$ is its degree of degeneracy.

The issues come up when the number of linearly independent eigenvectors corresponding to degenerate eigenvalue is smaller than its degree of degeneracy (the two numbers are always equal for hermitian, anti-hermitian, and, in general, normal matrices). If this is the case, then there is not enough (linearly independent) eigenvectors to form a basis and we call the matrix defective.

The simplest example of defective matrix is

\begin{align*}
\bm{A} = \left[\begin{array}{cc} a & 1 \\ 0 & a \\ \end{array}\right].
\end{align*}

Its degenerate eigenvalue is easy to find,

\begin{align*}
0 &= \operatorname{det}\left(\bm{A}-\mu\bm{I}\right) = (a-\mu)^2\ \Rightarrow\ \mu_1=a,\ n_1=2;
\end{align*}

and the corresponding eigenvector,

\begin{align*}
\left[\begin{array}{c} 0 \\ 0 \\ \end{array}\right] &= \left(\bm{A}-a\bm{I}\right)\bm{r}(a)
\left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \\ \end{array}\right]
\left[\begin{array}{c} r_1 \\ r_2 \\\end{array}\right]
\left[\begin{array}{c} r_2 \\ 0 \\\end{array}\right]
\ \Rightarrow\
\bm{r}(a) = \left[\begin{array}{c} 1 \\ 0\\ \end{array}\right].
\end{align*}

And so, there is only one eigenvector—not enough to form a basis.

When the matrix $\bm{B}$ is not defective, one can switch to the basis of its eigenvectors to diagonalize it,

\begin{align*}
\bm{S}^{-1}\bm{B}\bm{S} = \left[\begin{array}{ccccc}
\lambda_1 & & & & \\
& \lambda_2 & & & \\
& & \lambda_2 & & \\
& &  & \lambda_2 & \\
& & & & \ddots\\
\end{array}\right],
\end{align*}

where $\bm{S}$ is a transformation to the basis of eigenvectors, and $\lambda_i$ are the eigenvalues (visible here, $\lambda_2$ is degenerate and $\lambda_1$ is not degenerate). Once diagonalized, it is a trivial matter to calculate any function of the matrix, e.g., the exp,

\begin{align*}
\exp(t\bm{B}) &= \sum_{k=0}^\infty\frac{1}{k!}t^k\bm{B}^k = \sum_{k=0}^\infty\frac{1}{k!}t^k \bm{S}(\bm{S}^{-1}\bm{B}\bm{S})^k\bm{S}^{-1}\\
\left[\begin{array}{ccccc}
e^{\lambda_1t} & & & & \\
& e^{\lambda_2 t} & & & \\
& & e^{\lambda_2 t} & & \\
& &  & e^{\lambda_2 t} & \\
& & & & \ddots\\
\end{array}\right]
\bm{S}^{-1}
\end{align*}

Of course, it is impossible to diagonalize a defective matrix. Nevertheless, there is an algorithm for constructing the missing basis elements—the so-called generalized eigenvectors—in such a way that the matrix in this new basis is closest to the diagonal form as possible. Let $\bm{C}$ be defective matrix with a defective eigenvalue $\mu$ of some degeneration degree $n>1$, then it is possible to find a transformation $\bm{N}$ such that

\begin{align*}
\bm{N}^{-1}\bm{C}\bm{N} &=
\left[\begin{array}{ccccccc}
\ddots & & & & & &\\
& \mu & 1 & & & &\\
& & \mu & 1 & & &\\
& & & \ddots & \ddots & &\\
& & & & \mu & 1 & \\
& & & & & \mu & \\
& & & &  & &\ddots\\
\end{array}\right].
\end{align*}

The $n\times n$ block with $\mu$ on diagonal, $1$'s on super-diagonal and $0$'s everywhere else, is called the Jordan block. Using the basis completed with generalized eigenvectors one can now calculate functions of $\bm{C}$, in particular, the exp,

\begin{align*}
\exp(t\bm{C}) &= \bm{N}\exp(t\bm{N}^{-1}\bm{C}\bm{N})\bm{N}^{-1}\\
\ddots & & & & & &\\
& e^{\mu t} & t e^{\mu t} & \frac{t^2}{2!}e^{\mu t}&  \cdots &\frac{t^n}{n!}e^{\mu t} &\\
& & e^{\mu t} & t e^{\mu t} & \cdots & \frac{t^{n-1}}{(n-1)!}e^{\mu t} &\\
& & & \ddots & \ddots & \vdots &\\
& & & & e^{\mu t} & t e^{\mu t} & \\
& & & & & e^{\mu t} & \\
& & & &  & &\ddots\\
\end{array}\right]\bm{N}^{-1}.
\end{align*}

As we can see, in the case of exp, each non-zero matrix element within the Jordan block scales with $t$ as a polynomial times $\exp(\mu t)$—this is what we were looking for.

Now let us come back to the generator $i\mathcal{J}+\mathcal{G}$. To make sure that the map it generates leads to the total dephasing, we need to show that the real part of its eigenvalues, $\operatorname{Re}\mu$ ($\mu$  are complex, in general), are strictly negative; $\operatorname{Re}\mu = 0$ would not do, because then there would be some blocks that do not decay to zero over time.

As it turns out, we just need to prove that

\begin{align}\label{eq:G_neg}
\forall_{\hat A}\ \langle\!\langle \hat A|\mathcal{G} \hat A\rangle\!\rangle < 0,
\end{align}

i.e., that the dissipative part of the generator is a negative matrix. Let us proceed then,

\begin{align}
\nonumber
\langle\!\langle \hat A|\mathcal{G}\hat A\rangle\!\rangle &=
\sum_{a\neq b}\sum_{c\neq d}\delta(\omega_{ab}-\omega_{cd})2e^{\beta\epsilon_a}A_{ab}^*A_{cd}V_{ac}^*V_{bd}\frac{S(\omega_{ac})}{e^{\beta\omega_{ac}}+1}\\
-\sum_{a\neq b}e^{\beta\epsilon_a}|A_{ab}|^2\sum_{e}
\left(\frac{|V_{ea}|^2S(\omega_{ea})}{e^{\beta\omega_{ea}}+1}+\frac{|V_{eb}|^2S(\omega_{eb})}{e^{\beta\omega_{eb}}+1}\right).
\end{align}

Take two index pairs $(a,b)$ and $(c,d)$ belonging to one energy shell $I(\omega_{ab})$ (i.e., $\omega_{ab}=\omega_{cd}$) and focus on those terms only,

\begin{align}
\nonumber
\label{eq:pt1}&\phantom{=}
= \frac{e^{\beta\epsilon_a}S(\omega_{ac})}{e^{\beta\omega_{ac}}+1}\left(
\right),
\end{align}

where we have used the symmetries $S(\omega_{ac}) = S(-\omega_{ca}) =S(\omega_{ca})$, $V_{ca}^* = V_{ac}$. The same index pairs in the second sum gives us

\begin{align}
\nonumber
&-|A_{ab}|^2e^{\beta\epsilon_a}\sum_{e}\left(\frac{|V_{ea}|^2 S(\omega_{ea})}{e^{\beta\omega_{ea}}+1}+\frac{|V_{eb}|^2 S(\omega_{eb})}{e^{\beta\omega_{eb}}+1}\right)\\
\nonumber
-|A_{cd}|^2e^{\beta\epsilon_c}\sum_{e}\left(\frac{|V_{ec}|^2 S(\omega_{ec})}{e^{\beta\omega_{ec}}+1}+\frac{|V_{ed}|^2 S(\omega_{ed})}{e^{\beta\omega_{ed}}+1}\right)\\
-|A_{ab}|^2 e^{\beta\epsilon_a}R_{ab}-|A_{cd}|^2e^{\beta\epsilon_c}R_{cd}.
\end{align}

Note that $R_{ab}>0$ because we are assuming that $V_{ab}\neq 0$. The “mixed” terms (<ref>) are problematic because their sign is unspecified. Fortunately, we can take care of them by absorbing them into full moduli squares; to do this we need to draw a few terms from $R_{ab}$ and $R_{cd}$ in the second term, namely

\begin{align}
\nonumber
&-|A_{ab}|^2e^{\beta\epsilon_a}\left(\frac{|V_{ca}|^2 S(\omega_{ca})}{e^{\beta\omega_{ca}}+1}+\frac{|V_{db}|^2 S(\omega_{db})}{e^{\beta\omega_{db}}+1}\right)
-|A_{cd}|^2e^{\beta\epsilon_c}\left(\frac{|V_{ac}|^2 S(\omega_{ac})}{e^{\beta\omega_{ac}}+1}+\frac{|V_{bd}|^2 S(\omega_{bd})}{e^{\beta\omega_{bd}}+1}\right)\\
\label{eq:pt2}&\phantom{=}
\left(e^{\beta\omega_{ac}}|A_{ab}|^2+e^{-\beta\omega_{ac}}|A_{cd}|^2\right)\left(|V_{ac}|^2+|V_{bd}|^2\right),
\end{align}

Now we can combine it with the mixed terms to get

\begin{align}
\nonumber
+ |e^{\frac{\beta\omega_{ac}}{2}}A_{ab}V_{bd}^*-e^{-\frac{\beta\omega_{ac}}{2}}A_{cd}V_{ac}^*|^2
\right) \leqslant 0.
\end{align}

Given that $\hat H_S'$ is non-degenerate, there will always be at least two strictly positive terms left in $R_{ab}$ after we draw from it all the terms needed to balance the mixed terms within each energy shell $I(\omega_{ab})$—this means that $\langle\!\langle \hat A|\mathcal{G}\hat A\rangle\!\rangle < 0$ for any non-zero $\hat A$.

Of course, the inequality (<ref>) is also satisfied for all eigenoperators $\hat r(\mu)$ of the total generator, so that we can write,

\begin{align}
(\operatorname{Re}\mu + i \operatorname{Im}\mu)\langle\!\langle\hat r(\mu)|\hat r(\mu)\rangle\!\rangle &=\langle\!\langle \hat r(\mu)|(i\mathcal{J}+\mathcal{G})\hat r(\mu)\rangle\!\rangle
= i\langle\!\langle \hat r(\mu)|\mathcal{J}\hat r(\mu)\rangle\!\rangle + \langle\!\langle \hat r(\mu)|\mathcal{G}\hat r(\mu)\rangle\!\rangle,
\end{align}

then, since $\langle\!\langle \hat A|\hat A\rangle\!\rangle > 0$ and $\langle\!\langle \hat A|\mathcal{J}\hat A\rangle\!\rangle,\langle\!\langle \hat A|\mathcal{G}\hat A\rangle\!\rangle\in\mathbb{R}$ ($\mathcal{J}$ is hermitian because $i\mathcal{J}$ is anti-hermitian) for any non-zero $\hat A$, we get

\begin{align}
\operatorname{Re}\mu &= \frac{\langle\!\langle \hat r(\mu)|\mathcal{G}\hat r(\mu)\rangle\!\rangle}{\langle\!\langle \hat r(\mu)|\hat r(\mu)\rangle\!\rangle} < 0.
\end{align}

Therefore, when $\hat H_S'$ is non-degenerate and all $V_{ab}\neq 0$, we have the total dephasing,

\begin{align}
\nonumber
&\sum_{c\neq d}[e^{\lambda^2t(i\mathcal{J}+\mathcal{G})}]_{ab,cd}\rho_{cd}(0)\\
\nonumber
= \sum_{c\neq d}[\bm{N}
\left[\begin{array}{ccc}
e^{-\lambda^2t|\operatorname{Re}\mu_1|}e^{i\lambda^2t\operatorname{Im}\mu_1}\bm{B}_1(t) & &   \\
& e^{-\lambda^2t|\operatorname{Re}\mu_2|}e^{i\lambda^2t\operatorname{Im}\mu_2}\bm{B}_2(t) &  \\
& &  \ddots \\
\end{array}\right]
\bm{N}^{-1}]_{ab,cd}\rho_{cd}(0)\\
&\phantom{=}\xrightarrow{t\to\infty} 0,
\end{align}

\begin{align}
\bm{B}_i(t) &=\left[\begin{array}{cccc}
1 & t & \cdots & \frac{1}{n_i!}t^{n_i} \\
& 1 & t & \cdots \\
& &  \ddots & \\
& & & 1 \\
\end{array}\right],
\end{align}

if the eigenvalue $\mu_i$ is defective, and $\bm{B}_i(t) = \bm{I}$ when it is not.

§.§ Physics of thermalisation

Here is a bold but surprisingly apt analogy: the process of thermalisation can be compared to the spread of infectious disease. When a “diseasesed” system $E$, that manifests symptoms in a form of $\hat\rho_E\propto\exp(-\beta \hat H_E)$, comes into contact with other system $S$, the “pathogen” is passed on through the coupling $\lambda\hat V\otimes\hat F$. Over time, the infection takes hold and the state of $S$ starts to display symptoms, $\hat\rho_S(t)\to \hat p = \exp(-\beta \hat H_S')/Z$. Then, $S$ becomes a potential vector for a further spread of the infection.

If we are willing to fully embrace this analogy, then we could ask what is this “pathogen” that causes the thermalisation to spread from system to system? For example, the thermal state itself, $\hat\rho_E\propto\exp(-\beta\hat H_E)$, is not it; like we said, the final form of the density matrix is more of a symptom rather than the cause. Indeed, the dynamical map we have been investigating in this chapter does not contain any explicit reference to the initial state of $E$. Therefore, the process of thermalisation is not directly influenced by particular form of the thermal state.

Similarly, the process of thermalisation is largely independent of a detailed form of the coupling operator $\hat V$. However, this does not mean that $\hat V$ has no impact at all. For example, we had to constraint the $S$-side of the coupling with the condition $V_{ab}\neq 0$ to ensure the total dephasing, and thus, to complete the thermalisation. If any $V_{ab} = 0$, then it would open a possibility that $\operatorname{Re}\mu = 0$ for some of the eigenvalues $\mu$ of generator $i\mathcal{J}+\mathcal{G}$, and that would lead to pockets of perpetually oscillating off-diagonal elements of the density matrix. So, within our analogy, the shape of $\hat V$ would be one of the factors affecting how resistant the system is to the transmission of the pathogen. If all energy levels are directly coupled (i.e., $V_{ab}\neq 0$), the pathogen spreads freely and the diseases of thermalisation overcomes the system totally. If some energy levels are unreachable (a corresponding $V_{ab} = 0$), then parts of the system will remain “healthy” and it will display less sever symptoms (i.e., some non-zero off-diagonal elements of density matrix). This is, most likely, an artifact of Born approximation: since we have restricted the map to the lowest order in $\lambda\tau_c$, only the direct transmission $V_{ab}$ is available. Although, at the moment, it is only a speculation on our part, it seems highly probable that going beyond Born would enable indirect transmission paths, e.g., $V_{ac}V_{cb}$ to connect $a$ with $b$ through $c$ when the direct path is blocked, $V_{ab}=0$. If this intuition is correct, then the disease could defeat system's resistance in the end.

Finally, we have the qumulant $\langle\!\langle F(t)F\rangle\!\rangle$, or rather, the real part of its Fourier transform [e.g., see Eq. (<ref>)],

\begin{align*}
\int_{-\infty}^\infty\operatorname{Re}\{\langle\!\langle F(\tau)F\rangle\!\rangle\}e^{-i\omega\tau}d\tau = S(\omega) - i \kappa(\omega).
\end{align*}

As we have already noted, there are no direct references to the thermal state in the Davies map. If so, then how does $S$ know what is the temperature it is supposed to thermalise to? This is where qumulant comes in: in Davies map, it is the only “channel” able to transmit the information about $\beta$ from $E$ to $S$. To be more precise, this information is not passed on by the spectral density $S(\omega)$ or the susceptibility $\kappa(\omega)$ individually (see Sec. <ref> and assignment <ref>), but rather it is encoded in the very specific relation between them—the fluctuation-dissipation theorem (FDT)—which gives us the omnipresent temperature-dependent factor

\begin{align*}
S(\omega) - i \kappa(\omega) = \frac{2S(\omega)}{e^{\beta\omega}+1}.
\end{align*}

The spectral density $S(\omega)$, which, essentially, contains in itself $\hat F$ and $\hat H_E$, does not really influence the process of thermalisation in any qualitative way; it mostly (together with $\hat V$) determines the magnitude of the eigenvalues of generator $\mathcal{L}_D$, which in turn determine the speed of the process. The thermal factor $(e^{\beta\omega}+1)^{-1}$ is the key element here: it is the reason for the symmetries of $\gamma_{a,b}$ and $\mathcal{L}_D$ [see, Eq. (<ref>)]. Precisely because of these symmetries, the thermal state $\hat p\propto \exp(-\beta\hat H_S')$ is the steady-state solution of Davies equation (or, $\mu=0$ eigenoperator of $\mathcal{L}_D$) and we have the H-theorem. Therefore, if we go back to our analogy, the thermal factor can be pointed out as a cause of the disease. Still, we think that the factor is not the pathogen analogue, rather it is only an expression of the actual pathogen—the fluctuation-dissipation theorem. But how a concept like FDT can be analogous to something like virus? Consider this: a typical virus consists of a protein and lipid structure and the genetic material enclosed inside. The purpose of the protein structure is to transport the genetic material and help with depositing it into its target—the host's cell. So, in a sense, the structure is just a part of the transmission mechanism, along with aerosol droplets, bodily fluids, or any other intermediary between the genetic material and the cell's nucleus. The part of the virus that is indispensable for it to be a pathogen is the genetic material (in fact, a primitive kind of viruses called viroids are made of genetic material only). Even the genetic material—a strand of DNA or RNA—can be though of as merely a vessel for the information it carries. It is this information that is the essence of virus as a pathogen: it contains the instructions how to hijack the cell, how to use it to replicate itself, and how to facilitate its further spread. The FDT is like this abstract viral information carried in the genetic material; it is behind the transmission mechanism in the form of the thermal factor, and it takes care of self-replication by leading to thermal state, which is the necessary and sufficient condition for FDT (see assignment <ref>).

We are not done with our analogy yet. We now see that the thermal equilibrium spreads like a disease from infected system to healthy systems (i.e., system not in thermal state). But how this epidemic started? How the “patient zero” become infected? In other words, how the thermal equilibrium gets established when the system is closed and insulated? It is a very interesting and yet unsolved problem. Fortunately, this subject received a lot of attention over the years; if you wish to read more about it, we recommend looking up Eigenstate Thermalization Hypothesis as a starting point.

§ THE ORIGIN OF EXTERNAL FIELDS IN QUANTUM MECHANICS

The origin of external fields

§.§ Where do external fields in quantum mechanics come from?

In classical mechanics, a natural way to simplify a many-body problem is to “replace” some of the elements of the composite system with a surrogate force fields. Take for example a motion of a planet orbiting around its star. Formally, it is a two-body problem and to find the exact motion of the planet one has to solve the coupled equations for the trajectories of both involved parties. However, the formalism of the classical theory allows, through straightforward change of coordinate system, to effectively eliminate the massive star from the picture and, instead, depict the planet as if it moved through space filled with an external (i.e., independent of the planet) field of force. In fact, it is always possible to transform the coupled equations of motion into a coordinate system where the trajectories are decoupled and constrained by external (sometimes time-dependent) fields.

In general, tricks like that do not work in quantum mechanics: due to the product structure of the composite quantum systems' Hilbert spaces and the state-focused description, it is virtually impossible to decouple dynamical equations with a simple change of the coordinate system (although, there are specific contexts where it can be done). Nevertheless, the external fields are routinely used to accurately model the dynamics of quantum systems, presumably in contact with their quantum environments, as systems that are closed but not insulated. Examples range from simple Pauli model of a spin-$1/2$ precssing in an external magnetic field,

\begin{align}\label{eq:pauli_model}
\hat H = \sum_{\alpha=x,y,z}\frac{\mu}{2}\hat\sigma_\alpha B_\alpha,
\end{align}

to stochastic models we have discussed in detail in Chapter <ref>.

Here is the conundrum. Given the general incompatibility of the formalism of quantum mechanics with replacing-subsystems-with-force-fields, how is it possible that the surrogate field models work at all? What is the mechanism behind it? What kind of conditions the composite system must satisfy for such models to be valid? Since the point of surrogate field models is to “eliminate” a coupled subsystem from the description (and replace it with an external field), it seems obvious that the open system theory should be able to give us answers; after all, “eliminating” a subsystem (the environment) and accounting for its influence on what is left is the thing the theory is all about.

§.§ Surrogate field representation

The surrogate field model, where the subsystem $E$ of a composite total system $SE$ is superseded with an external surrogate field $F(t)$, is formally implemented by replacing the full $SE$ Hamiltonian with a $S$-only (stochastic) generator,

\begin{align}
\hat H_{SE} = \hat H_S\otimes\hat 1+\hat 1\otimes\hat H_E + \lambda\hat V\otimes\hat F \quad\to\quad \hat H_F[F](t) = \hat H_S + \lambda F(t)\hat V.
\end{align}

The field $F(t)$ is, in general, a stochastic process—remember that deterministic field [like the magnetic field in the Pauli model (<ref>)] is also a stochastic process but with a trivial probability distribution, $P_F[f] = \delta[f-f_0]$.

The model, or the surrogate field representation [19] as we will call it, is valid when the stochastic dynamics generated by $\hat H_F$ is a good approximation (or an accurate simulation) of the dynamical map originating from the actual $SE$ Hamiltonian, i.e.,

\begin{align}
\left\langle T_\mathcal{W}e^{-i\lambda\int_0^t\mathcal{W}(F(s),\bar{F}(s),s)ds}\right\rangle
\approx \overline{T_\mathcal{W}e^{-i\lambda\int_0^t\mathcal{W}(F(s),F(s),s)ds}}.
\end{align}

In other words, the surrogate field representation is valid when the quasi-average over joint quasiprobability distributions of the coupling operator $\hat F$, $\{Q_F^{(k)}\}_{k=1}^\infty$, can be well approximated with a stochastic average over surrogate field $F(t)$.