Witnessing non-Markovianity by quantum quasi-probability distributions

A topical issue in the quantum theory of open systems [1] is the definition, detection and quantification of non-Markovian dynamics, i.e. dynamical processes determined by memory effects [2–4]. One approach to the definition of non-Markovianity in the quantum regime is based on the concept of an information flow between the open system and its environment [5]. According to this concept Markovian dynamics is characterized by a continuous flow of information from the open system to its environment, while non-Markovian behavior features a flow of information from the environment back to the open system. To formulate these ideas in mathematical terms one has to introduce a suitable measure for the information inside the open system. A natural such measure is given by the trace distance [6] between pairs of open system states, since this distance has a direct interpretation in terms of the distinguishability of the quantum states [7]. Thus, a quantum process is said to be Markovian if the trace distance and, hence, the distinguishablility decreases monotonically in time and, vice versa, it is non-Markovian if the trace distance shows a non-monotonic behavior, implying a temporary increase of the distinguishability which characterizes the information backflow. We note that there are also entropic measures for this information with similar properties [8, 9]. As a further concept one introduces the dynamical maps $\Lambda_t$, $t \geqslant 0$, which describe the evolution in time of the density matrix ρ of the open system as $\boldsymbol{\rho}(t) = \Lambda_t[\boldsymbol{\rho}(0)]$. This map is known to be completely positive and trace preserving (CPTP) and has the important property to be a contraction for the trace distance [6, 10]. This leads to the following measure for the degree of memory effects of a certain pair of initial states $\boldsymbol{\rho}_1$ and $\boldsymbol{\rho}_2$ of the open system [5, 11],

$$\begin{align} \mathcal{N}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, \{\Lambda_t\}) &:= \underset{\sigma \geqslant 0}{\int} \mathrm{d}t \ \sigma(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, t) \end{align} \tag{ 1 }$$

$$\begin{align} \mathrm{with} \quad \sigma(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, t) &:= \frac{\mathrm{d}}{\mathrm{d}t} \mathrm{d}_\mathrm{tr}\big(\Lambda_t[\boldsymbol{\rho}_1], \Lambda_t[\boldsymbol{\rho}_2] \big), \end{align} \tag{ 2 }$$

quantifying the total backflow of information form the environment into the system. Obviously, this quantity depends on the chosen pair of initial states. To make this quantity independent of the initial state pair and to obtain a function which only depends on the family of dynamical maps, one can take the maximum over all initial state pair leading to the non-Markovianity measure

$$\begin{align} \mathcal{N}(\{\Lambda_t \}) \ := \ \underset{\boldsymbol{\rho}_1, \boldsymbol{\rho}_2}{\mathrm{max}}\, \mathcal{N}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, \{\Lambda_t \}). \end{align} \tag{ 3 }$$

This information flow approach to quantum non-Markovianity has been applied to many theoretical models and experimental systems (see, e.g. the reviews [2, 3] and references therein). However, most applications deal with relatively simple open systems, describing for example a qubit or a small number of interacting qubits (e.g. a spin-chain) linearly coupled to bosonic or fermionic reservoirs with a structured spectral density. For more complicated systems not only the determination of the dynamical evolution is much more involved, even if it can be represented in terms of an exact or approximate master equation, but also the calculation of the trace distance between the quantum states can become highly demanding, especially for infinite dimensional continuous variable (CV) systems. In this paper we want to tackle these difficulties by investigating non-Markovianity measures or at least non-Markovianity witnesses which are more suitable for CV systems, employing representations of quantum states by means of quasi-probability distributions.

To this end, we start in section 2 by introducing a vector-like representation of quantum states based on generalized bases on bounded operators consisting of a set of fixed but not mutually orthogonal pure states which we will address as quantum frames. We show that for any quantum state its decomposition into such a quantum frame can be understood as a quasi-probabilistic mixture, and we will recapitulate how one can connect these quantum frames to so called informationally complete positive operator valued measure (IC-POVM) [12, 13], a class of quantum measurements especially useful in quantum tomography and even experientially realizable [14, 15]. The frame based and the IC-POVM based representation of a quantum state can then be used in section 3 to define distance measures for quantum states by applying the Kolmogorov distance for probability distribution. We then will show how to use these measures in order to approximate the trace distance.

While section 2 considers finite dimensional quantum systems only, we will reformulate in section 4 the idea for the continuous quantum frame of coherent states in CV systems and its connection to the Glauber–Sudarshan P-function and the Husimi Q-function [16]. Again we apply the Kolmogorov distance to these quasi-probability distributions in order to approximate the trace distance between given quantum states. We also discuss the performance of this approximation with the help of a class of randomly generated Gaussian states [17]. On the basis of this approximation we will construct in section 5 a suitable witness for the non-Markovianity in CV systems using P- and Q-functions section 5.1 and illustrate this by means of the example of a non-Markovian quantum oscillator section 5.2. In section 5.3 we focus on a special class of dynamical maps which preserve the Gaussianity of an initial quantum state and discuss the time evolution of the Kolmogorov distance between P-functions under those maps. This will lead us to a novel measure for non-Markovianity for the case of Gaussian dynamics, since the Kolmogorov distance turns out to be contracting under Gaussian CPTP maps. Finally we draw our conclusions in section 6 and give an outlook on possible future studies, where one could apply the witnesses and measures of non-Markovianity developed here.

We begin with a brief introduction to frames in general Hilbert spaces in order to clarify the terminology and to sketch basic concepts. For more details, the reader is referred to reference [18]. A frame $\mathcal{F} = \{\vert f_i\rangle\} \subset \mathcal{H}$ in a Hilbert space $\mathcal{H}$ is a set of vectors satisfying $\mathrm{span} \left\lbrace \mathcal{F}\right\rbrace = \mathcal{H}$. The immediate question now is: what is the difference to a basis in $\mathcal{H}$? The answer is the following:

• Firstly, a basis has to be a minimal set spanning $\mathcal{H}$ while a frame might be overcomplete. In this sense a frame is just a generalized notion of a basis.
• Secondly, when speaking of a basis we usually mean an orthonormal basis (ONB) for which $\langle f_i|f_j\rangle = \delta_{ij}$. Although mathematically a basis is neither required to be orthogonal nor normalized, we will use the term frame to stress that orthonormality of its elements is not assumed.

The possible non-orthogonality causes the major technical differences between ONBs and frames. An important object in frame theory is the frame operator $S: \mathcal{H} \to \mathcal{H}$ defined by

$$\begin{align} S \vert v\rangle = \sum_i \langle f_i|v\rangle\vert f_i\rangle. \end{align} \tag{ 4 }$$

In case of an ONB the frame operator is simply the identity operator while for general frames it is not. Later on, this frame operator will provide the link between a frame induced decomposition and a frame induced measurement of a state. Furthermore, we can define a canonical dual frame via the frame operator as

$$\begin{align} \tilde{\mathcal{F}} := \left\{\vert \tilde{f}_i\rangle = S^{-1} \vert f_i\rangle \right\} \qquad \mathrm{with} \quad \vert v\rangle = \sum_i \langle \tilde{f}_i| v\rangle \vert f_i\rangle. \end{align} \tag{ 5 }$$

Again, in the case of an ONB we find $\tilde{\mathcal{F}} = \mathcal{F}$ and, hence, an ONB can be addressed as a minimal and self-dual frame. However, although orthogonality simplifies decomposition and spanning significantly, it might be of advantage to use non-orthogonal bases with certain other useful properties instead as we will do in this paper (see below). Even further, it might be desirable to have a spanning set, i.e. a frame, which covers some part of the vector space at hand more densely such that vectors one wants to decompose are already closer (in a certain sense) to some frame elements and have decompositions more pronounced on a single such element.

Let us now consider a quantum system with finite-dimensional Hilbert space $\mathcal{H}$. We denote the space of bounded operators on $\mathcal{H}$ by $\mathcal{B}(\mathcal{H})$, which becomes a Hilbert spaces on its own via the Hilbert–Schmidt scalar product

$$\begin{align} \langle A,B\rangle_\mathrm{HS} := \mathrm{tr}\left[A^\dagger B\right]. \end{align} \tag{ 6 }$$

The physical states of the quantum system are represented by density matrices ρ which are Hermitian and positive operators with trace one [6]. Pure states of the system are represented by rank-one projections $\boldsymbol{\rho} = |\varphi\rangle\langle\varphi|$, where $\vert \varphi\rangle \in \mathcal{H}$ is a normalized state vector. Any density operator has a spectral decomposition of the form $\boldsymbol{\rho} = \sum_i p_i |i\rangle\langle i|$ with $\vec{p} = ({\ldots}, p_i, {\ldots})^T$ a probability distribution and $\{\vert i\rangle\}$ an ONB in $\mathcal{H}$. Note that two different mixed quantum states not only differ in their probability distribution $\vec{p}$ but typically have different ONBs in their spectral decomposition, too. By contrast, the concept of frames allows the following definition which will give a decomposition into a fixed set of pure states.

Definition 1. A frame $\mathcal{F} : = \left\{|\psi_i\rangle\langle\psi_i| \right\} \subset \mathcal{B}(\mathcal{H})$ consisting of pure quantum states is called a quantum frame.

In the following we will assume $\mathcal{F}$ to be minimal (i.e. the $|\psi_i\rangle\langle\psi_i|$ form a basis in $\mathcal{B}(\mathcal{H})$ which, however, cannot be orthogonal). Accordingly, any density matrix can be decomposed as

$$\begin{align} \boldsymbol{\rho} = \sum_i f_i |\psi_i\rangle\langle\psi_i|, \end{align} \tag{ 7 }$$

where we have

$$\begin{align} \boldsymbol{\rho}^\dagger = \boldsymbol{\rho} &\Leftrightarrow \ f_i \in \mathbb{R} \ \forall i \end{align} \tag{ 8 }$$

$$\begin{align} \mathrm{tr}\left[\boldsymbol{\rho}\right] = 1 \ &\Leftrightarrow \ \sum_i f_i = 1 \end{align} \tag{ 9 }$$

$$\begin{align} \langle \varphi|\boldsymbol{\rho}|\varphi\rangle \geqslant 0 \ &\Leftrightarrow \ \sum_i f_i |\langle \psi_i|\varphi\rangle|^2 \geqslant 0. \end{align} \tag{ 10 }$$

