Lindblad Dynamics and Disentanglement in Multi-Mode Bosonic Systems

Abstract

In this paper, we consider the thermal bath Lindblad master equation to describe the quantum nonunitary dynamics of quantum states in a multi-mode bosonic system. For the two-mode bosonic system interacting with an environment, we analyse how both the coupling between the modes and the coupling with the environment characterised by the frequency and the relaxation rate vectors affect dynamics of the entanglement. We discuss how the revivals of entanglement can be induced by the dynamic coupling between the different modes. For the system, initially prepared in a two-mode squeezed state, we find the logarithmic negativity as defined by the magnitude and orientation of the frequency and the relaxation rate vectors. We show that, in the regime of finite-time disentanglement, reorientation of the relaxation rate vector may significantly increase the time of disentanglement.

1. Introduction

The theory of open quantum systems is very important for describing the transfer and storage of quantum information. In quantum information theory, these processes are generally described in terms of completely positive trace-preserving maps known as the quantum channels (see, e.g., [1,2,3,4,5] for analysis of the mathematical structures related to the quantum channels). There is also a variety of master equations for the reduced density matrix derived using different assumptions and approximations [6,7,8,9,10,11].

In particular, within the Markov approximation, master equations can often be cast into the well-known Lindblad form [12,13,14] which preserves complete positivity of the dynamics. This equation is also sometimes referred to as the Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) equation.

Though the general physics behind Lindblad-type master equations has been extensively discussed (see, e.g., References [15,16,17,18]), the mathematical techniques developed for the analysis of bosonic systems are mostly applicable to the single-mode Lindblad equation [19,20,21,22,23,24,25] and cannot be employed to treat its multi-mode generalizations (the most general form of the multi-mode Lindblad equation is described, e.g., in [26]).

In a recent paper [27], Fock-like eigenstates of a Lindbladian are constructed using Lie algebras induced by the master equation for a linear chain of coupled harmonic oscillators. The special case of two coupled oscillators presents a family of models that has been the subject of intense studies [28,29,30,31,32,33,34,35,36,37] dealing with the dynamics of entanglement in open continuous-variable systems (see, e.g., [38] for a review).

Typically, this model assumes that there is no interaction between the centre-of-mass and relative-distance modes in the course of relaxation, so that the relaxation part of the Lindbladian is a sum of two commuting relaxation superoperators. By contrast to mechanical systems, in photonic systems with quantised polarisation modes interacting with an optically anisotropic environment, both the dynamical and environment-induced intermode couplings can be important as they manifest themselves in the birefringence and dichroism of absorption [39,40,41]. In this paper, we shall relax this assumption to explore the effects induced by intermode couplings arising from the interaction between the system and the environment. The effects related to the dynamics of entanglement will be our primary concern.

Entanglement is known to be a resource of vital importance for rapidly developing quantum technologies, such as quantum communications and quantum computations [42,43], and its dynamics have been extensively studied during the last two decades (a recent review can be found in [44]). For the cases of discrete and continuous variables, which are typically represented by two interacting oscillators and qubits in the pioneering papers [45,46], it was found that the decay time of entanglement (the time of disentanglement) may be shorter than the time of decoherence and, under certain conditions, revivals of entanglement may occur. Finite-time disentanglement, commonly known as the “sudden death of entanglement”, is another important effect [47,48,49,50,51,52].

The structure of this paper is as follows: In Section 2, after introducing the GKSL master equation for the multi-mode bosonic system interacting with the thermal bath, in Section 2.1, we derive dynamical equations for the mean values of operators that preserve their normal ordering. In Section 2.2, the method of characteristics is utilised to solve the dynamical equation for the normally ordered characteristic function, ${\chi }_N$. In Section 3, the general theoretical results are applied to the important special case of a two-mode bosonic system in the thermal bath. This system can be regarded as a model describing the propagation of quantised polarisation modes in an optical fibre [41]. In Section 3.1, we derive the analytical results needed to evaluate time dependence of the averages that entre the elements of the covariance matrix. The logarithmic negativity of the Gaussian states introduced in Section 3.2 is numerically studied in Section 3.3. Finally, in Section 4, we discuss the results and make some concluding remarks.

2. Lindblad Dynamics

2.1. Master Equation

The starting point of our analysis is the Markovian thermal bath version of the Lindblad equation for the density matrix of an N-mode bosonic system:

$$\begin{array}{cc}\hfill & \frac{\partial \widehat{\rho }}{\partial t}=\mathcal{L}\widehat{\rho }=-i\sum _{n,m=1}^N{\Omega }_{nm}{\mathcal{C}}_{{\widehat{a}}_n^{†}{\widehat{a}}_m}\widehat{\rho }\hfill \\ \hfill & +\sum _{n,m=1}^N{K}_{nm}\left({\mathcal{D}}_{{\widehat{a}}_m{\widehat{a}}_n^{†}}\widehat{\rho }+{\mathrm{e}}^{-z_T}{\mathcal{D}}_{{\widehat{a}}_n^{†}{\widehat{a}}_m}\widehat{\rho }\right)\hfill \end{array}$$

