A Analytical expression for the covariance matrix

# Gaussian tripartite entanglement out of equilibrium

## Abstract

The stationary multipartite entanglement between three interacting harmonic oscillators subjected to decoherence is analyzed in the largely unexplored non-equilibrium strong dissipation regime. We compute the exact asymptotic Gaussian state of the system and elucidate its separability properties, qualitatively assessing the regions of the space of parameters in which fully inseparable states are generated. Interestingly, the sharing structure of bipartite entanglement is seen to degrade as dissipation increases even for very low temperatures, at which the system approaches its ground state. We also find that establishing stationary energy currents across the harmonic chain does not correspond with the build-up of biseparable steady states, which relates instead just to the relative intensity of thermal fluctuations.

###### pacs:
03.65.Yz, 03.67.Mn, 03.67.Bg, 42.50.Lc

## I Introduction

Entangled states of continuous variable (CV) systems have come to occupy a prominent position in quantum technologies Braunstein and van Loock (2005) for both experimental and theoretical convenience. On the experimental side, the high degree of control in the preparation, manipulation and measurement of Gaussian CV states Weedbrook et al. (2012) in a range of quantum physical supports including optical cavities, trapped ions (Leibfried et al., 2003) or nanomechanical devices (Ekinci and Roukes, 2005), makes them ideal for the efficient implementation of quantum information protocols. In particular, entangled CV multipartite Gaussian states are a valuable resource for communication schemes involving many parties (5); (6); (7), whose quantum-enhanced performance has been already demonstrated in experiments (8); (9).

This outperformance over classical protocols crucially relies on the amount and distribution of the entanglement shared by the multiple ‘modes’, which makes the precise quantification of multipartite entanglement a matter of paramount importance. The general assessment of entanglement even in low-dimensional quantum systems remains an open and challenging problem to date Horodecki et al. (2009); Bennett et al. (2011) and yet tremendous progress has been made towards its characterization in the CV Gaussian multipartite scenario Giedke et al. (2001); Adesso and Illuminati (2007); Adesso et al. (2012). This fact, combined with the simple mathematical description that CV multi-mode Gaussian states enjoy, further highlights their practical convenience.

Unfortunately, entanglement is very fragile to the unavoidable decorrelating external environments and therefore, the successful implementation of quantum technologies with CVs should start with a complete understanding of noise and dissipation, so that they may be avoided or eventually engineered to protect quantum coherences. In this line, a number of recent works have extensively analyzed the dynamics and asymptotic properties of bimodal entanglement in CV Gaussian states under realistic models of noise and dissipation Eisert et al. (2004); Plenio et al. (2004); Paz and Roncaglia (2008); (19); (20); (21); (22); (23); (24); Ludwig et al. (2010); Correa et al. (2012). Concretely, the stationary two-mode entanglement under weak correlated and uncorrelated local noise was addressed in Paz and Roncaglia (2008); (20); (21) for identical oscillators, and in (19); (22) for the non-resonant case. Moreover, the problem may be solved exactly once one abandons the assumption of weak interaction between system and environment, thus making it possible to probe into the strongly non-Markovian and non-equilibrium regimes Ludwig et al. (2010); Correa et al. (2012). In contrast, much less is known about noise and dissipation in the CV Gaussian multipartite scenario (27); Ferraro and Paris (2005); Adesso et al. (2006); Xiang et al. (2009); xiang Li et al. (2010); Li et al. (2011); (33) where, to our knowledge, all available results are limited by either the weak dissipation or equilibration assumptions.

The present paper aims to study multipartite stationary entanglement in the little-studied non-equilibrium strongly dissipative regime, through the extension of the exact techniques of Correa et al. (2012). We focus on the stationary Gaussian states that result from the contact of an interacting three-mode CV system with three local structured heat baths. A rich physical picture is gained by preparing the baths at generally different equilibrium temperatures, thus inducing steady-state energy transport. Endowed with all the versatility of an exact unconstrained stationary solution, we address the question whether robust tripartite entangled states may be generated out of equilibrium. As we shall see below, we can answer in the positive.

More precisely, we take three (generally non-resonant) modes arranged in an open chain with linear nearest-neighbour interactions and locally dissipating into uncorrelated Ohmic baths. We are then able to compute the exact Gaussian steady state of the system, under the sole assumption of initially uncorrelated system and environmental degrees of freedom Weiss (1999). Our model is particularly suited for the theoretical description of tripartite CV systems in which thermal relaxation is the main source of decoherence, as it may occur, for instance, to trapped ions in a Paul trap Brown et al. (2011) or clamped interacting nanomechanical oscillators Buks and Roukes (2002); Cleland and Roukes (2002).

Taking the exact steady state as starting point, we issue a comprehensive study of the tripartite entanglement distribution according to the classification introduced in Giedke et al. (2001). When the three equilibrium temperatures of the reservoirs are set to the same value and identical oscillators are considered, we observe the expected competition between decoherence and inter-oscillator coupling in the build-up of stationary tripartite entanglement. Most interestingly, we find limiting dissipation rates above which the ground state of the interacting oscillators switches from the ‘weak dissipation’ fully inseparable phase into a ‘strong dissipation’ bound entangled phase, passing through an intermediate two-mode biseparable stage. As we shall see, these changes in the entanglement-sharing structure occur as a consequence of the non-negligible renormalization effects introduced by the system-bath interaction, in spite of the vanishing thermal fluctuations.

Imposing a temperature gradient across the chain proves detrimental to the formation of robust fully inseparable states unless the system is set up in an asymmetrical configuration. Nevertheless, the resulting separability structure does not seem to depend on the stationary energy currents induced across the system, but rather, with the relative intensity of thermal fluctuations on each of the modes.

Finally, we discuss how the asymptotic tripartite entanglement may be enhanced with a suitable choice of parameters leading to well separated time scales for the thermal fluctuations and the free dynamics of the interacting modes.

This paper is organized as follows: We start by introducing the microscopic model for the system, the baths and their dissipative interaction in Sec. II. The reduced dynamics of the oscillators is tackled via the generalized quantum Langevin equation, introduced in Sec. III.1, and solved in the stationary regime in Secs. III.2 and III.3. For a detailed derivation of the closed formula of the exact steady state, the interested reader is directed to Appendix A. We then briefly review the classification criteria for tripartite entanglement in CV Gaussian states in Sec. IV, and apply them to the steady states of our system in Sec. V: The separability properties in the case of identical equilibrium temperatures are discussed in Sec. V.1, and the results on the steady-state entanglement under a temperature gradient are presented in Sec. V.2. Finally, in Sec. VI, we summarize and draw our conclusions.

## Ii The system

As already mentioned, our system consist of three quantum harmonic oscillators, labelled by after ‘left’, ‘center’ and ‘right’, respectively. They have bare oscillation frequencies and equal mass is assumed:

 HS0=∑αp2α2m+12mω2αx2α. (1)

Here and stand for the corresponding position and momentum operators. We connect the oscillators through a generic quadratic interaction term of the form

 HSI=12∑αβxαVαβ xβ, (2)

where are the entries of an Hermitian interaction matrix . In particular, we shall arrange the oscillators in an open chain with nearest neighbour interactions of strength connecting and , that is (see Fig. 1 below)

 \tensorsymV=⎛⎜⎝k−k0−k2k−k0−kk⎞⎟⎠. (3)

We address the local dissipation mechanism with the paradigmatic Caldeira-Legget model Caldeira and Leggett (1983); Weiss (1999). Therefore, three independent bosonic reservoirs are introduced, also labeled , comprised of non-interacting modes linearly coupled to their local oscillator with strength :

 HSB=∑αμ\mathpzcp2αμ2mαμ+12mαμω2αμ(\mathpzcqαμ−gαμmαμω2αμxα)2. (4)

Apart from the free Hamiltonian of the reservoirs and their linear interaction with the system (i.e. the terms of the form ), Eq. (4) also explicitly includes the renormalization term

 HR=∑αμg2αμ2mαμω2αμx2α, (5)

which is necessary in order to compensate the distortion exerted by the system-bath coupling on Weiss (1999). The effects of this term only start to become relevant as the system-baths interaction grows stronger Correa et al. (2012). The coupling constants define the spectral densities

 Jα(ω)≡π∑μg2αμ2mαμωαμ δ(ω−ωαμ), (6)

which receive a phenomenological functional form suitable for a correct description of dissipation. In particular, in Sec. III.3, we shall consider Ohmic spectral densities with Lorentz-Drude high frequency cutoff

 Jα(ω)=mγαω1+ω2/ω2c, (7)

where stands for the dissipation rate, and carries the order of magnitude of the system-bath interaction, and is the cutoff frequency, that places a lower bound in the characteristic time scale of the thermal fluctuations of the baths (40).

We initialize system and environment as , where is any state of the three oscillators, is a (Gaussian) thermal equilibrium state of reservoir at temperature , and where denotes the Boltzmann constant. The normalization factors are and stands for the free Hamiltonian of the corresponding reservoir. The linearity of the system’s effective dynamics, guaranteed by the overall linear Hamiltonian and the ‘Gaussianity’ of the baths, leads to Gaussian reduced stationary states (41).

Any Gaussian three-mode state is fully determined (up to local displacements) by its second order moments, arranged in the covariance matrix

 σ≡(\tensorsymC\vectorsymX\vectorsymX(0)\tensorsymC\vectorsymX\vectorsymP(0)\tensorsymC\vectorsymP\vectorsymX(0)\tensorsymC\vectorsymP\vectorsymP(0)). (8)

The blocks are defined as

 \tensorsymC\vectorsymA\vectorsymB(t−t′)≡12⟨\vectorsymA(t)\vectorsymBT(t′)+\vectorsymB(t′)\vectorsymAT(t)⟩ρ0, (9)

where and , are column vectors collecting position and momentum operators of the modes.

## Iii Exact stationary states

### iii.1 Generalized quantum Langevin equation

We shall now calculate the stationary matrices and thus, the steady state of the system, by making use of the generalized quantum Langevin equation (QLE) formalism Weiss (1999), which is widespreadly used in the study of quantum Brownian motion Hanggi and Ingold (2005). The QLE follows from the elimination of the environment in the Heisenberg equations of motion for and , and may be compactly written as

 \tensorsymM¨\vectorsymX+\tensorsymϕ\vectorsymX=\vectorsymη(t)+1ℏt∫−∞dτ \tensorsymχ(t−τ)\vectorsymX(τ). (10)

Note that this equation does not rely on any approximations and therefore, it remains valid in all regimes of parameters. We remark as well that we took the initial condition at so that for any finite , it already describes the asymptotic properties of the system.

The matrix is diagonal and carries the masses of the oscillators , where stands for Kronecker delta. The effective potential is encoded in , where the frequency shift

 m ΔΩα≡1π∫∞0dω Jα(ω)ω, (11)

directly follows from the renormalization term of Eq. (5).

In addition to the free dynamics of the interacting oscillators, Eq. (10) also accounts for decoherence: On the one hand, the oscillators are locally driven by the stochastic quantum forces that enclose the effects of the thermal noise. These form the column vector . On the other hand, the last term on the right-hand side stands for a ‘friction memory kernel’ or ‘generalized susceptibility’ and describes dissipation. Since the three baths are uncorrelated, the susceptibility matrix has elements

 χαβ(t)≡δαβ Θ(t) 2ℏπ∫∞0dω Jα(ω)sinωt, (12)

where stands for the Heaviside step function. Thermal noise and friction are connected via the Kubo relation

 \tensorsymχ(t−t′)=−iΘ(t−t′)⟨\vectorsymη(t)\vectorsymηT(t′)−\vectorsymη(t′)\vectorsymηT(t)⟩B, (13)

where denotes an average over the environmental degrees of freedom.

### iii.2 Formal stationary solution

Quite generically, the matrices may be extracted from Eq. (10) by taking its Fourier transform . One thus arrives to the linear expression

 ~\vectorsymX(ω)=\tensorsymα(ω)~\vectorsymη(ω), (14)

where the complex matrix is defined as

 \tensorsymα(ω)≡−(ω2\tensorsymM−\tensorsymϕ+1ℏ ~\tensorsymχ(ω))−1, (15)

and the Fourier transform of the generalized susceptibility matrix has elements such that

 −Im ~χαα(ω)ℏ=Jα(ω)Θ(ω)−Jα(−ω)Θ(−ω). (16)

The causality argument that renders also ensures that is analytic in the upper-half plane of complex frequencies Weiss (1999). By virtue of the Kramers-Kronig relations we then have

 Re ~χαα(ω)=P∫∞−∞dω′π Im ~χαα(ω′)ω′−ω, (17)

where stands for the principal value of the integral. Let us now introduce the notation

 Γα(ω)≡−Im ~χαα(ω)ℏcothℏω2kBTα, (18)

for the symmetrized power spectrum of the quantum stochastic force Ludwig et al. (2010); Correa et al. (2012), and the vector . Then, the matrix writes as

 \tensorsymC\vectorsymX\vectorsymX(t)=ℏ∫dω2πe−iωt\tensorsymα(ω)\tensorsymΓ(ω)\tensorsymα(−ω)T, (19)

while the remaining correlations are:

 \tensorsymC\vectorsymP\vectorsymP(t)=ℏm2∫dω2π ω2e−iωt\tensorsymα(ω)\tensorsymΓ(ω)\tensorsymα(−ω)T, (20)

and

 \tensorsymC\vectorsymX\vectorsymP(t)=iℏm∫dω2π ωe−iωt\tensorsymα(ω)\tensorsymΓ(ω)\tensorsymα(−ω)T. (21)

Eqs. (15)-(21) thus formally provide the desired exact stationary states of the system for arbitrary spectral densities .

### iii.3 Stationary solution for Ohmic baths

As already anticipated, in order to compute the steady state from Eqs. (15)-(21), we will restrict ourselves to the Ohmic spectral densities of Eq. (7) and further assume symmetric dissipation rates . In this case, reduces to

 ~χαβ(ω)=δαβ mℏγω2ciω−ωc, (22)

which gives and by immediate substitution into Eqs. (15) and (18). Note that the frequency shift of Eq. (11) is now .

It is indeed possible to carry out the integration in Eqs. (19)-(21) and get closed formulas for the exact correlations by means of contour integration in the plane of complex frequencies, as in (42). Unfortunately, little can be gained from the cumbersome expressions that result, neither from the physical, nor from the practical point of view. Their discussion is hence postponed until Appendix A, and in what follows, we shall evaluate of Eqs. (19)-(21) numerically.

In the next section, we briefly review the basic tools to be employed in the characterization of the entanglement distribution in the stationary states of our system.

## Iv Gaussian tripartite entanglement

As already mentioned, the precise quantification of genuine multipartite entanglement in general mixed states, still proves challenging Bennett et al. (2011); Horodecki et al. (2009) even in the simplest case of tripartite systems. For instance, when dealing with qubits, quantities that prove to be bona fide measures in the bipartite scenario, such as the concurrence (43) or the negativity Vidal and Werner (2002), have to be replaced with a suitable entanglement mononotone that additionally satisfies the Coffman-Kundu-Wootters (CKW) monogamy inequality, like the residual tangle, computed from the convex roof of the squared concurrence (45).

In complete analogy, a continuous variable residual tangle, or (Gaussian) cotangle, was introduced in Adesso et al. (2006) that satisfies the CKW inequality for all three-mode Gaussian states. It follows from the infimum of the squared logarithmic negativity Vidal and Werner (2002) taken over all possible (Gaussian) pure-state decompositions of . Alternatively, a monogamous Gaussian entanglement measure may also be defined in terms of the Rényi-2 entropy Adesso et al. (2012).

As a bipartite entanglement measure, the (logarithmic) negativity exploits the positivity-of-the-partial-transpose (PPT) separability criterion (46); (47) which turns out to be not only necessary, but also sufficient for all multi-mode Gaussian states (48). Therefore, even if the (logarithmic) negativity fails to faithfully account for genuine multipartite correlations, the PPT criterion does allow for a qualitative description of the distribution of Gaussian entanglement in a three-mode CV system, according to the number of non-separable bipartitions out of the three possible. We shall denote them as , and . This entails the following classification for tripartite Gaussian states, as introduced in Giedke et al. (2001):

1. Fully inseparable states, that are not separable in any of the bipartitions.

2. One-mode biseparable states, which are separable only in one out of the three possible bipartitions.

3. Two-mode biseparable states, for which now two of the bipartitions are separable.

4. Three-mode biseparable or bound entangled states, which are separable under all bipartitions, but cannot be written as a mixture of product states only.

5. Fully separable states, that unlike those of C4, can be written as a mixture of product states.

In order to distinguish between the PPT-equivalent classes C4 and C5, we make use of the criterion for full separability of Giedke et al. (2001). In what follows, rather than attempting to quantify genuine tripartite entanglement, we resort to the previous qualitative characterization and apply it to the exact stationary states of our system.

## V Results and discussion

Finally, we are in a position to analyze the distribution of the stationary tripartite entanglement classes in the space of parameters of the system. Even if Eqs. (15)-(21) are not underpinned by any restrictive assumptions, we shall focus on the low temperature regime, which is optimal for the build up of entanglement, and exploit our steady-state solution to probe into the strongly dissipative regime.

We shall also restrict to low effective inter-oscillator coupling strengths , as strong couplings are rather unrealistic in experiments. This translates into , where . Indeed, by noting that is a characteristic time for energy transport across the system when isolated from the environment, it becomes clear that the condition amounts to a separation of time scales that renders transport inefficient. Consequently, the typical time scale governing the closed evolution of the whole interacting system may be approximated as .

In the study of quantum Brownian motion, one usually assumes fast thermal fluctuations (, ) as compared with the free evolution and the dissipation time (40). On the contrary, we shall work with relatively low temperatures and strong dissipation rates (, ) so that the system is much more insensitive to noise. In this regime, picking a cutoff frequency of the order of gives rise to non-perturbative renormalization frequency shifts that should be expected to become relevant. It is also important to note that under strong dissipation, the stationary states of the system are generally not of thermal equilibrium (Gibbs states) (41); (42); Haake and Reibold (1985), even when the temperatures of the local baths coincide and no steady-state energy transport is established.

Under these conditions, the stationary tripartite entanglement is studied in absence of energy currents through the system (Sec. V.1), and when the equilibrium temperatures of the baths are arranged in a gradient (Sec. V.2).

### v.1 Identical equilibrium temperatures

We shall start by taking resonant frequencies and (see caption of Fig. 1). The tripartite entanglement class of the resulting stationary states is plotted in Fig. 2 as a function of the coupling strength and the equilibrium temperatures of the baths. Not surprisingly, the higher the temperatures, the higher the corresponding coupling that is required to keep the system in a fully inseparable state (C1). Note as well that one mode biseparable states (C2) do not build up asymptotically in this configuration.

In fact, the stationary entanglement in the bipartition proves more resilient to noise than in either or . This is obviously due to our choice of potential in Eq. (3), that only puts mode in direct interaction with the remaining two. Now, given that in this configuration the system is invariant under the exchange , its stationary states must be bisymmetric and, therefore, as the temperatures increase, steady-state entanglement in bipartitions and must disappear jointly, which entails a direct transition from C1 to C3. Increasing the temperatures further, the system also becomes separable with respect to , thus giving rise to stationary bound entangled states (C4). Even though class C5 only appears for extremely low coupling in Fig. 2, at any given there exist a temperature above which the steady states become fully separable Anders (2008).

Most interestingly, in the inset of Fig. 2 we can see how the separability properties of the ground state (GS) of the chain depend on and : For any above a temperature-dependent threshold (in the figure ), there exist dissipation rates at which the GS undergoes transitions C1C3 and C3C4. On the contrary, for , it remains in the fully inseparable phase C5 regardless of the dissipation strength. The sharing structure of bipartite entanglement in the GS of a harmonic chain thus depends on when decohering far from the Born-Markov regime.

This can be, at least, qualitatively understood by recalling that the system Hamiltonian includes the renormalization term of Eq. (4), that amounts to a shift on the frequencies . Hence, one may argue that the effective coupling strength decreases as the dissipation rate grows, thus potentially downgrading the GS to an entanglement class of higher separability.

### v.2 Temperature gradient across the system

We now arrange the baths in a temperature gradient by allowing for (see Fig. 1) so that stationary energy transport may be established across the harmonic chain. Let us first consider , and . This configuration is invariant with respect to the combined exchange of and and thus, the distribution of entanglement phases must be symmetric about , as seen in Figs. 3 and 3.

In Fig. 3 we fix and plot the entanglement classes as a function of and . First, notice that one-mode biseparable stationary states (C2) do build up, now that the symmetry argument invoked in Sec. V.1 is not applicable.

One sees as well that in general, whenever increases, the free dynamics of the central mode becomes more insensitive to noise since decreases. This helps to reduce the stationary biseparability and eventually yields fully inseparable states (C1). However, as illustrated in the inset, very large values of may also cause an effective decoupling of the central mode from the rest as becomes smaller. In other words, given a fixed interaction , fully inseparable stationary states may be generated by tuning the frequencies to a compromise between shielding the system from thermal noise and keeping the effective interaction between its modes sufficiently strong.

Finally, note that arranging the baths in a temperature gradient proves detrimental to the asymptotic formation of states in any of the bipartite entangled classes (C1C3). This seems to occur due to the intensification of thermal noise at the hot end of the chain rather than as a consequence of the stationary energy currents established across the system. We illustrate this point further in Figs. 3 and 3, where and are taken as the free parameters.

In Fig. 3 we consider resonant modes (), while in Fig. 3 the oscillators are set up in the asymmetrical configuration: , and . In the first case, keeping the steady state within the fully inseparable class requires stronger couplings as the temperature gradient increases in either direction. On the contrary, the asymmetric setting of Fig. 3 favors the formation of class C1 steady states at moderate negative temperature gradients, as these provide the low frequency mode with the lowest temperature () and the high frequency mode with the highest one (), which optimally shields the system from thermal noise.

It is also noticeable how the one-mode biseparable class (C2) takes over bound entangled steady states (C4) in Fig. 3 as contrasted with Fig. 3, even though it may be seen that the magnitude of the stationary energy currents (54) are comparable in either case. This observation further suggests that the build-up of steady-state quantum correlations might indeed not share a causal relation with the efficient transport of energy at microscopic scale, as already pointed out in different contexts such as excitation transfer in biological systems (51), thermal conduction in spin chains Wu and Segal (2011) or the optimized performance of quantum refrigerators (53).

## Vi Conclusions

We have addressed the qualitative classification of the bipartite entanglement distribution across three linearly coupled harmonic oscillators dissipating into independent structured baths. By making use of the quantum Langevin equation formalism, we were able to compute their exact stationary Gaussian states and then, issue a comprehensive analysis of the different entanglement classes that build-up asymptotically in terms of the parameters of the system. It is important to remark that this approach is not limited by the customary assumptions of equilibrium and/or weak-memoryless system-bath interactions, so that it allows to probe into the largely unexplored non-equilibrium strong dissipation regime.

Interestingly, we saw how the ground state of the harmonic chain undergoes structural transitions between different schemes of entanglement sharing, increasing its bipartite separability as the dissipation grows stronger. This is a direct consequence of the non-negligible back action of the system-bath coupling on the system itself.

It was also noted that inducing stationary energy transport by means of a temperature gradient is generally detrimental to the formation of fully inseparable steady states due to the more intense thermal fluctuations at the hot end of the system. The resulting stationary energy currents do not seem to correlate to the asymptotic formation of biseparable states.

We finally discussed how a suitable choice of frequencies may shield the system from thermal noise while keeping the effective inter-oscillator coupling strong enough, so that potentially useful fully inseparable states may build up asymptotically in spite of the strong decoherence.

As it was already pointed out, our model is appropriate for the theoretical description of a range of systems of interest in quantum technologies, especially arrays of interacting nanomechanical resonators. Indeed, considering typical frequencies in the range of MHz and masses around kg, the region of the space of parameters probed in our numerics may be achieved in present-day experiments.

One could also think of applying the powerful exact techniques illustrated here to the study of steady-state multipartite entanglement under the action of correlated thermal noise in a more realistic structured bath of spatial dimension greater than one. This problem is worthy of detailed study and will be considered elsewhere.

## Acknowledgments

The authors warmly thank A. Ruiz for reading and extensively commenting on the manuscript, and J. P. Palao, G. Adesso, D. Girolami, G. De Chiara, S. Kohler and N. García Marco for fruitful discussions and helpful criticism. L.A.C. wants to thank N. Ramos García in memoriam for his unconditional support through the years. This project was funded by the Spanish MICINN (Grant No. FIS2010-19998) and the European Union (FEDER). A.A.V. and L.A.C. acknowledge the Canary Islands Government for financial support through the ACIISI fellowships (85% cofunded by European Social Fund).

## Appendix A Analytical expression for the covariance matrix

As already mentioned in Sec. III.3, in order to get an analytical expressions for e.g., Eq. (19), one can use the customary toolbox of complex analysis to explicitly carry out the integration. Therefore, complete knowledge about the roots of the denominator of the integrand is required. Let us start by alternatively writing as

where the matrix is defined as . The notation stands for the adjugate matrix of , and the asterisk represents conjugate transposition. Note that from Eq. (22) it follows that .

The denominator of Eq. (23) is a real polynomial of degree eighteen comprised of the determinants and , which are complex polynomials of degree nine. Provided that is diagonalizable, may be written as the product of three polynomials of degree three, and therefore, its roots can be analytically worked out, even if the resulting expressions are rather involved. When it comes to the multiplicity of those complex roots, it can be checked that they are all simple for our choice of interaction potential in Eq. (3). We shall label them so that lie in the lower half plane of complex frequencies (and are their corresponding complex conjugates).

We may now decompose the integrand of Eq. (23) into partial fractions as

 Missing or unrecognized delimiter for \right (24)

with . We shall also make use of the identity

 cothx=1x+1iπ[ψ(1+ixπ)−ψ(1−ixπ)], (25)

where stands for the digamma or psi-fuction, i.e. the logarithmic derivative of Euler’s gamma function (55).

Combining Eq. (25) with (24), Eq. (23) may be evaluated by making the analytical continuation of the integrand into the plane of complex frequencies and calculating residues. Notice that the extended function has simple poles along the entire positive (negative) imaginary axis. We shall choose integration contours either in the lower or upper plane for each of the resulting terms in Eq. (23), such that those non-analyticities are avoided. The elements of the correlation thus result in

Similarly, may be computed from Eq. (20) to yield

 [\tensorsymC\vectorsymP\vectorsymP(0)]αδ=ℏγω2cm3∑β9∑j=1[kBTβ Re z2j2ℏ Im zj−2Re z3jπ Im zj Im ψ(1+iℏzj2πkBTβ)] Re adj[\tensorsymF(zj)]αβadj[\tensorsymF(zj)∗]βδ∏k≠j(zj−zk)∏k≠j(zj−¯¯¯zk), (27)

and finally, Eq. (21) translates into

 [\tensorsymC\vectorsymX\vectorsymP(0)]αδ=−ℏγω2cm4∑β9∑j=1[kBTβ Re zj2ℏ Im zj−2Re%  z2jπ Im zj Im ψ(1+iℏzj2πkBTβ)] Im adj[\tensorsymF(zj)]αβadj[\tensorsymF(zj)∗]βδ∏k≠j(zj−zk)∏k≠j(zj−¯¯¯zk), (28)

which provides us with the desired explicit formulas for the exact stationary Gaussian state of the system.

### References

1. S. L. Braunstein and P. van Loock,  Rev. Mod. Phys. 77, 513 (2005).
2. C. Weedbrook, S. Pirandola, R. García Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro,  and S. Lloyd,  Rev. Mod. Phys. 84, 621 (2012).
3. D. Leibfried, R. Blatt, C. Monroe,  and D. Wineland, Rev. Mod. Phys. 75, 281 (2003).
4. K. Ekinci and M. Roukes, Rev. Sci. Instrum. 76, 061101 (2005).
5. P. van Loock and S. L. Braunstein, Phys. Rev. Lett. 84, 3482 (2000).
6. P. van Loock and S. L. Braunstein, Phys. Rev. Lett. 87, 247901 (2001).
7. G. Adesso, A. Serafini, and Fabrizio Illuminati, New J. Phys. 9, 60 (2006).
8. H. Yonezawa, T. Aoki, and A. Furusawa, Nature 431, 430 (2004).
9. S. Koike, H. Takahashi, H. Yonezawa, N. Takei, S. L. Braunstein, T. Aoki, and A. Furusawa, Phys. Rev. Lett. 96, 060504 (2006).
10. A. Ferraro, S. Olivares,  and M. Paris, Gaussian states in quantum information (Bibliopolis, 2005).
11. R. Horodecki, P. Horodecki, M. Horodecki,  and K. Horodecki,  Rev. Mod. Phys. 81, 865 (2009).
12. C. H. Bennett, A. Grudka, M. Horodecki, P. Horodecki,  and R. Horodecki,  Phys. Rev. A 83, 012312 (2011).
13. G. Giedke, B. Kraus, M. Lewenstein,  and J. I. Cirac,  Phys. Rev. A 64, 052303 (2001).
14. G. Adesso and F. Illuminati, J. Phys. A 40, 7821 (2007).
15. G. Adesso, D. Girolami,  and A. Serafini,  Phys. Rev. Lett. 109, 190502 (2012).
16. J. Eisert, M. Plenio, S. Bose,  and J. Hartley, Phys. Rev. Lett. 93, 190402 (2004).
17. M. Plenio, J. Hartley,  and J. Eisert, New J. Phys. 6, 36 (2004).
18. J. P. Paz and A. J. Roncaglia,  Phys. Rev. Lett. 100, 220401 (2008).
19. J. P. Paz and A. J. Roncaglia, Phys. Rev. A 79, 032102 (2009).
20. R. Vasile, S. Olivares, M. G. A. Paris, and S. Maniscalco, Phys. Rev. A 80, 062324 (2009).
21. R. Vasile, P. Giorda, S. Olivares, M. G. A. Paris, and S. Maniscalco, Phys. Rev. A 82, 012313 (2010).
22. F. Galve, G. L. Giorgi, R. Zambrini, Phys. Rev. A 81, 062117 (2010).
23. F. Galve, L. A. Pachón, and D. Zueco, Phys. Rev. Lett. 105, 180501 (2010).
24. A. Wolf, G. De Chiara, E. Kajari, E. Lutz, and G. Morigi, Europhys. Lett. 95, 60008 (2011).
25. M. Ludwig, K. Hammerer,  and F. Marquardt,  Phys. Rev. A 82, 012333 (2010).
26. L. A. Correa, A. A. Valido,  and D. Alonso,  Phys. Rev. A 86, 012110 (2012).
27. M. M. Cola, M. G. A. Paris, and N. Piovella, Phys. Rev. A 70, 043809 (2004).
28. A. Ferraro and M. G. A. Paris,  Phys. Rev. A 72, 032312 (2005).
29. G. Adesso, A. Serafini,  and F. Illuminati,  Phys. Rev. A 73, 032345 (2006).
30. S. Xiang, B. Shao, K. Song,  and J. Zou, Phys. Rev. A 79, 032333 (2009).
31. G. xiang Li, L. hui Sun,  and Z. Ficek,  J. Phys. B 43, 135501 (2010).
32. J. Li, T. Fogarty, C. Cormick, J. Goold, T. Busch,  and M. Paternostro,  Phys. Rev. A 84, 022321 (2011).
33. G. Manzano, F. Galve, and R. Zambrini, Phys. Rev. A 87, 032114 (2013).
34. U. Weiss, Quantum dissipative systems, edited by W. Scientific, Series in Modern Condensed Matter Physics (World Scientific, 1999).
35. K. Brown, C. Ospelkaus, Y. Colombe, A. Wilson, D. Leibfried,  and D. Wineland, Nature 471, 196 (2011).
36. E. Buks and M. Roukes, J. Microelectromech. S. 11, 802 (2002).
37. A. Cleland and M. Roukes, J. Appl. Phys. 92, 2758 (2002).
38. A. O. Caldeira and A. J. Leggett, Ann. Phys. 149, 374 (1983).
39. P. Hanggi and G.-L. Ingold, Chaos 15, 026105 (2005).
40. H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, USA 2012).
41. H. Garbert, U. Weiss, and P. Talkner, Z. Phys. B 55, 87 (1984).
42. P. Riseborough, P. Hanggi, and U. Weiss, Phys. Rev. A 31, 471 (1985).
43. C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Phys. Rev. A 54, 3824â3851 (1996).
44. G. Vidal and R. F. Werner,  Phys. Rev. A 65, 032314 (2002).
45. V. Coffman, J. Kundu, and W. K. Wootters, Phys. Rev. A 61, 052306 (2000).
46. A. Peres, Phys. Rev. Lett. 77, 1413 (1996).
47. M. Horodecki, P. Horodecki, and R. Horodecki, Phys. Lett. A 223, 1 (1996).
48. R. F. Werner and M. M. Wolf, Phys. Rev. Lett. 86, 3658 (2001).
49. F. Haake and R. Reibold,  Phys. Rev. A 32, 2462 (1985).
50. J. Anders,  Phys. Rev. A 77, 062102 (2008).
51. M. B. Plenio and S. F. Huelga, New. J. Phys. 10, 113019 (2008).
52. L.-A. Wu and D. Segal,  Phys. Rev. A 84, 012319 (2011).
53. L. A. Correa, J. P. Palao, G. Adesso, and D. Alonso, Phys. Rev. E 87, 042131 (2013).
54. A. Dhar and D. Roy, J. Stat. Phys. 127, 801 (2006).
55. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover, New York 1972).
72407