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

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 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.


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 The Gaussian states can be determined by their covariance matrix σ and mean values q j and p j . 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.

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 operatorsp andq.
with δ 1 = δ 2 = 0. To show this, one can see that the differential equations for λ 1,2 correspond to the classical equations with the substitution λ 1 → q and λ 2 → −p; in other words, they correspond to the time inversion of the classical equations. In the case of the differential equation for λ 3 , one can show, in view of the Hamilton equations, that it corresponds to the classical Lagrangian (with δ 1,2 = 0) plus the time variation of the function pq of the system, that is, where L = pq − 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 , t) = x|Û(t)|x using each one of the solutions of the classical problem [35].

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: λ 1 (0) = 1, λ 2 (0) = 0, λ 3 (0) = 0 and λ 1 (0) = 0, λ 2 (0) = 1, λ 3 (0) = 0. These two sets define two different invariants calledP andQ, respectively, which can be written as: where λ 4,5,6 satisfy the same differential equations as λ 1,2,3 , respectively, with the different sets of initial conditions mentioned above. The operatorsP andQ fulfill the commutation relation [Q,P] = i, implying the relation λ 1 λ 5 − λ 2 λ 4 = 1 and the fact that the matrix Λ is symplectic. Furthermore, they can be related to operatorsp andq through the evolution operator as follows: it is not difficult to show [33,34] that these expressions can be used to obtain the propagator of the system G(x, x , t) = x|Û(t)|x , which reads: 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 ρ(x, x , t = 0) = x|ρ|x equal to: with complex parameters a 1 = a 1R + ia 1I and b 1 = b 1R + ib 1I and the real parameter a 12 ∈ R, which also satisfy the integrability conditions a 1R > a 12 /2 ≥ 0 and have a normalization constant N expressed as: .
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 operatorsp andq are: and the initial covariance matrix of the system reads: Here, the covariance between arbitrary operatorsx andŷ is given in terms of the expectation value of the anticommutator, i.e., σ xy = 1 2 Tr (ρ(xŷ +ŷx)) − Tr (ρx) Tr (ρŷ). 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, 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: as a result, it provides a Gaussian state with the same purity as the original state (Trρ 2 (t) = Trρ 2 (0)), since unitary transforms do not change purity, which then can infer the determinant invariance of the covariance matrix det σ(0) = det σ(t).
In a similar way, we can write the final state in an analogous way as the initial one, that is: 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: 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: and their corresponding complex conjugates. It is worth mentioning that these equations can be corroborated by the use of the von Neumann equation for ρ(x, x , t) given in Equation (11), namely 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 δ 1 = δ 2 = 0, which can be also written in terms of the symplectic transformation Λ of Equation (5), i.e., σ(t) = Λ −1 σ(0)Λ −1 . Then, the covariance matrix evolution σ(t) can be obtained as: After differentiating each covarianceσ pp (t),σ qq (t), andσ pq (t), by using Equation (2), the inverse expression of Equation (14), the purity conservation condition pq (t), and the condition λ 1 λ 5 − λ 2 λ 4 = 1, we arrive at the following differential equations for the covariances:σ σ pq (t) = 2 (ω 1 (t)σ pp (t) − ω 3 (t)σ qq (t)) .
One can also check that these differential equations imply that the derivative of the determinant of σ(t) is equal to zero, i.e., d dt σ p p(t)σ qq (t) − σ 2 pq = 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: where the matrix B(t) contains the Hamiltonian coefficients, while Σ is a symplectic matrix, i.e., On the other hand, one can also check, using the inverse expression of Equation (5), that the mean values ofp andq follow the classical Equation (3), i.e., d dt 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 σ(0) and mean values p (0) and q (0).