Note that the last line has to be fulfilled for all $\vert \varphi\rangle\in\mathcal{H}$ in order to ensure the positivity of ρ . The coefficients f_i are referred to as frame decomposition coefficients (FDCs). We see that using quantum frames we can decompose any mixed quantum state quasi-probabilistically (i.e. all f_i are real and sum up to one but might be negative) into a fixed set of pure states. Furthermore, $\vec{f} = (\dots, f_i, \dots )^T$ is a vector like representation of the mixed quantum state at hand which we will call the frame vector and indicate by means of

$$\begin{align} \vec{f} \simeq \boldsymbol{\rho} \end{align} \tag{ 11 }$$

the one-to-one correspondence between the frame vector and the density matrix. The frame operator is represented by

$$\begin{align} S_{ij} = \mathrm{tr}\left[|\psi_i\rangle\langle\psi_i| |\psi_j\rangle\langle\psi_j| \right] = |\langle \psi_i|\psi_j\rangle|^2, \end{align} \tag{ 12 }$$

where $\boldsymbol{S} \vec{f}$ is the frame vector of $S(\boldsymbol{\rho})$.

A given quantum frame $\mathcal{F} = \{|\psi_i\rangle\langle\psi_i| \}$ can also be associated to a rank-one or IC-POVM given by $\{\mathcal{E}_i : = \frac{1}{c_i} |\psi_i\rangle\langle\psi_i| \}$ such that $\sum_i \mathcal{E}_i = \unicode{x1D7D9}$ [12, 13]. Such an IC-POVM is a mathematical expression of a tomographic quantum measurement mapping a quantum state ρ to a probability distribution $\vec{p}$ with [6]

$$\begin{align} p_i = \mathrm{tr}\left[\mathcal{E}_i \boldsymbol{\rho}\right] = \sum_j f_j \frac{1}{c_i} \mathrm{tr}\left[|\psi_i\rangle\langle\psi_i| |\psi_j\rangle\langle\psi_j|\right] = \sum_j \frac{1}{c_i} S_{ij} f_j =: \sum_j M_{ij} f_j, \end{align} \tag{ 13 }$$

where we have introduced the matrix M with elements $M_{ij} = \frac{1}{c_i} S_{ij}$. This means that $\mathcal{E}_i$ is the effect operator of outcome i while p_i is the probability to measure this outcome if the system at hand is prepared in state ρ . The term informationally complete refers to the fact that due to the underlying frame structure different quantum states always have different IC-POVM probability distributions as well, i.e. the vector $\vec{p}$ encodes complete information about its quantum state ρ and hence is a vector like representation of it, too,

$$\begin{align} \vec{p} \simeq \boldsymbol{\rho}, \end{align} \tag{ 14 }$$

as used e.g. in [19, 20] (indeed, it is the decomposition into the canonical dual frame of $\mathcal{F}$ which does not consist of pure quantum states anymore), or similarly in [21]. Assuming a minimal quantum frame the matrix M is bijective and we find

$$\begin{align} \boldsymbol{M}\vec{f} = \vec{p} \qquad \mathrm{and} \qquad \vec{f} = \boldsymbol{M}^{-1} \vec{p}. \end{align} \tag{ 15 }$$

One can understand the frame vector $\vec{f}$ as a decomposition based representation of ρ since it tells us how to construct the quantum state at hand using the given frame. Yet the IC-POVM probability vector $\vec{p}$ is a measurement based representation and encodes the results of measuring ρ using the frame-corresponding IC-POVM and the link between both representations is the frame operator S reshaped as the matrix M , as can be seen from equation (13). Again, the difference between both representations is caused by the non-orthogonality of quantum frames.

We also want to point out once more that, although $\vec{f}$ might have single negative entries, $\vec{p}$ will always have non-negative entries since it captures the actual probabilities corresponding to the tomographic measurement with the IC-POVM at hand. Accordingly, $\vec{f}$ can only be guarantied to be quasi-probabilistic while $\vec{p}$ is always properly probabilistic. However, we will use the term quasi-probability distributions further on as the more general class which of course includes proper probability distributions as well.

A prominent and even experimentally realized example (see, e.g. [14, 15]) for quantum frame based IC-POVMs are symmetric-IC-POVMs (or SIC-POVMs) which are associated to a minimal frame with $|\langle \psi_i|\psi_j\rangle|^2 = \frac{1 + n \cdot \delta_{ij}}{1+n}$ and $c_i = n \ \forall i$, where $\dim \mathcal{H} = n$ [13], see figure 1.

Zoom In Zoom Out Reset image size

In section 2 we have introduced a frame based quasi-probabilistic decomposition of density matrices into a fixed set of pure states forming a so-called quantum frame. In this section we use these representations to introduce distance measures for quantum states and compare them with the trace distance which is known as a measure for the distinguishability of quantum states.

The trace distance between two quantum states ρ and σ defined by [6, 7]

$$\begin{align} \mathrm{d}_\mathrm{tr} (\boldsymbol{\rho}, \boldsymbol{\sigma}) := \frac{1}{2} ||\boldsymbol{\rho} - \boldsymbol{\sigma} || \equiv \frac{1}{2}\mathrm{tr} |\boldsymbol{\rho} - \boldsymbol{\sigma}|, \end{align} \tag{ 16 }$$

where $||A||$ denotes the trace norm of an operator A, while the modulus is given by $|A| = \sqrt{A^\dagger A}$ [22]. Defining the selfadjoint operator $\boldsymbol{\Delta} = \boldsymbol{\rho} - \boldsymbol{\sigma}$ and introducing its spectral decomposition $\boldsymbol{\Delta} = \sum_i \lambda_i |i\rangle\langle i|$, we can thus write the trace distance as a sum over the moduli of the eigenvalues λ_i of Δ,

$$\begin{align} \mathrm{d}_\mathrm{tr} (\boldsymbol{\rho}, \boldsymbol{\sigma}) = \frac{1}{2} \sum_i |\lambda_i|. \end{align} \tag{ 17 }$$

The trace distance between two quantum states ρ and σ can be shown to be a measure for the distinguishability of these states [6, 7] and has the following useful properties:

• $0 \leqslant \mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}, \boldsymbol{\sigma}) \leqslant 1$ and $\mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}, \boldsymbol{\sigma}) = 1 \ \Leftrightarrow \ \boldsymbol{\rho} \perp \boldsymbol{\sigma}$
• $\mathrm{d}_\mathrm{tr}(\mathcal{V}[\boldsymbol{\rho}], \mathcal{V}[\boldsymbol{\sigma}]) = \mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}, \boldsymbol{\sigma})$ for any unitary transformation $\mathcal{V}$
• $\mathrm{d}_\mathrm{tr}(\Lambda[\boldsymbol{\rho}], \Lambda[\boldsymbol{\sigma}]) \leqslant \mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}, \boldsymbol{\sigma})$ for any positive trace preserving (PTP) map Λ

We will use the trace distance later in section 5 to quantify the information flow between an open system and its environment and, hence, to quantify non-Markovianity.

According to equation (17) we have to determine the spectrum of $\boldsymbol{\Delta} = \boldsymbol{\rho} - \boldsymbol{\sigma}$ in order to compute the trace distance between two quantum states ρ and σ . However, in many practical applications this can be a demanding problem, in particular in the case of high-dimensional or even infinite-dimensional Hilbert spaces. If, for example, the quantum states at hand represent mixtures of many basis states of a CV system (see below), it can be a hard task to compute a full matrix representation of Δ and evaluate its eigenvalues.

One possible strategy to avoid these difficulties — which we will investigate in this paper — is to use the vector like frame and IC-POVM based representations of the quantum states under consideration. Since both are quasi-probabilistic one can apply distance measures for classical probability distributions to them, e.g. the Kolmogorov distance [6]

$$\begin{align} \mathrm{d}_\mathrm{kol}(\vec{p}, \vec{q}) := \frac{1}{2}\sum_i |p_i - q_i|, \end{align} \tag{ 18 }$$

where $\vec{p}$ and $\vec{q}$ are two probability distributions. Please note that, although originally developed for proper probability distributions only, the Kolmogorov distance applies as a distance measure for the even larger class of quasi-probability distributions, too, since positivity, symmetry and triangular inequality also hold. Though, in contrast to proper probability distributions the Kolmogorov distance between two quasi-probability distributions might take values greater than one.

With these definitions one finds the following chain of inequalities,

$$\begin{align} \mathrm{d}_\mathrm{kol}(\vec{p}_1, \vec{p}_2) \leqslant \mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2) \leqslant \mathrm{d}_\mathrm{kol}(\vec{f}_1, \vec{f}_2), \end{align} \tag{ 19 }$$

which states that the trace distance between any pair of quantum states is bounded from below by the Kolmogorov distance of the corresponding quasi-probability vectors, and from above by the Kolmogorov distance between the corresponding frame vectors. The proof is rather straight forward and given in appendix A. Furthermore, please note that the left inequality can also be deduced from the more general statement [6]

$$\begin{align} \mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2) \ = \ \underset{\left\{\mathcal{E}_i\right\} \in \mathrm{POVMs}}{\mathrm{max}} \ \mathrm{d}_\mathrm{kol}\left(\vec{p}_{\left\{\mathcal{E}_i\right\}}(\boldsymbol{\rho}_1), \vec{p}_{\left\{\mathcal{E}_i\right\}}(\boldsymbol{\rho}_2)\right) \end{align} \tag{ 20 }$$

where the maximum is taken over all possible POVMs (not just IC-POVMs) and $\vec{p}_{\{\mathcal{E}_i\}}(\boldsymbol{\rho})$ means the probability distribution of a POVM $\{\mathcal{E}_i\}$ applied to a quantum state ρ .

