Differential Parametric Formalism for the Evolution of Gaussian States: Nonunitary Evolution and Invariant States

Abstract

In the differential approach elaborated, we study the evolution of the parameters of Gaussian, mixed, continuous variable density matrices, whose dynamics are given by Hermitian Hamiltonians expressed as quadratic forms of the position and momentum operators or quadrature components. Specifically, we obtain in generic form the differential equations for the covariance matrix, the mean values, and the density matrix parameters of a multipartite Gaussian state, unitarily evolving according to a Hamiltonian $\widehat{H}$. We also present the corresponding differential equations, which describe the nonunitary evolution of the subsystems. The resulting nonlinear equations are used to solve the dynamics of the system instead of the Schrödinger equation. The formalism elaborated allows us to define new specific invariant and quasi-invariant states, as well as states with invariant covariance matrices, i.e., states were only the mean values evolve according to the classical Hamilton equations. By using density matrices in the position and in the tomographic-probability representations, we study examples of these properties. As examples, we present novel invariant states for the two-mode frequency converter and quasi-invariant states for the bipartite parametric amplifier.

1. Introduction

The study of Gaussian states has been of essential interest in the last few decades. These types of states, associated with classical random fields, were considered as a possibility to connect covariance matrices of the states as quantum density matrices and, with this definition, to study the quantum–classical relation of randomness with the quantization procedure [1,2]. The problems of the new developments of the foundations of quantum mechanics and applications of new results in quantum information and quantum probabilities, as well as in areas like mathematical finance and economics have attracted the attention of the researchers; they are intensely discussed in the literature [3,4,5]. An important role in this development is played by discussing the problems that appeared from the very beginning of quantum mechanics, like the notion of quantum system states and the interpretation of the states associated in the conventional formulation of quantum mechanics with Hilbert space vectors and density operators, using the quasiprobability distributions and the probability distributions containing the complete information on quantum states. There exists increasing interest in quantum foundations since a deeper understanding of the essence and formalism of quantum theory is needed for the development of quantum technologies and the possibilities to extend the applications of quantum formalism in physics to all other areas of science like the economy, finance, and social disciplines.

Some examples of the Gaussian states of quantum fields as the coherent, squeezed, and thermal light states are regularly used in the theoretical and experimental framework of quantum mechanics, optics, information, and computing. The use of these states in quantum information has been of particular importance [6,7,8]. One can list some of the most recent applications of the use of Gaussian systems: it has been demonstrated [9] that it is not possible to distill more entanglement from a bipartite Gaussian state, using local Gaussian transformations. In [10], several properties of the purity of Gaussian states were found. The connection between the symplectic invariants of bipartite Gaussian states, the von Neumann entropy, and the mutual information was established in [11]. The extremality of entanglement measures and secret key rates for Gaussian states was observed in [12]. It was shown [13] that Gaussian attacks are characterized by an optimum efficiency against eavesdropping protocols. Quantum illumination of a target using Gaussian light states was explored by Tan et al. [14]. A quantum discord for systems of continuous variables, such as Gaussian states, was implemented in [15]. In [16], an invariant describing the nonclassicality in a two-mode Gaussian state was reported. The entanglement of m modes with other n modes of a Gaussian multipartite system was treated in [17]. The linear response for systems close to steady-states under Gaussian processes was obtained in [18]. The optimal measurement of the fidelity of multimode Gaussian states was studied in [19]. On the other hand, the study of Gaussian wave packets by nonlinear differential equations, as the Riccati equation, was done in [20,21,22]. Several coherent states have been defined by the use of quadratic operators [23]. The behavior of different quantities as covariances in thermal relaxation phenomena was also studied in [24].

The aim of this work is to present a new way to characterize the dynamics of Gaussian states using the differential equations for the parameters, which determine their continuous variable density matrices. The proposed method makes use of the integrals of the motion of such systems, and it can be used to clarify new aspects of multimode Gaussian quantum states, such as an explicit form of the nonunitary evolution of the states of subsystems and the existence of invariant states with constant covariance matrices and mean values.

The time evolution of a quantum system was first established by Schrödinger [25]. The dynamics of the system given by a Hamiltonian operator $\widehat{H}$ for a pure state $\left|\psi \right(t)\rangle$ must follow the Schrödinger equation:

$$\widehat{H}|\psi \left(t\right)\rangle =i\hslash \frac{\partial }{\partial t}|\psi \left(t\right)\rangle ;$$

this expression corresponds to a second-order differential equation in the position representation. In the case of an arbitrary state represented by the density matrix $\widehat{\rho }\left(t\right)$, which cannot be pure [26,27], the evolution is determined by the von Neumann equation:

$$i\hslash \frac{\partial }{\partial t}\widehat{\rho }\left(t\right)=[\widehat{H},\widehat{\rho }\left(t\right)];$$

the general solution of this equation is given by the unitary transform $\widehat{U}\left(t\right)$, i.e., $|\psi \left(t\right)\rangle =\widehat{U}\left(t\right)|\psi \left(0\right)\rangle$ or $\widehat{\rho }\left(t\right)=\widehat{U}\left(t\right)\widehat{\rho }\left(0\right){\widehat{U}}^{†}\left(t\right)$, where $\left|\psi \right(0)\rangle$ and $\widehat{\rho }\left(0\right)$ describe the system at time $t=0$. Furthermore, it is a common knowledge that, when the system interacts with an environment, its dynamics is described by the master equation [28,29,30].

The Gaussian states can be determined by their covariance matrix $\sigma$ and mean values $\langle q_j\rangle$ and $\langle p_j\rangle$. This property also implies that the evolution of a Gaussian state can be obtained, if the time dependence of these parameters is known. In this work, we review the differential equations that the covariance matrix and the mean values satisfy [31,32]; employing these results, we can define differential equations for the density matrix parameters of a general multimode state satisfying these equations and then use the equations to discuss some physical characteristics of the unitary and nonunitary evolutions of Gaussian states.

This paper is organized as follows.

In Section 2, the evolution of non-pure Gaussian states for a one-dimensional quadratic Hamiltonian is presented. To obtain this evolution, we make use of the derivatives of the covariance matrix, the mean values, and the parameters of the density operator; also, we define and obtain invariant states for this system. The generalization of these results to the case of a multidimensional quadratic system is explored in Section 3. Examples of the application of the general results to the nonunitary evolution of the subsystems of a two-mode state, as well as the definition of invariant and quasi-invariant states are given in Section 4. Furthermore, in Section 5, we obtain new invariant states for the frequency converter and quasi-invariant states for the parametric amplifier. The detection of these invariant states using the quantum tomographic representation of the states is discussed for single-mode Gaussian states in Section 5 and for the bipartite system in Section 6; in these sections, the correspondence between the time-independent states and thermal density matrices is mentioned. Finally, we give our conclusions.

2. One-Dimensional Quantum Quadratic Hamiltonian and Its Linear Invariant Operators

In this section, we analyze some properties of the one-dimensional quadratic Hamiltonian. In particular, we are interested in the invariant operators, which in the quadratic case happen to be linear in the quadrature operators $\widehat{p}$ and $\widehat{q}$.

The most general (in a unit system where $\hslash =m=1$), one-dimensional quantum quadratic Hamiltonian can be obtained in terms of the quadrature operators $\widehat{p}$ and $\widehat{q}$ as follows:

$$\widehat{H}=(\widehat{p},\widehat{q})\left(\begin{array}{cc}{\omega }_1\left(t\right)& {\omega }_2\left(t\right)\\ {\omega }_2\left(t\right)& {\omega }_3\left(t\right)\end{array}\right)\left(\begin{array}{c}\widehat{p}\\ \widehat{q}\end{array}\right)+(\widehat{p},\widehat{q})\left(\begin{array}{c}{\delta }_1\\ {\delta }_2\end{array}\right),$$

(1)

where the parameters ${\omega }_1\left(t\right)$, ${\omega }_2\left(t\right)$, ${\omega }_3\left(t\right)$, and ${\delta }_{1,2}$ are real functions of time. The dynamics associated with this Hamiltonian can be solved by different methods. One of them is the method of time-dependent invariants (integrals of motion) [33,34]. These invariants are quantum operators $\widehat{R}\left(t\right)$, whose total time derivative is equal to zero ${\displaystyle \frac{d\widehat{R}\left(t\right)}{dt}}=0$. In the quadratic case, it is known that there exist invariants linearly depending on the quadrature operators $\widehat{p}$ and $\widehat{q}$, i.e., $\widehat{R}\left(t\right)={\lambda }_1\left(t\right)\widehat{p}+{\lambda }_2\left(t\right)\widehat{q}+{\lambda }_3\left(t\right)$.

By substituting this expression into the von Neumann equation, which determines the dynamics of $\widehat{R}$, i.e., ${\displaystyle \frac{d\widehat{R}\left(t\right)}{dt}}={\displaystyle \frac{i}{\hslash }}[\widehat{H}\left(t\right),\widehat{R}\left(t\right)]+{\displaystyle \frac{\partial \widehat{R}\left(t\right)}{\partial t}}=0$, one can show that $\widehat{R}\left(t\right)$ is an invariant operator, if the following differential equations are satisfied:

$$\begin{array}{c}\hfill {\dot{\lambda }}_1=2({\omega }_2{\lambda }_1-{\omega }_1{\lambda }_2),{\dot{\lambda }}_2=2(-{\omega }_2{\lambda }_2+{\omega }_3{\lambda }_1),{\dot{\lambda }}_3={\delta }_2{\lambda }_1-{\delta }_1{\lambda }_2.\end{array}$$

(2)

We point out that parameters ${\lambda }_{1,2}$ can be obtained by solving the classical Hamilton equations:

$$\begin{array}{c}\hfill \dot{p}=-2({\omega }_{2}p+{\omega }_{3}q)-{\delta }_2,\dot{q}=2({\omega }_{2}q+{\omega }_{1}p)+{\delta }_1.\end{array}$$

(3)

with ${\delta }_1={\delta }_2=0$. To show this, one can see that the differential equations for ${\lambda }_{1,2}$ correspond to the classical equations with the substitution ${\lambda }_1\to q$ and ${\lambda }_2\to -p$; in other words, they correspond to the time inversion of the classical equations. In the case of the differential equation for ${\lambda }_3$, one can show, in view of the Hamilton equations, that it corresponds to the classical Lagrangian (with ${\delta }_{1,2}\ne 0$) plus the time variation of the function $pq$ of the system, that is,

$${\dot{\lambda }}_3=-2\mathcal{L}+\dot{L},$$

(4)