(1)

written in terms of two superoperators given by

$$\begin{array}{cc}\hfill & {\mathcal{C}}_{\widehat{A}\widehat{B}}:\widehat{\rho }\mapsto {\mathcal{C}}_{\widehat{A}\widehat{B}}\widehat{\rho }=[\widehat{A}\widehat{B},\widehat{\rho }],\hfill \\ \hfill & {\mathcal{D}}_{\widehat{A}\widehat{B}}:\widehat{\rho }\mapsto {\mathcal{D}}_{\widehat{A}\widehat{B}}\widehat{\rho }=2\widehat{A}\widehat{\rho }\widehat{B}-\widehat{B}\widehat{A}\widehat{\rho }-\widehat{\rho }\widehat{B}\widehat{A}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & =[\widehat{A},\widehat{\rho }\widehat{B}]+[\widehat{A}\widehat{\rho },\widehat{B}],\hfill \end{array}$$

(3)

where the dagger denotes Hermitian conjugation, $\widehat{\rho }$ is the density matrix representing the quantum state; ${\widehat{a}}_n^{†}$ (${\widehat{a}}_n$) is the creation (annihilation) operator of the nth mode; $[\widehat{A},\widehat{B}]=\widehat{A}\widehat{B}-\widehat{B}\widehat{A}$ stands for the commutator; ${\Omega }_{nm}$ ($K_{nm}$) is the element of the frequency (relaxation) matrix, $\mathsf{\Omega }$ ($\mathbf{K}$); $z_T$ is the dimensionless inverse temperature parameter given by

$$\begin{array}{c}\hfill z_T=\frac{\hslash {\Omega }_0}{k_{B}T},\end{array}$$

(4)

where ${\Omega }_0$ is the bare frequency, ℏ is the reduced Planck constant, $k_B$ is the Boltzmann constant and T is the temperature of the environment. The frequency and relaxation matrices are both Hermitian: $\mathsf{\Omega }={\mathsf{\Omega }}^{†}$ and $\mathbf{K}={\mathbf{K}}^{†}$. The relaxation matrix $\mathbf{K}$ with elements giving the rates of thermalization is also positive definite: $\mathbf{K}>0$.

Note that the system (1) is conveniently rewritten in the following alternative form:

$$\begin{array}{cc}\hfill & \frac{\partial \widehat{\rho }}{\partial t}=-i[\widehat{\Omega },\widehat{\rho }]-(1+{\mathrm{e}}^{-z_T})[\widehat{K},\widehat{\rho }]\hfill \\ \hfill & +2\sum _{n,m}{K}_{nm}\left([{\widehat{a}}_m,\widehat{\rho }{\widehat{a}}_n^{†}]+{\mathrm{e}}^{-z_T}[{\widehat{a}}_n^{†},\widehat{\rho }{\widehat{a}}_m]\right),\hfill \end{array}$$

(5)

where the operators

$$\begin{array}{cc}\hfill & \widehat{\Omega }=\sum _{n,m}{\Omega }_{nm}{\widehat{a}}_n^{†}{\widehat{a}}_m,\widehat{K}=\sum _{n,m}{K}_{nm}{\widehat{a}}_n^{†}{\widehat{a}}_m\hfill \end{array}$$

(6)

are determined by the frequency and relaxation matrices, $\Omega$ and $\Gamma$, through the Jordan mapping that maps a Hermitian matrix J to the quadratic boson operator $\widehat{J}$ as follows:

$$\begin{array}{c}\hfill J\mapsto \widehat{J}=\sum _{n,m}{J}_{nm}{\widehat{a}}_n^{†}{\widehat{a}}_m.\end{array}$$

(7)

Our next step is to deduce the dynamic equation for the mean value of an operator $\widehat{S}$: $\langle \widehat{S}\rangle =\mathrm{Tr}\left(\widehat{S}\widehat{\rho }\right)$. From Equation (1) combined with the algebraic identities

$$\begin{array}{cc}\hfill & \mathrm{Tr}\left(\widehat{S}{\mathcal{C}}_{\widehat{A}\widehat{B}}\widehat{\rho }\right)=\langle [\widehat{S},\widehat{A}]\widehat{B}\rangle +\langle \widehat{A}[\widehat{S},\widehat{B}]\rangle ,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & \mathrm{Tr}\left(\widehat{S}{\mathcal{D}}_{\widehat{A}\widehat{B}}\widehat{\rho }\right)=\langle \widehat{B}[\widehat{S},\widehat{A}]\rangle -\langle [\widehat{S},\widehat{B}]\widehat{A}\rangle \hfill \end{array}$$

(9)

we have