Example
As an example, we consider the following Hamiltonian: In view of Equations (16) and (17), the matrix B can be identified as: and one can show that, in the case of constant frequencies (ω, ν), the evolution is determined by the same differential equations for all the covariances: which at ω 2 > ν 2 describes an oscillating motion with parameter f 2 = 4(ω 2 − ν 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: , and z 0 = σ pq (0). The solution for the classical equations of motion for the mean values reads: 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: one can also obtain the solution for the nonlinear Equation In both cases, we took frequencies ω = 2 and ν = 1.

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 conditioṅ σ(t) = 0, then one needs to obtain all the covariance matrices, which satisfy the condition σ(t)B(t)Σ − ΣB(t)σ(t) = 0. By taking the vectorṽ = (σ pp , σ pq , σ qq ), the equations forσ(t) = 0 can be written as follows: As the matrix M has rank R = 2, one can conclude that there is one nontrivial vector satisfying Equation (22). Exploring the null-space of M, one can check that the vector: with C being a constant, is the solution to Equation (22). We infer that all the states with a covariance matrix, given by: 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: with S = C 2 (ω 3 /ω 1 − ω 2 2 /ω 2 1 ) being the determinant of the invariant covariance matrix. In the case of the Hamiltonian (18), we have S = C 2 (ω 2 − ν 2 ), ω 3 /ω 1 = ω 2 , and ω 2 /ω 1 = ν, 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, ω = 2, and ν = 1, we have that the states with: are covariance invariant. The states that satisfyσ = 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 p (0) = q (0) = 0, the resulting states will be invariant, i.e., they will not change any of their properties over time (for δ 1 = δ 2 = 0). An example of such states for the Hamiltonian (18) are the ones in Equation (11), with parameters given by (25) and b(t) = b * (t) = 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.

Multidimensional Quadratic System
In this section, we review the equations determining the evolution of the covariance matrix and mean values p j and q j 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:Ĥ =rB(t)r +∆(t)r , where the tilde corresponds to the transposition operation and the vectorr = (p 1 ,q 1 ,p 2 ,q 2 , . . . ,p N ,q N ) corresponds to the vector of quadrature operators. The time dependence of this Hamiltonian is contained in the matrices: where B(t) is a real and symmetric matrix and ∆(t) is a real vector. As in the one-dimensional case, there exist 2N linear time-dependent operatorsR j (j = 1, . . . , 2N), whose time derivatives are equal to zero dR j dt = 0. . These operators can be arranged on a vector as follows: with the matrix Λ(t) and the vectorΓ(t) = (γ 1 (t), γ 2 (t), . . . , γ 2N (t)). By taking the time-derivative of the operator R j =R j and equating it to zero dR j dt = 0. , one can demonstrate that Λ(t) and Γ(t) must satisfy the following differential equations: The solution to these differential equations with the initial conditionsR j (0) = r j provides, as a result, the invariant operatorsR = P 1 ,Q 1 ,P 2 ,Q 2 , . . . ,P N ,Q N , satisfying the standard commutation rules for the operators r at time equal to zero: i.e., [R j , R k ] = [r j , r k ] = D j,k . This property leads us to the conclusion that the matrix Λ(t) must be symplectic and satisfies the equation: This relation can be then used to obtain the inverse of Λ(t), which results in the expression: The other important property of these invariant operators is that they correspond to the inverse evolution of the original operators, in other words, 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Λ(t) in Equation (29) satisfy the classical Hamilton equations after the time reversal operation, that is after the change p i → −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 r ( r (t)) and their covariances σ j,k = 1 2 {r j , r k } − r j r k . From the inverse of Equation (28), one can demonstrate that: as the invariant operators in R have a time derivative equal to zero, and they are equal to the standard operators r at zero time, then one can conclude that: From an analogous argument, one can see that the covariance matrix reads: Then, to obtain the expression for the time-derivative of the mean values r (t) and the covariance matrix σ(t), we make use of Equations (29) and (31)- (33) and arrive to the expressions: 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 × 2N symplectic matrix D contains in its diagonal blocks made of the 2 × 2 symplectic matrix Σ, that is, 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: where v is an N(2N + 1)-dimensional vector, which contains all the independent covariances in the N-partite system, and the matrix M is a square matrix of the same dimension that contains the Hamiltonian coefficients.

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 byp j andq 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: wherex = (x 1 , x 2 , . . . , x N ),x = (x 1 , x 2 , . . . , x N ), andỹ = (x 1 , x 2 , . . . , x N , x 1 , x 2 , . . . , x N ) are real vectors. Furthermore, we define the vectorb = (b 1 , b 2 , . . . , b N , b * 1 , b * 2 , . . . , b * N ) and the matrix: where the block matrices u and v can be written as: 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 x|ρ(t)|x , 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 × 2 covariance matrix and the mean values ( p (t) and q (t)), 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.

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: where ω j,k and γ 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 B(t), D, and σ(t) in the 2 × 2 block representation; in such a case, we have: where σ 1 (t) and σ 2 (t) are the covariance matrices for Subsystems 1 and 2, respectively, and σ 1,2 (t) 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., B(t) where the block matrices B 1 (t) and B 2 (t) are associated with Subsystems 1 and 2, respectively, while 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: Then, one can recognize the term 2(σ j B j Σ − ΣB j σ j ) for j = 1, 2, as the term corresponding to a unitary evolution of each subsystem (16). The extra term 2(σ 1,2B1,2 Σ − ΣB 1,2σ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.

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σ = 0, which can be expressed in terms of the following equation regarding the covariance and the Hamiltonian matrices σB(t)D − DB(t)σ = 0. As discussed before, this system of differential equations can be replaced byv = Mv (35) with the following correspondences: v = (σ p 1 p 1 , σ p 1 q 1 , σ p 1 p 2 , σ p 1 q 2 , σ q 1 q 1 , σ q 1 p 2 , σ q 1 q 2 , σ p 2 p 2 , σ p 2 q 2 , σ q 2 q 2 ) , and the matrix M containing the Hamiltonian coefficients is presented in Equation (A4) of Appendix B. It is possible to see that the matrix M has a determinant det M = 0 and a rank R = 8. From these properties, one can see that the systemσ(t) =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.

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: where the frequencies ω 1,2 are the input and output frequencies, respectively, ω is the pump field frequency, and the bosonic operatorsâ 1,2 are the annihilation operators of the input and output fields, respectively.
These beams interact with an intensity κ in a nonlinear medium as, e.g., a nonlinear crystal. In this case, the Hamiltonian matrix B(t) from (26) (in a unit system whereh = 1) reads: 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σ(t) = 0. As previously discussed,σ(t) = 0 can be expressed in vector form as: where v is defined as: v = σ p 1 p 1 , σ p 1 q 1 , σ p 1 p 2 , σ p 1 q 2 , σ q 1 q 1 , σ q 1 p 2 , σ q 1 q 2 , σ p 2 p 2 , σ p 2 q 2 , σ q 2 q 2 , and the matrix M is a 10 × 10 matrix of rank R = 8, which contains the Hamiltonian parameters of B(t).
To solve Equation (40), one can explore the null space of matrix M. The resulting null space contains two different vectors; one contains a non-physical solution. In this solution, σ p 1 q 2 = −ω 1 tan(ωt), σ q 1 q 1 = ω 1/2 2 (ω 2 − ω 1 ) sec(ωt)/(κω 1/2 1 ), σ q 1 p 2 = ω 2 tan(ωt), σ q 1 q 2 = 1, while all the other covariances are equal to zero. In particular, it contains the nonphysical terms σ p 2 p 2 = σ 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: 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 ω, ω 1,2 , and κ). 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: 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 σ p i p i σ q i q i − σ 2 p i q i ≥ 1/4 (i = 1, 2) and det σ ≥ 1/16.

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: where the frequencies ω 1,2 are the frequencies of the input and output beam channels and ω is the frequency of a pump field, which allows the amplification of the input channel. Then, the Hamiltonian matrix B(t) is: Following an analogous procedure to obtain the solutions of the equationσ = 0, one can show that the null space of the corresponding matrix 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 ω 1,2 > 0 can be written as: while all the other covariances are equal to zero. The other vector on the null space is: As the condition ω 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 Ĥ (t). These type of states can be defined by considering the initial covariances equal to C = 1/(ω 1 ω 2 ) times the ones presented in Equation (41); in other words, σ p 1 ,p 1 = 1, σ q 1 ,q 1 = 1/ω 2 1 , σ p 2 ,p 2 = ω 2 /ω 1 , σ q 2 ,q 2 = 1/(ω 1 ω 2 ) .
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 Ĥ (t), providing the strong coupling between the subsystems. In this particular example, one can see that σ p 2 p 2 (t) = σ q 2 q 2 (t) and σ p 1 q 1 (t) = σ p 2 q 2 (t) = 0.

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 | µ, ν) of the photon quadrature −∞ < X < ∞, to be measured in a reference frame with parameters µ = s cos θ and ν = s −1 sin θ. Here, −∞ < µ, ν < ∞, s is a time scaling parameter, and θ is the local oscillator phase, which is used in experiments [46] to obtain the Wigner function of the photon state.
The symplectic tomogram w ρ (X | µ, ν) is determined by the photon density operatorρ [45] as: whereq andp 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: and inversely, it determines the density operatorρ of the photon state: The optical tomogram of the photon state w opt (X | θ) ≡ w(X | µ = cos θ, ν = sin θ) 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: For Gaussian states (7), the tomographic-probability distribution of random photon quadrature X has the conventional form of a normal distribution: In view of (45), one has the mean value of the quadrature: and the covariance of the quadrature σ(µ, ν) reads: For the measured optical tomogram, the dispersion σ(θ) is: In the quantum system with Hamiltonian (18), the tomographic quadrature dispersion (51) evolves according to the formula: where σ qq (t), σ pp (t), and σ pq (t) are provided by the explicit expressions (19) and the parameters p(t) and q(t) 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 ν = 0. Since the dispersion matrix for the quadrature X is the linear combination of covariances σ qq (t), σ pp (t), and σ pq (t), 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 | E n E n | have the oscillator tomograms obtained from energy states | E n , whereĤ | E n = E n | E n . Tomograms of invariant Gaussian states satisfy the equality: where the parameter P n G is the probability to detect the properties of the stationary state | E n with energy value E n in the Gaussian state with the time-dependent tomogram w G (X | µ, ν). This state also does not depend on time.
Any convex sum of states | E n E n | is a density operator. One can conjecture that there is the decomposition of normal distribution w G (X | µ, ν) corresponding, e.g., to a thermal state withρ = exp(−Ĥ/(kT))/Tr(exp(−Ĥ/(kT))) (T being the temperature and k the Boltzmann constant), which can be presented as: 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. . The density matrix of the thermal equilibrium state with temperature T = β −1 in the position representation reads (here, we assumeh = ω = m = k = 1): The Green function of the oscillator has the Gaussian form: Since the density matrix (56) is determined by the Green function (57), i.e., with the partition function Z(β) given by the formula: Tr exp −βĤ | n n | = 1 2 sinh(β/2) ; here, we use the propertyĤ | n = (n + 1/2) | n . 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 | n n | do not depend on time.
In the position representation, the density matrix of Fock state | n n | reads: 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: where Z(β) 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, 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: The dispersion of quadrature X 2 = σ(µ, ν) given by Equation (63) reads: where the state with density matrix (56): Thus, the symplectic tomogram (63) is given by an invariant normal probability distribution: Since for optical tomogram µ = cos θ, ν = sin θ, and µ 2 + ν 2 = 1, in the case of the thermal single-mode photon state, its optical tomogram is: 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σ = 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.
The inverse transform provides the density operatorρ(1, 2) expressed in terms of the tomographic-probability distribution: The subsystem tomogram given by the partial trace of the density operatorρ(1) = Tr 2ρ (1, 2) reads: 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 M shown in Appendix B.

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.

