Solutions of generic bilinear master equations for a quantum oscillator -- positive and factorized conditions on stationary states

TL;DR

This paper derives solutions for a broad class of bilinear master equations governing quantum oscillators, characterizing stationary states' positivity and factorization, and identifying conditions for positive evolution even with non-completely positive generators.

Contribution

It provides a unified framework for solving bilinear master equations, characterizes stationary states via Minkowski vectors, and explores positivity conditions beyond complete positivity.

Findings

Solutions in Gaussian form for generic bilinear master equations.

Characterization of stationary states using Minkowski vectors.

Conditions for positive evolution with non-completely positive generators.

Abstract

We obtain the solutions of the generic bilinear master equation for a quantum oscillator with constant coefficients in the Gaussian form. The well-behavedness and positive semidefiniteness of the stationary states could be characterized by a three-dimensional Minkowski vector. By requiring the stationary states to satisfy a factorized condition, we obtain a generic class of master equations that includes the well-known ones and their generalizations, some of which are completely positive. A further subset of the master equations with the Gibbs states as stationary states is also obtained. For master equations with not completely positive generators, an analysis on the stationary states suggests conditions on the coefficients of the master equations that generate positive evolution for a given initial state.

Tables

Table 1. Table 1: Different classes of master equations with stationary states that satisfy the factorized condition ( 78 ). In all the classes, γ > 0 𝛾 0 \gamma>0 , θ 0 = 2 ​ ω 0 subscript 𝜃 0 2 subscript 𝜔 0 \theta_{0}=2\omega_{0} and η 0 < 0 subscript 𝜂 0 0 \eta_{0}<0 . (IIB) and (IIC) reduce to the KL equation at θ 1 = 0 subscript 𝜃 1 0 \theta_{1}=0 .

Class	${\eta }_2$	${\eta }_1$	${\theta }_2$	$\left|{\theta }_1\right|$	Completely positive
HPZ (IA)	$\ne 0$	${\eta }_0$	$-\gamma$		no
Conjugate HPZ (IB)	$\ne 0$	$-{\eta }_0$	$\gamma$		no
CL (IIAa)	0	${\eta }_0$	$-\gamma$		no
Conjugate CL (IIAb)	0	$-{\eta }_0$	$\gamma$		no
Generalized CL (IIAc)	0	$-{\eta }_0\frac{{\theta }_2}{\gamma }$	$\left|{\theta }_2\right|\le \gamma$	0	$\frac{\gamma }{\sqrt{1-{\theta }_2^2/{\gamma }^2}}\le \left|{\eta }_0\right|$
Generalized KL 1 (IIB)	0	${\eta }_0\frac{{\theta }_1}{{\theta }_0}$	0		yes
Generalized KL 2 (IIC)	0	0	$\gamma \frac{{\theta }_1}{{\theta }_0}$		yes

Equations

$$\frac{\partial \rho }{\partial t}=-K\rho ,$$

$$K$$

$$K_0$$

$$K_1$$

$$i{L}_0=\frac{i}{2}\left(-\frac{{\partial }^2}{\partial Q\partial r}+Qr\right),i{M}_1=\frac{i}{2}\left(\frac{{\partial }^2}{\partial Q\partial r}+Qr\right),i{M}_2=-\frac{1}{2}\left(Q\frac{\partial }{\partial Q}+r\frac{\partial }{\partial r}+1\right),O_0=-\frac{1}{2}\left(Q\frac{\partial }{\partial Q}-r\frac{\partial }{\partial r}\right),$$

$$O_{+}=\frac{1}{4}\left(\frac{{\partial }^2}{\partial Q^2}-r^2\right),L_{1+}=-\frac{1}{4}\left(\frac{{\partial }^2}{\partial Q^2}+r^2\right),L_{2+}=-\frac{i}{2}r\frac{\partial }{\partial Q},$$

$$Q$$

$$\rho (Q,r;t)$$

$$\frac{d\mu }{dt}$$

$$\frac{d\kappa }{dt}$$

$$\frac{d(\mu +\nu )}{dt}$$

$${\int }_{-\infty }^{\infty }\rho (Q,0;t)dQ=1,$$

$$\rho (Q,r=0;t)=\sqrt{\frac{2\mu (t)}{\pi }}{e}^{-2\mu (t)Q^2}$$

$$\mu (t)>0,\text{and}\nu (t)\ge 0,$$

$$\displaystyle\langle\hat{O}\rangle_{t}=\text{tr}\big{[}\hat{O}\hat{\rho}(t)\big{]}=\int_{-\infty}^{\infty}O(Q,r)\rho(Q,r;t)\big{|}_{r=0}dQ\,,$$

$$x(Q,r)=Q+\frac{r}{2},p(Q,r)=-i\frac{\partial }{\partial x}=-i\left(\frac{1}{2}\frac{\partial }{\partial Q}+\frac{\partial }{\partial r}\right),$$

$$⟨{\hat{x}}^2{⟩}_t$$

$$\underline{X}$$

$$\displaystyle\underline{\beta}\equiv\left(\begin{array}[]{cc}0&\underline{I}\\ -\underline{I}&0\end{array}\right),\qquad\underline{I}=\left(\begin{array}[]{cc}1&0\\ 0&1\end{array}\right).$$

$$S(J)\equiv \exp (\theta J),$$

$$X_i^{\prime }$$

$$\underline{S}(J)\equiv \exp (\theta \underline{J}),$$

$$[J,X_i]=-\sum _j{J}_{ij}{X}_j.$$

$$\displaystyle i\underline{L}_{0}=\frac{i}{2}\left(\begin{array}[]{cc}0&\underline{\sigma}_{1}\\ \underline{\sigma}_{1}&0\end{array}\right),\quad i\underline{M}_{1}=\frac{i}{2}\left(\begin{array}[]{cc}0&\underline{\sigma}_{1}\\ -\underline{\sigma}_{1}&0\end{array}\right),\quad i\underline{M}_{2}=\frac{1}{2}\left(\begin{array}[]{cc}-\underline{I}&0\\ 0&\underline{I}\end{array}\right),\quad\underline{O}_{0}=\frac{1}{2}\left(\begin{array}[]{cc}-\underline{\sigma}_{3}&0\\ 0&\underline{\sigma}_{3}\end{array}\right),$$

$$\displaystyle\underline{O}_{+}=-\frac{1}{2}\left(\begin{array}[]{cc}0&\underline{\sigma}_{d}\\ \underline{\sigma}_{u}&0\end{array}\right),\quad\underline{L}_{1+}=\frac{1}{2}\left(\begin{array}[]{cc}0&-\underline{\sigma}_{d}\\ \underline{\sigma}_{u}&0\end{array}\right),\quad\underline{L}_{2+}=\frac{i}{2}\left(\begin{array}[]{cc}-\underline{\sigma}_{-}&0\\ 0&\underline{\sigma}_{+}\end{array}\right),$$

$$\displaystyle\underline{\sigma}_{1}=\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right),\quad\underline{\sigma}_{+}=\left(\begin{array}[]{cc}0&1\\ 0&0\end{array}\right),\quad\underline{\sigma}_{-}=\left(\begin{array}[]{cc}0&0\\ 1&0\end{array}\right),\quad\underline{\sigma}_{u}=\left(\begin{array}[]{cc}1&0\\ 0&0\end{array}\right),\quad\underline{\sigma}_{d}=\left(\begin{array}[]{cc}0&0\\ 0&1\end{array}\right).$$

$${\underline{S}}^{\text{T}}(J)\underline{\beta }\underline{S}(J)=\underline{\beta },$$

$$e^{-t\underline{K}}$$

$$\displaystyle\qquad\quad-2\cosh(\omega t/2)\sinh(\gamma t/2)\underline{O}_{0}-\Sigma_{0}\underline{O}_{+}-e^{\gamma t/2}\Sigma_{1}\underline{L}_{1+}-e^{\gamma t/2}\Sigma_{2}\underline{L}_{2+}\bigg{]}\,,$$

$$\Sigma$$