$$\begin{array}{cc}\hfill & \frac{\partial \langle \widehat{S}\rangle }{\partial t}=\sum _{n,m}(i{\Omega }_{nm}-{\Gamma }_{nm})\langle {\widehat{a}}_n^{†}[{\widehat{a}}_m,\widehat{S}]\rangle \hfill \\ \hfill & -\sum _{n,m}(i{\Omega }_{nm}+{\Gamma }_{nm})\langle [\widehat{S},{\widehat{a}}_n^{†}]{\widehat{a}}_m\rangle \hfill \\ \hfill & +2\sum _{n,m}{\Gamma }_{nm}{n}_T\langle [{\widehat{a}}_m,[\widehat{S},{\widehat{a}}_n^{†}]]\rangle ,\hfill \end{array}$$

(10)

where ${\Gamma }_{nm}=(1-{\mathrm{e}}^{-z_T})K_{nm}$ and $n_T={({\mathrm{e}}^{z_T}-1)}^{-1}$ is the mean number of thermal photons.

An important point is that, for a normally ordered operator $\widehat{S}$ with $\widehat{S}=:\widehat{S}:$, the algebraic operations that entre the averages on the right hand side of Equation (10) preserve normal ordering.

2.2. Characteristic Function

In this section, we concentrate on the normally ordered characteristic function given by

$$\begin{array}{c}\hfill {\chi }_N\left(\mathit{\alpha }\right)=\langle {\mathrm{e}}^{(\mathit{\alpha },{\widehat{\mathbf{a}}}^{†})}{\mathrm{e}}^{-({\mathit{\alpha }}^{*},\widehat{\mathbf{a}})}\rangle ,\end{array}$$

(11)

where

$$\begin{array}{c}\hfill ({\mathit{\alpha }}^{*},\widehat{\mathbf{a}})\equiv \sum _{i=1}^N{\alpha }_i^{*}{\widehat{a}}_i,(\mathit{\alpha },{\widehat{\mathbf{a}}}^{†})\equiv \sum _{i=1}^N{\alpha }_i{\widehat{a}}_i^{†},\end{array}$$

(12)

$f\left(\mathit{\alpha }\right)\equiv f({\alpha }_1,\cdots ,{\alpha }_N,{\alpha }_1^{*},\cdots ,{\alpha }_N^{*})$ and an asterisk will indicate complex conjugation.

By using Equation (10), it is rather straightforward to deduce the following equation for ${\chi }_N$:

$$\begin{array}{cc}\hfill & \frac{\partial {\chi }_N(\mathit{\alpha },t)}{\partial t}=\widehat{L}{\chi }_N(\mathit{\alpha },t)\hfill \\ \hfill & \widehat{L}=i\sum _{n,m=1}^N{\Omega }_{nm}{\widehat{D}}_{mn}^{(-)}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & -\sum _{n,m=1}^N{\Gamma }_{nm}({\widehat{D}}_{mn}^{(+)}+2n_T{\alpha }_m{\alpha }_n^{*}),\hfill \end{array}$$

(14)

where

$$\begin{array}{c}\hfill {\widehat{D}}_{mn}^{(\pm )}={\alpha }_m\frac{\partial }{\partial {\alpha }_n}\pm {\alpha }_n^{*}\frac{\partial }{\partial {\alpha }_m^{*}}.\end{array}$$

(15)

The temporal evolution of the characteristic function ${\chi }_N$ is governed by the dynamical Equation (13) supplemented with the initial condition

$$\begin{array}{cc}\hfill & {\chi }_N(\mathit{\alpha },0)\equiv {\chi }_{\mathrm{ini}}\left(\mathit{\alpha }\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\int \mathrm{d}\mu \left(\mathit{\beta }\right){\mathrm{e}}^{[({\mathit{\alpha }}^{*},\mathit{\beta })-(\mathit{\alpha },{\mathit{\beta }}^{*})]}P(\mathit{\beta },0),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & \mathrm{d}\mu \left(\mathit{\beta }\right)=\prod _{n=1}^N\mathrm{d}\mu \left({\beta }_i\right),\mathrm{d}\mu \left({\beta }_i\right)=\frac{{\mathrm{d}}^2{\beta }_i}{\pi },\hfill \end{array}$$

(17)

where

$$\begin{array}{c}\hfill ({\mathit{\alpha }}^{*},\mathit{\beta })\equiv \sum _{i=1}^N{\alpha }_i^{*}{\beta }_i\end{array}$$

(18)

and ${\chi }_N$ is expressed in terms of the Glauber–Sudarshan P function (quasidistribution), $P(\beta ,t)$, related to the P-representation of the density matrix as follows:

$$\begin{array}{cc}\hfill & \widehat{\rho }\left(t\right)=\int \mathrm{d}\mu \left(\mathit{\beta }\right)P(\beta ,t)\left|\mathit{\beta }\right.〉\left.〈\mathit{\beta }\right|.\hfill \end{array}$$

(19)

We can now employ the method of characteristics [53] to solve the above initial value problem. According to this method, we begin with the system of characteristic equations

$$\begin{array}{c}\hfill \frac{\partial {\alpha }_n}{\partial t}=-\sum _{m=1}^N{Q}_{nm}{\alpha }_m,\mathbf{Q}=i\mathsf{\Omega }-\mathsf{\Gamma }\end{array}$$

(20)

and its solution written in the matrix form as follows:

$$\begin{array}{c}\hfill \mathit{\alpha }={\mathrm{e}}^{-\mathbf{Q}t}{\mathit{\alpha }}_0\equiv \mathbf{A}(-t){\mathit{\alpha }}_0,\mathbf{A}\left(t\right)={\mathrm{e}}^{(i\mathsf{\Omega }-\mathsf{\Gamma })t}\equiv {\mathrm{e}}^{\mathbf{Q}t},\end{array}$$

(21)

where ${\mathit{\alpha }}_0=\mathit{\alpha }\left(0\right)$. It is not difficult to obtain the solution along the characteristic curves (21) given by

$$\begin{array}{cc}\hfill & {\chi }_N({\mathit{\alpha }}_0,t)=exp\left[-2n_T{\int }_0^t({\mathit{\alpha }}_0^{*},{\mathbf{A}}^{†}(-\tau )\mathsf{\Gamma }\mathbf{A}(-\tau ){\mathit{\alpha }}_0)\mathrm{d}\tau \right]\hfill \\ \hfill & \times {\chi }_{\mathrm{ini}}\left({\mathit{\alpha }}_0\right).\hfill \end{array}$$

(22)

We can now express ${\mathit{\alpha }}_0$ in terms of $\mathit{\alpha }$ with the help of Equation (21): ${\mathit{\alpha }}_0=\mathbf{A}\left(t\right)\mathit{\alpha }$ and use the identity

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\tau }{\mathbf{A}}^{†}(-\tau )\mathbf{A}(-\tau )=2{\mathbf{A}}^{†}(-\tau )\mathsf{\Gamma }\mathbf{A}(-\tau ),\end{array}$$

(23)

to transform Formula (22) into the final expression for the characteristic function:

$$\begin{array}{cc}\hfill & {\mathit{\chi }}_N(\alpha ,t)={\mathrm{e}}^{-({\mathit{\alpha }}^{*},\mathbf{B}\left(t\right)\mathit{\alpha })}{\chi }_{\mathrm{ini}}\left(\mathbf{A}\left(t\right)\mathit{\alpha }\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & \mathbf{B}\left(t\right)=n_T({\mathbf{I}}_N-{\mathbf{A}}^{†}\left(t\right)\mathbf{A}\left(t\right)),\hfill \end{array}$$

(25)

where ${\mathbf{I}}_N$ is the identity $N\times N$ matrix. This formula is our central analytical result that will be used in the subsequent sections.

3. Two-Mode System

In this section, we focus our attention on the task of computing the time dependence of the averages that can be regarded as one-point correlation functions. One of the approaches to this important problem is to derive and solve dynamical equations for the averaged operators. For example, Equation (10) and the algebraic identities for the bosonic bilinear forms (7)

$$\begin{array}{cc}\hfill & \sum _{n,m}{\Gamma }_{nm}\langle {\widehat{a}}_n^{†}[{\widehat{a}}_m,\widehat{J}]\rangle =\sum _{n,m}{(J\Gamma )}_{nm}{\widehat{a}}_n^{†}{\widehat{a}}_m\equiv \widehat{J\Gamma },\hfill \\ \hfill & \sum _{n,m}{\Gamma }_{nm}\langle [\widehat{J},{\widehat{a}}_n^{†}]{\widehat{a}}_m\rangle =\sum _{n,m}{(\Gamma J)}_{nm}{\widehat{a}}_n^{†}{\widehat{a}}_m\equiv \widehat{\Gamma J}\hfill \end{array}$$

(26)

can be utilised to deduce the equation for the mean value $\langle \widehat{J}\rangle$ given by

$$\begin{array}{cc}\hfill & \frac{\partial \langle \widehat{J}\rangle }{\partial t}=-i\langle \widehat{[J,\Omega ]}\rangle -\langle \widehat{\{J,\Gamma \}}\rangle +2n_T\mathrm{Tr}(J\Gamma ),\hfill \end{array}$$

(27)

where $\{A,B\}=AB+BA$ denotes the anticommutator.

For illustrative purposes, in what follows, we consider a photonic system with two orthogonally polarised quantised modes [40,41,54], thus restricting our analysis to the special case of the two-mode bosonic system with $N=2$. Then, the above general result can be applied to the so-called Stokes operators that can be expressed in terms of the Pauli matrices

$$\begin{array}{cc}\hfill & {\sigma }_1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\sigma }_2=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),{\sigma }_3=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\hfill \end{array}$$

(28)

as follows [55]:

$$\begin{array}{c}\hfill {\widehat{S}}_i=\sum _{n,m=1}^2{\sigma }_{nm}^{\left(i\right)}{\widehat{a}}_n^{†}{\widehat{a}}_m,\end{array}$$

(29)

where $0\le i\le 3$ and ${\sigma }_0={\mathbf{I}}_2$ is the $2\times 2$ identity matrix. The mean values of these operators are known as the Stokes parameters and describe the state of polarization of the photonic system (dynamical regimes of the Stokes parameters are studied in [54]). These parameters, however, appear to be insufficient for complete characterization of the quantum states that generally requires the knowledge of higher order moments of the Stokes operators [55]. For such moments, the approach based on dynamical equations quickly becomes rather involved and unnecessarily complicated.

3.1. Exact Dynamics of Averages and Covariance Matrix

An alternative method is to use a formula for the characteristic Function (24). The derivatives of this function can be easily evaluated, giving the expressions for the mean values of normally ordered operators. In particular, the second order moments are computed as follows:

$$\begin{array}{cc}\hfill & \langle {\widehat{a}}_i^{†}{\widehat{a}}_j\rangle \left(t\right)=-\frac{{\partial }^2{\chi }_N(\mathit{\alpha },t)}{\partial {\alpha }_i\partial {\alpha }_j^{*}}{|}_{\mathit{\alpha }=0}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =B_{ji}\left(t\right)+\sum _{n,m}{A}_{ni}\left(t\right)A_{mj}^{*}\left(t\right)\langle {\widehat{a}}_n^{†}{\widehat{a}}_m\rangle \left(0\right)\hfill \\ \hfill & \langle {\widehat{a}}_i^{†}{\widehat{a}}_j^{†}\rangle \left(t\right)=\frac{{\partial }^2{\chi }_N(\mathit{\alpha },t)}{\partial {\alpha }_i\partial {\alpha }_j}{|}_{\mathit{\alpha }=0}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & =\sum _{n,m}{A}_{ni}\left(t\right)A_{mj}\left(t\right)\langle {\widehat{a}}_n^{†}{\widehat{a}}_m^{†}\rangle \left(0\right).\hfill \end{array}$$

(31)

These moments determine the elements of the covariance matrix of our two-mode photonic system [56] and the relations (30) and (31) yield the starting point of our analysis in the subsequent section.

For the two-mode system, the frequency and relaxation matrices can be written as a linear combination of the Pauli matrices

$$\begin{array}{cc}\hfill & i\mathsf{\Omega }-\mathsf{\Gamma }=\sum _{k=0}^3(i{\omega }_k-{\gamma }_k){\sigma }_k=(i{\omega }_0-{\gamma }_0){\mathit{\sigma }}_0\hfill \\ \hfill & +(i\mathit{\omega }-\mathit{\gamma },\mathit{\sigma }),\hfill \end{array}$$

(32)

where $\mathit{\sigma }\equiv ({\mathit{\sigma }}_1,{\mathit{\sigma }}_2,{\sigma }_3)$; $\mathit{\omega }\equiv ({\omega }_1,{\omega }_2,{\omega }_3)$ and $\mathit{\gamma }\equiv ({\gamma }_1,{\gamma }_2,{\gamma }_3)$ are the frequency and the relaxation rate vector, respectively; and it is rather straightforward to obtain the matrix exponential for $\mathbf{A}\left(t\right)$ in the following explicit form:

$$\begin{array}{cc}\hfill & \mathbf{A}\left(t\right)={\mathrm{e}}^{(i{\omega }_0-{\gamma }_0)t}\{cosh\left((\gamma -i\omega )t\right){\mathit{\sigma }}_0\hfill \\ \hfill & +\frac{sinh\left(\right(\gamma -i\omega \left)t\right)}{\gamma -i\omega }(i\mathit{\omega }-\mathit{\gamma },\mathit{\sigma })\},\hfill \end{array}$$

(33)

where

$$\begin{array}{c}\hfill \gamma -i\omega =\sqrt{(\mathit{\gamma }-i\mathit{\omega },\mathit{\gamma }-i\mathit{\omega })}=\sqrt{{\sum }_{k=1}^3{({\gamma }_k-i{\omega }_k)}^2}.\end{array}$$

(34)

We can now use Formulas (30)–(34) to evaluate the covariance matrix, $\mathsf{\Sigma }$, of our two-mode bosonic system. This matrix can be defined in terms of the quadrature operators (quadratures), ${\widehat{x}}_i$ and ${\widehat{p}}_i$, expressed in terms of the annihilation and creation operators, ${\widehat{a}}_i$ and ${\widehat{a}}_i^{†}$, as follows (see, e.g., references [38,56]):

$$\begin{array}{c}\hfill {\widehat{\mathbf{r}}}_i\equiv \left(\begin{array}{c}{\widehat{x}}_i\\ {\widehat{p}}_i\end{array}\right)=\mathbf{C}\left(\begin{array}{c}{\widehat{a}}_i\\ {\widehat{a}}_i^{†}\end{array}\right),\mathbf{C}=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ -i& i\end{array}\right).\end{array}$$

(35)

In our case, the block structure of the covariance matrix is given by

$$\begin{array}{cc}\hfill & \mathsf{\Sigma }=\left(\begin{array}{cc}{\mathsf{\Sigma }}_{11}& {\mathsf{\Sigma }}_{12}\\ {\mathsf{\Sigma }}_{21}& {\mathsf{\Sigma }}_{22}\end{array}\right),{\mathsf{\Sigma }}_{ij}={\mathsf{\Sigma }}_{ij}^{\left(0\right)}-2\langle {\widehat{\mathbf{r}}}_i\rangle \otimes \langle {\widehat{\mathbf{r}}}_j\rangle ,\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & {\mathsf{\Sigma }}_{ij}^{\left(0\right)}=\left(\begin{array}{cc}\langle \{{\widehat{x}}_i,{\widehat{x}}_j\}\rangle & \langle \{{\widehat{x}}_i,{\widehat{p}}_j\}\rangle \\ \langle \{{\widehat{p}}_i,{\widehat{x}}_j\}\rangle & \langle \{{\widehat{p}}_i,{\widehat{p}}_j\}\rangle \end{array}\right),\hfill \end{array}$$

(37)

where the diagonal and off-diagonal block matrices

$$\begin{array}{cc}\hfill & {\mathsf{\Sigma }}_{ii}^{\left(0\right)}=\left(\begin{array}{cc}2\langle {\widehat{a}}_i^{†}{\widehat{a}}_i\rangle +1+2\mathrm{Re}\langle {\widehat{a}}_i^2\rangle & 2\mathrm{Im}\langle {\widehat{a}}_i^2\rangle \\ 2\mathrm{Im}\langle {\widehat{a}}_i^2\rangle & 2\langle {\widehat{a}}_i^{†}{\widehat{a}}_i\rangle +1-2\mathrm{Re}\langle {\widehat{a}}_i^2\rangle \end{array}\right),\hfill \\ \hfill & {\mathsf{\Sigma }}_{12}^{\left(0\right)}=2\mathbf{C}\left(\begin{array}{cc}\langle {\widehat{a}}_1{\widehat{a}}_2^{†}\rangle & \langle {\widehat{a}}_1{\widehat{a}}_2\rangle \\ \langle {\widehat{a}}_1^{†}{\widehat{a}}_2^{†}\rangle & \langle {\widehat{a}}_1^{†}{\widehat{a}}_2\rangle \end{array}\right){\mathbf{C}}^{†}=\left(\begin{array}{cc}{x}_1+x_2& y_1+y_2\\ -y_1+y_2& x_1-x_2\end{array}\right),\hfill \\ \hfill & {\mathsf{\Sigma }}_{21}={\mathsf{\Sigma }}_{12}^T,\langle {\widehat{a}}_1^{†}{\widehat{a}}_2\rangle =\frac{x_1+i{y}_1}{2},\langle {\widehat{a}}_1{\widehat{a}}_2\rangle =\frac{x_2+i{y}_2}{2}\hfill \end{array}$$

(38)

are expressed in terms of the averages: $\langle {\widehat{a}}_i^{†}{\widehat{a}}_j\rangle$ and $2\langle {\widehat{a}}_i{\widehat{a}}_j\rangle$.

3.2. Symplectic Eigenvalues and Logarithmic Negativity

In what follows, we shall restrict our analysis to the important special case of the Gaussian states. These states are characterised by the symplectic eigenvalues of the covariance matrix (36) that can be computed using the relations [38,56]:

$$\begin{array}{cc}\hfill & 2{\mu }_{\pm }^2=\Delta \pm \sqrt{{\Delta }^2-4det\mathsf{\Sigma }},\hfill \end{array}$$

(39)

where $\Delta =det{\mathsf{\Sigma }}_{11}+det{\mathsf{\Sigma }}_{22}+2det{\mathsf{\Sigma }}_{12}$ is known as the “seralian” invariant.

The separability and entanglement criteria for bipartite continuous variable systems formulated in terms of the covariance matrix [29,57,58] and more general criteria [59,60] involving higher-order correlators provide a number of entanglement witnesses. For example, such witnesses are used for the analysis of the experimental data presented in [61] and the coherent state quantum key distribution suggested in [62,63].

For two-mode Gaussian states, different measures of entanglement have been proposed. These include the entanglement of formation, the Bures distance and the Gaussian measures of entanglement [64,65,66]. In this work, we deploy the logarithmic negativity [56] as a useful quantifier of bipartite entanglement in Gaussian states given by

$$\begin{array}{c}\hfill E_{\mathcal{N}}\left(\widehat{\rho }\right)={log}_2\left|\right|{\widehat{\rho }}^{PT}{\left|\right|}_1=max\{0,-{log}_2{\lambda }_{-}\},\end{array}$$