Figure 2 illustrates the inequalities (19) for a SIC-POVM and its minimal and symmetric frame $\mathcal{F}_\mathrm{SIC}$ for a qubit. For randomly generated pairs of quantum states the different distance measures are plotted over the trace distance. First we see that in accordance with (19) the Kolmogorov distance between the frame vectors and between the SIC-POVM probability distributions — encoded in their probability vectors $\vec{p}$ — is always above and below the trace distance, respectively. Furthermore, we find that with increasing trace distance the spreading of the Kolmogorov distances between frame vectors and between probability vectors increases as well. This is due to the fact that almost orthogonal state pairs allow for a larger variety of positioning relatively to the convex hull of the quantum frame $\mathcal{F}$ at hand. Either one of the two states is within the convex hull and its orthogonal complement is not, or both states are outside the convex hull. By definition states within the convex hull have frame vectors consisting of entries $f_i \in [0,1]$ only, while states outside the convex hull must have at least one negative entry and the sum of the other entries must be greater than one to compensate this negativity. Accordingly, the absolute values $|f_i|$ of the FDCs of a state inside the convex hull are on average smaller than in case of states outside the convex hull. Thus, if at least one state lies within the convex hull the Kolmogorov distance between the frame vectors will be smaller than in the case when both states are outside the convex hull. Further, this also explains that there is no general functional relation between trace distance and Kolmogorov distance between frame vectors since two pairs of quantum states might have the same trace distance but the Kolmogorov distance between their frame vectors might be different anyway. The same holds for probability vectors, though less pronounced since their Kolmogorov distances have to be smaller in general compared the distances between their corresponding frame vectors and hence, the probability vectors are closer to each other by construction. Interestingly, we find that the average of both Kolmogorov distances, i.e. $\frac{1}{2}\big(\mathrm{d}_\mathrm{kol}(\vec{f}_1, \vec{f}_2) + \mathrm{d}_\mathrm{kol}(\vec{p}_1, \vec{p}_2)\big)$, spreads as well with increasing trace distance but almost symmetrically around the trace distance.

Zoom In Zoom Out Reset image size

We extend the concepts and investigations of section 3 to CV systems, employing coherent states often used, e.g. in quantum optical systems, which actually form a continuous quantum frame [23] as will be explained below. We start with a brief summary of the features which are relevant for our purposes, further details may be found in references [16, 17, 24, 25].

Let us consider a CV system, describing a set of continuous degrees of freedom of a quantum mechanical system, for example a number k of independent harmonic oscillators or bosonic field modes. As usual we denote the vacuum or ground state by $\vert 0\rangle$ and the annihilation and creation operators of the ith mode by a _i and $\boldsymbol{a}_i^\dagger$, respectively. The Weyl displacement operator is given by $\boldsymbol{D}(\vec{\alpha}) = \bigotimes_{i = 1}^k \exp[\alpha_i \boldsymbol{a}_i^\dagger - \alpha_i^* \boldsymbol{a}_i]$ with $\vec{\alpha} \in \mathbb{C}^n$ and $\vec{\alpha}^*$ its complex conjugate. A coherent state of a single mode is defined as [16]

$$\begin{align} \vert \alpha\rangle := \boldsymbol{D}(\alpha)\vert 0\rangle = {\mathbb{e}}^{-|\alpha|^2/2} \sum_{n=0}^\infty \frac{\alpha^n}{\sqrt{n!}} \vert n\rangle \end{align} \tag{ 21 }$$

with $\vert n\rangle$ the number states (Fock basis) of this single mode. For a multi-mode CV system a coherent state is accordingly a tensor product of single mode coherent states

$$\begin{align} \vert \vec{\alpha}\rangle := \vert \alpha_1\rangle\vert \alpha_2\rangle \dots \vert \alpha_k\rangle \end{align} \tag{ 22 }$$

The coherent state are eigenstates of the annihilation operator, $\boldsymbol{a} \vert \alpha\rangle = \alpha \vert \alpha\rangle$. Furthermore, the real and imaginary part of the complex number α describes the expectation values of position q and of momentum p of the coherent state. In fact, choosing appropriate units for position and momentum we have

$$\begin{align} q = \sqrt{2} \ \Re(\alpha) \quad \mathrm{and} \quad p = \sqrt{2} \ \Im(\alpha). \end{align} \tag{ 23 }$$

Thus, we can regard the $\mathbb{C}^k$ as a representation of the phase space of the multi-mode system at hand where the ith complex plane represents the phase space of mode i.

Using (23) we can express the Weyl displacement operator and coherent states in proper phase space coordinates $\vec{x} : = (q_1,p_1, \dots, q_k, p_k)^T$. Let q _i and p _i be position and momentum operator of mode i, respectively, which are connected to creation and annihilation operators via the relations

$$\begin{align} \boldsymbol{q}_i = \sqrt{2}(\boldsymbol{a}_i^\dagger + \boldsymbol{a}_i) \qquad \mathrm{and} \qquad \boldsymbol{p}_i = \sqrt{2}(\boldsymbol{a}_i^\dagger - \boldsymbol{a}_i). \end{align} \tag{ 24 }$$

Using the symplectic form of the phase space

$$\begin{align} \boldsymbol{\Omega} = \bigoplus_{i = 1}^k \boldsymbol{\omega} \qquad \mathrm{with} \qquad \boldsymbol{\omega} = \left(\begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array}\right) \end{align} \tag{ 25 }$$

and $\vec{\boldsymbol{X}} : = (\boldsymbol{q}_1, \boldsymbol{p}_1, \dots, \boldsymbol{q}_k, \boldsymbol{p}_k)^T$ we can finally write [16]

$$\begin{align} \boldsymbol{D}(-\vec{x}) := \boldsymbol{D}\big(\vec{\alpha}(\vec{x})\big) = \exp\left[-i \vec{x}^T \ \boldsymbol{\Omega} \ \vec{\boldsymbol{X}}\right], \end{align} \tag{ 26 }$$

where $\vec{\alpha}(\vec{x})$ simply means the complex vector representing the point $\vec{x}$ in phase space as described by (23). Accordingly, one can equivalently define a single mode i coherent states as

$$\begin{align} \vert \vec{x}_i\rangle := \vert \alpha_i(\vec{x}_i)\rangle = \mathbb{e}^{-|\vec{x}_i|^2/4} \sum_{n=0}^\infty \frac{(q_i + ip_i)^n}{\sqrt{2^n n!}} \vert n\rangle. \end{align} \tag{ 27 }$$

Although coherent states and characteristic functions (see below) are usually introduced using the complex plane representation, we will continue using proper phase space vectors $\vec{x}$ since this vector representation simplifies the notation of Gaussian states (see below). For simplicity we also assume only a single mode further on but the concepts can easily be extended to multi-mode systems.

For mixed quantum states in CV systems there is a useful representation by characteristic functions and quasi-probability distributions [16, 17]. The characteristic function of a mixed quantum state ρ is defined by

$$\begin{align} \chi_{\boldsymbol{\rho}}^s (\vec{x}) = \mathrm{tr}\left[\boldsymbol{\rho} \boldsymbol{D}(-\vec{x})\right] \mathbb{e}^{\frac{s}{4} |\vec{x}|^2}, \end{align} \tag{ 28 }$$

where $s \in [-1, 1]$. Its Fourier transform is the corresponding quasi-probability distribution

$$\begin{align} W_{\boldsymbol{\rho}}^{\,s} (\vec{x}) = \frac{1}{(2\pi)^2} \int_{\mathbb{R}^2} \mathrm{d} \vec{y} \ \chi_{\boldsymbol{\rho}}^s (\vec{y}) \mathbb{e}^{-i \vec{y}^T \boldsymbol{\Omega}^T \vec{x}}. \end{align} \tag{ 29 }$$

Note that these distributions are real-valued, integrate to one but might be locally negative. For s = 0 one finds the so called Wigner function $W_{\boldsymbol{\rho}}^{\,0} = W_{\boldsymbol{\rho}}$ [26], in the case $s = -1$ we obtain the Husimi Q-function $W_{\boldsymbol{\rho}}^{-1} = Q_{\boldsymbol{\rho}}$ [27], while s = 1 yields the Glauber–Sudarshan P-function $W_{\boldsymbol{\rho}}^{\,1} = P_{\boldsymbol{\rho}}$ [28, 29].

As already mentioned, coherent states also form a continuous quantum frame closely related to the P- and Q-functions. It can be shown [16] that

$$\begin{align} \boldsymbol{\rho} = \int_{\mathbb{R}^2} \mathrm{d}^2x \ P_{\boldsymbol{\rho}}(\vec{x}) |\vec{x}\rangle\langle\vec{x}| \end{align} \tag{ 30 }$$

and, invoking the completeness relation $\unicode{x1D7D9} = \frac{1}{\pi} \int_{\mathbb{R}^2} \mathrm{d}^2x |\vec{x}\rangle\langle\vec{x}|$, one also finds

$$\begin{align} Q_{\boldsymbol{\rho}}(\vec{x}) \ = \ \frac{1}{\pi} \langle \vec{x}|\boldsymbol{\rho}|\vec{x}\rangle = \mathrm{tr}\left[\mathcal{E}(\vec{x}) \boldsymbol{\rho}\right], \end{align} \tag{ 31 }$$

where

$$\begin{align} \mathcal{E}(\vec{x}) = \frac{1}{\pi} |\vec{x}\rangle\langle\vec{x}| \end{align} \tag{ 32 }$$

are the elements of the continuous IC-POVM associated with the frame of coherent states. Thus, in the context of coherent states as quantum frame we find that the P-function and the Q-function are analogous to the frame vector $\vec{f}$ and to the IC-POVM probability distribution $\vec{p}$ of section 2, respectively. This analogy also explains why the Q-function is a real probability distribution for any ρ , i.e. takes only non-negative values for any quantum state. Furthermore, we can use (30) and (31) to express the transformation from one distribution to the other as

$$\begin{align} Q(\vec{x}) = \mathcal{M}[P](\vec{x}) = \int_{\mathbb{R}^2} \mathrm{d}^2y \ m(\vec{x}, \vec{y}) P(\vec{y}), \end{align} \tag{ 33 }$$

where

$$\begin{align} m(\vec{x}, \vec{y}) = \frac{1}{\pi} |\langle \vec{x}|\vec{y}\rangle|^2 = \frac{1}{\pi} \mathbb{e}^{-|\vec{x} - \vec{y}|^2/2}. \end{align} \tag{ 34 }$$

The transition kernel $m(\vec{x}, \vec{y})$ represents the continuous version of the matrix M_ij in (13) and again basically consists of the kernel expression of the frame operator $s(\vec{x}, \vec{y}) = |\langle \vec{x}|\vec{y}\rangle|^2 = \mathbb{e}^{-|\vec{x} -\vec{y}|^2/2}$. And again one can derive the same chain of inequalities between the trace distance of two quantum states and the Kolmogorov distance between their Q- and P-functions

$$\begin{align} \mathrm{d}_\mathrm{kol}(Q_1, Q_2) \ \leqslant \ \mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2) \ \leqslant \ \mathrm{d}_\mathrm{kol}(P_1, P_2). \end{align} \tag{ 35 }$$