$$\displaystyle\mbox{\boldmath$\eta$}=\eta_{0}\mbox{\boldmath$e$}_{0}+\eta_{1}\mbox{\boldmath$e$}_{1}+\eta_{2}\mbox{\boldmath$e$}_{2}\,,$$

$$\displaystyle\mbox{\boldmath$e$}_{0}\cdot\mbox{\boldmath$e$}_{0}$$

$$\displaystyle\mbox{\boldmath$e$}_{0}\wedge\mbox{\boldmath$e$}_{1}$$

$$\displaystyle\mbox{\boldmath$a$}\cdot(\mbox{\boldmath$b$}\wedge\mbox{\boldmath$c$})$$

Full text

††journal: Physica A

1 Introduction

When we obtain the master equations of a system in contact with an environment, there are often assumptions made to facilitate the derivation of the effective influence of the environment on the system [1, 2]. For instance, we assume that the initial density matrix of the system and environment is factorized, the coupling between the system and its environment is weak, the rotating wave approximation is valid, the memory effects are negligible, or the generators of the master equation should be completely positive, and etc. Appropriate derivations should produce master equations with well-behaved solutions.

However, anomaly could arise if we carelessly apply the master equations to situations that are not consistent with the assumptions made in their derivations, such as studying low temperature behavior of a system with the Caldeira-Leggett equation [3, 4, 5], obtaining an Markovian equation by ignoring the memory effects or using the rotating wave approximation [6, 7, 8, 9], or starting with initial conditions inconsistent with the factorized assumptions [10, 11].

An issue that has received wide attention is the completely positive [12, 13] nature of the generators of the master equations. Initial states that are positive could evolve outside their positive domain during some interval of the time evolution when the master equations do not have completely positive generators. However, there is no a priori reason why master equations must be completely positive [10]. Not completely positive generators with appropriately chosen coefficients could maintain the positivity of the states for a given initial state for all time [14]. The systems are then properly behaved as far as positivity is concerned.

It is in this spirit that we start with a generic master equation that satisfies the essential hermitian and trace preserving requirements of a reduced dynamics [15, 16, 17]. We obtain its explicit solutions in Gaussian form and analyzed the stationary states. We further identify a factorized condition on the stationary states that is strong enough to produce a generic class of master equations. This class includes the well-known master equations as special cases, i.e., the Kossakowski-Lindblad (KL) equation for quantum optical systems, the Caldeira-Leggett equation (CL) and the Hu-Paz-Zhang (HPZ) equation for quantum Brownian motion, as well as their generalizations, some of which are found to be completely positive. In this respect, Ref. [18] discussed the role of the Gibbs states, which are examples of states satisfying the factorized condition, and other requirements in deciding the form of quantum master equations. By analyzing the positivity of the stationary states of the master equation, we could infer the coefficients of the master equation that yield positive evolution for a given initial state.

Solutions to various master equations had been obtained, such as for a family of equations whose forms are closely related to the Fokker-Planck equations [19, 20], the KL equation [21], the HPZ equation [22], the generic equations [23] that include the quantum Brownian motion, as well as solutions in terms of the second moments [24]. The uncertainty principle [25, 18] and the Schwartz inequality [16] were invoked to give constraints on the second moments of the observables of the equations. Numerical solutions could be obtained using the Gaussian ansatz [26].

Our method is closely related to Ref. [23], where a different basis of generators was used to derive the general solutions to master equations for quantum oscillators with time-dependent coefficients. Generic conditions on the positivity of the equations were also given. However, owing to the complexity of the exact time-dependent solutions, the effects of the coefficients on the behaviours of the solutions are not transparent, and it is difficult to analyze them. Therefore, an analysis on the solutions to the generic master equation with constant coefficients could provide insights on the structure of the master equations and their roles in deciding the positivity of the time evolution.

We first provide the most general form of the master equations in Section 2, and discuss the necessary and sufficient conditions for positive semidefinite density matrices in the Gaussian form. We then obtain the four-dimensional (4D) matrix representation of the time evolution operator in Section 3, and decompose it into a Baker-Campbell-Hausdorff (BCH) formula in Section 4. The solutions to the generic master equation are then obtained in Section 5. This is followed by a study on the properties of the generic stationary states in Section 6. In Section 7, we make use of the factorized condition to obtain the known master equations and their generalizations. This is followed by discussions in Section 8. Some of the details of the calculations are presented in the appendices.

2 The generic master equation

2.1 Generator of time evolution

The density matrices of quantum oscillators evolve according to the equation

[TABLE]

where we separate the generator of time evolution $K$ into two parts

[TABLE]

in which $\gamma$, $\theta_{i},\eta_{i},i=0,1,2$ are real coefficients. $K$ is the most general form of the generators that preserves the hermiticity of density matrices and conserves the probability of reduced dynamics [17]. We consider only constant coefficients here. In the position coordinates, the operators are

[TABLE]

where the center and relative coordinates are defined as

[TABLE]

respectively, in which $\widetilde{x}$ and $x$ are the position coordinates in the bra- $|\widetilde{x}\rangle$ and ket-space $\langle x|$, respectively. The commutation relations between these generators can be found in Ref. [17].

Let us now discuss the effect of each generator. $iL_{0}$ is the free Liouvillean of the harmonic oscillator. Its coefficient is related to the natural frequency of the oscillator by $\theta_{0}=2\omega_{0}$. The operator $iM_{1}$ renormalizes or shifts the natural frequency to $\omega_{\text{renorm}}\equiv(\omega_{0}^{2}-\theta_{1}^{2}/4)^{1/2}$. This renormalization occurs, for example, when the oscillator is coupled to a field, which in the open quantum system context eventually reduces to the environment.

The $iM_{2}$ introduces damping to the oscillator. From the solution obtained later in Section 5.3, we find that $iM_{2}$ modifies the frequency of the damped oscillator to $\omega_{\text{damp}}\equiv(\omega_{\text{renorm}}^{2}-\theta_{2}^{2}/4)^{1/2}$. The relative magnitude of the coefficient of $iM_{2}$ decides whether the oscillator is underdamped, critically damped, or overdamped, when $\omega_{\text{renorm}}$ is greater than, equal to, or smaller than $\theta_{2}/2$, respectively. In the overdamped situation, $\omega_{\text{damp}}$ turns imaginary. The $O_{0}-1/2$ term contributes to the relaxation and excitation of the oscillator.

The coefficients of $O_{+}$ and $L_{1+}$ is a function of the temperature of the environment. The $\partial^{2}/\partial Q^{2}$ term is analogous to the diffusion term in the classical diffusion equation. It originates from the random motion of the environment degrees of freedom. It tends to spread the probability distribution function along the $Q$ or diagonal direction. On the other hand, the $r^{2}$ term is responsible for destroying the off-diagonal components of the density matrices, or decoherence [27].

The generator $L_{2+}$ is a diffusion term that has no classical counterpart in one space dimension. We will find that it leads to negative probability when its effects are dominant over the $O_{+}$ and $L_{1+}$ terms.

2.2 Gaussian ansatz

We consider solutions in the Gaussian form,

[TABLE]

where $\mu,\kappa$ and $\nu$ are real functions of time [26]. Density matrices of this form have zero first moments. Inserting Eq. (5) into the equation of motion (1), we find that the coefficients evolve according to the following set of equations,

[TABLE]

Eqs. (6a)-(6c) are coupled nonlinear equations of $\mu,\kappa$ and $\nu$. If we add kernels linear in the position coordinates, i.e. $\exp[-\xi(t)Q-\chi(t)ir]$, to $\rho$, the evolutions of $\mu,\kappa,\nu$ (6a)-(6c) are unaffected, though the first moments are now nonzero. The kernels linear in the coordinates can be transformed away by displacement operators [28, 2]. Hence, $\mu,\kappa$ and $\nu$ alone are sufficient to determine the behaviours of the dissipative dynamics and without loss of generality we will not consider kernels linear in the coordinates here.

We normalize the density matrices (5) according to

[TABLE]

where

[TABLE]

is a probability distribution function. The necessary and sufficient conditions for the positive semidefiniteness of $\rho$ is [29]

[TABLE]

for all time.

Our objective is to find the solutions to the set of equations of motion (6a)-(6c) that yield positive semidefinite $\rho$.

2.3 Second moments

Before we present the solutions to the equations, we divert to study the relationships between the second moments with the coefficients $\mu,\kappa,\nu$, and the positive semidefiniteness conditions (9). The expectation value of an operator $\hat{O}$ is defined as

[TABLE]

where the subscript $t$ denotes its time-dependence. In the $Q,r$ coordinates, the position and the momentum operators are respectively given by

[TABLE]

where we use the units $\hbar=1$.

For density matrices of the Gaussian form (5), the first moment of the position and momentum operators vanish, $\langle\hat{x}\rangle_{t}=0=\langle\hat{p}\rangle_{t}$. The second moments are

[TABLE]

Substituting Eqs. (12) to the Schwartz inequality $\langle\hat{x}^{2}\rangle_{t}\langle\hat{p}^{2}\rangle_{t}\geq\langle\hat{x}\hat{p}\rangle_{t}\langle\hat{p}\hat{x}\rangle_{t}$, or the generalized uncertainty relation [25], \big{[}\langle\hat{p}^{2}\rangle_{t}-\langle\hat{p}\rangle_{t}^{2}\big{]}\big{[}\langle\hat{x}^{2}\rangle_{t}-\langle\hat{x}\rangle_{t}^{2}\big{]}-\big{[}\frac{1}{2}(\langle\hat{x}\hat{p}\rangle_{t}+\langle\hat{p}\hat{x}\rangle_{t})-\langle\hat{p}\rangle_{t}\langle\hat{x}\rangle_{t}\big{]}^{2}\geq 1/4, where we use the units $\hbar=1$, the condition $\nu(t)\geq 0$ is reproduced. Together with $\mu>0$, we recover the necessary and sufficient conditions for the positive semidefiniteness of the time evolution (9). We arrive at the same conclusion when we apply the second moments of the creation and annihilation operators to the Schwartz inequality \langle a^{\dagger}a\rangle_{t}\big{(}\langle a^{\dagger}a\rangle_{t}+1\big{)}\geq\langle{a^{\dagger}}^{2}\rangle_{t}\langle a^{2}\rangle_{t} [16].

3 Four-dimensional (4D) matrix representation of time evolution operator

We will obtain the solutions to Eqs. (6a)-(6c) indirectly through group theory method. As a preparation, we first obtain the four-dimensional (4D) matrix representation of the time evolution operator. The matrix representation of the generators (3a)-(3b) can be obtained as follows. We introduce the column matrix

[TABLE]

where underlined symbol denotes matrix. Its components are labeled by $X_{i}$, $i=1,2,3,4$. We find that $[X_{i},X_{j}]=\beta_{ij}$, where $\beta_{ij}$ are the components of the skew-symmetric matrix

[TABLE]

Under a similarity transformation by the operator

[TABLE]

in which $J$ denotes the generic generators in Eqs. (3a)-(3b), $X_{i}$ transforms as [29]

[TABLE]

$S_{ij}$ are components of the matrix

[TABLE]

whereas $\underline{J}$ is the 4D matrix representation of the generator $J$. For infinitesimal transformations, we expand both sides of Eq. (24) to first order in $\theta$ using Eqs. (23) and (25) to obtain

[TABLE]

Applying Eq. (26) to each of the generators (3a)-(3b), we can extract the components $J_{ij}$ of the matrix representation of the generators. They are

[TABLE]

where $\underline{\sigma}_{i}$ denote the matrices

[TABLE]

The generators (27i)-(27p) satisfy $\underline{\beta}\underline{J}=(\underline{\beta}\underline{J})^{\text{T}}$, so that the commutation relation $[X^{\prime}_{i},X^{\prime}_{j}]=\beta_{ij}$ is preserved. Moreover, they obey the quadratic condition,

[TABLE]

which means that they belong to the sympletic group in 4D [30].

Using the matrix representation of the generators, we show in A that the time evolution operator can be cast into the following form,

[TABLE]

Notice that the prefactor $\exp(\gamma t/2)$ in Eq. (30) arises from the $-\gamma/2$ term in $K_{1}$ (2c).

In Eq. (30b), we have introduced three-dimensional Minkowski space vectors with metric signature $(-,+,+)$, labeled by boldfaced letters such as

[TABLE]

and similarly for $\theta$ and $\Sigma$. The basis vectors have the scalar products [31]

[TABLE]

and zero otherwise. The cross products are defined in close analogy to the commutation relations of the su(1,1) algebra

[TABLE]

and zero otherwise. We note that the cross products differ by an overall sign compared to those introduced in Ref. [31]. Some useful identities of the vectors are

[TABLE]

4 Baker-Campbell-Hausdorff (BCH) formula of time evolution operator

In the next step, we propose the following Baker-Campbell-Hausdorff (BCH) formula for the time evolution operator,

[TABLE]

where it is more convenient to use a different combination of generators for $K_{0}$, namely,

[TABLE]

because in the 4D matrix representation, we find that $\underline{M}_{\pm}^{2}=\underline{0}$.

Using the Wei-Norman method [32], we find that the coefficients satisfy the following equations,

[TABLE]

Eqs. (37a)-(37e) still hold if the coefficients in the original generator, $\gamma,\theta_{i},\eta_{i},$ have time-dependence. Since we consider only constant coefficients, the generic solution is simpler and can be worked out explicitly.

Instead of solving the set of differential equations directly, we solve for the coefficients by group theory method. We work out the coefficients in the BCH formula in terms of the coefficients $\gamma,\theta_{i},\eta_{i}$. After writing out the right hand side (RHS) of the BCH formula (35) in the matrix representation, it is equated to the corresponding matrix representation on the RHS of Eq. (30) to yield the desired results.

Noticing that $\underline{J}^{2}=\underline{0}$ for the generators except $\underline{O}_{0}$ and $\underline{M}_{0}$, the BCH formula (35) takes the form

[TABLE]

We can then carry out the matrix multiplication straight-forwardly. As shown in B, after setting the coefficients of the matrices equal to those of Eq. (30), we obtain the following results,

[TABLE]

with the initial conditions $h(0)=0,m_{0}(0)=1,m_{\pm}(0)=0,$ and \mbox{\boldmathg}(0)=\mbox{\boldmath0}. We have defined a unit vector $\hat{\theta}$ with components

[TABLE]

and

[TABLE]

We note that $\omega=4i\omega_{\text{damp}}$. As the coefficients vary, $\Gamma$ probes the Minkowski space spanned by the three linearly independent vectors, \mbox{\boldmath\eta},\mbox{\boldmath\theta}, and \mbox{\boldmath\theta}\wedge\mbox{\boldmath\eta}. \mbox{\boldmath\Pi}_{\pm} projects $\eta$ into a light-like vector,

[TABLE]

which lies on a plane Lorentz orthogonal to $\hat{\theta}$. The projectors have the following properties,

[TABLE]

We can verify straight-forwardly that Eqs. (39a)-(39d) are indeed the solutions to Eqs. (37a)-(37e).

5 Solutions of generic master equation

5.1 Matrix representation of the density matrices

Let us start with an initially normalized density matrix,

[TABLE]

It evolves into

[TABLE]

where we have extracted the factor $\exp(\gamma t/2)$ from $\exp(-Kt)$ by defining

[TABLE]

The matrix representation of $\exp(-K^{\prime}t)$ can be inferred from the first equality of Eq. (4). On the other hand, the matrix representation of $\rho^{\prime}(0)$ can be obtained by using another set of generators

[TABLE]

whose matrix representations are given in C. We find that

[TABLE]

where

[TABLE]

We can then multiply the matrices to obtain

[TABLE]

where $\underline{D}_{ij}$ are matrices that are functions of $\mu_{0},\nu_{0},\kappa_{0}$ and the coefficients from the BCH formula (4). The $\underline{D}_{ij}$ can be worked out explicitly.

Since $e^{-\underline{K}^{\prime}t}\underline{\rho}^{\prime}(0)$ belongs to the complex symplectic group [33], it satisfies the quadratic condition (29) [30, 29], which gives rise to the following matrix relations [29],

[TABLE]

where the superscript T denotes matrix transposition.