(40)

where ${\left|\right|.\left|\right|}_1$ stands for the trace norm and ${\widehat{\rho }}^{PT}$ is the partial transpose of $\widehat{\rho }$. The right-hand side of Equation (40) gives the expression for $E_{\mathcal{N}}\left(\widehat{\rho }\right)$ in terms of the lowest symplectic eigenvalue ${\lambda }_{-}$ of the partially transposed density matrix, ${\widehat{\rho }}^{PT}$, with the symplectic eigenvalues given by

$$\begin{array}{cc}\hfill & 2{\lambda }_{\pm }^2={\Delta }_{-}\pm \sqrt{{\Delta }_{-}^2-4det\mathsf{\Sigma }},\hfill \end{array}$$

(41)

where ${\Delta }_{-}=det{\mathsf{\Sigma }}_{11}+det{\mathsf{\Sigma }}_{22}-2det{\mathsf{\Sigma }}_{12}$. The logarithmic negativity being an entanglement monotone (a quantity which cannot be increased using local operations and classical communication) is known to bound the distillable entanglement contained in $\widehat{\rho }$ [67]. Note that, recently, the link between the logarithmic negativity and the quadrature coherence scale introduced as a nonclassicality measure was studied in [68].

3.3. Numerical Results

We assume that the system is initially prepared in the two-mode squeezed state

$$\begin{array}{c}\hfill {\left|\eta \right.〉}_{\mathrm{sq}}=exp[\eta {\widehat{K}}_{+}-{\eta }^{*}{\widehat{K}}_{-}]\left|0\right.〉,\eta =r{\mathrm{e}}^{i\theta },\end{array}$$

(42)

where r is the squeezing parameter, ${\widehat{K}}_{+}={\widehat{a}}_1^{†}{\widehat{a}}_2^{†}$ (${\widehat{K}}_{-}={\widehat{a}}_1{\widehat{a}}_2$) is the raising (lowering) generator of the $\mathfrak{su}(1,1)$ Lie algebra, with the covariance matrix given by

$$\begin{array}{cc}\hfill & {\mathsf{\Sigma }}_{11}\left(0\right)={\mathsf{\Sigma }}_{22}\left(0\right)=cosh\left(2r\right){\mathit{\sigma }}_0,\hfill \\ \hfill & {\mathsf{\Sigma }}_{12}\left(0\right)=sinh\left(2r\right)[cos\left(\theta \right){\mathit{\sigma }}_3+sin\left(\theta \right){\mathit{\sigma }}_1].\hfill \end{array}$$

(43)

According to Equation (33), the frequency and the relaxation rate vectors, $\omega$ and $\gamma$, describe the dynamical and environment mediated intermode couplings, respectively. These vectors determine the regime of the dissipative dynamics and will be parametrised using the angular representation:

$$\begin{array}{cc}\hfill & \mathit{\gamma }=\Gamma (sin{\theta }_{\Gamma }cos{\varphi }_{\Gamma },sin{\theta }_{\Gamma }sin{\varphi }_{\Gamma },cos{\theta }_{\Gamma }),\hfill \\ \hfill & \mathit{\omega }=\Omega (sin{\theta }_{\Omega }cos{\varphi }_{\Omega },sin{\theta }_{\Omega }sin{\varphi }_{\Omega },cos{\theta }_{\Omega }).\hfill \end{array}$$

(44)

In order to investigate the effects of the intermode couplings in the dynamics of entanglement, we have computed the time dependencies of the logarithmic negativity for differently oriented vectors $\mathit{\gamma }$ and $\mathit{\omega }$. The list of parameters used in our calculations is as follows: $n_T=0.1$, $r=2.0$, $\Gamma /{\gamma }_0=0.9$ and $\theta ={\varphi }_{\Gamma }={\varphi }_{\Omega }=0$${}^{\circ }$. At these parameters, we have the initial values of the symplectic eigenvalues ${\mu }_{\pm }\left(0\right)=1$, ${\mu }_{1,2}\left(0\right)=cosh\left(2r\right)$ and ${\lambda }_{-}\left(0\right)={\mathrm{e}}^{-2r}$ giving the corresponding value for the logarithmic negativity, $E_{\mathcal{N}}\left(\widehat{\rho }\right){|}_{t=0}\approx 5.77$.

For the dynamics governed by the thermal bath Lindblad Equation (1), in the long time limit with $t\to \infty$, the density matrix approaches the equilibrium state ${\widehat{\rho }}_{\mathrm{eq}}=\widehat{\rho }(\infty )$ with the covariance matrix given by

$$\begin{array}{c}\hfill \mathsf{\Sigma }\left(t\right)\stackrel{t\to \infty }{\to }{\Sigma }_{\mathrm{eq}}=(1+2n_T){\mathbf{I}}_4.\end{array}$$