Here, the Kolmogorov distance of two continuous distributions R₁ and R₂ over phase space is defined by

$$\begin{align} \mathrm{d}_\mathrm{kol}(R_1, R_2) = \frac{1}{2} \int_{\mathbb{R}^2} \mathrm{d}^2x \ |R_1(\vec{x}) - R_2(\vec{x})|. \end{align} \tag{ 36 }$$

Again one can find the proof for both inequalities in appendix B, though the second inequality can also been found in [30].

We also want to stress the point of P- and Q- function both playing two rather independent roles. On the one hand they represent the continuous version of frame vector and IC-POVM probability vector for the continuous quantum frame of coherent states and on the other hand they are also respectively the symmetrically (s = 1) and anti-symmetrically ($s = -1$) ordered quasi-probability distributions over the phase space of the CV system at hand. By contrast the Wigner function is just the normal ordered (s = 0) quasi-probability distribution over phase space and has no meaning so far in context of coherent states as continuous quantum frame. Accordingly, the Wigner function is not considered further on within this article.

A special class of quantum states in CV systems are so called Gaussian states which are usually defined by having a characteristic function for s = 0 that takes a Gaussian form [17]

$$\begin{align} \chi^0_{\boldsymbol{\rho}}(\vec{x}) \ = \ \exp \left[ -\frac{1}{4} \vec{x}^T \boldsymbol{\Omega}^T \boldsymbol{\sigma} \boldsymbol{\Omega} \vec{x} \ + \ i \vec{x}^T \boldsymbol{\Omega}^T \vec{d} \right], \end{align} \tag{ 37 }$$

where σ captures the statistical moments of second order (covariances) and $\vec{d}$ the statistical moments of first order (mean values) of the operators $\boldsymbol{q} = \boldsymbol{X}_1$ and $\boldsymbol{p} = \boldsymbol{X}_2$

$$\begin{align} \sigma_{ij} = \langle \boldsymbol{X}_i \boldsymbol{X}_j + \boldsymbol{X}_j \boldsymbol{X}_i\rangle - 2\langle \boldsymbol{X}_i\rangle\langle \boldsymbol{X}_j\rangle \quad \mathrm{and} \quad d_i = \langle \boldsymbol{X}_i\rangle. \end{align} \tag{ 38 }$$

In case of a Gaussian product state in a multi-mode system the characteristic function can be separated as a product of independent single-mode Gaussian states

$$\begin{align} \chi^0_{\boldsymbol{\rho}}(\vec{x}) \ = \ \prod_{i=1}^k \exp \left[ -\frac{1}{4} \vec{x}_i^T \boldsymbol{\omega}^T \boldsymbol{\sigma}_i \boldsymbol{\omega} \vec{x}_i \ + \ i \vec{x}_i^T \boldsymbol{\Omega}^T \vec{d}_i \right] \end{align} \tag{ 39 }$$

where $\boldsymbol{\sigma}_i$ and $\vec{d}_i$ is the covariance matrix and displacement of the ith mode and the complete covariance matrix and displacement is given as direct sum $\boldsymbol{\sigma} = \bigoplus_{i = 1}^k \boldsymbol{\sigma}_i$ and $\vec{d} = \bigoplus_{i = 1}^k \vec{d}_i$.

With

$$\begin{align} \chi_{\boldsymbol{\rho}}^s(\vec{x}) \ = \ \chi_{\boldsymbol{\rho}}^0(\vec{x}) \ \mathbb{e}^{\frac{s}{4}|\vec{x}|^2} \end{align} \tag{ 40 }$$

we see immediately that for any value of s the characteristic function will be Gaussian as well and since Gaussianity is preserved under Fourier transformation same is true for any quasi-probability distribution. Using $|\vec{x}|^2 = \vec{x}^T \vec{x} = \vec{x}^T \boldsymbol{\Omega}^T \boldsymbol{\Omega} \vec{x}$ one finds

$$\begin{align} \chi_{\boldsymbol{\rho}}^s(\vec{x}) \ &= \ \exp \left[ -\frac{1}{4} \vec{x}^T \boldsymbol{\Omega}^T (\boldsymbol{\sigma} - s\unicode{x1D7D9}) \boldsymbol{\Omega} \vec{x} \ + \ i \vec{x}^T \boldsymbol{\Omega}^T \vec{d} \right] \end{align} \tag{ 41 }$$

$$\begin{align} &= \ \exp \left[ -\frac{1}{4} \vec{x}^T \boldsymbol{\Omega}^T \boldsymbol{\sigma}_s \boldsymbol{\Omega} \vec{x} \ + \ i \vec{x}^T \boldsymbol{\Omega}^T \vec{d} \right] \end{align} \tag{ 42 }$$

with $\boldsymbol{\sigma}_s : = \boldsymbol{\sigma} - s\unicode{x1D7D9}$ and by this

$$\begin{align} W_{\boldsymbol{\rho}}^{\,s}(\vec{x}) \ = \ \frac{1}{\pi \sqrt{\det \boldsymbol{\sigma}_s}} \exp \left[ -\left(\vec{x} - \vec{d}\right)^T \boldsymbol{\sigma}_s^{-1} \left(\vec{x} - \vec{d}\right) \right]. \end{align} \tag{ 43 }$$

Note that for coherent states $|\vec{x}\rangle\langle\vec{x}|$ we have $\boldsymbol{\sigma} = \unicode{x1D7D9}$ and $\vec{d} = \vec{x}$ such that the P function reduces to $P_{|\vec{x}\rangle\langle\vec{x}|}(\vec{y}) = W_{|\vec{x}\rangle\langle\vec{x}|}^{\,s = 1}(\vec{y}) = \delta(\vec{y} - \vec{x})$ as expected. Further on we define $\boldsymbol{\sigma}_P : = \boldsymbol{\sigma}_1$ and $\boldsymbol{\sigma}_Q : = \boldsymbol{\sigma}_{-1}$ as the covariance matrices for P and Q function over phase space. These definitions together with $\boldsymbol{\sigma}_s : = \boldsymbol{\sigma} - s\unicode{x1D7D9}$ lead to

$$\begin{align} \boldsymbol{\sigma}_Q = \boldsymbol{\sigma}_P + 2\cdot\unicode{x1D7D9}, \end{align} \tag{ 44 }$$

i.e. for a given Gaussian quantum state its P- and Q-function are both Gaussian shaped, centered around the same mean value $\vec{d}$ and the Q-function is broadened compared to the P-function by an offset of $2\cdot\unicode{x1D7D9}$ in the covariance matrix.

Figure 3 shows a sample plot for a CV-system similar to figure 2 for the qubit. For simplicity we restrict ourselves to Gaussian quantum states as introduced above with their P-function given as

$$\begin{align} P(\vec{x}) \ = \ \frac{1}{\pi\sqrt{\det \boldsymbol{\sigma}_P}}\exp\left[-\left(\vec{x}-\vec{d}\right)^{\!\!T}\boldsymbol{\sigma}_P^{-1}\left(\vec{x}-\vec{d}\right)\right]. \end{align} \tag{ 45 }$$

Zoom In Zoom Out Reset image size

We plot the Kolmogorov distances $\mathrm{d}_\mathrm{kol}(P_1,P_2)$, $\mathrm{d}_\mathrm{kol}(Q_1,Q_2)$ and their difference $\mathrm{d}_\mathrm{kol}(P_1,P_2) -$ $\mathrm{d}_\mathrm{kol}(Q_1,Q_2)$ over $\mathrm{d}_\mathrm{kol}(P_1,P_2)$ for different scenarios of state pairs. As stated above each Gaussian state is uniquely defined by the displacement of its center $\vec{d}$ and its covariance matrix $\boldsymbol{\sigma}_P$. The displacement and the covariance matrix is expressed as

$$\begin{align} \vec{d} = d \left(\begin{array}{cc} \cos(\varphi_d) \\ \sin(\varphi_d) \end{array}\right) \qquad \mathrm{and} \qquad \boldsymbol{\sigma}_P = S(\varphi) \left(\begin{array}{cc} \sigma_{00} & 0 \\ 0 & \sigma_{11} \end{array}\right) S(\varphi)^T, \end{align} \tag{ 46 }$$

where $S(\varphi)$ is a rotation by ϕ in the phase plane. Thus we can characterize any Gaussian state by its randomly drawn parameters d, ϕ_d , σ_ii and ϕ.

In figure 3 we observe a transition of the differences between the two distances for increasing values of σ_ii . If on the one hand the σ_ii are small the state $\boldsymbol{\rho} = \int_{\mathbb{R}^2}\mathrm{d}^2x P(\vec{x}) |\vec{x}\rangle\langle\vec{x}|$ has a low von Neumann entropy (see below) and we can approximate it by means of the matrix representation $\rho_{mn} = \langle m|\boldsymbol{\rho}|n\rangle$ using only a small number of Fock basis states $\vert m\rangle$. This is the scenario plotted in figure 3(a) from which one can immediately see that the difference between $\mathrm{d}_\mathrm{kol}(P_1,P_2)$ and $\mathrm{d}_\mathrm{kol}(Q_1,Q_2)$ is rather large and, hence, both distances do not yield a good approximation of the trace distance $\mathrm{d}_\mathrm{tr}$. This effect can be explained by the relation (44) between the covariance matrices of the P-function and the Q-function. In case of small diagonal entries in $\boldsymbol{\sigma}_P$ the offset of $2\cdot\unicode{x1D7D9}$ widens the distribution significantly. As a consequence, although two close but very thin P-functions might have a Kolmogorov distance of almost one, their corresponding Q-functions will show significant overlap and hence a much smaller Kolmogorov distance. However, note that in this low entropy regime a direct evaluation of the trace distance by means of a truncated Fock space representation is typically feasible and much more efficient.

On the other hand, if the σ_ii are large this results in a large von Neumann entropy and due to the broad distribution a large number of relevant entries in the Fock-state representation. Figure 3(d) shows this case and we see that $\mathrm{d}_\mathrm{kol}(P_1,P_2)$ and $\mathrm{d}_\mathrm{kol}(Q_1,Q_2)$ are always very close to each other. Again, we can understand this effect by using (44) since now the width $\boldsymbol{\sigma}_P$ is already so large that the offset of $2\cdot\unicode{x1D7D9}$ does not cause a significant change. Accordingly, in this large entropy regime the Kolmogorov distance approximates the trace distance quite well. This is of particular interest since the Fock-state matrix representation will include such a huge number of relevant entries that a direct computation of the trace distance would turn out to be highly challenging. The plots in figures 3(b) and (c) indicate the transition between the low and high entropy regimes, which might also be characterized as quantum and classical like regimes.