5.2 Baker-Campbell-Hausdorff (BCH) formula of density matrix

As the final step before we arrive at the solution, we propose another form of BCH formula for $\exp(-K^{\prime}t)\rho^{\prime}(0)$,

[TABLE]

The matrix representation of the RHS of Eq. (59) are respectively given by

[TABLE]

After carrying out the matrix multiplications on the RHS of Eqs. (59) using Eqs. (60e)-(60t), we equate the resulting matrix to Eq. (57). This procedure yields four matrix equations,

[TABLE]

The unknown matrices on the RHS can be solved in terms of the $\underline{D}_{ij}$ matrices,

[TABLE]

which give the ten coefficients in the BCH formula (59). Eq. (61a) reproduces Eq. (58c), so it does not provide new information. Since the identity (59) is obtained in the defining 4D matrix representation, it is also valid in all representations [30], including the infinite dimensional representation we start with.

Lastly, we operate Eq. (59) on constant so that the exponentials containing differential operators on the RHS of Eq. (59) reduce to identity. Substituting the results back into Eq. (46), we are left with a simple expression,

[TABLE]

where $\mu(t),\kappa(t)$ and $\nu(t)$ can be worked out directly from Eq. (62a) by matrix multiplication, whereas $\exp[u(t)+v(t)]$ can be solved explicitly from Eq. (62c) using Eqs. (60j) and (60o) to yield $\exp[u(t)+v(t)]=1/\text{det}\underline{D}_{22}$. By inspecting the first equation in Eq. (64), Eqs. (65a) and (66a), we then realize that the factor $\mu_{0}e^{\gamma t}e^{u(t)+v(t)}$ under the square root is none other than $\mu(t)$. In this way we recover $\rho(Q,r;t)$ (5), which shows the consistency of this method. We note that a comparison of the RHS of Eqs. (59) and (63) reveals the fact that the density matrices has an infinite degeneracy in the group space.

5.3 Generic solution

The generic solution can be summarised as

[TABLE]

with the denominator

[TABLE]

The numerators are

[TABLE]

In D, we discuss how we verify that Eqs. (64) with Eqs. (65a)-(66e) indeed solve the nonlinear equations (6a)-(6c).

6 Properties of stationary states

Although the exact solution enables us to follow the time evolution of the density matrices, the expressions are very complicated to analyze. We instead discuss the effects of the coefficients on the stationary states by requiring their existence and positive semidefiniteness. We will also require the stationary states to satisfy a factorized condition in their coordinates. The latter requirement is strong enough to produce the known master equations and their generalizations.

6.1 Existence of stationary states

The existence of the stationary states depends on the sign of the damping constant $\gamma$. When $\gamma$ takes on negative value, from Eqs. (39d), (64), (65a), and (66a), we deduce that in the $t\rightarrow\infty$ limit, the dominant terms that govern the following expressions are

[TABLE]

regardless of real or imaginary $\omega$. As a result, the probability distribution function of the oscillator (8) vanishes in the limit $t\rightarrow\infty$ and the solution does not exist.

For positive $\gamma$, we recall from the discussion of Section 2.1 that both the underdamped and critically damped oscillators have imaginary $\omega$. Therefore, in the limit $t\rightarrow\infty$, we obtain the following behaviors,

[TABLE]

As a result, Eqs. (64)-(66e) give

[TABLE]

for the stationary states.

On the other hand, for an overdamped oscillator with real $\omega$, only $\omega<\gamma$ produces meaningful stationary states. In this case, a similar analysis shows that it behaves like Eq. (68), and hence its solution exhibits similar behaviours to the underdamped and critically damped oscillator (69a)-(69c). When an overdamped oscillator has $\omega>\gamma$ or $\omega=\gamma$, its $\mu(t)$ behaves like $e^{-(\omega-\gamma)t}$ and $e^{-\gamma t}$ in the limit $t\rightarrow\infty$, respectively. Hence the solution does not exist.

In summary, the stationary states of the density matrices of the underdamped, critically damped, and overdamped oscillator have similar behaviours. Stable solutions exist only when $\gamma>0$, and in the case of the overdamped oscillator, we need to limit the frequency to $\omega<\gamma$ for stable solutions to exist.

6.2 Positive semidefinite stationary states

The coefficient $\mu(t)$ characterizes the probability distribution function (8). A well-behaved probability distribution function requires $\mu(t)>0$ (9), or equivalently,

[TABLE]

As discussed in the previous section, for the overdamped oscillator we need to consider only $\omega<\gamma$. Consequently, in all types of oscillator, the prefactor on the RHS of Eq. (42) for $\Gamma$ can be written as a positive quantity, 1/\big{[}\gamma(\gamma^{2}-\mbox{\boldmath\theta}^{2})\big{]}, for real or imaginary $\omega$, multiplied by -\gamma^{2}\mbox{\boldmath\eta}+(\mbox{\boldmath\theta}\cdot\mbox{\boldmath\eta})\mbox{\boldmath\theta}+\gamma\mbox{\boldmath\theta}\wedge\mbox{\boldmath\eta}. After rewriting this expression in component form, Eq. (70) is equivalent to

[TABLE]

for both real or imaginary $\omega$.

The sufficient condition for the positive semidefiniteness of the stationary states is provided by $\nu(t)\geq 0$ (9), or from Eq. (69b),

[TABLE]

Similar to the analysis of $\Gamma$, -\mbox{\boldmath\Gamma}^{2} can be written as a positive quantity, 1/\big{[}\gamma^{2}(\gamma^{2}-\mbox{\boldmath\theta}^{2})\big{]}, for real or imaginary $\omega$, multiplied by (\mbox{\boldmath\theta}\cdot\mbox{\boldmath\eta})^{2}-\gamma^{2}\mbox{\boldmath\eta}^{2}. This means that we require

[TABLE]

Although the $\theta_{i}$ form the unitary part of the reduced dynamics, they affect the positivity of the stationary states through the \mbox{\boldmath\theta}\cdot\mbox{\boldmath\eta} term.

6.3 Factorized condition

The Wigner functions [2] of the Gaussian density matrix (5) takes the form,

[TABLE]

The coefficient $\kappa(t)$ determines whether or not the Wigner function can be factorized into two separate functions of $Q$ and $P$. If $\kappa(t)$ vanishes, we have

[TABLE]

We note that $g(P)$ is the momentum space distribution function conjugates to the position space distribution function $f(Q)$ only for pure states. If stationary states satisfy

[TABLE]

then their density matrices are factorized in the $Q,r$ coordinates, so do their Wigner function (75). It is in this sense that we call Eq. (78) a factorized condition. In terms of the coefficients of the master equation, Eq. (78) can be written as

[TABLE]

The Gibbs states,

[TABLE]

which are the stationary states of a system in thermal equilibrium with a thermal reservoir, are examples of stationary states obeying the factorized condition. This can be seen by comparing Eqs. (69a)-(69c) with $\rho_{\text{Gibbs}}$ to yield

[TABLE]

Examples of master equations with factorized stationary states are the Kossakowski-Lindblad (KL) equation [12, 13], and the Caldeira-Leggett (CL) equation [34], with generators given by

[TABLE]

respectively. The KL equation is found in quantum optical systems, whereas the CL equation is used to study quantum Brownian motion. We find that

[TABLE]

where

[TABLE]

in which $b_{\text{CL}}$ is the high temperature limit of $b$, and we use the units $\hbar=1$. Moreover, for $\theta_{1}=0$, the stationary states of the CL equation are the Gibbs states.

When the system is in contact with non-standard reservoir, or when the system is strongly correlated to the reservoir, the stationary states might not be the Gibbs states. For instance, the Hu-Paz-Zhang (HPZ) equation [35] with the generator

[TABLE]

is obtained under quite general conditions, including strong coupling to the reservoir. Though its stationary states are not the Gibbs states, they satisfy the factorized condition,

[TABLE]

In the following discussion, we are going to impose the factorized condition on the generic stationary state. We find that it provides a systematic way to classify the master equations.

6.4 General requirements on the coefficients

Before imposing the factorized condition on the stationary states, we first discuss some general requirements on the coefficients of the master equations. As discussed in Section 2.1, we observe that $\theta_{0}>0$ as it is related to the natural frequency of the oscillator. We assume that the coefficient $\theta_{1}$ renormalizes or shifts the natural frequency in a small amount, so that $|\theta_{1}|<|\theta_{0}|$, or $\theta_{0}-\theta_{1}>0$.

From the discussion in Section 2.1, the $\partial^{2}/\partial Q^{2}$-term causes diffusion-like effect on the probability distribution functions. Its coefficient $-\eta_{0}+\eta_{1}$ is required to be positive semidefinite,

[TABLE]

On the other hand, the $r^{2}$-term causes decoherence [27]. When this term becomes dominant at large $r$, it introduces to the density matrices a factor $\exp[(\eta_{0}+\eta_{1})r^{2}t]$ that decoheres their off-diagonal components. We then require

[TABLE]

From Eqs. (90) and (91) we conclude that $\eta_{0}$ and $\eta_{1}$ satisfy the condition,

[TABLE]

Due to Eq. (92) and the assumption that $|\theta_{1}|<|\theta_{0}|$, the effect of $\eta_{1}$ is usually less prominent than $\eta_{0}$.

6.5 Comparison with positive conditions obtained from the uncertainty principle

In Section 2.3 we found that by requiring second moments to satisfy the uncertainty principle, we arrive at the same necessary and sufficient condition (9) for the positive semidefiniteness of $\rho(t)$. On the other hand, Ref. [25] considered the time evolution of damped oscillator under the generator (in our notations)

[TABLE]

where $D_{pq}=D_{qp}$ and for simplicity, we assume that that there is no frequency shift to the oscillator. It was deduced that the uncertainty principle then imposes the following constraints on the diffusion coefficients [25]

[TABLE]

Consequently, it was concluded that the CL and HPZ equation must be rejected because they do not satisfy Eq. (94).

Comparing Eq. (93) with Eqs. (2a)-(3b), we find that

[TABLE]

From Eq. (95a), $\theta_{2}=-\gamma$ shows that the master equation considered is closely related to the class of CL and HPZ equations, cf. Table 1 at the end of Section 7.

The constraints on $D_{qq}$ and $D_{pp}$ in Eq. (94) are much stronger than Eqs. (90) and (91). The KL, HPZ equation and their conjugates could allow equalities in these conditions, cf. Table 1. We will also find that the equalities in Eqs. (90) and (91) also appear in other master equations that satisfy the factorized condition. Moreover, from the relation $16(D_{qq}D_{pp}-D_{qp}D_{pq})=\eta_{0}^{2}-\eta_{1}^{2}-\eta_{2}^{2}$, the third constraint in Eq. (94) can be written as

[TABLE]

On the other hand, in Section 2.3 we have shown that by requiring the second moments to satisfy the uncertainty principle at all time, we recover the positive condition (9). When this condition is applied to the stationary state of generic master equation, we obtain condition (73) which can also be written as

[TABLE]

Eq. (97) reveals that on the left hand side of the inequality, Eq. (97) contains two extra terms involving coefficients from the unitary part of the reduced dynamics compared to Eq. (96). Though the CL and HPZ equation do not satisfy Eq. (96), there are coefficients of the CL and HPZ equation that could satisfy Eq. (97). Furthermore, in the numerical analysis in Section 7, we show that there exist initial state and coefficients of the CL and HPZ equation that could generate positive evolution for all time.

The discrepancy between Eqs. (96) and (97) arises because the third inequality of Eq. (94) or (96) is too restrictive. Even if $0\leq D_{qq}D_{pp}-D_{qp}D_{pq}<\lambda^{2}/4$, there are coefficients of the master equation for which the sum of the first two positive terms in Eq. (10) of Ref. [25] are larger than its third negative term, to ensure that Eq. (10) of Ref. [25] is satisfied. Therefore, the positive semidefinite conditions we obtain are more general.

7 Master equations with factorized stationary states

Let us now impose the factorized condition (78) on generic stationary state. Keeping in mind that $\eta_{i}$ are independent coefficients, the master equations can be divided into two classes, according to either $\eta_{2}\neq 0$ or $\eta_{2}=0$.

$\eta_{2}\neq 0$.

We require the coefficient of $\eta_{2}$ in Eq. (79) to vanish, which divides the master equations into two classes, depending on whether $\theta_{2}=-\gamma$ or $\gamma$. In both cases, a zero $\eta_{0}$ of the master equation is not consistent the assumption that the frequency shift is small compared to the natural frequency. For, if $\eta_{0}=0$, Eq. (92) requires $\eta_{1}=0$. As a result, the positive condition -\mbox{\boldmath\Gamma}^{2}\geq 1 or Eq. (73) requires $\theta_{1}=\pm\theta_{0}$, which contradicts the assumption. As a result, we consider only $\eta_{0}\neq 0$.

1. (a)

$\theta_{2}=-\gamma$. HPZ equation.

The factorized condition further gives

[TABLE]

The equation does not contain the diffusion-like $\partial^{2}/\partial Q^{2}$-term. It produces the HPZ equation (88). In the Wigner representation, it is

[TABLE]

Furthermore, the conditions (71) and (73) give the following constraints on the coefficients of stationary states that are positive semidefinite,

[TABLE]

The requirement on the decoherence $r^{2}$ term, or the $-\partial^{2}/\partial P^{2}$ term in the Wigner representation, $\eta_{0}+\eta_{1}\leq 0$, then gives $\eta_{0}\leq 0$. Since $\eta_{0}=0$ is not permitted as discussed at the end of Section 6.4, $\eta_{0}<0$ for the HPZ equation. The HPZ equation with $\eta_{2}=0$ coincides with the CL equation (111). 2. (b)

$\theta_{2}=\gamma$. Conjugate of HPZ equation.

The factorized condition (79) reduces to $(-\theta_{0}+\theta_{1})(\eta_{0}+\eta_{1})=0$. The independence of $\theta_{0},\theta_{1}$ then requires

[TABLE]

which means that the equation does not contain the $r^{2}$ decoherence term. In the Wigner representation, this class of equation has the generator

[TABLE]

We shall call $K_{\text{cHPZ}}$ the conjugate of the HPZ equation. The positive semidefinite conditions (71) and (73) then give

[TABLE]

The conjugate equation also satisfies the requirements on the diffusion coefficients, Eqs. (90)-(92).

Let us define a unitary transformation that maps the coordinates $(Q,P)$ into $(P,-Q)$ (this is a canonical transformation if $Q,P$ are classical phase space coordinates [36]), together with inversions in the space components $i=1,2,$ of the vectors \mbox{\boldmath\theta},\mbox{\boldmath\eta}. A generator K(Q,P,\gamma,\mbox{\boldmath\theta},\mbox{\boldmath\eta}) that undergoes this transformation is denoted by

[TABLE]

where \mbox{\boldmath\theta}^{\prime}=(\theta_{0},-\theta_{1},-\theta_{2}), and etc. We find that

[TABLE]

Hence, $K_{\text{cHPZ}}$ is the image of $K_{\text{HPZ}}$ under this transformation in the Wigner representation. However, they are not physically equivalent because an exchange between $P$ and $Q$ in the Wigner representation will result in an exchange between $Q$ and $r$ in the space coordinates, which exchanges the roles of the probability component and the correlation component of the density matrices.

In Fig. 1, we illustrate the behaviors of the HPZ equation and its conjugate. The left plot of Fig. 1 depicts the evolution of $\nu(t)$ for $K_{\text{HPZ}}$ with the initial conditions $\mu_{0}=\nu_{0}=1,\kappa_{0}=1$ and the coefficients $\theta_{0}=2,\theta_{2}=-\gamma,\eta_{0}=\eta_{1}=-2\gamma b,\gamma=1,b=1$. The dotted, solid, long-dashed and dot-dashed curves label $K_{\text{HPZ}}$ with $\theta_{1}=0.5$, and $\eta_{2}=1,0,-1,-1.5$, respectively. We find that the positivity of the stationary states improves as $\eta_{2}$ reduces from 1 across the threshold of 0.875 down to $-1.5$. At the threshold, $\nu(t)$ approaches 0 in the limit $t\rightarrow\infty$, cf. Eq. (105). This behavior is consistent with Eq. (101). Though a large magnitude of $\eta_{2}$ is favourable for the positivity of the stationary states, however, the positivity in the early period of the evolution deteriorates. This is illustrated by the dot-dashed curve. A large $|\eta_{2}|$ also tend to violate the uncertainty principle [25]. On the other hand, the short-dashed curve shows that a negative $\theta_{1}=-0.5$ for $\eta_{2}=-1.5$ saves the positivity of its evolution. We conclude that large $|\eta_{2}|$ leads to negative evolution.

In the right plot of Fig. 1, we simulate the evolution of $\nu(t)$ for $K_{\text{cHPZ}}$ with initial conditions $\mu_{0}=\nu_{0}=1,\kappa_{0}=1$ and coefficients $\theta_{0}=2,\theta_{1}=0.5,\theta_{2}=\gamma,\eta_{0}=-\eta_{1}=-2\gamma b,\gamma=1,b=1$. The dotted, dot-dashed, solid, and long-dashed curves have $\eta_{2}=-3,-2,3,$ and 7, respectively. The short-dashed curve refers to $\theta_{1}=-0.5,\eta_{2}=7$, which is now positive. The threshold occurs at $\theta_{2}=-2.125$. The improvement of positivity occurs in the opposite direction in the value of $\eta_{2}$ compared to $K_{\text{HPZ}}$. Large values of $\eta_{2}$ again tend to destroy the positivity in the time evolution. Close to the threshold, the stationary state is positive but part of the early period of the evolution could be negative, as illustrated by the dot-dashed curve.

$\eta_{2}=0$.

If $\eta_{0}=0$, then Eq. (92) requires $\eta_{1}=0$. We find that well-behaved probability distribution functions cannot exist since Eq. (70) or (71) cannot be fulfilled. Consequently, we must consider $\eta_{0}\neq 0$.

There are three classes to consider. For $\eta_{1}\neq 0$, we can consider class (IIA) with $\theta_{2}\neq 0$ and class (IIB) with $\theta_{2}=0$, whereas for $\eta_{1}=0$, we consider class (IIC).

1. (a)

$\eta_{1}\neq 0$ and $\theta_{2}\neq 0$.

The factorized condition (78) becomes

[TABLE]

Since $\theta_{0},\theta_{1}$ are independent, we could have either $\theta_{1}\neq 0$ or $\theta_{1}=0$. For $\theta_{1}\neq 0$, the factorized condition (108) gives $\eta_{1}=-\eta_{0}\theta_{2}/\gamma=-\eta_{0}\gamma/\theta_{2}$, which permits only the solutions $\theta_{2}=-\gamma$ or $\gamma$, whereas for $\theta_{1}=0$, the factorized condition yields $\eta_{1}=-\eta_{0}\theta_{2}/\gamma$. There are altogether three possibilities.

1. i.

$\theta_{1}\neq 0,\theta_{2}=-\gamma,\eta_{1}=\eta_{0}$. CL equation.

The positive semidefinite conditions (71) and (73) yield

[TABLE]

This case gives $\eta_{0}+\eta_{1}=2\eta_{0}<0$ and $-\eta_{0}+\eta_{1}=0$, i.e., the $\partial^{2}/\partial Q^{2}$ diffusion-like term does not appear in the master equation. This gives the CL equation for quantum Brownian motion,

[TABLE] 2. ii.

$\theta_{1}\neq 0,\theta_{2}=\gamma,\eta_{1}=-\eta_{0}$. Conjugate of CL equation.

The positive semidefinite conditions (71) and (73) require

[TABLE]

This case gives $-\eta_{0}+\eta_{1}=-2\eta_{0}>0$ and $\eta_{0}+\eta_{1}=0$, i.e., the decoherence $r^{2}$ term does not appear in the master equation. In the Wigner representation, we have

[TABLE]

$K_{\text{cCL}}$ is the image of $K_{\text{CL}}$ under the unitary transformation defined in Eq. (106) in the Wigner representation,

[TABLE]

3. iii.

$\theta_{1}=0,\theta_{2}\neq 0,\eta_{1}=-\eta_{0}\theta_{2}/\gamma$. Generalized CL equation.

Applying the positive semidefinite conditions (71) and (73), we obtain

[TABLE]

At this stage, we need to impose $-\eta_{0}+\eta_{1}\geq 0$ (90) and $\eta_{0}+\eta_{1}\geq 0$ (91) to yield

[TABLE]

As a result, the Wigner representation of its dissipative part interpolates between the two extremes, i.e., the CL equation (111) with $\theta_{1}=0,\theta_{2}=-\gamma,\eta_{1}=\eta_{0}$, to its conjugate (114) with $\theta_{2}=\gamma,\eta_{1}=-\eta_{0}$,

[TABLE]

This equation generalizes the CL equation. It is the image of itself under the unitary transformation (106) in the Wigner representation,

[TABLE]

In Fig. 2, we illustrate the behaviors of the three types of CL equations in this class. The left plot depicts $K_{\text{CL}}$ and $K_{\text{cCL}}$ with initial conditions $4\mu_{0}=1/b_{0},\mu_{0}+\nu_{0}=b_{0}=0.6,\kappa_{0}=1$ and coefficients $\theta_{0}=2,\eta_{2}=0,\gamma=1$. The solid, dotted, and dot-dashed curves are $K_{\text{CL}}$ with $\theta_{2}=-\gamma,\eta_{0}=\eta_{1}=-2\gamma b$, and $(\theta_{1},b)=(0.2,1),(0.2,0.6)$, and $(-0.2,0.6)$, respectively. When this equation is applied to low temperature, it tends to give rise to unphysical behaviour, because the contribution from the vacuum fluctuation, i.e., the 1/2 term on the RHS of Eq. (86) that is neglected in $b_{\text{CL}}$ (87), becomes important [5]. This fact is illustrated in the solid and dotted curves when the temperature reduces from $b=1$ to 0.6. The low temperature curve gives negative evolution in the early interval of its evolution, even though its stationary state is positive. This occurs because of an inconsistent application of the CL equation to the low temperature environment. However, a negative $\theta_{1}=-0.2$ saves the positivity of the evolution in the dotted-dashed curve. The long-dashed and short-dashed curves are $K_{\text{cCL}}$ with $\theta_{2}=\gamma,\eta_{0}=-\eta_{1}=-2\gamma b$, and $(\theta_{1},b_{\text{CL}})=(0.2,0.6)$ and $(-0.2,0.6)$, respectively. Both curves show positive evolution.

The right plot of Fig. 2 refers to $K_{\text{gCL}}$ with initial conditions $4\mu_{0}=1/b_{0},\mu_{0}+\nu_{0}=b_{0}=0.6,\kappa_{0}=1$ and coefficients $\theta_{0}=2,\theta_{1}=0,\gamma=1,b=0.6,\eta_{0}=-2\gamma b,\eta_{1}=-\eta_{0}\theta_{2}/\gamma$. The curves interpolates from $\theta_{2}=-\gamma,-0.86\gamma,-0.553\gamma,0$, to $\gamma$, labeled by the dotted, dot-dashed, solid, long-dashed and short-dashed curves, respectively. The dotted and solid curves give the CL and the KL equation, respectively. The plots show that an increase in the value of $\theta_{2}$ improves the positivity of the evolution. In Section 8.2, we show that $K_{\text{gCL}}$ is completely positive when the constraint (135) is satisfied.

$\eta_{2}=0$.

1. (b)

$\eta_{1}\neq 0$ and $\theta_{2}=0$. Generalized KL equation type 1.

The factorized condition (78) yields $\eta_{1}=\eta_{0}\theta_{1}/\theta_{0}$. The case of $\theta_{1}=0$ overlaps with the situation $\eta_{1}=0$ considered in class (IIC). The positive semidefinite conditions (71) and (73) require

[TABLE]

This class generalizes the KL equation,

[TABLE]

$K_{\text{gKL1}}$ is the image of itself under the unitary transformation (106) in the Wigner representation,

[TABLE]

In the left plot of Fig. 3, we illustrate the behaviors of the generalized KL type 1 class of equations, with initial conditions $4\mu_{0}=1/b_{0},\mu_{0}+\nu_{0}=b_{0}=0.6,\kappa_{0}=1$ and coefficients $\theta_{0}=2,\theta_{2}=0,\eta_{0}=-2\gamma b,\eta_{1}=\eta_{0}\theta_{1}/\theta_{0},\gamma=1,b=0.6$. $K_{\text{gKL1}}$ interpolates between the threshold values $\theta_{1}=\pm 1.106$, cf. (122), when $\nu(t)$ approaches 0 in the limit $t\rightarrow\infty$. Outside the threshold values the stationary states are negative. The dotted, dot-dashed, solid, long-dashed and short-dashed curves label $\theta_{1}=-1.106,-0.5,0,0.5$, and 1.106, respectively. The solid curve with $\theta_{1}=0$ refers to the KL equation. The plots show that all the evolutions are positive.

$\eta_{2}=0$.

1. (c)

$\eta_{1}=0$. Generalized KL equation type 2.

The factorized condition (79) gives $\theta_{2}=\gamma\theta_{1}/\theta_{0}$. When $\theta_{1}=0$, we have $\theta_{2}=0$ to recover the KL equation (82). The positive semidefinite conditions (71) and (73) yield

[TABLE]

If we parameterize $\eta_{0}$ as $\eta_{0}=-2\gamma b$, then both conditions are always satisfied. In contrast to $K_{\text{gKL1}}$ in which the diffusion terms are interpolated, now it is the drift terms of $K_{\text{gKL2}}$ that are interpolated,

[TABLE]

It is the image of itself under the transformation (106)

[TABLE]

The right plot of Fig. 3 shows $K_{\text{gKL2}}$ with the initial conditions $4\mu_{0}=1/b_{0},\mu_{0}+\nu_{0}=b_{0}=0.6,\kappa_{0}=1$ and coefficients $\theta_{0}=2,\theta_{2}=\gamma\theta_{1}/\theta_{0},\gamma=1,\eta_{0}=-2\gamma b,\eta_{1}=0$. When $|\theta_{1}|$ is greater than the threshold value 1.789, the oscillator becomes overdamped. The plots shows the behaviors of $\nu(t)$ when $\theta_{1}$ interpolates from the overdamped region, -1.8 and 1.8 (dotted and short-dashed curves, respectively), into the underdamped region, -1 and 1 (dot-dashed and long-dashed curves, respectively). The solid curve with $\theta_{1}=0$ labels the KL equation. The plots of $\nu(t)$ are all positive as expected for completely positive $K_{\text{gKL2}}$, a fact which will be shown in Section 8.2

We observe that for all the equations from the class (II) and class (IB), $\eta_{0}<0$ guarantees that $\mu_{\text{st}}>0$. The only exceptions is the HPZ equation from class (IA), in which $\eta_{2}$ takes part in ensuring that $\mu_{\text{st}}>0$. In all the plots, the oscillators are in the underdamped region. For larger damping $\gamma$, $\theta_{2}=\pm\gamma$ could bring $\omega$ into the overdamped region. Due to the way we parameterized $\eta_{0}=-2\gamma b$, the positivity of the evolution is not much affected.

In summary, the HPZ equation and its conjugates, the CL equation, its conjugates and generalizations, and the two types of generalized KL equations (including the KL equation), exhaust the list of equations that are consistent with the reduced dynamics and simultaneously possess stationary states with factorized Wigner functions. Whether the conjugate and generalized equations have any realizations in real system is not known. We summarized the results of this section in Table 1.

8 Discussions

8.1 Gibbs states as stationary states

After imposing the factorized condition, if we further require $\Gamma_{1}=0$, then we obtain an even smaller subset of master equations whose stationary states are the Gibbs states. All of them will have $\Gamma_{0}=-\eta_{0}/\gamma$, or $\Gamma_{0}=2b$ if we parameterize $\eta_{0}=-2\gamma b$. We list the corresponding master equations below.

Both the HPZ equations (IA) and its conjugate (IB) will be subjected to the additional constraint

[TABLE]

If $\theta_{1}=0$, then it reduces to the CL equation (83) with $\theta_{1}=0$.

Both the CL equations (IIAa) and its conjugate (IIAb) should have $\theta_{1}=0$ to ensure that $\Gamma_{1}=0$. In contrast, the generalized CL equations (IIAc) automatically satisfies $\Gamma_{1}=0$. Hence, all the stationary states of the generalized CL equations are the Gibbs states.

The generalized KL equations of type 1 (IIB) and type 2 (IIC) also should have $\theta_{1}=0$. As a result, both of them reduce to the KL equation (82).

8.2 Completely positive generators

The plots in Fig. 3 suggest that $K_{\text{gKL1}}$ and $K_{\text{gKL2}}$ are positive operators. We can show that this is indeed true by writing $K_{\text{gKL1,2}}$ in the form of completely positive operators, which have the generic form

[TABLE]

where $a$ and $a^{\dagger}$ are the annihilation and creation operators of the oscillator, $c_{i}$ and $d_{i}$ are complex coefficients. Eq. (130a) can be written in the form $K_{\text{CP}}=\gamma(O_{0}-1/2)+\eta_{0}O_{+}+\eta_{1}L_{1+}+\eta_{2}L_{2+}$, provided

[TABLE]

A possible realization of $c_{i}$ and $d_{i}$ is

[TABLE]

where $c,s_{i}$ are real. Solving for $c$ and $s_{i}$ from Eqs. (131a)-(131d), we find that

[TABLE]

Since $c,s_{i}$ are real, the expressions under the square roots must be positive semidefinite, which lead to the inequalities,

[TABLE]

For $K_{\text{gKL1}}$, substituting $\eta_{1}=\eta_{0}\theta_{1}/\theta_{0}$ into the second inequality of Eq. (134) returns the positive condition in Eq. (122). Therefore, real solutions always exist and $K_{\text{gKL1}}$ is completely positive. For $K_{\text{gKL2}}$ with $\eta_{1}=0$, it is completely positive provided $\gamma\leq|\eta_{0}|$, which is always true if we parameterize $\eta_{0}=-2\gamma b$.

We can carry out the same consideration on $K_{\text{gCL}}$. We find that for$\eta_{1}=-\eta_{0}\theta_{2}/\gamma$, real solutions exist provided $|\theta_{2}|<\gamma$. If we parameterize $\eta_{0}=-2\gamma b$, then there always exist real solutions at high enough temperature, or large $b$, that fulfill the constraint,

[TABLE]

When $\theta_{2}=\pm\gamma$, Eq. (135) cannot be fulfilled with finite $b$. Hence, $K_{\text{gCL}}$ is completely positive within the constraint (135). We should emphasize that being a not completely positive generator does not rule out the possibility of the master equation to generate positive evolution for a given initial state provided the coefficients of the master equation are chosen appropriately. This is illustrated in the right plot of Fig. 2. The threshold value of $\theta_{2}$ within which $K_{\text{gCL}}$ is completely positive is $|\theta_{2}|\leq 0.553\gamma$. However, the dot-dashed curve with a value of $\theta_{2}=0.86\gamma$ lying outside the threshold still gives positive evolution, although it is not completely positive.

It is impossible for $K_{\text{CL}},K_{\text{cCL}},K_{\text{HPZ}}$ and $K_{\text{cHPZ}}$ to satisfy the second inequality of Eq. (134). For these generators, $\eta_{1}=\pm\eta_{0}$, Eq. (134) could be satisfied if $\gamma=0$, in which case we have the trivial case of a free oscillator.

8.3 Dissipation caused by unitary components of the generator

We note that even though $iM_{2}$ belongs to the unitary part of the dynamics, i.e. $\exp(-\theta_{2}iM_{2}t)$ is unitary, but it could affect the positivity of the time evolution through the term (\mbox{\boldmath\theta}\cdot\mbox{\boldmath\eta})\mbox{\boldmath\theta} and \mbox{\boldmath\theta}\wedge\mbox{\boldmath\eta} in \mbox{\boldmathg}(t), cf. Eq. (39d). Mathematically, this is because though the sets of operators $J_{0}=\{iL_{0},iM_{1},iM_{2}\}$ and $J_{+}=\{O_{+},L_{1+},L_{2+}\}$ are closed separately under the commutator brackets, but $[J_{0},J_{+}]\in J_{+}$ [17]. Hence, the coefficients of $J_{0}$ could influence the dissipative part of the dynamics even though they belong to the unitary part.

