No-activation theorem for Gaussian nonclassical correlations by Gaussian operations

We study general quantum correlations of continuous variable Gaussian states and their interplay with entanglement. Specifically, we investigate the existence of a quantum protocol activating all nonclassical correlations between the subsystems of an input bipartite continuous variable system, into output entanglement between the system and a set of ancillae. For input Gaussian states, we prove that such an activation protocol cannot be accomplished with Gaussian operations, as the latter are unable to create any output entanglement from an initial separable yet nonclassical state in a worst-case scenario. We then construct a faithful non-Gaussian activation protocol, encompassing infinite-dimensional generalizations of controlled-NOT gates to generate entanglement between system and ancillae, in direct analogy with the finite-dimensional case. We finally calculate the negativity of quantumness, an operational measure of nonclassical correlations defined in terms of the performance of the activation protocol, for relevant classes of two-mode Gaussian states.


I. INTRODUCTION
Quantum correlations in composite systems transcend entanglement [1]. A bipartite quantum state ρ AB can be defined as nonclassical or nonclassically correlated if it cannot be expressed as a convex mixture of local basis states of subsystems A and B [2]. Consequently, all inseparable (entangled) states as well as the majority of separable states are nonclassical.
General nonclassical correlations, however, can be mapped to entanglement in a very precise sense, which provides an insightful framework for their characterization and operational interpretation. Specifically, it was proven in [3][4][5] and very recently experimentally observed in [6] that all nonclassical states of a finite-dimensional system can be turned into states with distillable entanglement between the system and a set of ancillae by an activation protocol. Focusing on a bipartite setting, the protocol runs as follows. The subsystems A and B are first subject to arbitrary local unitary transformations U A,B ; then, each system j = A, B interacts via a controlled-NOT (CNOT) operation U CNOT j j (i.e. a so-called premeasurement interaction) with an auxiliary system j , j = A, B, initialized in a pure state |0 j .The activation protocol then possesses two key properties: i) for all classical states ρ AB at the input of the protocol, there exist local unitaries U A,B for which the output state ρ ABA B is separable across the AB|A B splitting, and ii) for all nonclassical states ρ AB and for all local unitaries, the output state is entangled across the AB|A B splitting.
Let us stress that both criteria i) and ii) must be met by any scheme in order to be a valid activation protocol. In particular, they allow us to define faithful measures of nonclassical correlations for the input state ρ AB in terms of the output AB|A B entanglement, minimized over U A,B . One such measure, when the output entanglement is quantified by the negativity [7], has been termed negativity of quantumness [3,8], and has been experimentally investigated in [6,9] In this paper we study activation of nonclassical correlations in multimode bipartite Gaussian states ρ AB of continuous variable systems [10]. Nonclassical correlations of Gaus-sian states have been studied extensively both theoretically and experimentally [11][12][13][14] but their interplay with entanglement has not been pinned down so far in terms of the activation framework. Attempts to devise activation-like protocols for Gaussian states have been explored [15]. However, these differed significantly from the original prescription in that nonunitary operations were employed between system and ancillae, so that the entanglement generation was obtained as a dynamical feature, and conditions i) and ii) were not generally verified.
Here we consider a general Gaussian activation protocol in which U A,B are Gaussian unitaries and the CNOT gates are replaced with a global Gaussian unitary on subsystems A, B, A , B . In Section II we then prove that any such protocol satisfying condition i) will unavoidably violate condition ii), which implies that activation of Gaussian nonclassical correlations by Gaussian operations is impossible. This fact establishes a new no-go theorem for Gaussian quantum information processing, which can be enlisted alongside other well known no-go results such as the the no-distillation theorem, according to which distilling entanglement from Gaussian states by using only Gaussian operations is impossible [16]. We then show in Section III how, by using non-Gaussian operations which properly extend the CNOT to infinite dimensions, one can construct the continuous variable counterpart of the activation protocol of [3], verifying criteria i) and ii). This allows us to define the negativity of quantumness for Gaussian states and to calculate it for relevant examples in Section IV. This work provides an operational setting to understand and manipulate nonclassical correlations in paradigmatic infinitedimensional systems. We draw our conclusions in Section V, while some technical derivations (which can be of independent interest) are deferred to the Appendices.