where $\mathcal{L}=p\dot{q}-H$ and $L=pq$. From these identifications, one can conclude that the classical dynamics given by the Hamilton equation or the equation of motion can lead to the solution of the quantum dynamics given by the Hamiltonian operator of Equation (2). For example, one can derive the propagator of the system $G(x,x^{\prime },t)=\langle x|\widehat{U}\left(t\right)|x^{\prime }\rangle$ using each one of the solutions of the classical problem [35].

2.1. Dynamics of Non-Pure States

Here, we demonstrate that the dynamics of a generic Gaussian state, which may not be pure, can be given by solving differential equations for the covariance matrix or the density matrix parameters. We show that these differential equations imply the invariance of the determinant of the covariance matrix, when the time evolution is unitary.

As discussed above, the propagator of the system can be obtained using the time-dependent invariants resulting from the solution to Equation (2) for two sets of initial conditions: ${\lambda }_1\left(0\right)=1$, ${\lambda }_2\left(0\right)=0$, ${\lambda }_3\left(0\right)=0$ and ${\lambda }_1\left(0\right)=0$, ${\lambda }_2\left(0\right)=1$, ${\lambda }_3\left(0\right)=0$. These two sets define two different invariants called $\widehat{P}$ and $\widehat{Q}$, respectively, which can be written as:

$$\left(\begin{array}{c}\widehat{P}\\ \widehat{Q}\end{array}\right)=\Lambda \left(\begin{array}{c}\widehat{p}\\ \widehat{q}\end{array}\right)+\left(\begin{array}{c}{\lambda }_3\\ {\lambda }_6\end{array}\right),\mathrm{with}\Lambda =\left(\begin{array}{cc}{\lambda }_1& {\lambda }_2\\ {\lambda }_4& {\lambda }_5\end{array}\right),$$

(5)

where ${\lambda }_{4,5,6}$ satisfy the same differential equations as ${\lambda }_{1,2,3}$, respectively, with the different sets of initial conditions mentioned above. The operators $\widehat{P}$ and $\widehat{Q}$ fulfill the commutation relation $[\widehat{Q},\widehat{P}]=i$, implying the relation ${\lambda }_1{\lambda }_5-{\lambda }_2{\lambda }_4=1$ and the fact that the matrix $\Lambda$ is symplectic.

Furthermore, they can be related to operators $\widehat{p}$ and $\widehat{q}$ through the evolution operator as follows:

$$\widehat{P}=\widehat{U}\widehat{p}{\widehat{U}}^{†},\widehat{Q}=\widehat{U}\widehat{q}{\widehat{U}}^{†};$$

it is not difficult to show [33,34] that these expressions can be used to obtain the propagator of the system $G(x,x^{\prime },t)=\langle x|\widehat{U}\left(t\right)|x^{\prime }\rangle$, which reads:

$$\begin{array}{lll}G(x,x^{\prime },t)& =& \frac{1}{\sqrt{-2\pi i{\lambda }_4}}exp\left\{-\frac{i}{2{\lambda }_4}[{\lambda }_5{x}^2-2x{x}^{\prime }+{\lambda }_1{x^{\prime }}^2+2x{\lambda }_6+2x^{\prime }({\lambda }_3{\lambda }_4-{\lambda }_1{\lambda }_6)\right.\\ & & \left.+{\lambda }_1{\lambda }_6^2-2{\lambda }_4{\int }_0^t{\dot{\lambda }}_3{\lambda }_{6}d\tau ]\right\},\end{array}$$

(6)

which we immediately identify as a Gaussian function.

In view of this propagator, the dynamics of any initial state in the position representation can be found by the integration of the propagator and the wave function of the initial state. In this work, we will suppose the case of the initial state given by a generic mixed Gaussian system with the density matrix in the position representation $\rho (x,x^{\prime },t=0)=\langle x|\widehat{\rho }|x^{\prime }\rangle$ equal to:

$$\rho (x,x^{\prime },t=0)=Nexp\left\{-a_1{x}^2+a_{12}x{x}^{\prime }-a_1^{*}{x^{\prime }}^2+b_{1}x+b_1^{*}{x}^{\prime }\right\},$$

(7)

with complex parameters $a_1=a_{1R}+i{a}_{1I}$ and $b_1=b_{1R}+i{b}_{1I}$ and the real parameter $a_{12}\in \mathbb{R}$, which also satisfy the integrability conditions $a_{1R}>a_{12}/2\ge 0$ and have a normalization constant N expressed as:

$$N={\left(\frac{a_1+a_1^{*}-a_{12}}{\pi }\right)}^{1/2}exp\left\{-\frac{{(b_1+b_1^{*})}^2}{4(a_1+a_1^{*}-a_{12})}\right\}.$$

As discussed above, the Gaussian states can be fully identified by their covariance matrix and the mean values of the quadrature components. In the case of the state (7), the mean values of the operators $\widehat{p}$ and $\widehat{q}$ are:

$$\langle \widehat{p}\rangle \left(0\right)=b_{1I}-\frac{2a_{1I}{b}_{1R}}{2a_{1R}-a_{12}},\langle \widehat{q}\rangle \left(0\right)=\frac{b_{1R}}{2a_{1R}-a_{12}},$$

(8)

and the initial covariance matrix of the system reads:

$$\sigma \left(0\right)=\left(\begin{array}{cc}{\sigma }_{pp}& {\sigma }_{pq}\\ {\sigma }_{pq}& {\sigma }_{qq}\end{array}\right)=\frac{1}{2(2a_{1R}-a_{12})}\left(\begin{array}{cc}4|a_1{|}^2-a_{12}^2& -2a_{1I}\\ -2a_{1I}& 1\end{array}\right).$$

(9)

Here, the covariance between arbitrary operators $\widehat{x}$ and $\widehat{y}$ is given in terms of the expectation value of the anticommutator, i.e., ${\sigma }_{xy}={\displaystyle \frac{1}{2}}\mathrm{Tr}\left(\widehat{\rho }(\widehat{x}\widehat{y}+\widehat{y}\widehat{x})\right)-\mathrm{Tr}\left(\widehat{\rho }\widehat{x}\right)\mathrm{Tr}\left(\widehat{\rho }\widehat{y}\right)$.

All properties of the Gaussian state can be obtained by the use of the covariance matrix and the mean values of the state. For example, the purity of the Gaussian state can be obtained by the determinant of its covariance matrix, that is,

$$\mathrm{Tr}{\widehat{\rho }}^2=\frac{1}{2\sqrt{det\sigma }}.$$

(10)

The unitary dynamics of the initial state of Equation (7) can be obtained by the integration of the propagators multiplied by the mixed-state density matrix:

$$\rho (x,x^{\prime },t)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }d{x}_{1}d{x}_2{G}^{*}(x_1,x,t)\rho (x_1,x_2,0)G(x_2,x^{\prime },t);$$

as a result, it provides a Gaussian state with the same purity as the original state ($\mathrm{Tr}{\widehat{\rho }}^2\left(t\right)=\mathrm{Tr}{\widehat{\rho }}^2\left(0\right)$), since unitary transforms do not change purity, which then can infer the determinant invariance of the covariance matrix $det\sigma \left(0\right)=det\sigma \left(t\right)$.

In a similar way, we can write the final state in an analogous way as the initial one, that is:

$$\rho (x,x^{\prime },t)=Nexp\left\{-a_1\left(t\right)x^2+a_{12}\left(t\right)x{x}^{\prime }-a_1^{*}\left(t\right){x^{\prime }}^2+b\left(t\right)x+b{\left(t\right)}^{*}{x}^{\prime }\right\},$$

(11)

where the Gaussian parameters are written in terms of the symplectic matrix associated with the invariants of Equation (5); thus, we arrive at the following expressions:

$$\begin{array}{ccc}\hfill a_1\left(t\right)& =& \frac{1}{2{\lambda }_4}\left(\frac{2a_1^{*}{\lambda }_4-i{\lambda }_1}{{\lambda }_1^2-2i(a_1-a_1^{*}){\lambda }_1{\lambda }_4+d{\lambda }_4^2}+i{\lambda }_5\right),\hfill \\ \hfill a_{12}\left(t\right)& =& \frac{a_{12}}{{\lambda }_1^2-2i(a_1-a_1^{*}){\lambda }_1{\lambda }_4+d{\lambda }_4^2},\hfill \\ \hfill b\left(t\right)& =& \frac{{\lambda }_1(b-i{\lambda }_3+(a_{12}-2a_1){\lambda }_6)+{\lambda }_4((2a_1^{*}-a_{12}){\lambda }_3+i(2a_1^{*}b+a_{12}{b}^{*}-\delta {\lambda }_6))}{{\lambda }_1^2-2i(a_1-a_1^{*}){\lambda }_1{\lambda }_4+d{\lambda }_4^2},\hfill \end{array}$$

(12)

with $d=4|a_1{|}^2-a_{12}^2$.

It is possible to obtain the differential equations, which these parameters must satisfy. This is done by taking the time derivative of the parameters, in view of Equation (2); after some algebra, we obtain:

$$\begin{array}{ccc}\hfill {\dot{a}}_1\left(t\right)& =& i(a_{12}^2\left(t\right)-4a_1^2\left(t\right)){\omega }_1\left(t\right)-4a_1\left(t\right){\omega }_2\left(t\right)+i{\omega }_3\left(t\right),\hfill \\ \hfill {\dot{a}}_{12}\left(t\right)& =& 4a_{12}\left(t\right)(i(a_1^{*}\left(t\right)-a_1\left(t\right)){\omega }_1\left(t\right)-{\omega }_2\left(t\right)),\hfill \\ \hfill \dot{b}\left(t\right)& =& (2a_1\left(t\right)-a_{12}\left(t\right)){\delta }_1\left(t\right)-i{\delta }_2\left(t\right)-2i{a}_{12}\left(t\right)b^{*}\left(t\right){\omega }_1\left(t\right)-2b\left(t\right)({\omega }_2\left(t\right)+2i{a}_1\left(t\right){\omega }_1\left(t\right)),\hfill \end{array}$$

(13)

and their corresponding complex conjugates. It is worth mentioning that these equations can be corroborated by the use of the von Neumann equation for $\rho (x,x^{\prime },t)$ given in Equation (11), namely $i{\displaystyle \frac{d\rho (x,x^{\prime },t)}{dt}}=\langle x|[\widehat{H},\widehat{\rho }]|x^{\prime }\rangle$. On the other hand, it is known that the covariance matrix of the system can be obtained, in view of the quantum solutions of Equation (2); as we have pointed out, this corresponds to the classical solutions (3) with ${\delta }_1={\delta }_2=0$, which can be also written in terms of the symplectic transformation $\Lambda$ of Equation (5), i.e., $\sigma \left(t\right)={\Lambda }^{-1}\sigma \left(0\right){\tilde{\Lambda }}^{-1}$. Then, the covariance matrix evolution $\sigma \left(t\right)$ can be obtained as:

$$\sigma \left(t\right)=\left(\begin{array}{cc}{\sigma }_{pp}\left(t\right)& {\sigma }_{pq}\left(t\right)\\ {\sigma }_{pq}\left(t\right)& {\sigma }_{qq}\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{\lambda }_5& -{\lambda }_2\\ -{\lambda }_4& {\lambda }_1\end{array}\right)\left(\begin{array}{cc}{\sigma }_{pp}\left(0\right)& {\sigma }_{pq}\left(0\right)\\ {\sigma }_{pq}\left(0\right)& {\sigma }_{qq}\left(0\right)\end{array}\right)\left(\begin{array}{cc}{\lambda }_5& -{\lambda }_4\\ -{\lambda }_2& {\lambda }_1\end{array}\right).$$

(14)

After differentiating each covariance ${\dot{\sigma }}_{pp}\left(t\right)$, ${\dot{\sigma }}_{qq}\left(t\right)$, and ${\dot{\sigma }}_{pq}\left(t\right)$, by using Equation (2), the inverse expression of Equation (14), the purity conservation condition ${\sigma }_{pp}\left(0\right){\sigma }_{qq}\left(0\right)-{\sigma }_{pq}^2\left(0\right)={\sigma }_{pp}\left(t\right){\sigma }_{qq}\left(t\right)-{\sigma }_{pq}^2\left(t\right)$, and the condition ${\lambda }_1{\lambda }_5-{\lambda }_2{\lambda }_4=1$, we arrive at the following differential equations for the covariances:

$$\begin{array}{c}{\dot{\sigma }}_{pp}\left(t\right)=-4({\omega }_2\left(t\right){\sigma }_{pp}\left(t\right)+{\omega }_3\left(t\right){\sigma }_{pq}\left(t\right)),\hfill \\ {\dot{\sigma }}_{qq}\left(t\right)=4({\omega }_2\left(t\right){\sigma }_{qq}\left(t\right)+{\omega }_1\left(t\right){\sigma }_{pq}\left(t\right)),\hfill \\ {\dot{\sigma }}_{pq}\left(t\right)=2({\omega }_1\left(t\right){\sigma }_{pp}\left(t\right)-{\omega }_3\left(t\right){\sigma }_{qq}\left(t\right)).\hfill \end{array}$$

(15)

One can also check that these differential equations imply that the derivative of the determinant of $\sigma \left(t\right)$ is equal to zero, i.e., ${\displaystyle \frac{d}{dt}}\left[{\sigma }_{p}p\left(t\right){\sigma }_{qq}\left(t\right)-{\sigma }_{pq}^2\right]=0$, which also implies that the purity of the state (11) is time invariant. It is noteworthy that the time-derivative expressions for the covariance matrix can be expressed as follows:

$$\dot{\sigma }\left(t\right)=2\left[\sigma \left(t\right)\mathbf{B}\left(t\right)\Sigma -\Sigma \mathbf{B}\left(t\right)\sigma \left(t\right)\right],$$

(16)

where the matrix $\mathbf{B}\left(t\right)$ contains the Hamiltonian coefficients, while $\Sigma$ is a symplectic matrix, i.e.,

$$\mathbf{B}\left(t\right)=\left(\begin{array}{cc}{\omega }_1\left(t\right)& {\omega }_2\left(t\right)\\ {\omega }_2\left(t\right)& {\omega }_3\left(t\right)\end{array}\right),\Sigma =\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right).$$

On the other hand, one can also check, using the inverse expression of Equation (5), that the mean values of $\widehat{p}$ and $\widehat{q}$ follow the classical Equation (3), i.e.,

$$\frac{d}{dt}\left(\begin{array}{c}\langle \widehat{p}\rangle \\ \langle \widehat{q}\rangle \end{array}\right)=2\left(\begin{array}{cc}-{\omega }_2& -{\omega }_3\\ {\omega }_1& {\omega }_2\end{array}\right)\left(\begin{array}{c}\langle \widehat{p}\rangle \\ \langle \widehat{q}\rangle \end{array}\right)+\left(\begin{array}{c}-{\delta }_1\\ {\delta }_2\end{array}\right).$$

(17)

All information regarding the evolution of the Gaussian state can then be obtained by solving the differential Equations (16) and (17). As an example, we can consider the evolution of a Gaussian state with the initial covariance matrix $\sigma \left(0\right)$ and mean values $\langle \widehat{p}\rangle \left(0\right)$ and $\langle \widehat{q}\rangle \left(0\right)$.

Example

As an example, we consider the following Hamiltonian:

$$\widehat{H}=\frac{1}{2}({\widehat{p}}^2+{\omega }^2{\widehat{q}}^2)+\frac{\nu }{2}(\widehat{p}\widehat{q}+\widehat{q}\widehat{p}).$$

(18)

In view of Equations (16) and (17), the matrix $\mathbf{B}$ can be identified as:

$$\mathbf{B}=\frac{1}{2}\left(\begin{array}{cc}1& \nu \\ \nu & {\omega }^2\end{array}\right),$$

and one can show that, in the case of constant frequencies ($\omega$, $\nu$), the evolution is determined by the same differential equations for all the covariances:

$${\stackrel{⃛}{\sigma }}_{pp}=-4({\omega }^2-{\nu }^2){\dot{\sigma }}_{pp},{\stackrel{⃛}{\sigma }}_{qq}=-4({\omega }^2-{\nu }^2){\dot{\sigma }}_{qq},{\stackrel{⃛}{\sigma }}_{pq}=-4({\omega }^2-{\nu }^2){\dot{\sigma }}_{pq},$$

which at ${\omega }^2>{\nu }^2$ describes an oscillating motion with parameter $f^2=4({\omega }^2-{\nu }^2)$. The solution to these equations, which satisfy the initial conditions for the derivatives and second derivatives at time $t=0$ implied by Equation (16), are:

$$\begin{array}{c}{\sigma }_{pp}\left(t\right)=-\frac{2\left(fsin\left(ft\right)\left(w^2{z}_0+\nu x_0\right)+cos\left(ft\right)\left({\omega }^4{y}_0-{\omega }^2(x_0-2\nu z_0)+2{\nu }^2{x}_0\right)-{\omega }^2\left({\omega }^2{y}_0+x_0+2\nu z_0\right)\right)}{f^2},\\ {\sigma }_{qq}\left(t\right)=\frac{2\left(-cos\left(ft\right)\left({\omega }^2(-y_0)+x_0+2\nu (\nu y_0+z_0)\right)+fsin\left(ft\right)(\nu y_0+z_0)+{\omega }^2{y}_0+x_0+2\nu z_0\right)}{f^2},\\ {\sigma }_{pq}\left(t\right)=\frac{2cos\left(ft\right)\left({\omega }^2(\nu y_0+2z_0)+\nu x_0\right)+fsin\left(ft\right)\left(x_0-{\omega }^2{y}_0\right)-2\nu \left({\omega }^2{y}_0+x_0+2\nu z_0\right)}{f^2},\end{array}$$

(19)

with $x_0={\sigma }_{pp}\left(0\right)$, $y_0={\sigma }_{qq}\left(0\right)$, and $z_0={\sigma }_{pq}\left(0\right)$. The solution for the classical equations of motion for the mean values reads:

$$\begin{array}{ccc}\hfill \langle \widehat{p}\rangle \left(t\right)& =& \langle \widehat{p}\rangle \left(0\right)cos\left(ft/2\right)-(2/f)\left[{\omega }^2\langle \widehat{q}\rangle \left(0\right)+\nu \langle \widehat{p}\rangle \left(0\right)\right]sin\left(ft/2\right),\hfill \\ \hfill \\ \hfill \langle \widehat{q}\rangle \left(t\right)& =& \langle \widehat{q}\rangle \left(0\right)cos\left(ft/2\right)+(2/f)\left[\nu \langle \widehat{q}\rangle \left(0\right)+\langle \widehat{p}\rangle \left(0\right)\right]sin\left(ft/2\right).\hfill \end{array}$$

(20)

With these solutions, one can characterize the state behavior. In Figure 1, we show the evolution of the mean values and covariance matrix given by Equations (19) and (20) for the Hamiltonian (18). Here, we observe the oscillating behavior of the system. We would like to remark that, using these solutions for the covariances and the correspondence between the density matrix parameters and the covariances:

$$\sigma \left(t\right)=\frac{1}{2(a_1\left(t\right)+a_1^{*}\left(t\right)-a_{12}\left(t\right))}\left(\begin{array}{cc}4a_1\left(t\right)a_1^{*}\left(t\right)-a_{12}^2\left(t\right)& i(a_1\left(t\right)-a_1^{*}\left(t\right))\\ i(a_1\left(t\right)-a_1^{*}\left(t\right))& 1\end{array}\right),$$

(21)

one can also obtain the solution for the nonlinear Equation (13).

2.2. Invariant States

After obtaining the differential equations determining the evolution of the covariance matrix (16), one can ask the question: Do invariant states exist under the evolution of the quadratic Hamiltonian? To answer this question, we should examine the properties of Equation (16). If we assume the condition $\dot{\sigma }\left(t\right)=0$, then one needs to obtain all the covariance matrices, which satisfy the condition $\sigma \left(t\right)\mathbf{B}\left(t\right)\Sigma -\Sigma \mathbf{B}\left(t\right)\sigma \left(t\right)=0$. By taking the vector $\tilde{\mathbf{v}}=({\sigma }_{pp},{\sigma }_{pq},{\sigma }_{qq})$, the equations for $\dot{\sigma }\left(t\right)=0$ can be written as follows:

$$\mathbf{Mv}=0,\mathrm{with}\mathbf{M}=\left(\begin{array}{ccc}-4{\omega }_2& -4{\omega }_3& 0\\ 2{\omega }_1& 0& -2{\omega }_3\\ 0& 4{\omega }_1& 4{\omega }_2\end{array}\right).$$

(22)

As the matrix $\mathbf{M}$ has rank $\mathcal{R}=2$, one can conclude that there is one nontrivial vector satisfying Equation (22). Exploring the null-space of $\mathbf{M}$, one can check that the vector:

$$\tilde{\mathbf{v}}=C\left(\frac{{\omega }_3}{{\omega }_1},-\frac{{\omega }_2}{{\omega }_1},1\right),$$

with C being a constant, is the solution to Equation (22). We infer that all the states with a covariance matrix, given by:

$$\sigma \left(0\right)=C\left(\begin{array}{cc}{\displaystyle \frac{{\omega }_3}{{\omega }_1}}& -{\displaystyle \frac{{\omega }_2}{{\omega }_1}}\\ -{\displaystyle \frac{{\omega }_2}{{\omega }_1}}& 1\end{array}\right),$$

(23)