The connection of the different scenarios of figure 3 to the von Neumann entropy is illustrated in figure 4, where one finds for the state pairs of figure 3 a plot of their differences of Kolmogorov distances between P- and Q-functions over their average von Neumann entropy $\frac{1}{2}\left(S_\mathrm{v.N.}(\boldsymbol{\rho}_1) + (S_\mathrm{v.N.}(\boldsymbol{\rho}_2) \right)$. The von Neumann entropy of a Gaussian state can be computed using the symplectic eigenvalues of its Wigner functions covariance matrix $\boldsymbol{\sigma}_0 = : \boldsymbol{\sigma}$ [16]. The symplectic eigenvalues ±v of σ are the eigenvalues of the matrix $i \boldsymbol{\Omega} \boldsymbol{\sigma}$ which come as a pair of two real numbers differing only in their sign [16], i.e. we can define the symplectic eigenvalue v as the positive one (in case of more than one mode we find one symplectic eigenvalue for each mode). The von Neumann entropy then reads [16]

$$\begin{align} S_\mathrm{v.N.}(\boldsymbol{\rho}) = -\mathrm{tr}\left\{\boldsymbol{\rho} \log_2 \boldsymbol{\rho} \right\} = \frac{v+1}{2} \log_2\left(\frac{v+1}{2} \right) - \frac{v-1}{2} \log_2\left(\frac{v-1}{2} \right). \end{align} \tag{ 47 }$$

Zoom In Zoom Out Reset image size

Of course the von Neumann entropy of two states in a pair might differ but since both states are generated using parameters from the same intervals, their von Neumann entropies are typically of the same order of magnitude. One finds a clear transition from large scattering in case of low entropies to very small scattering in case of high entropy. The four clusters seen in figure 4 correspond to the four scenarios shown in figure 3.

Finally, we have to make some general remarks concerning the plots of figure 3. Firstly, we see that even the Kolmogorov distance of the P functions — in contrast to the Kolmogorov distance of the frame vectors in the qubit case of figure 2 — is bounded by one. This is simply due to the fact that the plot is restricted to Gaussian states for which even the P-function and the W-function are proper probability distributions. Secondly, we have to note that due to our restriction to Gaussian states we have to be a bit careful when drawing general conclusions about the behavior of Kolmogorov distances between quasi-probability distributions in an entropic transition. Nevertheless, Gaussian states are often understood as prototypical states of CV systems and are widely used in many fields [17, 25, 31–33]. Thus, any knowledge of their behavior in an entropic transition is certainly of value and interest and might hint to more general statements of further investigations.

5.1. Witnessing non-Markovianity with P- and Q-functions

In this section we will use the distance measures introduced above in order to detect and quantify non-Markovianity, i.e. memory effects in the dynamics of open quantum systems described by a family of dynamical maps $\{\Lambda_t\}$. To this end, we employ the notion of quantum non-Markovianity based on the information flow between the open system and its environment. As explained in the introduction, this information flow is quantified by the trace distance between a pair of quantum states which represents a measure for the distinguishability of these states. However, as already mentioned in section 3 the calculation of the trace distance can become very demanding, in particular for CV systems. Though, based on the approximation of the trace distance by the Kolmogorov distances between P- and Q-functions as explained in section 4, we can define an easily accessible witness of non-Markovianity in such CV systems. In fact, if we find

$$\begin{align} \mathrm{d}_\mathrm{kol} \big(Q_1(t), Q_2(t)\big) > \mathrm{d}_\mathrm{kol}\big(P_1(s), P_2(s)\big) \quad {\mbox{for some}} \quad t > s \end{align} \tag{ 48 }$$

we conclude from the inequalities (35) that the trace distance $\mathrm{d}_\mathrm{tr}\big(\boldsymbol{\rho}_1(\tau), \boldsymbol{\rho}_2(\tau)\big)$ increases for some $\tau \in [s,t]$, i.e. non-Markovian effects emerge during this time interval. In a different context a similar strategy has been proposed in [34].

5.2. Example: non-Markovian damped oscillator

We discuss the criterion (48) with the help of the example of a damped harmonic oscillator given by the master equation

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{\rho}(t) \ = -i \omega [\boldsymbol{a}^\dagger\boldsymbol{a}, \boldsymbol{\rho}] + \gamma_-(t)\left( \boldsymbol{a}\boldsymbol{\rho}(t)\boldsymbol{a}^\dagger - \frac{1}{2} \left\{\boldsymbol{a}^\dagger \boldsymbol{a}, \boldsymbol{\rho}(t) \right\} \right) \nonumber\\ + \gamma_+(t)\left( \boldsymbol{a}^\dagger\boldsymbol{\rho}(t)\boldsymbol{a} - \frac{1}{2} \left\{\boldsymbol{a} \boldsymbol{a}^\dagger, \boldsymbol{\rho}(t) \right\} \right) \end{align} \tag{ 49 }$$

with constant frequency ω and time dependent emission and absorption rates $\gamma_-(t)$ and $\gamma_+(t)$, respectively. If $\gamma_\pm(t) \geqslant 0$ for all $t\geqslant 0$ the dynamical map is CP divisible and, hence, describes a Markovian dynamics [3]. We therefore allow temporarily negative rates setting $\gamma_\pm(t) : = \gamma_\pm \sin(\Omega t)$. Transforming the master equation to the interaction picture we thus have

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{\rho}(t) = \gamma_- \sin(\Omega t)\left( \boldsymbol{a}\boldsymbol{\rho}(t)\boldsymbol{a}^\dagger - \frac{1}{2} \left\{\boldsymbol{a}^\dagger \boldsymbol{a}, \boldsymbol{\rho}(t) \right\} \right) \nonumber\\ + \gamma_+ \sin(\Omega t)\left( \boldsymbol{a}^\dagger\boldsymbol{\rho}(t)\boldsymbol{a} - \frac{1}{2} \left\{\boldsymbol{a} \boldsymbol{a}^\dagger, \boldsymbol{\rho}(t) \right\} \right). \end{align} \tag{ 50 }$$

In the following we assume $\gamma_- \geqslant \gamma_+$ and set $\gamma_0 : = \gamma_- - \gamma_+ \geqslant 0$. As is shown in appendix C the master equation (50) can easily be solved for Gaussian initial states by determining the mean values $\langle \boldsymbol{a}\rangle(t)$, $\langle \boldsymbol{a}^2\rangle(t)$ and $\langle \boldsymbol{a}^\dagger \boldsymbol{a}\rangle(t)$ which directly yield the time evolution of $\vec{d}(t)$ and $\boldsymbol{\sigma}_P(t)$:

$$\begin{align} \vec{d}(t) \ &= \ \vec{d}(0) D(t) \end{align} \tag{ 51 }$$

$$\begin{align} \boldsymbol{\sigma}_P(t) \ &= \ \boldsymbol{\sigma}_P(0) D(t)^2 + \frac{2\gamma_+}{\gamma_0}\big(1-D(t)^2\big) \unicode{x1D7D9} , \end{align} \tag{ 52 }$$

where we have introduced

$$\begin{align} D(t) := \mathbb{e}^{-\frac{\gamma_0}{2\Omega}\big( 1 - \cos(\Omega t) \big)}. \end{align} \tag{ 53 }$$

The condition $\gamma_0 \geqslant 0$ ensures that $D(t)\leqslant 1$, and we see that $\gamma_\pm/\Omega$ can be understood as a measure for the size of the coupling (for emission and absorption, respectively).

Figure 5 shows an example of the time evolution of the Kolmogorov distances between P- and Q-functions for the non-Markovian damped oscillator. We can immediately see that the criterion (48) can be applied here to prove non-Markovianity of the dynamics. It can even be used to obtain a lower bound $\mathcal{N}_\mathrm{min}$ for the trace distance based non-Markovianity measure defined by (1). In order to approximate this measure for a time interval $[s,t]$ between a minimum of $\mathrm{d}_\mathrm{kol}(P_1,P_2)$ and a following maximum of $\mathrm{d}_\mathrm{kol}(Q_1,Q_2)$ we use inequalities (48) to get

$$\begin{align} \mathcal{N}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, \left\{\Lambda_\tau| \tau \in [s,t]\right\}) &\geqslant \mathrm{d}_\mathrm{tr}\big( \boldsymbol{\rho}_1(s), \boldsymbol{\rho}_2(s) \big) - \mathrm{d}_\mathrm{tr}\big( \boldsymbol{\rho}_1(t), \boldsymbol{\rho}_2(t) \big) \nonumber \\ &\geqslant \mathrm{d}_\mathrm{kol}\big(Q_1(t),Q_2(t)\big) - \mathrm{d}_\mathrm{kol}\big(P_1(s),P_2(s)\big) \nonumber \\ &=: \mathcal{N}_\mathrm{min}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, \left\{\Lambda_\tau| \tau \in [s,t]\right\}). \end{align} \tag{ 54 }$$

Zoom In Zoom Out Reset image size

For example, in the case shown in figure 5 we can use $[s,t] = [0.5 , 1.0]$ which gives $\mathcal{N}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, \{\Lambda_\tau|$ $\tau \in [0.5,1.0]\}) \ \geqslant \ 0.876$.

Figure 6 shows the lower bound $\mathcal{N}_\mathrm{min}$ of the non-Markovianity as a function of the parameters $\frac{\gamma_+}{\gamma_-}$ and $\frac{\gamma_-}{\Omega}$ for the same initial states as used in figure 5. As one can see for $\frac{\gamma_+}{\gamma_-} \to 1$ as well as for $\frac{\gamma_-}{\Omega} \to 0$ the witnessed non-Markovianity $\mathcal{N}_\mathrm{min}$ approaches zero. This can be understood if we rewrite expression (53) for D(t), which mainly controls the dynamics and the non-Markovianity of the process, as follows

$$\begin{align} D(t) = \exp \left[ -\frac{1}{2} \frac{\gamma_-}{\Omega} \left(1-\frac{\gamma_+}{\gamma_-}\right) \Big(1 - \cos(\Omega t)\Big) \right]. \end{align} \tag{ 55 }$$