II. GAUSSIAN NO-ACTIVATION THEOREM
Gaussian states are quantum states of systems with an infinite-dimensional Hilbert space (continuous variable sys-arXiv:1405.6859v2 [quant-ph] 5 Aug 2014 tems), e.g. a collection of harmonic oscillators, which possess a Gaussian-shaped Wigner function in phase space [10]. L modes are described by a vector r = (x 1 , p 1 , . . . , x L , p L ) T of quadrature operators x j , p j satisfying the canonical commutation rules expressible in terms of elements of the vector r as [r j , r k ] = iΩ jk , j, k = 1, . . . , L with Ω = ⊕ L j=1 iσ y , where σ y is the Pauli y-matrix. An L-mode Gaussian state ρ is fully characterized by a 2L × 1 vector r of the first moments with elements r i = Tr(ρr i ) and by its 2L × 2L covariance matrix (CM) γ with elements γ i j = ∆r i ∆r j + ∆r j ∆r i /2, i, j = 1, . . . , L, where ∆r i = r i − r i . Gaussian unitaries are generated by Hamiltonians that are quadratic in the quadrature operators and they preserve the Gaussian characteristic of quantum states. An L-mode Gaussian unitary U(S ) is represented in phase space by a 2L×2L real symplectic transformation S satisfying the condition S ΩS T = Ω, which transforms a CM γ to S γS T .
Here we are interested in the question of whether an activation protocol exists satisfying conditions i) and ii) which would rely solely on Gaussian states and Gaussian unitaries. We therefore assume the state ρ AB to be a Gaussian state of (N + M) modes with CM γ AB and the state ρ A B of the ancilla to be also a Gaussian state with CM γ A B . The local unitaries U A,B of the original discrete protocol are replaced with local Gaussian unitaries U A (S A ) and U B (S B ) represented by the symplectic matrices S A and S B , respectively. Likewise, the global operation U CNOT AA ⊗ U CNOT BB on the whole system ABA B is replaced with one global Gaussian unitary U(S ) represented by a symplectic matrix S .
Let us recall the definition of a fully classical state [2,3,5]. Suppose ρ AB is a bipartite state containing two subsystems A and B with N and M modes respectively, and let B j = {|B j (n j ) } be a basis of subsystem j, with n A = (n A 1 , . . . , n A N ), n B = (n B 1 , . . . , n B M ) and n j i ∈ N 0 . If there exists a basis B consisting of the tensor products of all elements of B A with all elements of B B , then ρ AB is a classical state if it is diagonal with respect to B. It has been shown in [11,17] that a twomode Gaussian state is classical if and only if it is a product state, i.e., its CM is represented by a direct sum γ A ⊕γ B of local CMs γ A,B . One can prove that this statement remains valid for the generic case of bipartite (N + M)-mode Gaussian states (see Appendix A for the proof). Therefore, all non-product bipartite Gaussian states (including separable ones) are nonclassical. According to condition i) in any Gaussian activation protocol with an input Gaussian product state there must exist local Gaussian unitaries U A,B for which one gets a separable state ρ ABA B across the AB|A B splitting at the output of the protocol. We will show however, that this implies that for all separable Gaussian states including nonclassical ones, there exist local Gaussian unitaries U A,B for which the output state is separable. That is, condition ii) is not satisfied. Thus, any Gaussian activation protocol described above cannot meet simultaneously criteria i) and ii), and hence does not exist.
The proof of this no-go theorem is depicted in Fig. 1. It follows from the decomposability of any Gaussian separable state into a product state and noise [18], and the linearity of symplectic transformations. Namely, for any separable Gaussian state with CM γ AB there exist local CMs γ A,B such that (2), cannot turn a separable state into an entangled state and therefore the protocol transforms the separable state ρ AB into a state which is separable across the AB|A B cut (thick dashed line). See text for details.
In other words, any separable Gaussian state with CM γ AB can be prepared from a suitable product state with CM γ A ⊕ γ B by the addition of noise, represented by a positivesemidefinite matrix P, i.e., γ AB = γ A ⊕ γ B + P. The noise can be created by displacing the vector of quadratures r = T of the product state as r → r + VR. Here V is a 2(N + M) × K matrix given by the first K columns of the matrix V bringing the matrix P to the diagonal form V T PV = diag (λ 1 , λ 2 , . . . , λ K , 0, 0, . . . , 0), where λ 1 , . . . , λ K denote K ≤ 2(N + M) strictly positive eigenvalues of the matrix (1), and R = (R 1 , . . . , R K ) T is the vector of classical displacements uncorrelated with the vector of quadratures r and distributed according to the Gaussian distribution with zero means and the diagonal correlation matrix diag (λ 1 , λ 2 , . . . , λ K ).
Let us now consider a separable state with CM γ AB at the input of a Gaussian activation protocol and let γ A ⊕γ B be a CM of the product state from which the state can be prepared using the aforementioned algorithm. Assume that the local symplectic matrices S A,B are chosen such that the CM γ (0) is the CM of the state of the ancilla, is separable across the AB|A B splitting. Hence, for the original separable state with CM γ AB , the output of the activation protocol is obtained by displacing the vector of quadratures r (0) for the state with CM γ (0) ABA B by where O is a 2T × 1 zero vector with T being the number of modes of the ancilla A B . However, for a separable state with CM γ (0) ABA B , where AB is separable from A B , the local displacements (2) cannot create a state in which the system AB is entangled with the system A B . Consequently, for any separable state (even nonclassical) it is always possible to find local Gaussian unitaries for which the output is separable, thus accomplishing the proof of the no-go theorem.
Therefore, Gaussian operations are unable to activate nonclassical correlations of Gaussian separable states into entanglement in the worst-case scenario: assuming condition i) holds, then for any Gaussian separable state there exist local Gaussian unitaries for which the output of the activation protocol remains a separable Gaussian state. This indicates that a non-Gaussian element, like a non-Gaussian global unitary U or a non-Gaussian state of the ancilla, is necessary for faithful activation of nonclassical correlations in Gaussian states. In the following we design such an activation protocol involving a non-Gaussian CNOT gate in the Fock basis and an ancillary system in a Gaussian state.

III. NON-GAUSSIAN ACTIVATION PROTOCOL
The main benefit of the activation protocol is that it allows one to quantify the amount of nonclassical correlations in a given quantum state as the potential to create entanglement in the activation protocol [3][4][5]. More precisely, if E AB|A B (ρ ABA B ) denotes an entanglement measure quantifying the amount of entanglement between systems AB and A B in a quantum state ρ ABA B , then we can define a measure of nonclassical correlations on the input state ρ AB as where the minimization is carried out over all local unitaries U A and U B on subsystems A and B. It has been proven in [5] that From now on we assume that systems A and B each contain one mode. The non-Gaussian activation protocol is obtained as a direct generalization of the finite-dimensional protocol [3]. At the input we allow for generally non-Gaussian states ρ AB of continuous variable systems, local unitaries U A and U B , and the global Gaussian unitary U(S ) of the preceding protocol is replaced with the tensor product V ≡ U CNOT AA ⊗ U CNOT BB of the infinite-dimensional generalizations of CNOT gates in the Fock basis, where |m, n j j ≡ |m j ⊗ |n j , m, n = 0, 1, . . ., and |k l is the kth Fock state of mode l. We also assume the initial state ρ A B of the ancilla A B to be the vacuum state |0 A 0| ⊗ |0 B 0|. Hence, the final output state can be expressed as whereρ By following arguments similar to the finite-dimensional case [3], one can show that the non-Gaussian activation protocol defined above satisfies both criteria i) and ii). For condition i) we assume that ρ AB is classically correlated and hence there exist local unitaries U A and U B such that the density matrixρ, Eq. (6), takes the formρ = ∞ n,m=0 p n,m |n, m AB n, m|. Making use of Eqs. (4) and (5) it then follows that the output state of the protocol is the following convex mixture of product states, ρ ABA B = ∞ n,m=0 p n,m |n, m AB n, m|⊗|n, m A B n, m|, and is thus a separable state across the AB|A B splitting as required.
For the proof of condition ii) we now suppose that the density matrix ρ AB is nonclassical and show that then the density matrix ρ ABA B given in Eq. (5) is entangled across the AB|A B cut for all local unitaries U A and U B . To prove the presence of entanglement in ρ ABA B we will use the negativity N defined in [7] as Here . 1 denotes the trace norm, ρ T AB ABA B is the partial transpose [19] of the state ρ ABA B with respect to subsystem AB, and a strictly positive value of negativity implies that the state ρ ABA B is (distillable) entangled with respect to the AB|A B splitting. The specific feature of the present activation protocol is that the output state ρ ABA B is a so-called maximally correlated state and therefore, following results in [3,8], the output negativity can be expressed as where m = (m 1 , m 2 ), n = (n 1 , n 2 ) andρ m,n = AB m 1 m 2 |ρ|n 1 n 2 AB are elements of the density matrixρ, Eq. (6), in the Fock basis. Since our input state ρ AB is nonclassical, the stateρ is also nonclassical for any choice of unitaries U A and U B . Thus, there must be at least one non-zero off-diagonal elementρ m,n for every choice of U A and U B . Hence, Eq. (8) implies N(ρ ABA B ) > 0 and the output state ρ ABA B is entangled for any nonclassical input state. This completes the proof of our non-Gaussian activation protocol.

IV. EXAMPLES
The optimization in Eq. (3) is generally carried out over all local unitary operations U A and U B , including non-Gaussian ones, which is not a tractable task. Here we consider input Gaussian states with CM in standard form [20], and consider the non-optimized output entanglement E AB|A B (ρ ABA B ) obtained when the local unitaries U A,B are selected to be identity matrices. Therefore, the state (6) remains a Gaussian state in standard form with the following CM where A = diag(a, a), B = diag(b, b) and C = diag(c 1 , c 2 ) are diagonal matrices. In what follows we determine the nonoptimized quantity for some classes of two-mode Gaussian states by considering the negativity (7) as an entanglement measure E, and using Eq. (8). The corresponding measure of nonclassical correlations Q N (ρ AB ) is called the negativity of quantumness [3] accordingly. Although our choice of local unitaries U A,B gives in general an upper bound on Q N , we find that it coincides with the true measure on pure states, leading us to conjecture that our choice is optimal for calculating the negativity of quantumness of all two-mode Gaussian states in standard form. Verifying this conjecture numerically is beyond the scope of this work.

A. Pure states
A closed form of the output negativity can be found for pure two-mode Gaussian states. The density matrixρ amounts to that of the two-mode squeezed vacuum state, withρ m,n = [1 − tanh 2 (r)](tanh r) m 1 +n 1 δ m 1 ,m 2 δ n 1 ,n 2 , where r ≥ 0 is the squeezing parameter. Hence, by a direct substitution into Eq. (8) we get Consequently, as the output negativity N(ρ ABA B ) is equal to the negativity of the input state ρ AB , it coincides with the true optimized negativity of quantumness Q N (ρ AB ) [5], and our choice of local unitaries is thus optimal for pure states. The negativity (10) is depicted by a solid red line in Fig. 2.

B. Unbiased mixtures of coherent states
These Gaussian states are of the form and can be prepared by splitting a thermal state with mean number of thermal photons 2σ 2 on a balanced beam splitter. Here α ∈ C, P(α) = exp(−|α| 2 /σ 2 )/(πσ 2 ) and d 2 α = d(Reα)d(Imα). The states are already in standard form with a CM (9) specified by a = b = σ 2 + 1/2 and c 1 = c 2 = σ 2 .
Making use of the components of a coherent state in Fock basis m|α = exp(−|α| 2 /2)α m / √ m! we get the following matrix elements of the state (11), where s( j) = σ 2 1/σ 2 + 2 j+1 . By substitution of the latter expression into Eq. (8) we get after some algebra The negativity (13) is depicted by a dashed blue line in Fig. 2, and is generally smaller than the one of pure states calculated in (10). Both classes of Gaussian states have a nonzero negativity of quantumness which increases with n > 0; this is in agreement with earlier studies of nonclassical correlations based on entropic measures of quantum discord [11,12].

C. Standard-form two-mode Gaussian states
In general we need the Fock basis elementsρ m,n for an arbitrary two-mode Gaussian state with zero first moments. Combining the results of Refs. [21][22][23] we can express them as whereγ is the CM of the state, 1 1 is the 4 × 4 identity matrix, and H (R) m 1 ,m 2 ,n 1 ,n 2 (0) is the four-dimensional Hermite polynomial [24] at the origin; see Appendix B for a complete derivation of Eq. (14). Here is the symmetric matrix defining the polynomial, where and For the standard-form CMγ, Eq. (9), we get in particular with 1 1 2 being the 2 × 2 identity matrix, and d j = (a + 1/2)(b + 1/2) − c 2 j ( j = 1, 2). One can then evaluate the negativity (8) by performing a numerical summation of the absolute values of the elements (14). The higherorder Hermite polynomials can be calculated from the lowerorder ones by using e.g. the recurrence formula derived in Appendix B.
We remark that the compact expression in equation (14) is of independent interest and can be useful for the characterization of hybrid information processing involving conversion between continuous and discrete variable entanglement [25], or particularly for studies of Bell nonlocality of arbitrary twomode Gaussian states by means of dichotomic pseudospin measurements [26], whose expectation value can be conveniently evaluated at the Fock space level.
In the context of the present paper, apart from the utility for numerical evaluation of the output negativity (8), equation (14) also enables us to derive a simple analytical lower bound on the output negativity. The bound results from the following chain of inequalities where the first inequality follows from the inequality 1/ √ n! ≥ 1/n! which holds for any n ≥ 0, the second inequality is a consequence of the triangular inequality for absolute values, and the last equation follows from the expression for the generating function of the four-dimensional Hermite polynomials at the origin [24], where h = (α * 1 , α * 2 , α 1 , α 2 ) T and R is the matrix (15). A comparison between the right-hand side (RHS) of the previous equation and the expression of the Husimi Q-quasiprobability distribution Φ A (α 1 , α 2 ) = α 1 α 2 |ρ|α 1 α 2 /π 2 in the Fock basis further yields as can be easily seen from the results of Appendix B. Therefore, the last expression in the chain of inequalities (20) can be written in the following compact form Now, making use of the inequalities (20) and equality (23) one finds that the sum in (8) is lower-bounded as which finally gives the following bound on the output negativity (8) The bound (25) can be evaluated for any zero-mean twomode Gaussian state with CMγ by calculating the matrix (15) and substituting it into the formula (23). To test the tightness of the bound we calculate it for the previous examples of pure states and mixtures of coherent states, and compare the obtained lower bounds with the exact values of the negativities (10) and (13), respectively. The CMγ is in the standard form (9) in both cases and therefore one can evaluate easily the matrix (15) using Eqs. (18) and (19) which gives, after substitution into Eq. (23), (πe) 2 Φ p A (1, 1) = e 2 tanh r cosh 2 (r) (26) for pure states, and for unbiased mixtures of coherent states. The corresponding negativities then satisfy N p ≥ 1 2 e 2 tanh r cosh 2 (r) − 1 ≡ L p (28) and The bounds L p and L m as well as the negativities N p , Eq. (10), and N m , Eq. (13), are depicted in Fig. 2. The figure shows that both bounds are tight in the region of small n (see the inset), which also proves that Eq. (13) amounts to the exact value of the negativity of quantumness for mixtures of coherent states with small mean number of thermal photons in each mode. Both lower bounds are then shown to increase with increasing n and the gap between the bounds L p,m and the numerically evaluated values of the output negativities N p,m gets larger. Further analysis reveals however that the lower bounds L p and L m are nonmonotonic for larger n ; they both attain a maximum at n ≈ 0.62 and n ≈ 0.52, respectively, and then both monotonically decrease for larger values of n ; eventually, both lower bounds become trivial as they enter the region of negative values, namely L p < 0 for n 5.26 and L m < 0 for n 1.97. As a final remark, note that the sum in negativity (8) just amounts to the so-called 1 -norm of the density matrixρ [8], i.e., ∞ m,n=0 |ρ m,n | = ρ 1 . The results of the present Section thus also describe how to calculate numerically the 1 -norm for an arbitrary two-mode Gaussian state with zero means and the inequality (24) gives a simple analytical lower bound ρ 1 ≥ (πe) 2 Φ A (1, 1) on such a norm.

V. CONCLUSIONS
We have shown that a protocol capable of activating nonclassical correlations in bipartite Gaussian states based solely on Gaussian operations cannot exist. We have also constructed a non-Gaussian activation protocol and we have investigated quantitatively its performance using the negativity of quantumness as a figure of merit. Our analysis suggests that optimal performance of the protocol is achieved if the input Gaussian state is in the standard form. Restricting to the local Gaussian unitaries the conjecture can be proved or disproved with the help of Eq. (14) by numerical minimization of the negativity (8) with respect to the unitaries, which is left for further research.
We believe that our results will stimulate further exploration of the negativity of quantumness and its interplay with other nonclassicality indicators [11,13] in the context of Gaussian states. This section is dedicated to the proof that a bipartite Gaussian state ρ AB of an N-mode subsystem A and an M-mode subsystem B is classically correlated across the A|B splitting if and only if it is a product state ρ A ⊗ ρ B .
The proof of the "only if" part is trivial because any product state is diagonal in the product of eigenbases of local states.
The "if" part can be proved using the necessary and sufficient condition for zero quantum discord [17]. Quantum discord D B (ρ AB ) of a quantum state ρ AB with a measurement on subsystem B is zero if an only if the state can be expressed as [27] where {|i B } is an orthonormal basis of subsystem B. The zerodiscord criterion [17] then says that a quantum state ρ AB can be expressed in the form (A1) if and only if for an informationally complete positive operator valued measurement (IC-POVM) on subsystem A, the conditional states ρ B|k of subsystem B corresponding to the measurement outcomes k, mutually commute, i.e., [ρ B|k , ρ B|k ] = 0, for all k and k .
We consider a Gaussian state ρ AB with zero means and covariance matrix (CM) γ AB . Modes A 1 , A 2 , . . . , A N comprising the subsystem A are subject to a Gaussian measurement characterized by a CM γ m and a vector of measurement outcomes k = (x A 1 , p A 1 , . . . , x A N , p A N ) T ∈ R 2N . If a measurement outcome k occurs then the state ρ AB collapses into the M-mode state ρ B|k of subsystem B with CM σ and vector of first moments d k of the form [28] where A, B and C are blocks of the CM γ AB expressed with respect to the A|B splitting, As in Ref. [17] we will now express criterion (A2) in terms of the characteristic function. For this purpose we will first use the fact that an M-mode quantum state ρ j can be expressed as [28] where C j (ξ) is the characteristic function of the state ρ j and W(ξ) = exp(−iξ T r) is the displacement operator with ξ = (ξ x 1 , ξ p 1 , . . . , ξ x M , ξ p M ) T ∈ R 2M and r = (x 1 , p 1 , . . . , x M , p M ) T is the vector of quadratures. Due to the validity of the relation Tr W † (ξ )W(ξ) = (2π) M δ(ξ−ξ ) we get from Eq. (A6) immediately the following expression for the characteristic function of the state ρ j : Making use of Eq. (A6) we can express the commutator on the left-hand side (LHS) of Eq. (A2) as where C k (ξ) and C k (ξ ) are the characteristic functions of the states ρ B|k and ρ B|k , respectively, and where we have used the relation with From Eqs. (A3) and (A4) it follows that the Gaussian states ρ B|k and ρ B|k possess the same CM σ and first moments d k and d k , respectively, and therefore their characteristic functions read Equation (A8) allows us to calculate the characteristic function of the commutator [ρ B|k , ρ B|k ] given by By inserting the RHS of the commutator from Eq. (A8) into Eq. (A12), using Eq. (A9) and carrying out the integration, we arrive at the characteristic function (A12) in the form From Eq. (A12) and the formula it follows that [ρ B|k , ρ B|k ] = 0 if and only if C kk (ξ) = 0 for all ξ. Assuming that the CM σ in Eq. (A3) has finite second moments and the measurement outcomes k and k and hence also the displacements d k and d k defined by Eq. (A4) are finite, the condition C kk (ξ) = 0 for all ξ is equivalent to the condition which can be rewritten using Eq. (A4) as Previous results allow us to rephrase the zero-discord criterion of Ref. [17] for bipartite Gaussian states and Gaussian IC-POVMs as follows. An N + M-mode Gaussian state ρ AB can be expressed in the form (A1) if and only if the condition (A16) is satisfied for all k, k , where k and k are measurement outcomes of an Gaussian IC-POVM on subsystem A characterized by the CM γ m . Condition (A16) is satisfied for all k, k (k k ) if and only if the matrix Consider now the heterodyne measurement which is an example of a Gaussian IC-POVM [29]. Then γ m = (1/2)1 1, the matrix 1 A+γ m is invertible and therefore condition (A17) is equivalent with the equation Cσ −1 Ω = 0. As both the matrices Ω and σ −1 are also invertible the latter condition is equivalent with the condition C = 0. For the heterodyne detection the condition (A17) is thus equivalent with the vanishment of the off-diagonal block C given in Eq. (A5), which carries intermodal correlations. This means in other words, that a bipartite (N + M)-mode Gaussian state can be expressed in the form (A1) if and only if it is a product state, i.e., ρ AB = ρ A ⊗ ρ B .
Let us now move to the necessary and sufficient condition for a bipartite (N + M)-mode Gaussian state to be a classically correlated state. A quantum state ρ AB is classically correlated if and only if D A (ρ AB ) = D B (ρ AB ) = 0 [30], where D A (ρ AB ) is the discord of ρ AB for measurement on subsystem A. A quantum state ρ AB is therefore classically correlated if and only if it can be expressed simultaneously in the form (A1) and in the form According to the criterion given in [17] a quantum state ρ AB can be expressed in the form (A18) if and only if for an IC-POVM on subsystem B the conditional states ρ A|k of subsystem A corresponding to the measurement outcomes k mutually commute, i.e., [ρ A|k , ρ A|k ] = 0, for all k and k .
Like in the previous case we can express the latter condition in terms of a characteristic function. We can proceed exactly along the same lines as in the case of the commutator (A8) with the only difference that now we consider measurement on the M-mode subsystem B. Consequently, the formulas which we get for the present case of the commutator (A19) are obtained from the formulas derived in the context of commutator (A8) by the replacements A ↔ B, C ↔ C T of the blocks of the matrix γ AB and by the replacement M → N. Thus we find that the commutator (A19) vanishes if and only if C T = 0. Therefore, the condition C = 0 is necessary and sufficient for an (N + M)-mode Gaussian state to be classical, which concludes our proof.

Appendix B: Matrix elements of a Gaussian state in Fock basis in terms of Hermite polynomials
Our aim is to express the elements of a density matrix of a Gaussian state ρ of two modes A and B in the Fock basis.
Here and in what follows we assume that the state has all first moments equal to zero. The present derivation combines the results obtained in Refs. [21][22][23]. Firstly we express the elements of the density matrix in the basis of coherent states as e |α 1 | 2 +|α 2 | 2 α 1 α 2 |ρ|α 1 α 2 (B1) The matrix element on the LHS of Eq. (B1) can be further expressed as where Φ A (α 1 , α 2 ) = 1 is the Husimi Q-quasiprobability distribution of the Gaussian state ρ [31]. Here, α = (α 1 , α * 1 , α 2 , α * 2 ) T and γ (c) A is the complex CM corresponding to antinormal ordering of the canonical operators. Substituting now from Eq. (B3) into the LHS of Eq. (B1) and making use of Eq. (B4) we arrive at the fol- The LHS of the latter equation can be expressed in terms of the multi-dimensional Hermite polynomials [24]. Specifically, the generating function of the four-dimensional Hermite polynomials is A can be expressed as where 1 1 is the 4 × 4 identity matrix, is a 4 × 4 unitary matrix, and γ is the standard real symmetrically ordered CM of the state ρ, with elements γ i j = r i r j + r j r i /2, i, j = 1, . . . , 4, where r i is the i-th component of the vector of quadratures r = (x A , p A , x B , p B ) T . Hence we get Furthermore, we can write where Consequently,