have an invariant covariance matrix. Using the inverse expressions of Equation (9), one can obtain an explicit form of the covariance invariant density matrix function. The parameters of the density operator (11) read:

$$a_1=\frac{4C^2{\omega }_1{\omega }_3+{({\omega }_1+2iC{\omega }_2)}^2}{8C{\omega }_1^2},a_{12}=\frac{4S-1}{4C},$$

(24)

with $S=C^2({\omega }_3/{\omega }_1-{\omega }_2^2/{\omega }_1^2)$ being the determinant of the invariant covariance matrix. In the case of the Hamiltonian (18), we have $S=C^2({\omega }^2-{\nu }^2)$, ${\omega }_3/{\omega }_1={\omega }^2$, and ${\omega }_2/{\omega }_1=\nu$, which lead us to realize that any state with parameters, as in Equation (24) with $S>1/4$, is a bonafide quantum state, which is covariance invariant. For $C=1$, $\omega =2$, and $\nu =1$, we have that the states with:

$$a_1=\frac{13}{8}+\frac{i}{2},a_{12}=\frac{11}{4}$$

(25)

are covariance invariant.

The states that satisfy $\dot{\sigma }=0$ or equivalently have parameters according to (24) are the states that do not change its shape on the phase space $q,p$; also, their mean values move following the classical equations of motion. If we assume these types of states with initial mean values $\langle \widehat{p}\rangle \left(0\right)=\langle \widehat{q}\rangle \left(0\right)=0$, the resulting states will be invariant, i.e., they will not change any of their properties over time (for ${\delta }_1={\delta }_2=0$). An example of such states for the Hamiltonian (18) are the ones in Equation (11), with parameters given by (25) and $b\left(t\right)=b^{*}\left(t\right)=0$. We would like to express that, in the case of an invariant system with vanishing mean values, the initial energy will be different from zero as the initial covariances are also different from zero.

This parametric formalism for the evolution of Gaussian states and the definition of invariant states can be generalized to any multidimensional quadratic system as seen in the following section.

3. Multidimensional Quadratic System

In this section, we review the equations determining the evolution of the covariance matrix and mean values $\langle p_j\rangle$ and $\langle q_j\rangle$ for an arbitrary system under the evolution of a quadratic Hamiltonian; also, we mention the connection and dynamics of the continuous density matrix parameters. To obtain these properties, we use, as in the one-dimensional case, the invariant operators defined in [33,34].

In the case of an N-dimensional quadratic system, the time evolution is characterized by the Hamiltonian:

$$\widehat{H}=\tilde{\mathbf{r}}\mathbf{B}\left(t\right)\mathbf{r}+\tilde{\Delta }\left(t\right)\mathbf{r},$$

(26)

where the tilde corresponds to the transposition operation and the vector $\tilde{\mathbf{r}}=\left({\widehat{p}}_1,{\widehat{q}}_1,{\widehat{p}}_2,{\widehat{q}}_2,\dots ,{\widehat{p}}_N,{\widehat{q}}_N\right)$ corresponds to the vector of quadrature operators. The time dependence of this Hamiltonian is contained in the matrices:

$$\mathbf{B}\left(t\right)=\left(\begin{array}{cccc}{\omega }_{1,1}\left(t\right)& {\omega }_{1,2}\left(t\right)& \cdots & {\omega }_{1,2N}\left(t\right)\\ {\omega }_{1,2}\left(t\right)& {\omega }_{2,2}\left(t\right)& \cdots & {\omega }_{2,2N}\left(t\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\omega }_{1,2N}\left(t\right)& {\omega }_{2,2N}\left(t\right)& \cdots & {\omega }_{2N,2N}\left(t\right)\end{array}\right),\Delta \left(t\right)=\left(\begin{array}{c}{\delta }_1\left(t\right)\\ {\delta }_2\left(t\right)\\ ⋮\\ {\delta }_{2N}\left(t\right)\end{array}\right),$$

(27)

where $\mathbf{B}\left(t\right)$ is a real and symmetric matrix and $\Delta \left(t\right)$ is a real vector. As in the one-dimensional case, there exist $2N$ linear time-dependent operators ${\widehat{R}}_j$ ($j=1,\dots ,2N$), whose time derivatives are equal to zero $\left({\displaystyle \frac{d{\widehat{R}}_j}{dt}}=0.\right)$. These operators can be arranged on a vector as follows:

$$\mathbf{R}=\Lambda \left(t\right)\mathbf{r}+\Gamma \left(t\right),$$

(28)

with the matrix $\Lambda \left(t\right)$ and the vector $\tilde{\Gamma }\left(t\right)=\left({\gamma }_1\left(t\right),{\gamma }_2\left(t\right),\dots ,{\gamma }_{2N}\left(t\right)\right)$. By taking the time-derivative of the operator ${\mathbf{R}}_j={\widehat{R}}_j$ and equating it to zero $\left({\displaystyle \frac{d{\widehat{R}}_j}{dt}}=0.\right)$, one can demonstrate that $\Lambda \left(t\right)$ and $\Gamma \left(t\right)$ must satisfy the following differential equations:

$$\dot{\Lambda }\left(t\right)=2i\Lambda \left(t\right)\mathbf{D}\mathbf{B}\left(t\right),\dot{\Gamma }\left(t\right)=i\Lambda \left(t\right)\mathbf{D}\Delta \left(t\right),\mathrm{with}{\mathbf{D}}_{j,k}=[{\mathbf{r}}_j,{\mathbf{r}}_k].$$

(29)

The solution to these differential equations with the initial conditions ${\widehat{R}}_j\left(0\right)={\mathbf{r}}_j$ provides, as a result, the invariant operators $\tilde{\mathbf{R}}=\left({\widehat{P}}_1,{\widehat{Q}}_1,{\widehat{P}}_2,{\widehat{Q}}_2,\dots ,{\widehat{P}}_N,{\widehat{Q}}_N\right)$, satisfying the standard commutation rules for the operators $\mathbf{r}$ at time equal to zero: i.e., $[{\mathbf{R}}_j,{\mathbf{R}}_k]=[{\mathbf{r}}_j,{\mathbf{r}}_k]={\mathbf{D}}_{j,k}$. This property leads us to the conclusion that the matrix $\Lambda \left(t\right)$ must be symplectic and satisfies the equation:

$$\Lambda \left(t\right)\mathbf{D}\tilde{\Lambda }\left(t\right)=\mathbf{D}.$$

(30)

This relation can be then used to obtain the inverse of $\Lambda \left(t\right)$, which results in the expression:

$${\Lambda }^{-1}\left(t\right)=\mathbf{D}\tilde{\Lambda }\left(t\right)\mathbf{D}.$$

(31)

The other important property of these invariant operators is that they correspond to the inverse evolution of the original operators, in other words,

$${\mathbf{R}}_j=\widehat{U}\left(t\right){\mathbf{r}}_j{\widehat{U}}^{†}\left(t\right),$$

which in the most cases can be obtained from the Heisenberg picture operators by assuming a time reversal operation. This property implies that the entries of $\dot{\Lambda }\left(t\right)$ in Equation (29) satisfy the classical Hamilton equations after the time reversal operation, that is after the change $p_i\to -p_i$ in the classical Hamilton equations.

By the use of these invariant operators, one can obtain the time dependence of the mean values of the operators in $\mathbf{r}$ ($\langle \mathbf{r}\rangle \left(t\right)$) and their covariances ${\sigma }_{j,k}={\displaystyle \frac{1}{2}}\langle \{{\mathbf{r}}_j,{\mathbf{r}}_k\}\rangle -\langle {\mathbf{r}}_j\rangle \langle {\mathbf{r}}_k\rangle$. From the inverse of Equation (28), one can demonstrate that:

$$\langle \mathbf{r}\rangle \left(t\right)={\Lambda }^{-1}\left(t\right)(\langle \mathbf{R}\rangle -\Gamma \left(t\right)),$$

as the invariant operators in $\mathbf{R}$ have a time derivative equal to zero, and they are equal to the standard operators $\mathbf{r}$ at zero time, then one can conclude that:

$$\langle \mathbf{r}\rangle \left(t\right)={\Lambda }^{-1}\left(t\right)(\langle \mathbf{r}\rangle \left(0\right)-\Gamma \left(t\right)).$$

(32)

From an analogous argument, one can see that the covariance matrix reads:

$$\sigma \left(t\right)={\Lambda }^{-1}\left(t\right)\sigma \left(0\right){\tilde{\Lambda }}^{-1}\left(t\right).$$

(33)

Then, to obtain the expression for the time-derivative of the mean values $\langle \mathbf{r}\rangle \left(t\right)$ and the covariance matrix $\sigma \left(t\right)$, we make use of Equations (29) and (31)–(33) and arrive to the expressions:

$$\frac{d}{dt}\langle \mathbf{r}\rangle \left(t\right)=-i\mathbf{D}(2\mathbf{B}\left(t\right)\langle \mathbf{r}\rangle \left(t\right)+\Delta \left(t\right)),\frac{d}{dt}\sigma \left(t\right)=2i(\sigma \left(t\right)\mathbf{B}\left(t\right)\mathbf{D}-\mathbf{D}\mathbf{B}\left(t\right)\sigma \left(t\right)).$$

(34)

These differential equations, being first obtained in [31,32], are the generalization of the one-dimensional case discussed in the previous section. In our case, the nonlinear differential equations for the density matrix parameters can be obtained by explicit calculation of the covariances at time t. The resulting equations can then be solved by the substitution of the solution of Equation (34) or by direct integration.

To make the relation easier to see, we point out that the $2N\times 2N$ symplectic matrix $\mathbf{D}$ contains in its diagonal blocks made of the $2\times 2$ symplectic matrix $\Sigma$, that is,

$$\mathbf{D}=-i\left(\begin{array}{ccccc}\Sigma & 0& 0& \cdots & 0\\ 0& \Sigma & 0& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 0& \Sigma & 0\\ 0& 0& \cdots & 0& \Sigma \end{array}\right).$$

One can also notice that the differential equations for the entries of the covariance matrix are linear and can be expressed in the following matrix form:

$$\frac{d}{dt}\mathbf{v}=\mathbf{M}\mathbf{v},$$

(35)

where $\mathbf{v}$ is an $N(2N+1)$-dimensional vector, which contains all the independent covariances in the N-partite system, and the matrix $\mathbf{M}$ is a square matrix of the same dimension that contains the Hamiltonian coefficients.

4. Nonunitary Evolution for Gaussian Subsystems

Assume that the operators in the Hamiltonian of Equation (26) correspond to the ones of a multipartite system, where the position and momentum for the jth part are given by ${\widehat{p}}_j$ and ${\widehat{q}}_j$, respectively. Given this, one can see that the evolution of the complete system is unitary, but each one of its parts evolves in a nonunitary way due to the correlations between these parts. When the complete system is Gaussian, each one of its parts is also Gaussian. To show this property, lets assume that the N-partite system can be determined by the following density matrix at time $t=0$:

$$\langle \mathbf{x}|\widehat{\rho }\left(0\right)|{\mathbf{x}}^{\prime }\rangle =Nexp\left\{-\frac{1}{2}\tilde{\mathbf{y}}\mathbf{A}\mathbf{y}+\tilde{\mathbf{b}}\mathbf{y}\right\},$$

(36)

where $\tilde{\mathbf{x}}=(x_1,x_2,\dots ,x_N)$, ${\tilde{\mathbf{x}}}^{\prime }=(x_1^{\prime },x_2^{\prime },\dots ,x_N^{\prime })$, and $\tilde{\mathbf{y}}=(x_1,x_2,\dots ,x_N,x_1^{\prime },x_2^{\prime },\dots ,x_N^{\prime })$ are real vectors. Furthermore, we define the vector $\tilde{\mathbf{b}}=(b_1,b_2,\dots ,b_N,b_1^{*},b_2^{*},\dots ,b_N^{*})$ and the matrix:

$$\mathbf{A}=\left(\begin{array}{cc}\mathbf{u}& -\mathbf{v}\\ -\tilde{\mathbf{v}}& {\mathbf{u}}^{*}\end{array}\right),$$

where the block matrices $\mathbf{u}$ and $\mathbf{v}$ can be written as:

$$\mathbf{u}=\left(\begin{array}{cccc}2a_{1,1}& -a_{1,2}& \cdots & -a_{1,N}\\ -a_{1,2}& 2a_{2,2}& \cdots & -a_{2,N}\\ ⋮& ⋮& \ddots & ⋮\\ -a_{1,N}& -a_{2,N}& \cdots & 2a_{N,N}\end{array}\right),\mathbf{v}=\left(\begin{array}{ccccc}{a}_{1,N+1}& a_{1,N+2}& \cdots & a_{1,2N-1}& a_{1,2N}\\ a_{1,N+2}^{*}& a_{2,N+2}& \cdots & a_{2,2N-1}& a_{2,2N}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{1,2N-1}^{*}& a_{2,2N-1}^{*}& \cdots & a_{N-1,2N-1}& a_{N-1,2N}\\ a_{1,2N}^{*}& a_{2,2N}^{*}& \cdots & a_{N,2N-1}^{*}& a_{N,2N}\end{array}\right).$$

As stated earlier, the dynamics of the composite system is determined by the evolution of its covariance matrix and the mean values of the position and momentum operators, i.e., by the solution of Equation (34) with the initial state (36). The resulting state has the same purity as the initial state, since the evolution is unitary. However, there exists a nonunitary evolution of the parts, which compose the N-partite system.

To obtain the dynamic evolution of one of the parts, we can use the partial trace method. In other words, one should take the partial trace of all the subsystems in $\langle \mathbf{x}\left|\widehat{\rho }\left(t\right)\right|{\mathbf{x}}^{\prime }\rangle$, except the one we want to study. Nevertheless, as the system is Gaussian, the partial traces should give us also a Gaussian state for the density matrix under study.

As the most general one-dimensional Gaussian state can be obtained by the $2\times 2$ covariance matrix and the mean values ($\langle \widehat{p}\rangle \left(t\right)$ and $\langle \widehat{q}\rangle \left(t\right)$), we can obtain the result from the solutions to Equation (34) without the necessity of the partial trace operation.

On the other hand, once the time derivatives of these properties are established, one can derive the differential equation that the density matrix for the subsystem must satisfy. To show this procedure, we can take the bipartite system as an example.

4.1. Nonunitary Evolution on a Bipartite System

To exemplify the nonunitary evolution of a subsystem within a system, one can take a bipartite Gaussian state, which evolves on the Hamiltonian:

$$\begin{array}{cc}\hfill \widehat{H}\left(t\right)=\tilde{\mathbf{r}}\mathbf{B}\left(t\right)\mathbf{r}+\tilde{\Gamma }\left(t\right)\mathbf{r}& =({\widehat{p}}_1,{\widehat{q}}_1,{\widehat{p}}_2,{\widehat{q}}_2)\left(\begin{array}{cccc}{\omega }_{1,1}& {\omega }_{1,2}& {\omega }_{1,3}& {\omega }_{1,4}\\ {\omega }_{1,2}& {\omega }_{2,2}& {\omega }_{2,3}& {\omega }_{2,4}\\ {\omega }_{1,3}& {\omega }_{2,3}& {\omega }_{3,3}& {\omega }_{3,4}\\ {\omega }_{1,4}& {\omega }_{2,4}& {\omega }_{3,4}& {\omega }_{4,4}\end{array}\right)\left(\begin{array}{c}{\widehat{p}}_1\\ {\widehat{q}}_1\\ {\widehat{p}}_2\\ {\widehat{q}}_2\end{array}\right)\hfill \\ \hfill & +({\gamma }_1,{\gamma }_2,{\gamma }_3,{\gamma }_4)\left(\begin{array}{c}{\widehat{p}}_1\\ {\widehat{q}}_1\\ {\widehat{p}}_2\\ {\widehat{q}}_2\end{array}\right),\hfill \end{array}$$

(37)

where ${\omega }_{j,k}$ and ${\gamma }_j$ may be functions of time. In order to determine the time evolution of the system, one can solve the differential equations defined for the covariance matrix and the mean values of the position and momentum operators or, similarly to the one-dimensional case, one can solve the equations for the density matrix parameters given in Equation (A3) of Appendix A. The differential equations for the covariance matrix and mean values of the position and momentum operators for the subsystem can be obtained using Equation (34). To solve the time derivative equations, we express the matrices $\mathbf{B}\left(t\right)$, $\mathbf{D}$, and $\sigma \left(t\right)$ in the $2\times 2$ block representation; in such a case, we have:

$$\mathbf{B}\left(t\right)=\left(\begin{array}{cc}{\mathbf{B}}_1\left(t\right)& {\mathbf{B}}_{1,2}\left(t\right)\\ {\tilde{\mathbf{B}}}_{1,2}\left(t\right)& {\mathbf{B}}_2\left(t\right)\end{array}\right),\mathbf{D}=\left(\begin{array}{cc}-i\Sigma & 0\\ 0& -i\Sigma \end{array}\right),\sigma \left(t\right)=\left(\begin{array}{cc}{\sigma }_1\left(t\right)& {\sigma }_{1,2}\left(t\right)\\ {\tilde{\sigma }}_{1,2}\left(t\right)& {\sigma }_2\left(t\right)\end{array}\right),$$

where ${\sigma }_1\left(t\right)$ and ${\sigma }_2\left(t\right)$ are the covariance matrices for Subsystems 1 and 2, respectively, and ${\sigma }_{1,2}\left(t\right)$ is a matrix containing the covariances associated with the correlations between the two subsystems. The same can be said for the matrix linked to the Hamiltonian (37), i.e., $\mathbf{B}\left(t\right)$ where the block matrices ${\mathbf{B}}_1\left(t\right)$ and ${\mathbf{B}}_2\left(t\right)$ are associated with Subsystems 1 and 2, respectively, while ${\mathbf{B}}_{1,2}$ is associated with the interactions between these two subsystems.

After this identification, the expression for the covariance matrices of the subsystems and the correlations can be given as follows:

$$\begin{array}{c}\hfill {\dot{\sigma }}_1\left(t\right)=2\left(({\sigma }_1\left(t\right){\mathbf{B}}_1\left(t\right)+{\sigma }_{1,2}\left(t\right){\tilde{\mathbf{B}}}_{1,2})\Sigma -\Sigma ({\mathbf{B}}_1\left(t\right){\sigma }_1\left(t\right)+{\mathbf{B}}_{1,2}{\tilde{\sigma }}_{1,2}\left(t\right))\right),\\ \hfill {\dot{\sigma }}_2\left(t\right)=2\left(({\sigma }_2\left(t\right){\mathbf{B}}_2\left(t\right)+{\tilde{\sigma }}_{1,2}\left(t\right){\mathbf{B}}_{1,2})\Sigma -\Sigma ({\mathbf{B}}_2\left(t\right){\sigma }_2\left(t\right)+{\tilde{\mathbf{B}}}_{1,2}{\sigma }_{1,2}\left(t\right))\right),\\ \hfill {\dot{\sigma }}_{1,2}\left(t\right)=2\left(({\sigma }_1\left(t\right){\mathbf{B}}_{1,2}\left(t\right)+{\sigma }_{1,2}\left(t\right){\mathbf{B}}_2)\Sigma -\Sigma ({\mathbf{B}}_1\left(t\right){\sigma }_{1,2}\left(t\right)+{\mathbf{B}}_{1,2}{\sigma }_2\left(t\right))\right).\end{array}$$

(38)

Then, one can recognize the term $2({\sigma }_j{\mathbf{B}}_j\Sigma -\Sigma {\mathbf{B}}_j{\sigma }_j)$ for $j=1,2$, as the term corresponding to a unitary evolution of each subsystem (16). The extra term $2({\sigma }_{1,2}{\tilde{\mathbf{B}}}_{1,2}\Sigma -\Sigma {\mathbf{B}}_{1,2}{\tilde{\sigma }}_{1,2})$ is associated with the nonunitary evolution of the subsystems.

It is worth noting that these results are in accordance with the ones described by Sandulescu et al. [36] and Isar [37,38], where those results were obtained by solving the Gorini–Kossakowski–Sudarshan–Lindblad master equation [28,29,30] for two coupled oscillators. The main difference here is that our results were obtained exactly from the von Neumann equation without introducing a master equation.

4.2. Invariant and Quasi-Invariant States

The expression for the derivatives of the covariance matrix can lead to the definition of different Gaussian states, which do not evolve in the Hamiltonian dynamics. These types of states can be found as solutions to the equation $\dot{\sigma }=0$, which can be expressed in terms of the following equation regarding the covariance and the Hamiltonian matrices $\sigma \mathbf{B}\left(t\right)\mathbf{D}-\mathbf{D}\mathbf{B}\left(t\right)\sigma =0$. As discussed before, this system of differential equations can be replaced by $\dot{\mathbf{v}}=\mathbf{Mv}$ (35) with the following correspondences:

$$\tilde{\mathbf{v}}=({\sigma }_{p_1{p}_1},{\sigma }_{p_1{q}_1},{\sigma }_{p_1{p}_2},{\sigma }_{p_1{q}_2},{\sigma }_{q_1{q}_1},{\sigma }_{q_1{p}_2},{\sigma }_{q_1{q}_2},{\sigma }_{p_2{p}_2},{\sigma }_{p_2{q}_2},{\sigma }_{q_2{q}_2}),$$