Author Contributions:
The following statements The original idea was given by J.A.L.-S. All authors contributed equally to the conception, design, and methodology of this study. All authors contributed equally to the analysis of the results and the conclusions. All authors contributed equally to the final writing of the manuscript. All authors have read and agreed to the published version of the manuscript. where the values of the rest of the covariances read: σ q 1 ,q 1 = 2(a 22 + a * 22 − a 24 ) 4(a 11 + a * 11 − a 13 )(a 22 + a * 22 − a 24 ) − (a 12 + a * 12 + a 14 + a * 14 ) 2 , σ q 1 ,q 2 = a 12 + a * 12 + a 14 + a * 14 4(a 11 + a * 11 − a 13 )(a 22 + a * 22 − a 24 ) − (a 12 + a * 12 + a 14 + a * 14 ) 2 , (A2) σ q 2 ,q 2 = 2(a 11 + a * 11 − a 13 ) 4(a 11 + a * 11 − a 13 )(a 22 + a * 22 − a 24 ) − (a 12 + a * 12 + a 14 + a * 14 ) 2 , which can be used to determine the time evolution of the Gaussian system at any time.

Appendix B. Matrix M for Bipartite System
As discussed in Section 4.2, the evolution of the covariance matrix of a bipartite system can be written as the linear system of equations: where the vector containing the different covariances reads: v = σ p 1 p 1 , σ p 1 q 1 , σ p 1 p 2 , σ p 1 q 2 , σ q 1 q 1 , σ q 1 p 2 , σ q 1 q 2 , σ p 2 p 2 , σ p 2 q 2 , σ q 2 q 2 , and the matrix M is obtained by analyzing the derivatives of the covariance matrix; this results in an expression for M given by: