Gaussian States: Evaluation of the Covariance Matrix from the Implementation with Primitive Component

: Quantum Gaussian states play a fundamental role in quantum communications and in quantum information. This paper deals with the implementation of multimode, and particularly of two-mode Gaussian unitaries and Gaussian states with primitive components (phase shifters, single-mode real squeezers, displacements, and beam splitters). The architecture thus obtained allows one to obtain an insight into the physical meaning of each variable involved. Moreover, following the implementation architecture, it is possible to formulate an easy algebra (radical free) for the main operations and transformations of Gaussian states.


Introduction
In recent years, quantum systems based on continuous variables have attracted great attention for the development of quantum information transmission and processing [1][2][3]. A large amount of theoretical work has been dedicated to multimode quantum states and transformations, with the prominence of Gaussian states and Gaussian transformations which can be more easily implemented and manipulated. In particular, the multimode case is attracting a lot of interest, since it may exhibit entanglement, which is a key resource for several quantum protocols and applications, such as teleportation, computing, and cryptography.
An efficient approach dealing with Gaussian unitaries and Gaussian states is based on the decomposition of the so-called fundamental Gaussian unitaries (FGUs), which are given by displacement, rotation, and squeezing, as shown in Figure 1. Combination of these three FGUs allows for the implementation of the whole class of multimode Gaussian unitaries and the generation of the whole class of multimode Gaussian states, arriving at explicit closed-form results [4]. Another important approach is obtained resorting to the Bloch-Messiah reduction [1,5] in which the multimode squeezer is decomposed into the parallel of local and real single-mode squeezers, as shown in Figure 2. There are several tools to efficiently describe Gaussian unitaries/states. For Gaussian unitaries, the main tools are given by Bogoliubov and symplectic transformations; in the N-mode, the former are specified by two complex matrices of order N, while the latter are specified by a 2N × 2N real matrix. For Gaussian states, the main tool is given by the covariance matrix, a 2N × 2N symmetric matrix with real elements (the mean vector is often neglected in the analysis).
A key role is played by the specification of the variables involved, which determine the degree of freedom of Gaussian unitary/states. For the decomposition of Figure 1, the variables are: the rotation matrix φ, the squeeze matrix z, and the displacement vector α; this is the usual specification for a unitary matrix, which we call for convenience the algebraic specification. For the decomposition of Figure 2, the variables are: two rotation matrices, φ and γ, a real diagonal squeeze matrix r D , and the displacement vector α which we call structural or physical specification. Since a Gaussian state is obtained by processing a multimode thermal state through a Gaussian unitary, in order to characterize a Gaussian state, we have to add the N single mode thermal states in the specification.
However, with the algebraic specification, the first step is the evaluation of the polar decomposition of the squeeze matrix z. The formula contains radicals of radicals, which "propagate" on all the subsequent elaborations, leading to cumbersome formulas. This is completely avoided if we start from the structural specification.
The target of the paper is the implementation of Gaussian unitaries/states developing the cascade of Figure 2 with primitive components, starting from the structural specification and using the technique of Bloch-Messiah reduction. The primitive components are rotations, squeezers, displacements in the single mode, and beam splitters [6]. The architecture thus obtained allows one to obtain an insight into the physical meaning of each component involved. Moreover, following the implementation architecture, it is possible to formulate an (radical-free) easy algebra for the main operations and transformations of Gaussian states. This paper is organized as follows. In Section 2, we formulate the Gaussian unitaries and their decompositions in FUGs, according to the Bloch-Messiah (BM) reduction. In Section 3, we formulate Gaussian states according to Williamson's theorem. Section 4 deals with the implementation of Gaussian unitaries with primitive components, where the squeeze matrix is decomposed according to the Takagi factorization [7]. The implementation could be considered for multimode [6], but for simplicity, it is carried out for the two-mode one. Continuing with the two-mode, from the architecture with primitive components in Section 5, we evaluate the symplectic transformations and in Section 6, the covariance matrix. As we will realize, all the results are simple and radical free (for comparison, the same evaluation is carried out in Appendix A using the algebraic approach). In the final part of the paper, we focus our attention on the covariance matrix, in which the two-mode depends on 10 real variables. Examples of application are outlined by fixing the variables to specific values. Furthermore, in Section 8, we consider the so-called standard form of the covariance, which only depends on four real variables (symplectic invariants) and contains all the relevant information on the Gaussian state, particularly the entanglement property. We outline a simple method for the evaluation of the symplectic invariants of arbitrary Gaussian states in the two mode.

Gaussian Unitary in Terms of Fundamental Unitaries (FGUs)
We remind that a quantum transformation is Gaussian when it transforms Gaussian states into Gaussian states. When the Gaussian transformation is defined by a unitary map, it is called Gaussian unitary.
It is known [2] that Gaussian unitaries can be expressed in the form U = exp(−iH/2), where the Hamiltonian H is a second-order polynomial in the bosonic operators a i and a * i or in the field operators q i , p i . When these unitaries are applied to the annihilators a = [a 1 , . . . , a N ] T , we obtain a Bogoliubov transformation, and, in terms of field operators, a symplectic transformation. For the purpose of describing Gaussian unitaries, the two types of transformations are equivalent. Any Gaussian unitary can be expressed by FGUs. The fundamental unitaries are defined in terms of the column vectors a * and a of size N, collecting the bosonic operators a * m and a m , as: 1.
N-mode displacement operator where α is called the displacement vector.

2.
N-mode rotation operator where φ is called the rotation matrix. where z is called the squeeze matrix. The polar decomposition applied to z gives z = re iθ where r is positive semidefinite and θ is Hermitian symmetric. The matrix r is uniquely determined by z as r = √ z×z. If r is positive definite, the phase matrix θ associated with the squeeze operator is determined as e iθ = r −1 z. Then [4], Theorem 1. The most general Gaussian unitary is obtained by the cascade of the three fundamental Gaussian unitaries D(α), S(z), and R(φ), where the order is arbitrary, that is, S(z) D(α) R(φ), R(φ) D(α) S(z), etc.
One can pass from one cascade to any other one using the so-called switching rules, which provide a slight change in the parameters [4]. In the following, without loss of generality, we refer to the cascade which is illustrated in Figure 1.

Gaussian Unitaries Decomposed according to the Block-Messiah Reduction
The Block-Messiah reduction [1,5] was reconsidered in [7] in terms of the Takagi factorization [8].
Theorem 2. The most general Gaussian unitary can be decomposed into the cascade of a rotation operator R(ψ), a local squeeze operator S(r D ) with a real diagonal matrix r D , another rotation operator R(γ), and finally a displacement operator D(α).
The theorem, illustrated in Figure 2, is essentially based on the Takagi factorization [8] of the squeeze matrix z = U z r D U T z = e iβ r D e −iβ (5) which leads to the decomposition of the squeeze operator as S(z) = R(β)S(r D )R(−β). Then, the general unitary given by (4), becomes where R(ψ) = R(−β)R(φ) (with ψ uniquely determined by e iψ = e −β e iφ ). From (6), one can easily obtain the original form (4), in particular the squeeze matrix Starting from other cascades in the application of Theorem 1, the Bloch-Messiah reduction leads to other architectures, but all are equivalent in representing the most general Gaussian unitary.

Specification and Degree of Freedom
Theorem 1 states that the specification of any N-mode Gaussian unitary is given by the matrices α, z, φ. Considering their symmetry, the associated degrees of freedom are As a verification, from Theorem 2 and the BM architecture of Figure 2, for the twomode, one finds that the number of degrees of freedom is given by 14, which is to be expected, due to the equivalence between the two architectures.

Gaussian States
Gaussian states can be obtained from a Gaussian unitary driven by thermal states; in particular, pure Gaussian states are obtained when the thermal states degenerate into vacuum states. To formalize this statement, it is convenient to recall that a Gaussian state is completely described by the covariance matrix and the mean vector. In particular, the covariance matrix can always be written in the form (Williamson's theorem) where S is an N-mode symplectic matrix and corresponds to a tensor product of N thermal states with an average number of thermal photons The quantities {ν k } are called the symplectic eigenvalues of the CM V and the operation performed by the matrix S is said to be the symplectic diagonalization of V.
With reference to the previous decompositions of the Gaussian unitaries shown in Figures 1 and 2, we find that, when the input is driven by N thermal states, at the output, we obtain the most general N-mode Gaussian state.
Note that a decomposition into unitaries of the Lie group Sp(2N, R) has also been used in [9,10] for the dynamical symmetry group describing the vibronic transitions in polyatomic molecules.

Specification and Degree of Freedom
To the specification of the Gaussian unitary, one has to add the specification of the input thermal noise, given by N real variables, as in (9). On the other hand, we have a reduction of the degrees of freedom because a rotation operator is irrelevant for the input thermal noise, such as the rotation R(φ) of Figure 1, or the rotation R(ψ) of Figure 2. Then, the degree of freedom of an N-mode Gaussian state results in

Gaussian States in the Two-Mode
Since we focus on two-mode Gaussian states, we review in detail their specification, which is given by the complex matrices (see Figure 1) and two thermal noises N 1 , N 2 . From (8), a two-mode Gaussian unitary has 14 real variables as degrees of freedom, while from (12), a two-mode Gaussian state has 12 real variables degrees of freedom.
We want to evaluate the quantities involved in the two-mode in terms of the parameters given by (13).

Implementation with Primitive Components
The implementation of Gaussian unitaries and the generation of Gaussian states with primitive components could be carried out in the general N-mode [6], but hereafter, it will be developed in great detail in the two mode.

Primitive Components
They are: (1) single-mode displacement; (2) single-mode rotation operators, briefly shifters; (3) single-mode real squeezers; and (4) beam splitters (BSs). Figure 3 shows the graphical symbols we use for these components. A shifter is specified by a phase β ∈ [0, 2π), which leads to the 1 × 1 exponential matrix e iβ . A single-mode squeezing operator is specified by the squeeze factor r ∈ R. A free-phase BS is defined by a two-mode rotation operator corresponding to the phase matrix with antidiagonal entries where φ describe the reflectivity (s := sin φ) and the transmittance (c := cos φ) of the BS. An equivalent specification is given by the rotation matrix

Implementation in the Two Mode
The target is the implementation in the two-mode of the architecture of Figure 2 with primitive components. The squeezer is already decomposed into primitive components. The implementation of the two-mode displacement operator D(α) is trivial, since it is given by two parallel single-mode displacement operators D(α 1 ) and D(α 2 ). For the rotation operators, we have (for the proof, see the appendixes of [6]). Proposition 1. An arbitrary two-mode rotation operator defined by the unitary matrix U = e iψ := [ρ hk e iγ hk ] can be implemented by the cascade of (1) two phase shifters with phases ψ 11 and ψ 12 , (2) a BS with reflectivity s = ρ 12 , and (3) a phase shifter with phase µ = ψ 22 − ψ 12 , as illustrated as in Figure 4. (The scheme of Figure 4 is sometimes called beam splitter with phase [11,12].) The matrix U = e iψ can be written as and also The interpretation of the above factorization proves the proposition.
We are now ready to draw the Gaussian unitary of Theorem 1 with primitive components (see Figure 5). Note that this architecture can be used to specify the whole class of two-mode Gaussian unitaries. It is given by six shifters, two BSs, two real squeezers, and two displacements, corresponding to 14 real variables, in agreement with (8). The second target is the generation of two-mode Gaussian states. They are obtained when the architecture of Figure 5 is driven by two thermal noises ( Figure 6). Now, the shifters γ 11 and γ 12 are removed because they are irrelevant when driven by thermal states. Furthermore, the final displacements were removed because they do not give a contribution to the covariance matrix V (whose evaluation is the main target of the paper).

Strategy for the Evaluation
The above architectures (Figures 5 and 6) represent the key basis of our theory, which has the goal to understand the action of each primitive component as well as obtain simple and meaningful formulas for the main descriptions of Gaussian unitaries/states, in particular to be "radical free". This possibility relies ultimately on the Bloch-Messiah reduction.
In the next sections, following the architectures, we will evaluate the symplectic matrix and the covariance matrix. This is called for convenience the structural approach or physical approach, in competition with the algebraic approach more often considered in the literature (in special cases, e.g., with EPR states, the algebraic approach turns out to be easy, but not in the general case).
Here, we develop only the first step of the algebraic approach to see how the complication is originated. From the algebraic specification given by (13), to proceed, one has to obtain the polar decomposition of the squeeze matrix z and the singular values (SVs) of z, which are given by where P = 4 z 11 z * 12 + z 12 z * To proceed, we have to distinguish between the cases of distinct and coincident SVs because the matrix functions have different formulas. This is a further complication. The details of the algebraic approach are given in Appendix A. Example 1. We consider squeezing with matrix In the algebraic approach, the squeeze matrix is one of the data and the first step requires the evaluation of the polar decomposition z = r e iθ , which results in Additionally, the SVs of z must be evaluated. They result in σ 1 = r 0 / √ 2 and σ 2 = √ 2r 0 . With the structural approach, the squeeze matrix (20) is obtained with the data: r 1 = r 0 , r 2 = r 0 /2 , γ 11 = π/2 , γ 12 = 0 and balanced BSs, so that the free variables are only r 0 and . However, neither the polar decomposition nor the SVs are needed to proceed.

Evaluation of Symplectic Matrices The Symplectic Matrix (SM)
Here, we consider the real interlaced SM S, where the phase-space variables are arranged in the form X := [q 1 , p 1 , q 2 , p 2 ] T . A symplectic transformation is of the form where S is a 4 × 4 real matrix and d ∈ R 4 . The commutation relations require the condition A matrix S is called symplectic if it verifies this condition. Following the architecture of Figure 5, we find where we can use the general formulas We developed the structural approach for the description of Gaussian unitaries based on the primitive-components architecture of Figure 5. Clearly all the formulas are "radical free".

Evaluation of the Covariance Matrix (CM)
This is an important topic considered by several authors, e.g., [13][14][15][16][17]. In particular, [15] uses a similar approach to characterize Gaussian states, based on Williamson's theorem and Bloch-Messiah decomposition, describing the symplectic transformations that correspond to the fundamental operators of squeezing, rotation and beam splitters. However, the final results are only provided for the single-mode and a couple of specific examples of two-mode states. In [13], the integration within an ordered product (IWOP) is used to express the covariance matrix of N-mode in terms of the squeeze matrix. Note, however, that the hyperbolic tangent of the squeeze matrix must be evaluated. The final derivation is obtained in two specific examples for the two-mode case. Ref. [14] derives the elements of the covariance matrix of a two-mode Gaussian state, with a generalized squeezing, by using the technique of integration within an ordered product (IWOP). The objective is to characterize the entanglement and non-locality properties of the corresponding state.
In order to characterize the CM, it is more convenient to express the result partitioned into 2 × 2 blocks. Letting S = S 11 S 12 S 21 S 22 we obtain: where Considering that the CM refers to Gaussian states, as remarked above in Figure 6, the dependence is on the 10 real parameters ν 1 , ν 2 , s, µ, r 1 , r 2 , γ 11 , γ 12 , p, , notwithstanding that the S ij also depend on the phases ψ 11 , ψ 12 .

Standard Form of the Covariance Matrix
We remind that, for any two-mode Gaussian states, there exists a local symplectic operation S l that brings the covariance matrix to the standard form, Simon2000, [18] where the correlation terms a, b, c + and c − are obtained from the ordinary CM V by the four local symplectic invariants The importance of the standard form and of the symplectic invariants lies on the fact that they concentrate all the relevant information on the two-mode Gaussian states, particularly that concerning the entanglement [19]. This means that the essential degree of freedom is reduced to 4 real variables, instead of the 10 variables listed above.
Evaluation of "invariants" (a, b, c ± ) from the ordinary CM V. The invariants can be obtained from the blocks of the ordinary (nonstandard) covariance matrix (see (30)).

Proposition 2.
From the blocks of the ordinary covariance matrix V, evaluate the quantities Then, the invariants result in: a = √ det A, b = √ det B, and The proof is in Appendix B. We give a numerical example. With the data

Physical Analysis of the Global Architectures
The architectures of Figures 5 and 6 allow us to obtain a physical insight into the dependence of the global performance, as stated by the symplectic matrix and by the covariance matrix, from each primitive component. In the limit, we can remove a single component from the architecture to see the effect on the ordinary covariance matrix V. In a forthcoming paper [20], we will show that a quantum state, having as ordinary covariance just the standard covariance matrix V s f , can be obtained by removing all the shifters. This shows that the shifters have less relevance with respect to the other primitive components.

Examples
By choosing the 14 parameters listed above, one can evaluate the description of all Gaussian unitaries in the two-mode. Analogously, by choosing the 10 parameters of Figure 6, one can evaluate the description of all two-mode Gaussian states. Here, we outline a few cases of Gaussian states focusing the attention on the covariance matrix.

Case 1: All-zero phases
A first case is considered with all zero phases, with parameters as in Table 1. thermal states The squeeze matrix results in The blocks which constitute the covariance matrix V are given by where c 1 = (ce r 2 q − e r 1 ps)(ce r 2 p + e r 1 qs)ν 2 − (ce r 1 q − e r 2 ps)(ce r 1 p + e r 2 qs)ν 1 Case 2: All-zero phases, r 1 = r 2 = r 0 , and the second BS balanced We consider another case, with all zero phases, the same squeezing on each mode and balanced BS, as in Table 2.
The squeeze matrix results in The blocks of the covariance matrix V are given by

Case 3: EPR state with noise
The EPR unitary is a squeeze with the following matrix z = 0 e iθ r 0 e iθ r 0 0 To obtain this unitary, the primitive-component architecture must have the following data r 1 = r 2 = r 0 , balanced BSs p = q = s = c = 1 √ 2 and the phases indicated in Table 3. Table 3. Parameters of the architecture in Case 3. In red the free variables.

Separability and Entanglement
The separability and the entanglement properties of any two-mode Gaussian state can be handled from both the ordinary CM V and the standard CM V s f . Now, a Gaussian state is separable if its standard CM is diagonal. Hence, from (32), the separability condition is When the standard CM is evaluated from the ordinary CM using Proposition 2, Equation (35) is not useful because it leads to degeneracy. Instead of (33), namely det V = (ab − c 2 + +)(ab − c 2 − ) and det C = c + c − , one finds the conditions on the symplectic invariants On the other hand, the separability based on the ordinary CM is that V must be diagonal. Then, from (30) and (31), the conditions read Recalling that ν 1 ≥ 1, ν 2 ≥ 1, from (31), we find the conditions S 11 S T 11 , S 12 S T 12 , S 21 S T 21 , S 22 S T 22 diagonal S 11 S T 21 = 0 , S 12 S T 22 = 0 (37)

Conclusions
For the description of Gaussian unitaries/states, we considered the algebraic and structural approaches. We remark that the class of noisy Gaussian states generated by the algebraic specification and by the structural specification coincide, since both approaches generate the whole class of noisy Gaussian states. In other words, one can choose the first class as well as the second class as the definition of noisy Gaussian states. In terms of performance, both the methods and also alternative solutions proposed in [13][14][15] achieve the same covariance matrix, since no approximations are introduced, but at different costs, in terms of complexity and often only in some particular cases. In particular, in [15], the solution is provided only for a couple of specific examples. In [13], the covariance matrix is obtained considering only a squeezed state, not a general one, by means of integration within an ordered product (IWOP), starting from the knowledge of the correlation operator. Note also that the hyperbolic tangent of the squeeze matrix must be evaluated. Furthermore, in [14], only the squeezing is considered, again with the same complexity of IWOP, and the objective is to characterize the entanglement of the corresponding state.
The main target of the paper is to show the advantages of the structural approach in the description of the most generic Gaussian states, namely:

1.
The structural approach is completely radical free compared to the algebraic approach and requires several matrix operations leading to results that contain the radicals of radicals. Note that the key to avoiding radicals is the following: the quantities which exhibit radicals in the algebraic approach become independent variables (data) in the structural approach.

2.
The structural approach is completely general, while the algebraic approach exhibits several degeneracies (mainly coincident eigenvalues) concerning some very important cases (see EPR states). These cases should be treated separately with ad hoc procedures. Such a distinction is not required in the structural approach. 3.
In the structural approach, all the variables have a precise physical meaning, related to the corresponding components of the architecture, i.e., squeezers, beam splitters, phase shifters, and one can choose the specific variables to achieve the desired properties of the covariance matrix, for example, the entanglement.
Finally, we note that the theory, here developed in detail for the two-mode, could be extended to higher modes, of course with the penalty of complication at the increase in the order. We developed the three-mode (not reported here) with the structural approach and without difficulty, resulting, however, in long formulas (which were radical free). With the algebraic approach, the three mode results to be very complicated, mainly for the large number of particular cases to be treated separately. Hence, with G(x) = cosh(x) and G(x) = sinh(x), we can obtain the desired matrix functions. With coincident eigenvalues, one obtains cosh(r) = cosh(σ 1 ) I 2 , sinh(r) = sinh(σ 1 ) I 2 .

Appendix
we find four solutions given by