Zoom In Zoom Out Reset image size

In figure 7 one can find the evolution of the two Kolmogorov distances for different parameters reflecting four different areas in the plot of figure 6. Subfigures (b) and (c) are both cases in which we can witness almost no non-Markovianity, although for completely different reasons. In the case of (b) the approximation of the trace distance by the two Kolmogorov distances is for the whole time interval rather inaccurate since both Kolmogorov distances differ significantly. Hence, we just cannot witness a lot of non-Markovianity although there might be much more. Yet, in the case of (c) the approximation is sufficiently accurate but the dynamics shows just a small amount of non-Markovianity due to the weak coupling to the environment. The maximal non-Markovianity can be witnessed in subplot (a) reflecting the ideal parameter range for witnessing non-Markovianity, while in subplot (d) the witness works sufficiently well, too, yet the dynamics only leads to a medium amount of non-Markovianity.

Zoom In Zoom Out Reset image size

5.3. Gaussian evolution of P-function Kolmogorov distance

Gaussian states are of special interest e.g. in quantum optics where they are often seen as a typical class of almost classical states (since for all values of s one finds an actually positive function $W_{\boldsymbol{\rho}}^{\,s}$) as well as in CV quantum information since they are easily experimentally realizable and controllable [17, 35]. Accordingly, one often focuses on CPTP dynamical maps which preserve the Gaussian character of the input state and are hence addressed as Gaussian quantum channels [16, 17, 36].

For any CPTP dynamical map Λ on a CV system we can define its action on P-functions as

$$\begin{align} \Lambda[P_{\boldsymbol{\rho}}](\vec{x}) = P_{\Lambda(\boldsymbol{\rho})}(\vec{x}), \end{align} \tag{ 56 }$$

which yields the explicit representation

$$\begin{align} P_{\Lambda(\boldsymbol{\rho})}(\vec{x}) = \int_{\mathbb{R}^2} \mathrm{d}^2y \ \lambda(\vec{x}, \vec{y})P_{\boldsymbol{\rho}}(\vec{y}), \end{align} \tag{ 57 }$$

where

$$\begin{align} \lambda(\vec{x},\vec{y}) = P_{\Lambda(|\vec{y}\rangle\langle\vec{y}|)}(\vec{x}) \end{align} \tag{ 58 }$$

is the P-function of the image under Λ of the coherent state $|\vec{y}\rangle\langle\vec{y}|$. The kernel λ can be understood as continuous analog of the transition matrix mapping an input distribution to an output distribution for a discrete classical stochastic process [37]. This leads us to the following theorem.

Theorem 1. For any pair of quantum states $\boldsymbol{\rho}_1$ and $\boldsymbol{\rho}_2$ on a CV quantum system the Kolmogorov distance between their Glauber–Sudarshan P-functions P₁ and P₂ is contracting under any Gaussian quantum channel.

Proof. Let $P^{^{\prime}} = \Lambda[P]$ be the image of P under a Gaussian quantum channel Λ. We find

$$\begin{align} \mathrm{d}_\mathrm{kol}(P^{^{\prime}}_1, P^{^{\prime}}_2) &= \frac{1}{2} \int_{\mathbb{R}^2} \mathrm{d}^2x \left| \int_{\mathbb{R}^2} \mathrm{d}^2y \ \lambda(\vec{x}, \vec{y}) (P_1(\vec{y}) - P_2(\vec{y})) \right| \end{align} \tag{ 59 }$$

$$\begin{align} ~~~~~~~~~~~~~~~~~~~~~~&\leqslant \ \frac{1}{2} \int_{\mathbb{R}^2} \mathrm{d}^2x \int_{\mathbb{R}^2} \mathrm{d}^2y \ |\lambda(\vec{x}, \vec{y})| \cdot |P_1(\vec{y}) - P_2(\vec{y})| \end{align} \tag{ 60 }$$

$$\begin{align} ~~~~~~~~~~~~~~~~~~&= \ \frac{1}{2} \int_{\mathbb{R}^2} \mathrm{d}^2x \int_{\mathbb{R}^2} \mathrm{d}^2y \ \lambda(\vec{x}, \vec{y}) |P_1(\vec{y}) - P_2(\vec{y})| \end{align} \tag{ 61 }$$

$$\begin{align} ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~&= \ \frac{1}{2} \int_{\mathbb{R}^2} \mathrm{d}^2y \ \underset{=1}{\underbrace{\left(\int_{\mathbb{R}^2} \mathrm{d}^2x \ P_{\Lambda(|\vec{y}\rangle\langle\vec{y}|)}(\vec{x})\right)}} |P_1(\vec{y}) - P_2(\vec{y})| \end{align} \tag{ 62 }$$

$$\begin{align} &= \ \mathrm{d}_\mathrm{kol}(P_1, P_2),~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \end{align} \tag{ 63 }$$

where we have used the fact that coherent states are Gaussian and, hence, $\lambda(\vec{x},\vec{y}) = P_{\Lambda(|\vec{y}\rangle\langle\vec{y}|)}(\vec{x})$ has to be again a Gaussian state for a Gaussian channel Λ, and is therefore a positive function.

This theorem has several interesting consequences. First, it implies that any increase of the Kolmogorov distance between two P-functions is due to a violation of the CPTP property or to a violation of Gaussianity preservation. Let us first suppose we have a family of dynamical maps $\{\Lambda_t\}$ which is known to preserve Gaussianity for all t. We can then use the Kolmogorov distance of the P-functions instead of the trace distance in the definition of non-Markovianity [3, 5] which leads to the novel non-Markovianity measure

$$\begin{align} \mathcal{N}_P(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, \{\Lambda_t\}) := \underset{\delta \geqslant 0}{\int} \mathrm{d}t \ \delta(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, t) \end{align} \tag{ 64 }$$

$$\begin{align} \mathrm{with} \quad \delta(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, t) := \frac{\mathrm{d}}{\mathrm{d}t} \mathrm{d}_\mathrm{kol}\big(\Lambda_t[P_1], \Lambda_t[P_2] \big). \end{align} \tag{ 65 }$$

This quantity measures the backflow of information for a certain pair of quantum states by means of the Kolmogorov distance of the corresponding P-functions P₁ and P₂. By maximizing over all possible pairs of initial states we obtain

$$\begin{align} \mathcal{N}_P(\{\Lambda_t \}) \ := \ \underset{\boldsymbol{\rho}_1, \boldsymbol{\rho}_2}{\mathrm{max}}\ \mathcal{N}_P(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2, \{\Lambda_t \}). \end{align} \tag{ 66 }$$

This is a measure of the non-Markovianity of a given family of Gaussian dynamical maps $\{\Lambda_t \}$ which is analogous to the trace distance based measure defined by equations (1)–(3), but which is significantly easier to use in the case of CV systems. An example is given by the non-Markovian damped oscillator discussed in the previous section 5.2 which is obviously described by a Gaussianity preserving family of dynamical maps. We emphasize again that this measure only requires the dynamical maps $\Lambda_t$ to be Gaussianity preserving, while the initial states are allowed to be arbitrary although we chose them to be Gaussian in figure 5 for reasons of computational simplicity. A further example is the Caldeira–Leggett model of quantum Brownian motion [38, 39], which may be treated analogously to the procedure used in [40].

Finally we suggest another application of theorem 1. If we know for sure that the dynamical map is CP divisible (see, e.g. [3]) any increase of the Kolmogorov distance between the P-functions provides a witness for the non-Gaussianity of the dynamical map. Note that this criterion of the increase of the Kolmogorov distance only represents a sufficient condition for non-Gaussianity since we might think of a CPTP map which maps a coherent state to a non-Gaussian state with nevertheless positive P function. Since the proof relies only on the positivity of $\lambda(\vec{y}, \vec{x})$ and not on its Gaussianity, this would result in a contraction of the Kolmogorov distance as well.

In this paper we used a vector like representations of quantum states by (quasi) probability distributions based on quantum frames and their induced IC-POVMs. The first we called the frame vector $\vec{f}$ and the latter the IC-POVM probability vector $\vec{p}$ of a quantum state ρ . To distinguish two quantum states $\boldsymbol{\rho}_1$ and $\boldsymbol{\rho}_2$ the Kolmogorov distance (as a distance measure for probability distributions) was applied to both $\vec{f}_{1,2}$ and $\vec{p}_{1,2}$ to find an upper bound (by $\mathrm{d}_\mathrm{kol}(\vec{f}_1, \vec{f}_2)$) and a lower bound (by $\mathrm{d}_\mathrm{kol}(\vec{p}_1, \vec{p}_2)$) of the trace distance between $\boldsymbol{\rho}_1$ and $\boldsymbol{\rho}_2$.

Having discussed these bounds for the case of a SIC-POVM and its symmetric quantum frame on a qubit, we applied this idea to the continuous quantum frame given by the coherent states in a CV system. We saw that the frame vector $\vec{f}$ has to be replaced by the Glauber–Sudarshan P-function and the probability vector $\vec{p}$ by the Husimi Q-function — both two widely used quasi-probability distributions for quantum states. As for the discrete case an inequality chain holds for these continuous quasi-probability distributions, providing again upper and lower bounds for the trace distance. We illustrated the performance of these bounds for a class of randomly generated Gaussian states and demonstrated that the bounds become steadily tighter for high entropy states. Thus, for large-entropy quantum states one can approximate successfully their trace distance by the Kolmogorov distances between their P- and Q-functions. This is of particular interest since the numerical computation of the trace distance of quantum states becomes highly demanding in the regime of large entropies, while on the other hand an efficient determination of the trace distance is relatively easy for low entropy states in a suitable basis of the state space. Thus, the method proposed here is in a certain sense complementary to the direct determination of the trace distance by diagonalization in low-dimensional Hilbert spaces.

Building on these result we constructed a witness for the non-Markovianity of the dynamics of open quantum systems which is based on the Kolmogorov distances between P- and Q-functions and applied it to the model of a non-Markovian oscillator. The witness was also shown to lead to a lower bound for the non-Markovianity measure based on the trace distance which is again particularly efficient in the regime of high entropy states. Furthermore, we showed that the Kolmogorov distance between P-functions has to contract under any Gaussian CPTP map (Gaussian quantum channel). Based on this we derived a proper non-Markovianity measure for Gaussian quantum dynamical maps, i.e. maps that preserve the Gaussianity of quantum states. Finally, we suggested another application of the contraction property, namely to use it as a witness for non-Gaussianity in case one knows for sure the dynamical map to be Markovian.

For future work we suggest to apply the non-Markovianity measure for Gaussian dynamics to more complex systems, e.g. the Caldeira–Leggett model of quantum Brownian motion, or to the dynamics of open systems described by more realistic non-Markovian quantum master equations. Beyond the regime of Gaussian dynamics one can of course still use our non-Markovianity witness, which might be a possible approach to investigate non-Markovian interacting many-body systems. To this end, it will be helpful to find a more detailed characterization of the scenarios in which this witness works especially efficient and to investigate if the observed quality of the approximation for high entropy quantum states actually holds in more general situations.

All data that support the findings of this study are included within the article (and any supplementary files).

We thank Jonathan Brugger, Christoph Dittel, Andreas Buchleitner and Tanja Schilling for fruitful discussions. This work has been supported by the German Research Foundation (DFG) through FOR 5099.

We start with the first inequality which states that the Kolmogorov distance between two IC-POVM probability vectors is smaller or equal to the trace distance between the quantum states at hand. Introducing the spectral decomposition of $\boldsymbol{\Delta} : = \boldsymbol{\rho}_1 - \boldsymbol{\rho}_2 = \sum_a \Delta_a |a\rangle\langle a|$ we obtain

$$\begin{align} \mathrm{d}_\mathrm{kol.}(\vec{p}_1, \vec{p}_2) &= \frac{1}{2} \sum_i \left|p_i^1 - p_i^2\right| = \frac{1}{2} \sum_i \left|\mathrm{tr}\left[\mathcal{E}_i \boldsymbol{\rho}_1\right] - \mathrm{tr}\left[\mathcal{E}_i \boldsymbol{\rho}_2\right]\right| \end{align} \tag{ A.1 }$$

$$\begin{align} ~~~~~~~~~~~~~~~~~~~~~~~~~&= \frac{1}{2} \sum_i \left|\mathrm{tr}\left[\mathcal{E}_i \boldsymbol{\Delta}\right]\right| = \frac{1}{2} \sum_i \left|\mathrm{tr}\left[\sum_a \Delta_a \mathcal{E}_i |a\rangle\langle a|\right]\right| \end{align} \tag{ A.2 }$$

$$\begin{align} ~~~~~~~~~~~~~~~~~~~~~~~~~&\leqslant \ \frac{1}{2} \sum_{i, a} |\Delta_a \underset{\geqslant 0}{\underbrace{\mathrm{tr}\left[\mathcal{E}_i |a\rangle\langle a|\right]}}| = \frac{1}{2} \sum_{i, a} |\Delta_a| \ \mathrm{tr}\left[\mathcal{E}_i |a\rangle\langle a|\right] \end{align} \tag{ A.3 }$$

$$\begin{align} ~~~~~~~~~~~~~~&= \frac{1}{2} \sum_a |\Delta_a| \ \underset{=\mathrm{tr}\left[\unicode{x1D7D9} |a\rangle\langle a|\right] = 1}{\underbrace{\mathrm{tr}\left[\sum_i \mathcal{E}_i |a\rangle\langle a|\right]}} = \frac{1}{2} \underset{=\mathrm{tr} |\boldsymbol{\Delta}|}{\underbrace{\sum_a |\Delta_a|}} \end{align} \tag{ A.4 }$$

$$\begin{align} &= \frac{1}{2} \mathrm{tr}|\boldsymbol{\rho}_1 - \boldsymbol{\rho}_2| = \mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2).~~~~~~ \end{align} \tag{ A.5 }$$

Note that we actually do not need $\{\mathcal{E}_i\}$ to be an IC-POVM associated to a quantum frame. In fact, the first inequality is actually true for any POVM.

For the second inequality stating that the trace distance between two quantum states is smaller or equal to the Kolmogorov distance between their quantum frame vectors we use $\boldsymbol{\rho}_1 - \boldsymbol{\rho}_2 = \sum_i (f_i^{\,1} - f_i^{\;2}) |\psi_i\rangle\langle\psi_i|$ to get

$$\begin{align} \mathrm{d}_\mathrm{tr} (\boldsymbol{\rho}_1, \boldsymbol{\rho}_2) = \frac{1}{2} \|\boldsymbol{\rho}_1 - \boldsymbol{\rho}_2 \| = \frac{1}{2} \left\| \sum_i (f_i^{\,1} - f_i^{\;2}) |\psi_i\rangle\langle\psi_i| \right\| \end{align} \tag{ A.6 }$$

$$\begin{align} \leqslant \frac{1}{2} \sum_i || (f_i^{\,1} - f_i^{\;2}) |\psi_i\rangle\langle\psi_i| || \end{align} \tag{ A.7 }$$

$$\begin{align} = \frac{1}{2} \sum_i |f_i^{\,1} - f_i^{\;2}| \cdot |||\psi_i\rangle\langle\psi_i||| \end{align} \tag{ A.8 }$$

$$\begin{align} = \frac{1}{2} \sum_i |f_i^{\,1} - f_i^{\;2}| = \mathrm{d}_\mathrm{kol}(\vec{f}_1, \vec{f}_2), \end{align} \tag{ A.9 }$$

where we have used that the trace norm satisfies the triangular inequality and that the trace norm of a rank-one projection is equal to one.                            $\quad \square$

The proof runs similar to the one of appendix A. Based on the completeness relation [16]

$$\begin{align} \frac{1}{\pi} \int \mathrm{d}^2\alpha \ |\alpha\rangle\langle\alpha| = \unicode{x1D7D9} \end{align} \tag{ B.1 }$$

we define for simplicity $\mathcal{E}(\alpha) = \frac{1}{\pi} |\alpha\rangle\langle\alpha|$, i.e. we have $Q_{\boldsymbol{\rho}}(\alpha) = \mathrm{tr}\left[\mathcal{E}(\alpha) \boldsymbol{\rho}\right]$, as well as again $\boldsymbol{\Delta} : = \boldsymbol{\rho}_1 - \boldsymbol{\rho}_2 = \sum_i \Delta_i |i\rangle\langle i|$. The first inequality is then obtained as follows,

$$\begin{align} &\mathrm{d}_\mathrm{kol}(Q_1, Q_2) = \frac{1}{2} \int \mathrm{d}^2\alpha \ |Q_1(\alpha) - Q_2(\alpha)| = \frac{1}{2} \int \mathrm{d}^2\alpha \ |\mathrm{tr}\left[\mathcal{E}(\alpha)\boldsymbol{\Delta}\right]| \end{align} \tag{ B.2 }$$

$$\begin{align} & = \frac{1}{2} \int \mathrm{d}^2\alpha \ \left|\sum_i \Delta_i \underset{\geqslant 0}{\underbrace{\mathrm{tr}\left[\mathcal{E}(\alpha) |i\rangle\langle i|\right]}}\right| \end{align} \tag{ B.3 }$$

$$\begin{align} & \leqslant \frac{1}{2} \int \mathrm{d}^2\alpha \ \sum_i |\Delta_i| \mathrm{tr}\left[\mathcal{E}(\alpha) |i\rangle\langle i|\right] \end{align} \tag{ B.4 }$$

$$\begin{align} ~~~~~~~~~~~~~~~~~~~~~~~~~~~&= \frac{1}{2} \sum_i |\Delta_i| \mathrm{tr}\left[\int \mathrm{d}^2\alpha \mathcal{E}(\alpha) |i\rangle\langle i|\right] = \frac{1}{2} \sum_i |\Delta_i| \underset{=1}{\underbrace{\mathrm{tr}\left[\unicode{x1D7D9} |i\rangle\langle i|\right]}} \end{align} \tag{ B.5 }$$

$$\begin{align} &= \mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2).~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \end{align} \tag{ B.6 }$$

To prove the second inequality we use $\boldsymbol{\rho}_1 - \boldsymbol{\rho}_2 = \int \mathrm{d}^2\alpha \ \big(P_1(\alpha) - P_2(\alpha)\big) |\alpha\rangle\langle\alpha|$ which yields

$$\begin{align} &\mathrm{d}_\mathrm{tr}(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2) = \frac{1}{2} \|\boldsymbol{\rho}_1 - \boldsymbol{\rho}_2\| = \frac{1}{2} \left\| \int \mathrm{d}^2\alpha \ \big(P_1(\alpha) - P_2(\alpha)\big) |\alpha\rangle\langle\alpha| \right\| \end{align} \tag{ B.7 }$$

$$\begin{align} ~~~~~~~~~~~~~~~~~~&\leqslant \frac{1}{2} \int \mathrm{d}^2 \alpha \ |P_1(\alpha) - P_2(\alpha)| \cdot \| |\alpha\rangle\langle\alpha| \| = \mathrm{d}_\mathrm{kol}(P_1, P_2), \end{align} \tag{ B.8 }$$

where we have again used properties of the trace norm.                    $\quad \square$

The master equation is given by

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{\rho}(t) \ = \ \gamma_- \sin(\Omega t)\big( \boldsymbol{a}\boldsymbol{\rho}(t)\boldsymbol{a}^\dagger - \frac{1}{2} \left\{\boldsymbol{a}^\dagger \boldsymbol{a}, \boldsymbol{\rho}(t) \right\} \big) & \nonumber\\ + \ \gamma_+ \sin(\Omega t)\big( \boldsymbol{a}^\dagger\boldsymbol{\rho}(t)\boldsymbol{a} - \frac{1}{2} \left\{\boldsymbol{a} \boldsymbol{a}^\dagger, \boldsymbol{\rho}(t) \right\} \big), \end{align} \tag{ C.1 }$$

and for any time independent operator $\hat{o}$ we have

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \langle \hat{o}\rangle = \frac{\mathrm{d}}{\mathrm{d}t}\mathrm{tr}\left[\hat{o}\boldsymbol{\rho}\right] = \mathrm{tr}\left[\hat{o}\frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{\rho}\right]. \end{align} \tag{ C.2 }$$

Furthermore, we know that for Gaussian states the Wigner function $W(q,p)$ over the phase space with is defined by its displacement vector and covariance matrix

$$\begin{align} d_i = \langle \boldsymbol{X}_i\rangle \quad \mathrm{and} \quad \sigma_{ij} = \langle \boldsymbol{X}_i \boldsymbol{X}_j + \boldsymbol{X}_j \boldsymbol{X}_i\rangle - 2\langle \boldsymbol{X}_i\rangle\langle \boldsymbol{X}_j\rangle \end{align} \tag{ C.3 }$$

with $\vec{\boldsymbol{X}} = (\boldsymbol{q}, \boldsymbol{p})^T = \frac{1}{\sqrt{2}}\left(\boldsymbol{a} + \boldsymbol{a}^\dagger, -i (\boldsymbol{a} - \boldsymbol{a}^\dagger)\right)^T$. Accordingly we find

$$\begin{align} d_0 = \frac{1}{\sqrt{2}} \left(\langle \boldsymbol{a}\rangle + \langle \boldsymbol{a}^\dagger\rangle\right) \quad \mathrm{and} \quad d_1 = \frac{-i}{\sqrt{2}} \left(\langle \boldsymbol{a}\rangle - \langle \boldsymbol{a}^\dagger\rangle\right) \end{align} \tag{ C.4 }$$

as well as

$$\begin{align} \sigma_{00} &= \mathrm{var}\left(\boldsymbol{a}\right) + \mathrm{var}\left(\boldsymbol{a}^\dagger\right) + 2\left( \langle \boldsymbol{a}^\dagger \boldsymbol{a}\rangle - \langle \boldsymbol{a}^\dagger\rangle\langle \boldsymbol{a}\rangle \right) + 1 \end{align} \tag{ C.5 }$$

$$\begin{align} \sigma_{11} &= -\mathrm{var}\left(\boldsymbol{a}\right) - \mathrm{var}\left(\boldsymbol{a}^\dagger\right) + 2\left( \langle \boldsymbol{a}^\dagger \boldsymbol{a}\rangle - \langle \boldsymbol{a}^\dagger\rangle\langle \boldsymbol{a}\rangle \right) + 1 \end{align} \tag{ C.6 }$$

$$\begin{align} \sigma_{01} &= \sigma_{10} = -i\left( \mathrm{var}\left(\boldsymbol{a}\right) - \mathrm{var}\left(\boldsymbol{a}^\dagger\right) \right) \end{align} \tag{ C.7 }$$

with $\mathrm{var}\left(\hat{o}\right) = \langle \hat{o}^2\rangle - \langle \hat{o}\rangle^2$. Using (C.2) we now find

$$\begin{align} & \frac{\mathrm{d}}{\mathrm{d}t} \langle \boldsymbol{a}\rangle(t) = \gamma_-(t)\mathrm{tr}\left[\boldsymbol{a}\boldsymbol{a}\boldsymbol{\rho}\boldsymbol{a}^\dagger - \frac{1}{2}(\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{\rho} + \boldsymbol{a}\boldsymbol{\rho}\boldsymbol{a}\boldsymbol{a}^\dagger)\right] + \gamma_+(t)\mathrm{tr}\left[\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{\rho}\boldsymbol{a} - \frac{1}{2}(\boldsymbol{a}\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{\rho} + \boldsymbol{a}\boldsymbol{\rho}\boldsymbol{a}\boldsymbol{a}^\dagger)\right] \end{align} \tag{ C.8 }$$

$$\begin{align} & = \ \gamma_-(t)\mathrm{tr}\left[\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{a}\boldsymbol{\rho} - \frac{1}{2}\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{\rho} - \frac{1}{2}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{a}\boldsymbol{\rho}\right] + \ \gamma_+(t)\mathrm{tr}\left[\boldsymbol{a}\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{\rho} - \frac{1}{2}\boldsymbol{a}\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{\rho} - \frac{1}{2}\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{\rho}\right] \end{align} \tag{ C.9 }$$

$$\begin{align} &= \ \frac{\gamma_-(t)}{2} \mathrm{tr}\left[-[\boldsymbol{a},\boldsymbol{a}^\dagger]\boldsymbol{a}\boldsymbol{\rho}\right] + \frac{\gamma_+(t)}{2}\mathrm{tr}\left[\boldsymbol{a}[\boldsymbol{a},\boldsymbol{a}^\dagger]\boldsymbol{\rho}\right] \end{align} \tag{ C.10 }$$

$$\begin{align} &= \ -\frac{\gamma_-(t)-\gamma_+(t)}{2}\langle \boldsymbol{a}\rangle(t) \end{align} \tag{ C.11 }$$

with $\gamma_\pm(t) = \gamma_\pm \sin(\Omega t)$ and in the last step we have used $[\boldsymbol{a},\boldsymbol{a}^\dagger] = \unicode{x1D7D9}$. Similar calculations also lead to

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \langle \boldsymbol{a}^\dagger\rangle(t) &= -\frac{\gamma_-(t)-\gamma_+(t)}{2}\langle \boldsymbol{a}^\dagger\rangle \end{align} \tag{ C.12 }$$

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \langle \boldsymbol{a}^2\rangle(t) &= -(\gamma_-(t)-\gamma_+(t))\langle \boldsymbol{a}^2\rangle \end{align} \tag{ C.13 }$$

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \langle (\boldsymbol{a}^\dagger)^2\rangle(t) &= -(\gamma_-(t)-\gamma_+(t))\langle (\boldsymbol{a}^\dagger)^2\rangle . \end{align} \tag{ C.14 }$$

To obtain $\langle \boldsymbol{a}^\dagger \boldsymbol{a}\rangle$ we proceed as follows,

$$\begin{align} & \frac{\mathrm{d}}{\mathrm{d}t}\langle \boldsymbol{a}^\dagger \boldsymbol{a}\rangle = \gamma_-(t) \mathrm{tr}\left[\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{a}\boldsymbol{\rho}\boldsymbol{a} - \frac{1}{2}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{\rho} - \frac{1}{2}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{\rho}\boldsymbol{a}^\dagger\boldsymbol{a}\right] + \gamma_+(t) \mathrm{tr}\left[\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{\rho}\boldsymbol{a} - \frac{1}{2}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{\rho} - \frac{1}{2}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{\rho}\boldsymbol{a}\boldsymbol{a}^\dagger\right] \end{align} \tag{ C.15 }$$

$$\begin{align} &= \ \gamma_-(t) \mathrm{tr}\left[\boldsymbol{a}^\dagger \ [\boldsymbol{a}^\dagger,\boldsymbol{a}] \ \boldsymbol{a}\boldsymbol{\rho}\right] \ + \ \gamma_+(t) \mathrm{tr}\big[\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{a}^\dagger \boldsymbol{\rho} - \frac{1}{2}(\boldsymbol{a}\boldsymbol{a}^\dagger + [\boldsymbol{a}^\dagger,\boldsymbol{a}])\boldsymbol{a}\boldsymbol{a}^\dagger\boldsymbol{\rho} - \frac{1}{2}\boldsymbol{a}\boldsymbol{a}^\dagger(\boldsymbol{a}\boldsymbol{a}^\dagger + [\boldsymbol{a}^\dagger,\boldsymbol{a}])\boldsymbol{\rho}\big]\nonumber\\ \end{align} \tag{ C.16 }$$

$$\begin{align} & = \ -\gamma_-(t)\mathrm{tr}\left[\boldsymbol{a}^\dagger\boldsymbol{a} \boldsymbol{\rho}\right] + \gamma_+(t)\mathrm{tr}\left[\boldsymbol{a}^\dagger\boldsymbol{a}\boldsymbol{\rho} + [\boldsymbol{a},\boldsymbol{a}^\dagger]\boldsymbol{\rho}\right] \end{align} \tag{ C.17 }$$

$$\begin{align} &= \ -\gamma_-(t)\langle \boldsymbol{a}^\dagger\boldsymbol{a}\rangle + \gamma_+(t)\big(\langle \boldsymbol{a}^\dagger\boldsymbol{a}\rangle + 1\big). \end{align} \tag{ C.18 }$$

With $\gamma_0 = \gamma_- - \gamma_+$ and $D(t) = \mathbb{e}^{\frac{\gamma_0}{\Omega}\big(\cos(\Omega t)-1\big)}$ these differential equations get solved by

$$\begin{align} \langle \boldsymbol{a}\rangle(t) &= \langle \boldsymbol{a}\rangle_0 D(t) \end{align} \tag{ C.19 }$$

$$\begin{align} \langle \boldsymbol{a}^\dagger\rangle(t) &= \langle \boldsymbol{a}^\dagger\rangle_0 D(t) \end{align} \tag{ C.20 }$$

$$\begin{align} \langle \boldsymbol{a}^2\rangle(t) &= \langle \boldsymbol{a}^2\rangle_0 D(t)^2 \end{align} \tag{ C.21 }$$

$$\begin{align} \langle (\boldsymbol{a}^\dagger)^2\rangle(t) &= \langle (\boldsymbol{a}^\dagger)^2\rangle_0 D(t)^2 \end{align} \tag{ C.22 }$$

$$\begin{align} \langle \boldsymbol{a}^\dagger\boldsymbol{a}\rangle(t) &= \left(\langle \boldsymbol{a}^\dagger\boldsymbol{a}\rangle_0 - \frac{\gamma_+}{\gamma_0}\right) D(t)^2 + \frac{\gamma_+}{\gamma_0}, \end{align} \tag{ C.23 }$$

which finally yields

$$\begin{align} \vec{d}(t) &= \vec{d}(0) D(t) \end{align} \tag{ C.24 }$$

$$\begin{align} \boldsymbol{\sigma}(t) &= \boldsymbol{\sigma}(0) D(t)^2 + \frac{2\gamma_+}{\gamma_0}\big(1-D(t)^2\big)\unicode{x1D7D9} + \unicode{x1D7D9}. \end{align} \tag{ C.25 }$$

Since $\boldsymbol{\sigma} = \boldsymbol{\sigma}_0$ is the covariance matrix of the Wigner function and $\boldsymbol{\sigma}_P = \boldsymbol{\sigma}_{1} = \boldsymbol{\sigma}_{0} - \unicode{x1D7D9}$ one finally gets the covariance matrix of the P-function

$$\begin{align} \boldsymbol{\sigma}_P(t) = \boldsymbol{\sigma}_P(0) D(t)^2 + \frac{2\gamma_+}{\gamma_0}\big(1-D(t)^2\big)\unicode{x1D7D9}, \end{align} \tag{ C.26 }$$

while the displacement is the same for Wigner function and P-function.