For a completely different reason, the operator $\exp(vO_{0})$ with real $v$ is unitary [37, 38, 39, 17] but it causes dissipation because it is not factorizable [40, 17]. As a consequence, it maps pure states into mixed states. For example, the oscillator ground state $|0\rangle\langle 0|$ is mapped by $\exp(vO_{0})$ into the Gibbs state with finite temperature, provided that $\tanh(v/2)=(2b-1)/(2b+1)$ [40].

9 Conclusion

We obtain the generic solution to the most general bilinear master equation with constant coefficients for a quantum oscillator in the form of Gaussian. The properties of the stationary states are determined by the components of a three-dimensional vector in the Minkowski space. We show that a factorized condition on the Wigner function of the stationary states is sufficient to generate a generic class of master equations that includes the well-known ones as special cases. In addition to this, it also generates their conjugates and generalizations. We also show that the generalized KL equations and some of the generalized CL equations are completely positive. For master equations that do not have completely positive generators, although positive semidefinite stationary states are not sufficient to warrant a positive evolution for a given initial state, they serve as references to identify the coefficients of the master equations that are able to generate positive evolution.

Acknowledgments

We thank Professor Sujin Suwanna, Dr. Fattah Sakuldee and the Department of Physics, Faculty of Science for hospitality and interesting discussions during our visit at the University of Mahidol, Bangkok, Thailand. This work is supported by the Ministry of Higher Education Malaysia (MOHE) under the Fundamental Research Grant Scheme (FRGS), Project No. FRGS/2/2014/ST02/UNIM/02/1.

Appendix A Matrix representation of the time evolution operator

We denote an anti-commutator bracket between $X,Y$ by the notation

[TABLE]

The identity operator that accompanies $O_{0}$ gives rise to an overall factor $\exp(\gamma t/2)$. For simplicity, we extract it from $K$ and consider $K^{\prime}=K+\gamma/2$ (47). We can then rewrite the time evolution operator as

[TABLE]

since

[TABLE]

In terms of the $M_{0}\equiv L_{0}$ and $M_{\pm}\equiv M_{1}\pm iM_{2}$ operators, where $\underline{M}_{\pm}^{2}=\underline{0}$, $K_{0}$ then takes the form

[TABLE]

Consequently, we deduced that

[TABLE]

and so on, where

[TABLE]

where we have made use of the identities,

[TABLE]

where $\Theta$ is defined by Eq. (31). Bold face letters denote three-dimensional Minskowski space vectors introduced in Section 3.

Eqs. (141a)-(141c) suggest that $A^{n}_{{\underline{K}}^{\prime}}(\underline{I})$ can be simplified by writing them in terms of $2\times 2$ matrices as follows,

[TABLE]

where $\underline{A}$ is a non-hermitian matrix

[TABLE]

and the column matrices denote

[TABLE]

As a result, the time evolution operator can be written as

[TABLE]

We can diagonalize this expression by using the biorthogonal basis of $\underline{A}$. The right and left eigenvalue problem of the non-hermitian matrix $\underline{A}$ are

[TABLE]

respectively. The eigenvalues are

[TABLE]

and the right and left eigenvectors are

[TABLE]

respectively. The eigenvectors form a pair biorthogonal basis

[TABLE]

We can form two matrices

[TABLE]

which satisfy the relations

[TABLE]

We can now diagonalize $\underline{A}$ through

[TABLE]

As a result,

[TABLE]

After simplifying this expression and making use of Eqs. (142c), (143c) and (143d), we finally obtain the matrix representation of the time evolution operator through

[TABLE]

which is given in Eqs. (30)-(30b). Here we obtain the BCH formula for real $\omega$ for the overdamped oscillator. After we obtain the solutions to the master equation, we analytically continue the solutions to imaginary $\omega$ for the underdamped oscillator.

Appendix B Coefficients of the BCH formula

Multiplying the matrix on the RHS of Eq. (4), we can simplify the BCH formula to yield

[TABLE]

where we have omitted the time-dependence of the coefficients to simplify the expression.

Comparing the coefficient of $\underline{I}$ in Eq. (B) and Eq. (30), we obtain Eq. (39a). Furthermore, a comparison of the coefficients of $\underline{O}_{0}$, $\underline{M}_{0}$ and $\underline{M}_{\pm}$ gives

[TABLE]

and Eq. (39c), respectively, where the unit vector $\hat{\theta}$ is defined in Eq. (40). Eqs. (206b) and (39c) are consistent with the coefficients obtained by comparing the coefficients of $\underline{M}_{0}\underline{O}_{0},\underline{M}_{+}\underline{O}_{0}$ and $\underline{M}_{-}\underline{O}_{0}$ in Eq. (B) and Eq. (30). Eqs. (206a)-(206b) can be solved for $m_{0}$ to give Eq. (39b). We note that the solutions to $m_{0},m_{\pm}$ are the same as those obtained in Ref. [41].

Substituting Eqs. (206b) and (39c) into the coefficients of $\underline{O}_{+},\underline{L}_{1+}$ and $\underline{L}_{2+}$ in Eq. (B), and equating them to the coefficients of Eq. (30), we obtain

[TABLE]

where $\Sigma$ is defined in Eq. (30b). Eq. (207) can be inverted for $g_{0},g_{1}$ and $g_{2}$ to yield Eq. (39d).

Comparing the coefficient of $\underline{M}_{0}\underline{O}_{+}$ in Eq. (B) and Eq. (30), we obtain

[TABLE]

If we take the dot product of Eqs. (207) and (30b) with $\theta$ and -\mbox{\boldmath\theta}, respectively, and set them equal, we obtain Eq. (208). Hence, Eq. (208) does not provide new information.

Appendix C 4D matrix representation of the generators

For completeness, we list three other generators that together with the set of generators in Eqs. (3a)-(3b) form a complete basis, cf. Ref. [33],

[TABLE]

The generators in Eq. (209) do not conserve probability of the reduced dynamics. Their 4D matrix representations are

[TABLE]

The matrix representation of the set of generators (48) can be extracted using Eq. (26) directly or by solving Eqs. (27i)-(27p) and (216) for them. The results are

[TABLE]

Appendix D Proof of Eqs. (64)-(66e) as solutions to equations of motion

Starting from the expression of $D(t)$ in Eq. (65a), we take the time derivative to obtain

[TABLE]

where we substitute the time derivative of $g_{i}$ using Eq. (37e). Then, we obtain

[TABLE]

Comparing Eq. (219) with Eq. (6a), we find that

[TABLE]

The proof is therefore reduced to showing that the coefficients of the various exponents, i.e., $1,\exp(\gamma t),\exp(\pm\omega t)$, on the RHS of Eqs. (D) and (220) are equal. In the proof, we convert the terms g_{0}+g_{1},\mbox{\hat{\theta}}_{0}+\mbox{\hat{\theta}}_{1} into g_{0}-g_{1},\theta_{0}-\theta_{1},\mbox{\hat{\theta}}_{2}g_{2},(\mbox{\boldmath\hat{\theta}}\wedge\mbox{\boldmathg})_{2},\mbox{\boldmath\hat{\theta}}\cdot\mbox{\boldmathg} by using the following identities,

[TABLE]

The equation of motion of $\kappa$ gives

[TABLE]

where we have used Eq. (220). Comparing Eq. (223) with Eq. (6b), we need to prove that

[TABLE]

by showing that the coefficients of the various exponents on both sides of Eq. (224) are equal. The following identities of the projectors \mbox{\boldmath\Pi}_{\pm}(\mbox{\boldmath\eta}) could be helpful in the proof,

[TABLE]

Lastly, the equation of motion of $\mu+\nu$ gives

[TABLE]

Comparing Eq. (226) with Eq. (6c), we need to prove that

[TABLE]