(45)

The equilibrium density matrix, ${\widehat{\rho }}_{\mathrm{eq}}\propto exp[-z_T({\widehat{a}}_1^{†}{\widehat{a}}_1+{\widehat{a}}_2^{†}{\widehat{a}}_2)]$, is the disentangled mixed Gaussian state with $E_{\mathcal{N}}\left({\widehat{\rho }}_{\mathrm{eq}}\right)=I(1;2){|}_{t\to \infty }=0$ and $S\left[{\widehat{\rho }}_{\mathrm{eq}}\right]\equiv S_{\mathrm{eq}}\approx 0.56$.

In Figure 1, we present the results on the temporal evolution of the logarithmic negativity computed at $\Omega =0$ for different values of the angle ${\theta }_{\Gamma }$. The curves clearly indicate the regime of finite-time disentanglement known as the “sudden death of entanglement”. In this regime, $E_{\mathcal{N}}$ vanishes at $t\ge t_{\mathrm{d}}$, where $t_{\mathrm{d}}$ is the time of disentanglement.

Figure 2 shows the disentanglement time calculated as a function of ${\theta }_{\Gamma }$. It can be seen that $t_{\mathrm{d}}$ is sensitive to orientation of the relaxation rate vector describing the intermode coupling induced by the interaction between the bosonic system and the environment.

The curves plotted in Figure 3 are computed at different values of the angle ${\theta }_{\Omega }$ that specify the orientation of the frequency vector with $\Omega /{\gamma }_0=\pi$ and illustrate the effect of the dynamical intermode coupling on entanglement dynamics. It is seen that oscillatory behaviour of the logarithmic negativity at ${\theta }_{\Omega }=90$${}^{\circ }$ translates into sudden revivals of entanglement as the angle ${\theta }_{\Omega }$ decreases.

4. Conclusions

In this paper, we have studied the dynamics of a multi-mode bosonic system governed by the thermal bath Lindblad master equation. Our general approach is based on an exact solution for the characteristic function obtained using the method of characteristics.

We have applied this approach to the special case of the two-mode bosonic system. In this case, the dynamics are determined by the intermode couplings that entre the dynamical and relaxation parts of the Lindblad superoperator $\mathcal{L}$ (see Equation (1)) and can be described in terms of the frequency and the relaxation rate vectors (see Equation (32)). We have focussed our attention on the effects of the intermode couplings in the dynamics of entanglement and have presented a number of the numerical results on time dependence of the logarithmic negativity.

It is found that, for a vanishing frequency vector with $\Omega =0$, the logarithmic negativity of the system initially prepared in the two-mode squeezed state monotonically decays, reaching zero at the disentanglement point in time, $t=t_{\mathrm{d}}$. In this regime of finite-time disentanglement, the disentanglement time appears to be dependent on the orientation of the relaxation rate vector. We have also shown that the presence of the dynamical intermode coupling with $\Omega >0$ complicates the dynamics of the logarithmic negativity, and its nonmonotonic behaviour results in revivals of entanglement.

Our theoretical considerations were motivated by the model describing mixed polarisation quantum states propagating in an optically anisotropic lossy environment [40,41]. Interestingly, our results are formulated in terms of the covariance matrix that can, in principle, be extracted from experimental data measured using either optical homodyne or heterodyne detection techniques [69]. Qualitatively, our results on the regime of finite-time disentanglement and revivals of entanglement are in agreement with the predictions for oscillators interacting with the thermal bath previously published in [28,29,32,33]. It turns out that the values of the disentaglement time reported in [29] and obtained in the limit of negligible relaxation anisotropy where $\Gamma =0$ are close to our estimate ${\gamma }_0{t}_{\mathrm{d}}\approx 0.6$ evaluated at ${\theta }_{\Gamma }=0$.

Our findings provide further insights about the protection techniques of entangled states from the detrimental effects of surrounding environments by suitably manipulating the intermode interaction. Generally, revivals of quantum correlations in composite quantum systems are a useful dynamical feature against these effects [70,71,72,73]. Different experimental methods and theoretical approaches to protect quantum resources have been put forward in [74,75,76,77,78,79,80,81].

We conclude with the remark that our analytical approach can be readily extended to study a number of problems, such as Gaussian Einstein–Podolsky–Rosen steering for two-mode squeezed states transmitted in lossy quantum communication channels [82] and the problems related to controlled quantum dynamics in a realistic setup involving open environments. Quantum navigation is an important class of controlled quantum dynamics whereby the objective is to transport one quantum state into another, or to generate quantum gates, in the shortest possible time under the influence of an uncontrollable external field. Problems of this kind can be thought of as representing the quantum counterpart of the classical Zermelo navigation problem of finding the time-optimal control that takes a ship from one location to another, under the influence of external wind or currents [83,84,85]. In a forthcoming publication, we will apply our results to the Zermelo navigation problem for open multi-mode bosonic systems.