and the matrix $\mathbf{M}$ containing the Hamiltonian coefficients is presented in Equation (A4) of Appendix B. It is possible to see that the matrix $\mathbf{M}$ has a determinant $det\mathbf{M}=0$ and a rank $\mathcal{R}=8$. From these properties, one can see that the system $\dot{\sigma }\left(t\right)=\dot{\mathbf{v}}=0$ has at most two different nontrivial solutions, which may be physical.

To exemplify the definition of bipartite states, which have a stationary behavior, we consider the frequency converter and the parametric amplifier. Both of these systems are quadratic and model the interaction between different electromagnetic fields in a nonlinear medium.

4.3. Frequency Converter

The quantum frequency converter is a device where two different unimodal electromagnetic fields, called the input and the output, interact with a semiclassical pump field on a nonlinear material. This interaction has the goal of interchanging the frequencies of the input and output beams at specific times. This behavior can be modeled using the following Hamiltonian:

$$\widehat{H}\left(t\right)=\hslash {\omega }_1\left({\widehat{a}}_1^{†}{\widehat{a}}_1+1/2\right)+\hslash {\omega }_2\left({\widehat{a}}_2^{†}{\widehat{a}}_2+1/2\right)-\hslash \kappa \left({\widehat{a}}_1^{†}{\widehat{a}}_2{e}^{-i\omega t}+{\widehat{a}}_1{\widehat{a}}_2^{†}{e}^{i\omega t}\right),$$

(39)

where the frequencies ${\omega }_{1,2}$ are the input and output frequencies, respectively, $\omega$ is the pump field frequency, and the bosonic operators ${\widehat{a}}_{1,2}$ are the annihilation operators of the input and output fields, respectively. These beams interact with an intensity $\kappa$ in a nonlinear medium as, e.g., a nonlinear crystal. In this case, the Hamiltonian matrix $\mathbf{B}\left(t\right)$ from (26) (in a unit system where $\hslash =1$) reads:

$$\mathbf{B}\left(t\right)=\frac{1}{2}\left(\begin{array}{cccc}1& 0& -{\displaystyle \frac{\kappa }{\sqrt{{\omega }_1{\omega }_2}}}cos\left(\omega t\right)& \kappa \sqrt{{\omega }_2/{\omega }_1}sin\left(\omega t\right)\\ 0& {\omega }_1^2& -\kappa \sqrt{{\omega }_1/{\omega }_2}sin\left(\omega t\right)& -\kappa \sqrt{{\omega }_1{\omega }_2}cos\left(\omega t\right)\\ -{\displaystyle \frac{\kappa }{\sqrt{{\omega }_1{\omega }_2}}}cos\left(\omega t\right)& -\kappa \sqrt{{\omega }_1/{\omega }_2}sin\left(\omega t\right)& 1& 0\\ \kappa \sqrt{{\omega }_2/{\omega }_1}sin\left(\omega t\right)& -\kappa \sqrt{{\omega }_1{\omega }_2}cos\left(\omega t\right)& 0& {\omega }_2^2\end{array}\right).$$

For this Hamiltonian, one can obtain different states that have dynamical equilibrium properties, i.e., states with a time derivative for the covariance matrix equal to zero. To characterize these types of states, one should solve Equation (34) with $\dot{\sigma }\left(t\right)=0$. As previously discussed, $\dot{\sigma }\left(t\right)=0$ can be expressed in vector form as:

$$\mathbf{Mv}=0,$$

(40)

where $\mathbf{v}$ is defined as:

$$\mathbf{v}=\left({\sigma }_{p_1{p}_1},{\sigma }_{p_1{q}_1},{\sigma }_{p_1{p}_2},{\sigma }_{p_1{q}_2},{\sigma }_{q_1{q}_1},{\sigma }_{q_1{p}_2},{\sigma }_{q_1{q}_2},{\sigma }_{p_2{p}_2},{\sigma }_{p_2{q}_2},{\sigma }_{q_2{q}_2}\right),$$

and the matrix $\mathbf{M}$ is a $10\times 10$ matrix of rank $\mathcal{R}=8$, which contains the Hamiltonian parameters of $\mathbf{B}\left(t\right)$. To solve Equation (40), one can explore the null space of matrix $\mathbf{M}$. The resulting null space contains two different vectors; one contains a non-physical solution. In this solution,

$$\begin{array}{ccc}\hfill {\sigma }_{p_1{p}_1}={\left({\omega }_1^3{\omega }_2\right)}^{1/2}({\omega }_2-{\omega }_1)sec\left(\omega t\right)/\kappa ,& & {\sigma }_{p_1{p}_2}={\omega }_1{\omega }_2,\hfill \\ \hfill {\sigma }_{p_1{q}_2}=-{\omega }_{1}tan\left(\omega t\right),& & {\sigma }_{q_1{q}_1}={\omega }_2^{1/2}({\omega }_2-{\omega }_1)sec\left(\omega t\right)/\left(\kappa {\omega }_1^{1/2}\right),\hfill \\ \hfill {\sigma }_{q_1{p}_2}={\omega }_{2}tan\left(\omega t\right),& & {\sigma }_{q_1{q}_2}=1,\hfill \end{array}$$

while all the other covariances are equal to zero. In particular, it contains the nonphysical terms ${\sigma }_{p_2{p}_2}={\sigma }_{q_2{q}_2}=0$ that resemble the case where one of the subsystem is classical, as in a classical system, the values of the covariances can be equal to zero. The null space also contains another vector, which has the following physical values:

$${\sigma }_{p_1{p}_1}={\omega }_1{\omega }_2,{\sigma }_{q_1{q}_1}={\omega }_2/{\omega }_1,{\sigma }_{p_2{p}_2}={\omega }_2^2,{\sigma }_{q_2{q}_2}=1,$$

(41)

with the remaining covariances equal to zero.

These results led us to the conclusion that a two-mode Gaussian state with initial covariances proportional to the ones established in (41) has the same covariances for any time $t>0$ (for time-independent parameters $\omega$, ${\omega }_{1,2}$, and $\kappa$). This property has several physical implications such as, for example, that the purity of the subsystems will always be the same regardless of the interaction between them and despite the interchange of their frequencies. The resulting states will only have different mean values of the quadrature components ($p_1$, $q_1$, $p_2$, and $q_2$), which evolve according to the classical Hamilton equations.

In the case where the mean values of the quadrature components $b_j$; $j=1,2$ are equal to zero in Equation (36), one can obtain different states, which do not change their properties over time. In this case, the entanglement of the system (which can be obtained by the logarithmic negativity of the covariance matrix) will also be an invariant of the system. These properties make the evolution of these types of states relevant to quantum computing and quantum information.

By the use of the inverse relations of Equations of the Appendix A: (A1) and (A2), one can then write a general state with an invariant covariance matrix, which only changes its mean values according to the classical motion equations. Such a state can be expressed as the one in (36), after making the identification:

$$\begin{array}{c}\hfill \mathbf{A}=\left(\begin{array}{cccc}{\displaystyle \frac{{\omega }_1}{4C{\omega }_2}}+C{\omega }_1{\omega }_2& 0& {\displaystyle \frac{{\omega }_1}{4C{\omega }_2}}-C{\omega }_1{\omega }_2& 0\\ 0& {\displaystyle \frac{1}{4C}}+C{\omega }_2^2& 0& {\displaystyle \frac{1}{4C}}-C{\omega }_2^2\\ {\displaystyle \frac{{\omega }_1}{4C{\omega }_2}}-C{\omega }_1{\omega }_2& 0& {\displaystyle \frac{{\omega }_1}{4C{\omega }_2}}+C{\omega }_1{\omega }_2\\ 0& {\displaystyle \frac{1}{4C}}-C{\omega }_2^2& 0& {\displaystyle \frac{1}{4C}}+C{\omega }_2^2\end{array}\right),\end{array}$$

(42)

with C being a constant, which needs to be chosen in order for the covariance matrix to be positive; in particular, to fulfill the Schrödinger–Robertson inequalities ${\sigma }_{p_i{p}_i}{\sigma }_{q_i{q}_i}-{\sigma }_{p_i{q}_i}^2\ge 1/4$ ($i=1,2$) and $det\sigma \ge 1/16$.

4.4. Parametric Amplifier

The other Hamiltonian, which can be taken as an example, is the nondegenerate parametric amplifier. This system also describes the interaction of the input and output beams with the pump field in a nonlinear medium. As a result of this interaction, one can obtain the amplification of the input beam. The Hamiltonian associated with the parametric amplifier reads:

$$\widehat{H}=\hslash {\omega }_1\left({\widehat{a}}_1^{†}{\widehat{a}}_1+1/2\right)+\hslash {\omega }_2\left({\widehat{a}}_2^{†}{\widehat{a}}_2+1/2\right)-\hslash \kappa \left({\widehat{a}}_1^{†}{\widehat{a}}_2^{†}{e}^{-i\omega t}+{\widehat{a}}_1{\widehat{a}}_2{e}^{i\omega t}\right),$$

(43)

where the frequencies ${\omega }_{1,2}$ are the frequencies of the input and output beam channels and $\omega$ is the frequency of a pump field, which allows the amplification of the input channel. Then, the Hamiltonian matrix $\mathbf{B}\left(t\right)$ is:

$$\mathbf{B}\left(t\right)=\frac{1}{2}\left(\begin{array}{cccc}1& 0& \kappa {\displaystyle \frac{cos\left(\omega t\right)}{\sqrt{{\omega }_1{\omega }_2}}}& \kappa \sqrt{{\omega }_2/{\omega }_1}sin\left(\omega t\right)\\ 0& {\omega }_1^2& \kappa \sqrt{{\omega }_1/{\omega }_2}sin\left(\omega t\right)& -\kappa \sqrt{{\omega }_1{\omega }_2}cos\left(\omega t\right)\\ \kappa {\displaystyle \frac{cos\left(\omega t\right)}{\sqrt{{\omega }_1{\omega }_2}}}& \kappa \sqrt{{\omega }_1/{\omega }_2}sin\left(\omega t\right)& 1& 0\\ \kappa \sqrt{{\omega }_2/{\omega }_1}sin\left(\omega t\right)& -\kappa \sqrt{{\omega }_1{\omega }_2}cos\left(\omega t\right)& 0& {\omega }_2^2\end{array}\right).$$

(44)

Following an analogous procedure to obtain the solutions of the equation $\dot{\sigma }=0$, one can show that the null space of the corresponding matrix $\mathbf{M}$ for this problem can lead us to nonphysical values for different covariances on the system. One of the vectors of the null space for the case ${\omega }_{1,2}>0$ can be written as:

$$\begin{array}{ccc}\hfill {\sigma }_{p_1{p}_1}& =& \frac{{\omega }_1\sqrt{{\omega }_1{\omega }_2}({\omega }_1+{\omega }_2)sec\left(\omega t\right)}{\kappa },{\sigma }_{p_1{p}_2}=-{\omega }_1{\omega }_2,{\sigma }_{p_1{q}_2}=-{\omega }_{1}tan\left(\omega t\right),\hfill \\ \hfill {\sigma }_{q_1{q}_1}& =& \frac{{\omega }_2({\omega }_1+{\omega }_2)sec\left(\omega t\right)}{\kappa \sqrt{{\omega }_1{\omega }_2}},{\sigma }_{p_2{q}_1}=-{\omega }_{2}tan\left(\omega t\right),{\sigma }_{q_1{q}_2}=1,\hfill \end{array}$$

while all the other covariances are equal to zero. The other vector on the null space is:

$${\sigma }_{p_1{p}_1}=-{\omega }_1{\omega }_2,{\sigma }_{q_1{q}_1}=-\frac{{\omega }_2}{{\omega }_1},{\sigma }_{p_2{p}_2}={\omega }_2^2,{\sigma }_{q_2{q}_2}=1.$$

As the condition ${\omega }_{1,2}>0$ was used to obtain these results, then both vectors lead to nonphysical covariances. Nevertheless, one can obtain states with a slow change ratio of the covariances compared with the change of the mean value of the system Hamiltonian $\langle \widehat{H}\rangle \left(t\right)$. These type of states can be defined by considering the initial covariances equal to $C=1/\left({\omega }_1{\omega }_2\right)$ times the ones presented in Equation (41); in other words,

$${\sigma }_{p_1,p_1}=1,{\sigma }_{q_1,q_1}=1/{\omega }_1^2,{\sigma }_{p_2,p_2}={\omega }_2/{\omega }_1,{\sigma }_{q_2,q_2}=1/\left({\omega }_1{\omega }_2\right).$$

The slow time dependence behavior of these covariances can be seen in Figure 2, where the time dependence of the covariances and the purity of the subsystems are illustrated. The evolution of the subsystems in the parametric amplifier normally varies very fast, as the photons from the pump field are transformed into the photons of both subsystems. Nevertheless, it can be seen in Figure 2 that the variation of the majority of the covariances is not as fast compared with the change of $\langle \widehat{H}\rangle \left(t\right)$, providing the strong coupling between the subsystems. In this particular example, one can see that ${\sigma }_{p_2{p}_2}\left(t\right)={\sigma }_{q_2{q}_2}\left(t\right)$ and ${\sigma }_{p_1{q}_1}\left(t\right)={\sigma }_{p_2{q}_2}\left(t\right)=0$.

The detection of these type of states can be done by the use of quantum tomography, as is discussed in the next section.

5. Gaussian States and Their Evolution in the Tomographic-Probability Representation

There exist different representations of quantum states [39,40,41,42,43,44], and among them, the probability tomographic representation is of particular interest. In this representation, e.g., one-mode photon states are identified with symplectic tomograms [45], which correspond to the conditional probability distribution $w(X\mid \mu ,\nu )$ of the photon quadrature $-\infty <X<\infty$, to be measured in a reference frame with parameters $\mu =scos\theta$ and $\nu =s^{-1}sin\theta$. Here, $-\infty <\mu ,\nu <\infty$, s is a time scaling parameter, and $\theta$ is the local oscillator phase, which is used in experiments [46] to obtain the Wigner function of the photon state.

The symplectic tomogram $w_{\rho }(X\mid \mu ,\nu )$ is determined by the photon density operator $\widehat{\rho }$ [45] as:

$$w_{\rho }(X\mid \mu ,\nu )=Tr\left[\widehat{\rho }\delta (X\widehat{1}-\mu \widehat{q}-\nu \widehat{p})\right],$$

(45)

where $\widehat{q}$ and $\widehat{p}$ are quadrature components—the position and momentum operators within the framework of the oscillator model of the one-mode electromagnetic-field photons. The symplectic tomogram satisfies the normalization condition:

$$\int w_{\rho }(X\mid \mu ,\nu )dX=1,$$

(46)

and inversely, it determines the density operator $\widehat{\rho }$ of the photon state:

$$\widehat{\rho }=\frac{1}{2\pi }\int w(X\mid \mu ,\nu )exp\left[i\delta (X\widehat{1}-\mu \widehat{q}-\nu \widehat{p})\right]dXd\mu d\nu .$$

(47)

The optical tomogram of the photon state $w_{\mathrm{opt}}(X\mid \theta )\equiv w(X\mid \mu =cos\theta ,\nu =sin\theta )$ is measured in experiments, and in view of the homogeneity property of the Dirac delta-function, the measured optical tomogram of the photon state determines the symplectic tomogram:

$$w(X\mid \mu ,\nu )=\frac{1}{\sqrt{{\mu }^2+{\nu }^2}}{w}_{\mathrm{opt}}\left[\frac{X}{\sqrt{{\mu }^2+{\nu }^2}}|arctan\frac{\nu }{\mu }\right].$$

(48)

For Gaussian states (7), the tomographic-probability distribution of random photon quadrature X has the conventional form of a normal distribution:

$$w(X\mid \mu ,\nu )=\frac{1}{\sqrt{2\pi \sigma (\mu ,\nu )}}exp\left[-\frac{{\left(X-\overline{X}(\mu ,\nu )\right)}^2}{2\sigma (\mu ,\nu )}\right].$$

(49)

In view of (45), one has the mean value of the quadrature:

$$\overline{X}(\mu ,\nu )=\mu \langle \widehat{q}\rangle +\nu \langle \widehat{p}\rangle$$

(50)

and the covariance of the quadrature $\sigma (\mu ,\nu )$ reads:

$$\sigma (\mu ,\nu )={\mu }^2{\sigma }_{qq}+{\nu }^2{\sigma }_{pp}+2\mu \nu {\sigma }_{pq}.$$

(51)

For the measured optical tomogram, the dispersion $\sigma \left(\theta \right)$ is:

$$\sigma \left(\theta \right)=\left({cos}^2\theta \right){\sigma }_{qq}+\left({sin}^2\theta \right){\sigma }_{pp}+(sin2\theta ){\sigma }_{pq}.$$

(52)

In the quantum system with Hamiltonian (18), the tomographic quadrature dispersion (51) evolves according to the formula:

$$\sigma (\mu ,\nu ,t)={\mu }^2{\sigma }_{qq}\left(t\right)+{\nu }^2{\sigma }_{pp}\left(t\right)+2\mu \nu {\sigma }_{pq}\left(t\right),$$

(53)

where ${\sigma }_{qq}\left(t\right)$, ${\sigma }_{pp}\left(t\right)$, and ${\sigma }_{pq}\left(t\right)$ are provided by the explicit expressions (19) and the parameters $\langle \widehat{p}\left(t\right)\rangle$ and $\langle \widehat{q}\left(t\right)\rangle$ are given by (20). Thus, the properties of the Gaussian states of the oscillator with time-dependent parameters described by the covariances of the position and momentum and their mean values can be checked by considering the covariance of the homodyne quadrature X, as well as the mean value evolution.

The properties of the invariant Gaussian states can be visualized by the properties of either the Wigner function or the tomographic-probability distributions. There are examples of the time-dependent Gaussian-packet solutions of the kinetic equation for the symplectic tomogram [47,48] in the case of the harmonic oscillator Hamiltonian (18) with $\nu =0$. Since the dispersion matrix for the quadrature X is the linear combination of covariances ${\sigma }_{qq}\left(t\right)$, ${\sigma }_{pp}\left(t\right)$, and ${\sigma }_{pq}\left(t\right)$, which in the case of invariant Gaussian states, do not depend on time, the state tomogram also does not depend on time. The invariant states with density operators $\mid E_n\rangle \langle E_n\mid$ have the oscillator tomograms obtained from energy states $\mid E_n\rangle$, where $\widehat{H}\mid E_n\rangle =E_n\mid E_n\rangle$. Tomograms of invariant Gaussian states satisfy the equality:

$${\mathcal{P}}_G^n=\frac{1}{2\pi }\int w_G(X\mid \mu ,\nu )w_{E_n}(Y\mid \mu ,\nu )e^{i(X+Y)}dXdYd\mu d\nu ,$$

(54)

where the parameter ${\mathcal{P}}_G^n$ is the probability to detect the properties of the stationary state $\mid E_n\rangle$ with energy value $E_n$ in the Gaussian state with the time-dependent tomogram $w_G(X\mid \mu ,\nu )$. This state also does not depend on time.

Any convex sum of states $\mid E_n\rangle \langle E_n\mid$ is a density operator. One can conjecture that there is the decomposition of normal distribution $w_G(X\mid \mu ,\nu )$ corresponding, e.g., to a thermal state with $\widehat{\rho }=exp(-\widehat{H}/\left(kT\right))/\mathrm{Tr}(exp(-\widehat{H}/\left(kT\right)))$ (T being the temperature and k the Boltzmann constant), which can be presented as:

$$w_G(X\mid \mu ,\nu )=\sum _n{\mathcal{P}}_G^n{w}_{E_n}(X\mid \mu ,\nu ),\sum _n{\mathcal{P}}_G^n=1.$$

(55)

An analogous relation can be then written also for the Wigner function of the invariant Gaussian state of the oscillator, as well as for the density matrix in the position representation.

Now, we consider a harmonic oscillator with the Hamiltonian $\widehat{H}={\displaystyle \frac{{\widehat{p}}^2}{2}}+{\displaystyle \frac{{\widehat{q}}^2}{2}}$. The density matrix of the thermal equilibrium state with temperature $T={\beta }^{-1}$ in the position representation reads (here, we assume $\hslash =\omega =m=k=1$):

$$\rho (x,x^{\prime },\beta )={\left[{\pi }^{-1}{tan}^2(\beta /2)\right]}^{1/2}exp\left[\frac{x{x}^{\prime }}{sinh\beta }-\frac{x^2+{x^{\prime }}^2}{2}coth\beta \right].$$

(56)

The Green function of the oscillator has the Gaussian form:

$$G(x,x^{\prime },t)=\langle x\mid e^{-it\widehat{H}}\mid x^{\prime }\rangle =\frac{1}{\sqrt{2\pi isint}}exp\left[\frac{i(x^2+{x^{\prime }}^2)}{2}cott-\frac{ix{x}^{\prime }}{sint}\right].$$

(57)

Since the density matrix (56) is determined by the Green function (57), i.e.,

$$\rho (x,x^{\prime },\beta )=\frac{G(x,x^{\prime },-i\beta )}{Z\left(\beta \right)},$$

(58)

with the partition function $Z\left(\beta \right)$ given by the formula:

$$Z\left(\beta \right)=\sum _{n=0}^{\infty }\mathrm{Tr}\left[exp\left(-\beta \widehat{H}\right)\mid n\rangle \langle n\mid \right]=\frac{1}{2sinh(\beta /2)};$$

(59)

here, we use the property $\widehat{H}\mid n\rangle =\left(n+1/2\right)\mid n\rangle$. The density matrix (56) does not depend on time; this means that in all other representations, as the Wigner function or tomographic-probability representation, it is time-invariant. The density matrices of Fock states $\mid n\rangle \langle n\mid$ do not depend on time.

In the position representation, the density matrix of Fock state $\mid n\rangle \langle n\mid$ reads:

$$\langle x\mid n\rangle \langle n\mid x^{\prime }\rangle =\frac{H_n\left(x\right)H_n\left(x^{\prime }\right)}{2^{n}n!\sqrt{\pi }}exp\left(-\frac{x^2}{2}-\frac{x^{\prime 2}}{2}\right),$$

(60)

and it does not depend on time.

The density matrix (56), being described by an invariant Gaussian function, is the convex sum of the Fock state density matrices. One has the relation:

$$\rho (x,x^{\prime },\beta )=\frac{1}{\sqrt{\pi }}exp\left(-\frac{x^2}{2}-\frac{x^{\prime 2}}{2}\right)\sum _{n=0}^{\infty }\frac{e^{-(n+1/2)\beta }}{Z\left(\beta \right)2^{n}n!}{H}_n\left(x\right)H_n\left(x^{\prime }\right),$$

(61)

where $Z\left(\beta \right)$ is given in (59).

In the tomographic-probability representation of the thermal equilibrium oscillator and Fock states, we can calculate the tomograms in explicit form. For Fock states,

$$w_n(X\mid \mu ,\nu ,\beta )={\left[\pi ({\mu }^2+{\nu }^2)\right]}^{-1/2}\frac{1}{2^{n}n!}exp\left(-\frac{X^2}{{\mu }^2+{\nu }^2}\right)H_n^2\left(\frac{X}{\sqrt{{\mu }^2+{\nu }^2}}\right).$$

(62)

With all these properties, one can check that the thermal equilibrium Gaussian state of Equation (58) has a symplectic tomogram in the form of the normal distribution:

$$w(X\mid \mu ,\nu ,\beta )=\frac{1}{\sqrt{2\pi \sigma (\mu ,\nu )}}exp\left(-\frac{X^2}{2\sigma (\mu ,\nu )}\right).$$

(63)

The dispersion of quadrature $\langle X^2\rangle =\sigma (\mu ,\nu )$ given by Equation (63) reads:

$$\sigma (\mu ,\nu )={\mu }^2\langle {\widehat{q}}^2\rangle +{\nu }^2\langle {\widehat{p}}^2\rangle ,$$

(64)

where the state with density matrix (56):

$$\langle {\widehat{q}}^2\rangle =\langle {\widehat{p}}^2\rangle =\frac{1}{2}{coth}^2(\beta /2).$$

(65)

Thus, the symplectic tomogram (63) is given by an invariant normal probability distribution:

$$w(X\mid \mu ,\nu ,\beta )=\frac{coth(\beta /2)}{\sqrt{\pi ({\mu }^2+{\nu }^2)}}exp\left(-\frac{X^2}{{\mu }^2+{\nu }^2}{coth}^2(\beta /2)\right).$$

(66)

Since for optical tomogram $\mu =cos\theta$, $\nu =sin\theta$, and ${\mu }^2+{\nu }^2=1$, in the case of the thermal single-mode photon state, its optical tomogram is:

$$w(X\mid \theta ,\beta )=\frac{coth(\beta /2)}{\sqrt{\pi }}exp\left(-X^2{coth}^2(\beta /2)\right);$$

(67)

it depends neither on the local oscillator phase, nor on time. These types of states are Gaussian and time-independent, so there is a connection between them and the invariant states discussed above. This connection can be checked by equaling the covariance matrices in both cases, which can be also checked experimentally, for example, by preparing an initial Gaussian state according to the invariance condition $\dot{\sigma }=0$. As seen in previous examples, this condition implies a value for the initial covariances in terms of the Hamiltonian parameters. Then, using the tomographic representation discussed here, the relation of these states with thermal states can be corroborated. As the result of this comparison, one can obtain certain thermodynamic properties such as the temperature of the system. This can also be extended for the bipartite harmonic oscillator, as seen in the next section.

6. Two-Mode Gaussian States in the Tomographic-Probability Representation

For a two-mode harmonic oscillator, the Gaussian-state tomogram is determined by the normal probability distribution of quadratures $X_1$ and $X_2$; it is expressed in terms of the state density operator $\widehat{\rho }(1,2)$ as follows:

$$w(X_1,X_2\mid {\mu }_1,{\nu }_1,{\mu }_2,{\nu }_2)=\mathrm{Tr}\left[\widehat{\rho }(1,2)\delta (X_1-{\mu }_1{\widehat{q}}_1-{\nu }_1{\widehat{p}}_1)\delta (X_2-{\mu }_2{\widehat{q}}_2-{\nu }_2{\widehat{p}}_2)\right];$$

(68)

in the case of $\langle {\widehat{q}}_1\rangle =\langle {\widehat{q}}_2\rangle =\langle {\widehat{p}}_1\rangle =\langle {\widehat{p}}_2\rangle =0$, it reads:

$$w(X_1,X_2\mid {\mu }_1,{\nu }_1,{\mu }_2,{\nu }_2)=\frac{1}{2\pi \sqrt{det\sigma ({\mu }_1,{\mu }_2,{\nu }_1,{\nu }_2)}}exp\left[-\frac{1}{2}\left(\tilde{\mathbf{X}}{\mathit{\sigma }}^{-1}({\mu }_1,{\nu }_1,{\mu }_2,{\nu }_2)\mathbf{X}\right)\right].$$

(69)

Here, $\tilde{\mathbf{X}}=(X_1,X_2)$ and $\mathit{\sigma }({\mu }_1,{\nu }_1,{\mu }_2,{\nu }_2)=\left(\begin{array}{cc}\langle X_1^2\rangle & \langle X_1{X}_2\rangle \\ \langle X_1{X}_2\rangle & \langle X_2^2\rangle \end{array}\right),$ with:

$$\begin{array}{c}\langle X_1^2\rangle ={\mu }_1^2\langle {\widehat{q}}_1^2\rangle +{\nu }_1^2\langle {\widehat{p}}_1^2\rangle +{\mu }_1{\nu }_1(\langle {\widehat{q}}_1{\widehat{p}}_1\rangle +\langle {\widehat{p}}_1{\widehat{q}}_1\rangle ),\langle X_2^2\rangle ={\mu }_2^2\langle {\widehat{q}}_2^2\rangle +{\nu }_2^2\langle {\widehat{p}}_2^2\rangle +{\mu }_2{\nu }_2(\langle {\widehat{q}}_2{\widehat{p}}_2\rangle +\langle {\widehat{p}}_2{\widehat{q}}_2\rangle ),\hfill \\ \langle X_1{X}_2\rangle ={\mu }_1{\mu }_2\langle {\widehat{q}}_1{\widehat{q}}_2\rangle +{\nu }_1{\nu }_2\langle {\widehat{p}}_1{\widehat{p}}_2\rangle +{\mu }_1{\nu }_2\langle {\widehat{q}}_1{\widehat{p}}_1\rangle +{\mu }_2{\nu }_1\langle {\widehat{q}}_2{\widehat{p}}_1\rangle .\hfill \end{array}$$

The inverse transform provides the density operator $\widehat{\rho }(1,2)$ expressed in terms of the tomographic-probability distribution:

$$\begin{array}{lll}\widehat{\rho }(1,2)& =& \frac{1}{4{\pi }^2}\int w\left(X_1,X_2\mid {\mu }_1,{\nu }_1,{\mu }_2,{\nu }_2\right)\\ & & \times exp\left[i\left(X_1+X_2-{\mu }_1{\widehat{q}}_1-{\nu }_1{\widehat{p}}_1-{\mu }_2{\widehat{q}}_2-{\nu }_2{\widehat{p}}_2\right)\right]d{X}_{1}d{X}_{2}d{\mu }_{1}d{\mu }_{2}d{\nu }_{1}d{\nu }_2.\end{array}$$

(70)

The subsystem tomogram given by the partial trace of the density operator $\widehat{\rho }\left(1\right)={\mathrm{Tr}}_2\widehat{\rho }(1,2)$ reads:

$$w_1(X_1\mid {\mu }_1,{\nu }_1)=\int w\left(X_1,X_2\mid {\mu }_1,{\nu }_1,{\mu }_2,{\nu }_2\right)d{X}_2;$$

(71)

it is also given by the normal distribution discussed in the previous section.

If the tomogram of the two-mode oscillator state corresponds to the solution of the time evolution equation with a quadratic Hamiltonian, the unitary evolution of the system can induce the nonunitary evolution of tomogram (71). These evolutions can be used to obtain the shape of the invariant states, which we discussed above using the matrix $\mathbf{M}$ shown in Appendix B.

7. Summary and Conclusions

A differential formalism to obtain the time evolution of a multidimensional, multipartite Gaussian state was defined and studied. This new formalism used the time derivative of the parameters of the continuous variable density matrix of the system. The general procedure to obtain the time evolution can be summarized as follows: using the derivative of the covariance matrix for the Gaussian state and the expressions for the covariances in terms of the parameters of the density matrix, the differential equations for the parameters of the density function of the system were obtained. The resulting nonlinear differential equations could be used to obtain new physical information of the state instead of the use of the Schrödinger equation.

This differential formalism can also be used to describe exactly the nonunitary evolution of the subsystems of a composite Gaussian state. As an example, we considered a two-mode Gaussian state and demonstrated that the resulting derivatives of the covariance matrices for the subsystems contained unitary and nonunitary terms.

This study also allowed us to define invariant states, i.e., states that do not change their properties over time. To show this, we considered unimodal and bipartite Gaussian systems with density matrices in the position representation and the corresponding tomographic-probability representation. As explicit examples, we presented the invariant states for the one-dimensional quadratic Hamiltonian and the invariant states for the two-mode frequency converter and mentioned the applicability of these type of states in quantum information and computing. Furthermore, quasi-invariant states for the two-mode parametric amplifier were presented. We pointed out that the discussed examples of studying parametric systems could be used to apply the results associated with the behavior of physical systems like photons in cavities with time-dependent locations of boundaries to the dynamical Casimir effect (see [49]) and its analog in superconducting circuits [50,51]. One can discuss the nonunitary evolution of systems, which have no subsystems, using hidden correlations [52], which are present in noncomposite systems.