Avoiding dissipation in a system of three quantum harmonic oscillators

We analyze the symmetries in an open quantum system composed by three coupled and detuned harmonic oscillators in the presence of a common heat bath. It is shown analytically how to engineer the couplings and frequencies of the system so as to have several degrees of freedom unaffected by decoherence, irrespective of the specific spectral density or initial state of the bath. This partial thermalization allows observing asymptotic entanglement at moderate temperatures, even in the nonresonant case. This latter feature cannot be seen in the simpler situation of only two oscillators, highlighting the richer structural variety of the three-body case. When departing from the strict conditions for partial thermalization, a hierarchical structure of dissipation rates for the normal modes is observed, leading to a long transient where quantum correlations such as the quantum discord are largely preserved, as well as to synchronous dynamics of the oscillators quadratures.


I. INTRODUCTION
Prevention of decoherence and dissipation in open quantum systems is a fundamental condition for the presence of quantum phenomena in warm macroscopic everyday world. While decoherence has been extensively studied from the early 1980s [1,2] to nowadays (for a list of reviews see for instance [3][4][5][6]), and different mechanisms to avoid it have been discussed in the last years, including strategies to engineer it for applications [7][8][9]. Furthermore, some macroscopic systems from marine algae [10] to metal carboxylates [11], can present quantum correlations at high temperatures. These recent experiments suggest that avoiding a complete quantum-to-classical transition induced by monitoring of the environment is not only artificially possible but can occur inherently in natural phenomena. The mechanisms that produce such survival or even construction of coherence at large time scales remain almost unclear, but different theoretical strategies have been proposed in order to predict it, mostly motivated in the context of quantum computation [12][13][14][15][16]. Indeed decoherence and the leak of information to the environment are the major obstacle in quantum processing of information and construction of quantum memories [17].
In this context, one of the strategies to bypass decoherence is exploiting dynamical symmetries in the systemenvironment interaction. In order to generate unitary evolution in a certain subspace of the Hilbert space of an open system, a common dissipation (where several units equally couple to the same environment) has been firstly used in a two-qubit system [18][19][20][21] and later extended to multiple qubits [22,23] and continuous variable systems [24-29, 31, 33]. A general framework has been developed with several contributions (see for example [15] and the references therein) agglutinating the main concepts of decoherence-free subspaces (DFSs), Noiseless subsystems (NSs) or more recently information-preserving structures (IPSs) [16] and proposing general methods to obtain them for arbitrary Hamiltonians. In our work we extend previous studies in the context of continuous variables exploring the vaster landscape offered when moving from two to three coupled harmonic oscillators in the search for NSs induced by non-trivial symmetries in the system.
Previous works on dissipative harmonic oscillators reported that in presence of identical frequencies and couplings between oscillators the symmetry of the collective motion and their interaction with the environment can lead to the effective decoupling of some of the normal modes of the system from the bath [24][25][26][27][28]. The cases of different frequencies [29,30] or couplings [31] is instead less studied and understood. The natural step of considering three harmonic oscillators beyond the symmetric configuration of identical oscillators already provides much more phenomenological richness while at the same time it allows for analytic treatment and gives valuable intuition when pursuing a further extension to the case of N oscillators.
We show in this work how to obtain NSs and bypass decoherence independently of the bath properties such as temperature or frequency distribution, considering different frequencies, couplings and boundary conditions, in the presence of a common bath. We also analyze how by using these NSs, quantum correlations like entanglement can persist (by two different mechanisms) in the asymptotic limit of the dynamical evolution, given in our case by a Markovian-type dynamics in the weak-coupling limit. Furthermore even for three different oscillators an homogeneous dynamics, such as collective quantum synchronization, is shown to be present in some parameters manifolds. Synchronization is a widespread phenomenon in physical, biological, chemical and social systems [32], only recently explored in the quantum regime [33,34]. Our results for the system of three oscillators in presence of common dissipation may be implemented with ions in linear Paul traps by following the proposal in [35]. Experimental realization of coupled harmonic oscillators appear in optical [36,37] and superconducting [38,39] cavities as well as trapped ions [40,41] or nanoelectromechanical resonators [42]. Three coupled elements architectures are also known to allow for isochronous synchronization of semiconductor lasers with delayed coupling or neuronal models [43].
The paper is organized as follows: in Sect. II we introduce the model for the system of three harmonic oscillators dissipating into a common bath, also in terms of the normal modes of the system. For certain particular values of the system's parameters (and independently of the bath characteristics), one of several normal modes can be protected from decoherence. We find analytically these conditions leading for NSs in Sect. III giving some specific cases that are analyzed in more detail. In the following we focus on the long time dynamics of the system considering in Sect. III A the creation and conservation of entanglement at asymptotic times. Section III B is dedicated to ability of the system to generate asymptotic synchronized states and exploring its connection to the robustness of quantum correlations. In section IV deviations of the NSs conditions are considered. This leads to dynamical relaxation of the system that converges towards a thermal state. We conclude with Sect. V summarizing our main results.

II. THE MODEL
We start with a Hamiltonian describing three coupled quantum harmonic oscillators with arbitrary frequencies and coupling constants between them. For simplicity we suppose unit masses and natural units ( = k B = 1): where p i and q i represents respectively the momentum and position operators of each harmonic oscillator ([q i , p j ] = iδ ij ) and N = 3 for our case. This equation is conveniently expressed in quadratic matrix form as: where I is the identity (N × N ) matrix, q T = (q 1 , ..., q N ) and H contains all the parameters of the system, i.e. the squared frequencies and couplings between oscillators. We will only consider H with positive eigenvalues, so as to have bounded states (attractive potential). The environment is introduced by equally coupling each oscillator of the system to the same thermal bath, which is described by an infinite collection of independent bosonic modes: where [X α , Π β ] = iδ αβ and masses are unitaries. We will use along the paper greek subindices to refer to bath modes while latin ones are reserved for system oscillators (i, j) and normal modes (k, n). The system-bath interaction reads: with a factorized form H I = S⊗B of an operator S acting only on the system's Hilbert space and B acting on the environment one. As usual, this type of interaction yields a renormalization of the frequencies that we may include directly in our model by changing The normal modes basis of the system (2) is obtained after a canonical transformation of the system operators through the orthogonal basis-change matrix F : which diagonalizes H (q T Hq = Q T ΩQ).
Here Ω = F T HF is a diagonal matrix containing the squared normal modes frequencies Ω n with n = 1, 2, ..., N . In this basis H S now represents the Hamiltonian for a N = 3 uncoupled harmonic oscillators, or normal modes, related with the original (natural) modes by F , so we can write: and the interaction Hamiltonian of Eq. (4) is now: We stress that even if the oscillators are coupled with the same strength to the bath, the couplings of the normal modes to the bath, κ n , are not homogeneous: The effective couplings κ n only depend on the canonical transformation, i.e. on the system's parameters and arrangement defined by H. This suggests a strategy to protect one or more normal modes from the environment action based on proper tuning of {ω i , λ ij }. Our analysis in section III addresses this point while deviations form the condition of vanishing effective coupling of a system normal mode and the environment are explored in section IV.
We mention that while this work focuses on the case of three coupled harmonic oscillators, the description in terms of effective couplings here provided is rather general and applies for arbitrary networks of N harmonic oscillators. The case of a common bath for all oscillators in the system corresponds to situations where the correlation length in the environment is larger than the system size. This assumption is not crucial for the framework we develop here, though any other choice would produce different specific analytic expressions. Finally, the case of a separate bath for each oscillator was shown in Ref. [33] to yield equal decoherence for all normal modes and therefore no NSs nor synchronization.

A. Unbalanced bath couplings
While the common bath situation is mostly a characteristic related to the coherence length of the bath as compared to the extension of the system, the equal coupling of each system oscillator to the latter might seem an arbitrary restriction. Imagine for example that each oscillator is at a different distance from the common heat bath, leading to an interaction where the different oscillators feel a coupling of strength γ i . The immediate consequence is that the effective couplings become Though in this paper we will consider γ i = 1, the unbalanced case would be solved following exactly the same procedure as we outline in the next section.

III. NOISELESS SUBSYSTEMS AND ASYMPTOTIC PROPERTIES
In this section we discuss the conditions to achieve noiseless subsystems where dissipation is avoided in one or two of the system's normal modes. The properties of our system are specified completely by the matrix H appearing in Eq. (2): and we aim to derive the set of conditions for the system parameters leading to one or two normal modes decoupled from the environment, i.e. whose effective coupling κ n is zero. Let us consider a normal mode δ with normal frequency Ω δ . The eigenvalue problem is expressed adequately by three equations: for the components of the normal vector C δ = (F 1δ , F 2δ , F 3δ ) T with F ij defined in Eqs. (5). The condition for normal mode δ to be non-dissipative (out of the bath' influence) leads to a constraint as follows From Eqs. (12), (13) and normalization condition we can obtain analytically C δ , Ω δ with a further constraint for the system parameters. In other words, not all parameter choices lead to NSs, but it is possible for some configurations of frequencies and couplings of the set of three oscillators (satisfying some constraint).
When Eq. (16) is fulfilled we obtain a NS composed by (at least) a single normal mode that is effectively uncoupled to the reservoir. This could be performed artificially by tuning one of the d = 6 parameters of H, such as, for instance, the natural frequency of one oscillator. In experiments where local addressing is possible, such as ions confined to individual traps, this modification should be rather straightforward (see e.g. [45]). It should be stressed that noise models for ion traps typically favor a SB interpretation in terms of fluctuating uncorrelated surface dipoles [46], though other microscopic models based on charge diffusion [47] in the electrode surface question whether the bath's correlation length could in fact be larger than the distance of the ion to the electrode. For the moment, this is an open problem.
Configurations in which a NS consisting of two normal modes is produced can also be obtained analytically in the specific case of three oscillators. Indeed, when two normal modes uncouple from the environmental action, the third one must necessarily coincide with the center of mass of the system. Explicitly, the condition for the center of mass being a normal mode is: (17) with C CM = (1, 1, 1) T / √ 3. The latter implication can be captured in the next two relations that have to be fulfilled simultaneously by the system parameters of the system: Since we want to remain in the domain of attractive potentials, we have to restrict ourselves to regions of the parameter space where Ω 2 CM = ω 2 1 + ω 2 3 − ω 2 2 + 2λ 13 > 0. In order to see the power of the conditions (16) and (18,19) we will give in the following some examples of configurations in which a NS of one or two modes is produced. We will consider for instance simpler situations in which the six-dimensional parameter space is reduced by assuming two of the three natural frequencies of oscillators to be equal, (ω ≡ ω 1 = ω 3 = ω 2 ) or in which two of the three couplings are equal (λ ≡ λ 12 = λ 23 = λ 13 ). This is sufficient to obtain some different scenarios appearing in open an closed chains configurations, as is schematically shown in Figure 1.
The presence of one or two non dissipating normal modes prevents full thermalization of the system and leads to an asymptotic state whose features are analyzed in the following, focusing on entanglement and on quantum synchronization between the oscillators.

A. Asymptotic Entanglement
When a NS is enabled, decoherence can be prevented in the system leading to asymptotic entanglement that would be absent in the thermal state. As a measure of entanglement, we will use the well known logarithmic negativity which is computable for bipartite Gaussian states [48,49] as is our case: where ν − is the minimum symplectic eigenvalue of the partial transposed covariance matrixṼ AB , corresponding to time reflection of one party. With the help of general expressions we can calculate analytically the asymptotic entanglement when a NS is produced. Here we present our results for the open chain with equal couplings (λ 12 = λ 23 ≡ λ) and same frequencies for the external oscillators (ω 1 = ω 3 ≡ ω), when only one of the three normal modes is not subjected to dissipation (Fig.  1a) and when only one of them is dissipating (Fig. 1e) by imposing λ =λ 0 . The details of the calculations are reported in Appendix A. As initial condition for the natural oscillators we choose a squeezed separable vacuum state given by: where any other first order or second order moments are zero.

One mode NS
As a paradigmatic example of the case in which there is one frozen normal mode, let us consider the configuration given in Fig. 1a. As for the initial condition, ω 1 = ω 3 ≡ ω in Eq. (24) and we will assume the same squeezing factor for the external pair, i.e. r 1 = r 3 ≡ r while the squeezing in the central oscillator r 2 will be irrelevant.
Normal modes coupled to the environment will reach an asymptotic thermal state whose variances are given in Eq. (A3) while the uncoupled one evolves freely. The asymptotic covariance matrix of the external oscillators is obtained by expressing the second order moments of natural oscillators in terms of the normal modes. Then we substitute respectively the asymptotic expressions corresponding to the frozen mode (not coupled to the bath) or the thermalized ones. This yields the following analytical expression for the entanglement: that is defined by a minimum value E 0 and an oscillatory term with amplitude ∆E and frequency 2ω: where r c = (r + 0 +r − 0 )/4 and the critical values are defined by the following expressions: Coefficients σ P and σ Q depend both on the bath's temperature and on the system parameters via the shapes and frequencies of the dissipative normal modes as can be seen in their definition in Eq. (A8). Note that while decoupling of normal modes from the bath is temperature independent, the amount of entanglement generated depends on it via the thermalized degrees of freedom. The presence of asymptotic entanglement between the external pair of oscillators in a symmetric chain (independent of the frequency of the central one, but depending on the temperature and initial squeezing) are illustrated in Figure 2. The minimum entanglement E 0 is plotted both for low (left panel) and high temperatures (right panel) in the relevant squeezing ranges. Different regions are distinguished in the map and are labeled following Paz and Roncaglia notation in Ref. [27]: sudden death is reached (SD), the asymptotic state consisting of an infinite sequence of sudden death and revivals (SDR) and finally when non-zero entanglement is present at all times (no sudden death, NSD).
An asymptotic entangled state with strictly E N > 0, can be generated both when r > r + 0 (> 2r c ) or equivalently when r < r − 0 (< 2r c ) with different origins. In the first case (r > r + 0 ) the entanglement oscillates between r − r + 0 and r + r − 0 and the initial squeezing in the natural oscillators is employed as a resource to generate an entangled state, while the bath contribution r + 0 is a source for its degradation. It's interesting to see that r + 0 is strongly dependent on the bath's temperature while the system parameters play a secondary role, only important at low temperatures. Indeed when temperature increases (T ≫ ω) sudden death of entanglement can be only avoided by increasing r to be greater On the other hand the amplitude of the oscillations in this case is ∆E = 2r c that has a very weak dependence on temperature quickly reaching a constant The second case (r < r − 0 ) is only present at low temperatures (of order 0.1ω). Here entanglement oscillates around r − 0 + r with amplitude 2r. This means that introducing no squeezing in the initial state leads to a constant entanglement at r − 0 , while adding squeezing (a re- source in the former case) makes entanglement tend to a SDR phase by widening its oscillatory amplitude. The fact that thermalization can lead to entanglement at low temperatures is well known [50].
Finally, we can relate critical values r ± 0 with the uncertainty induced by the environment in the virtual oscillatorq = (q 1 + q 3 )/ √ 2 position and momenta: this reveals that when r − 0 > 0 a squeezing in momentum is generated (∆p < ω 2 ) yielding entanglement as we have commented above. However note that we have never a minimum uncertainty state with r + 0 > r − 0 for all temperatures and physical regimes of system parameters. Indeed the uncertainty relation can be expressed for the virtual oscillatorq as ∆q∆p The quantity r c = (r + 0 +r − 0 )/4 can be also related with virtual oscillator uncertainties in position and momenta as e −rc = ∆p ω∆q < 1. In the left panel of Figure 2 we can see the two regions in which E 0 > 0 (NSD phases): the big one at the left top corner corresponding to entanglement generation by using the initial squeezing in the external oscillators as a resource (once r > r + 0 ) and the small left bottom island that represents the environment yielding entanglement via the squeezing generated in theq coordinate when r < r − 0 . The SDR phase is centered around 2r c (white dotted line) for low temperatures and their amplitude is given by the separation of the dashed colored line (r + 0 − r − 0 )/2 from the zero squeezing axis. For temperatures greater than that for which r − 0 = 0 (cross point between the dotted and dashed lines) they interchange their roles acting now (r + 0 − r − 0 )/2 (dashed colored line) as the center of the SDR region and 2r c as the amplitude. The SD phase is bounded by the quantity r + 0 − 4r c = −r − 0 corresponding to the case in which E 0 + 2∆E < 0 thus no entanglement is present in the asymptotic limit. For high temperatures (right panel) we can see how 2r c reaches a constant value while the pronounced curvature in the SDR region reveals that we can always obtain robust entanglement by increasing the squeezing parameter r logarithmically with temperature.
We have to point out that our results resemble those obtained in Ref. [27] for two resonant harmonic oscillators. There a similar entanglement phases diagram has been found and the same two different mechanism for entanglement generation appear. In that context, both oscillations and the appearance of the low temperatures NSD phase were attributed to non-markovian effects, while here follow by considering a final asymptotic (Gibbsian) state for the normal modes coupled to the bath (that can be reproduced by a Markovian Lindbland dynamics as is pointed in appendix B). Presence of a third oscillator in the system allows for manipulation of the width of entanglement phases at low temperatures (specially the low squeezings NSD one) by tuning the free system parameters ω 2 and λ.

Two modes NS
Let us now consider the case in which two modes become decoupled from the bath and in particular focusing on the configuration in Fig. 1e. This is indeed a symmetric open chain configuration as before, but now we have a special value of the couplings λ =λ 0 leading to a larger NS. The calculation is similar to the previous one while now we have that only one of the normal modes thermalizes and the other two have a free evolution decoupled from the bath. This leads to a less compact expression for asymptotic entanglement and more details are reported in Appendix A. Still, a similar phase diagram for entanglement can be found in this case by numerical evaluation of logarithmic negativity from Eq. (A9) given in the appendix A. Results are shown in Fig. 3 in the same range of squeezing and temperatures as in the previous (one mode NS) case. For low temperatures (left panel) the low temperatures low squeezing NSD island of Fig. 2 that corresponds to the environment acting as a resource for entanglement generation disappears, since the bigger one expands to low squeezing. Degradation of resources by the environmental action here is not sufficient to prevent entanglement production even in the non-squeezed (r = 0) case for T < T c since actually the mode ǫ is also contributing to entanglement generation. On the other hand the entanglement phases shows the same structure for high temperatures (right panel) where the only difference resides in the attenuated growth for entanglement when r increases (see color bars).
Notice that all the expressions have been calculated in the limit of small gamma assuming a final Gibbs state for decohered eigenmodes and free evolution for nondecohered ones. Of course this situation can only be perturbative, since for stronger coupling to the bath the eigenmodes become increasingly coupled among them, through second order processes mediated by the bath. This necessarily leads to decoherence of all eigenmodes, at a low rate though.

B. Quantum synchronization
In this section we analyze the dynamics of the system showing the existence of a parameter manifold where the oscillators oscillate in phase, synchronously, in spite of having different natural frequencies. Full dynamics for Gaussian states is characterized by mean values of position, momenta and variances (second order moments). The possibility to have synchronization in this system is important for two reasons: (i) this phenomenon has been largely studied in non-linear systems but we show that for dissipation in a common bath, it can arise even among harmonic oscillators; (ii) few attempts have been done to extend it to the quantum regime and we show here that this phenomenon witnesses the presence of robust quantum correlations, being actually a synchronous steady state accompanied by asymptotic entanglement.
Spontaneous synchronization -where no external forcing is considered-has been recently studied in some quantum systems like optomechanical cells or nanomechanical resonators [33,34,51,52]. Actually the phenomenon of mutual synchronization has been extended to the quantum regime considering a different scenario in Ref. [33] where two coupled harmonic oscillators with different frequencies and initially out of equilibrium are studied during their relaxation towards a thermal state. Synchronization was reported in first and second order moments during a long transient and accompanied by the robust preservation of quantum correlations (quantum discord [53,54]) between oscillators. Two oscillators dissipating in a common bath are actually preserving asymptotically their entanglement and retaining a larger energy than in the thermal state only if they are identical [29]. In this (symmetric) case they also evolve towards a synchronous asymptotic state [33]. When three elements are considered, we have shown above that the symmetric chain can reach an asymptotic regime with entanglement between the external oscillators, independently on the frequency of the central one.Then asymptotic synchronization between the external pair is also expected. Beyond this symmetric case, more interesting is the possibility offered by a chain to freeze all the oscillators out of the thermal state when their frequencies are all different, as discussed below.
The long time dynamics of our system can be straightforwardly calculated by assuming that normal modes which are coupled to the bath get thermalized, while uncoupled ones have a free evolution. This is sufficient to analyze the presence of synchronization in the asymptotic state. In general quantum mutual synchronization appears always in one-mode NSs among natural oscillators linked by the non-dissipative mode, as long as they have an asymptotic dynamics with only one oscillatory contribution. Phase or antiphase synchronization at the non-dissipating normal frequency is possible in first order moments depending on the sign of their F matrix coefficients, while only in-phase synchronization occurs for second order moments at twice the frequency of first order ones. Let us illustrate it in some situations and compare with the time evolution of q 2 i ∀i (Fig. 4) when considering a simple Markovian dynamics in the weak coupling limit (see Appendix B).
Consider firstly the specific case of an open chain with equal couplings and frequencies in the external oscillators (corresponding to situation in Fig. 1a). The form of the non-dissipative normal mode is C δ = (1, 0, −1) T / √ 2 so synchronization will emerge only between external oscillators in antiphase for the position and momenta at frequency Ω δ = ω (the normal mode frequency) and for the second order moments (necessarily in-phase and at  Fig. 1a where one mode NS is generated (ω = 1.3ω2 , λ = 0.4ω 2 2 ) synchronizing the external oscillators at 2ω, (b) configuration in Fig. 1c where a different mode NS is generated (ω1 = 1.2ω2 , ω3 = 1.3ω2 , λ = 0.4ω 2 2 ) producing synchronization in all pairs of oscillators at 2Ωǫ and (c) configuration in Fig. 1e when a two modes NS is generated (ω = 1.3ω2 , λ =λ0) and synchronization is lost. Bath parameters for the simulation are in all cases T = 10ω2, γ = 0.07ω2 and Λ = 50ω2. 2ω). The central oscillator instead decays into the thermal state, its initial oscillations being suppressed in the long time dynamics. This case is shown in Fig. 4a, where synchronization appears after a transient only for the external oscillators of the open chain while the central oscillator looses oscillation amplitude.
In the latter case synchronization appears between identical unlinked (λ 13 ) oscillators in a symmetric chain (Fig. 1a). More peculiar is the case in which all oscillators have different frequencies and eventually couplings. In the case of Fig. 1c we actually have that the non-dissipative mode involves all the three oscilla- This can actually give rise to synchronous dynamics of all the oscillators, in spite of the difference in their natural frequencies. Since one of the components has different sign than the other two in C ǫ two of the oscillators first moments will synchronize inphase between them and in antiphase with the third one. In Fig. 4b a total synchronization is produced involving all three (different) oscillators, consistently with the fact that the non dissipative normal mode ǫ involves all three oscillators.
A different situation is produced when we have a twomodes NSs since two oscillatory contributions are present in the asymptotic limit of the natural oscillators. Here synchronization is only possible when the two normal modes frequencies are the same. An example is the open chain with a two-modes NS (Fig. 1e) where apart from the previous non-dissipative mode C δ , actually the mode C ǫ = (1, −2, 1) T / √ 6 with frequency Ω ǫ = 2ω 2 2 − ω 2 does not dissipate either. In this case synchronization is destroyed by the presence of mode ǫ and it can be only recovered when Ω ǫ equals Ω δ , i.e. in the trivial case of independent (λ = 0) identical oscillators (ω 2 = ω). Lack of synchronization as well as a multimode oscillation are shown in Fig. 4c.
The initial state employed for simulations is a squeezed separable vacuum state, where the squeezing parameters have been chosen to be different (r 1 = 2, r 2 = 2.5 and r 3 = 3). In general, we have tried to avoid special initial conditions that could have filtered just one normal mode into the dynamics. What we discussed is therefore the emergence of synchronization as a dynamical process when considering more general initial states, leading to robust conclusions.
The scenarios here discussed allow to establish the effect of having one or two modes NS in the configurations of open chains (Fig. 1a and e). The same analysis can be extended to other cases where a different normal mode is uncoupled from the environment. For instance, the configurations in Fig. 1c and d admit only one nondissipative normal mode that involves the three oscillators, producing then a collective synchronization of the chain.

IV. THERMALIZATION AND ROBUSTNESS OF QUANTUM CORRELATIONS
Creation of NSs is a powerful tool to avoid decoherence and produce synchronized dynamics as we have seen in previous sections. However, the conditions leading to NSs are satisfied only in some parameter manifolds and it is relevant to see which is the effect of deviations from these couplings and detunings, that could also arise from the difficulty of experimental tuning. In this case, dissipation is present in all normal modes and the effective couplings of Eq. (7) are all different from zero. Thus a thermal state is finally reached in all the degrees of freedom.
In absence of NS entanglement is lost after a finite time. Although the asymptotic state is simply the thermal (Gibbs) state, damping dynamics of the normal modes with different decoherence and relaxation time scales is present, producing a rich behaviour in which synchronization or high quantum correlations may emerge during a large transient before the final thermalization of the system. These effects have been reported recently in the case of two harmonic oscillators, where disparate decay rates between the two normal modes is produced for small deviations from the resonat case [33].
A dynamical description of the system in terms of a Markovian Master Equation (Appendix B) reveals the central influence of the effective couplings in the relaxation time scales of the different normal modes. In this context, the ratio between the two smallest decay rates, defined as: provides important information about the dynamics of the system, being in fact one of the central figures in order to predict the robustness of correlations between oscillators or the emergence of synchronization as we will see in the following. In presence of disparate decay rates (R << 1), a large time interval appears between thermalization of the two modes with largest damping coefficients (strongly-damped modes) and thermalization of the mode with the smallest one (weakly-damped mode) . This produces the emergence, after a transient, of a long interval in which the weakly-damped mode is effectively the only one present in the dynamics, hence producing the synchronization between pairs of oscillators linked by this normal mode and the slow decay of quantum discord between these pairs. On the other hand, when the decay rates are similar (R ∼ 1) the different modes are present for all times inhibiting synchronization, and the survival of correlations associated to one of the modes for long times is lost. These phenomena will be exemplified in the scenario of an open chain with equal couplings (λ 12 = λ 23 ≡ λ).
In Figure 5 we represent R showing broad regions in which a weakly-damped mode exists (white regions) near the NSs manifolds corresponding to the configurations in

A. Quantum correlations
Even if out of the NS conditions entanglement suffers a sudden death, following Ref. [33] other quantumness indicators, such as quantum discord can remain robust in regions where disparate decay rates (R << 1) exist. By using an adapted measure of discord for Gaussian bipartite states [55,56] we observe the existence of a plateau in the dynamical evolution of discord between pairs of oscillators which are linked by a weakly-damped normal mode. More precisely, in the white region of Figure 5, close to the dashed-dotted hyperbola, the weakly-damped mode links the three natural oscillators, producing a plateau in the evolution of discord for all pairs. Moving to the tighter white region, close to the dashed diagonal line, the weakly-damped mode only involves the external oscillators pair of the open chain, leading to a slowly decaying discord only for this pair of oscillators. In blue regions no plateau is observed for discord, reaching for each pair the value corresponding to the thermal state in shorter times. Figure 6 shows time evolution of discord in logarithmic scale for the three pairs of oscillators (see colors) for a selection of parameters close and far away from the dashed-dotted hyperbola (Fig. 6a and 6b respectively). The initial condition has been taken to be a squeezed separable vacuum state with same squeezing parameters as in Fig. 4 and will be kept for further simulations. A Gaussian filter has been employed to eliminate fast oscillations (original quantities are plotted in gray), in order to make it easier to identify the plateau characterizing discord robustness.
As already seen for asymptotic entanglement, the effect of increasing the bath's temperature is, in general, a degradation of quantum effects. It is therefore important to see how robustness of discord is a feature present also in hotter environment. The main effect when increasing T is that the final thermal state (by virtue of Eq. B4) displays lower correlations so that the amount of discord that can be maintained in a robust way diminishes. In Figure 7 we show the evolution of discord for a pair of linked oscillators (1, 2) for increasing T by factors 3 and 6 (other parameters as in Fig. 6a): while the plateau is present for all temperatures and their (negative) slope is very similar, a lower amount of discord can be generated in the initial transient producing a shift of the curves to lower values (oscillations are increased by the logarithmic scale of the plot). This degradation by temperature effects can be avoided by increasing the squeezing in the initial separable vacuum state like in the case of entanglement presented in section III A.

B. Synchronous thermalization
With respect to the emergence of synchronization for pairs of oscillators when the NS is lost, we have to point out that when thermalizing the system reaches a station- ary state where oscillations are suppressed. We therefore restrict our analysis to a transient (which becomes longer the more we approach one of the NS conditions) where oscillations in the first and second order moments are still present in the dynamics. In this situation synchronization of first and second order moments can be estimated quantitatively by introducing a typical indicator, that is defined for two generic time-dependent functions h(t) and g(t) as: . When evolutions are phase or antiphase synchronized we will obtain |C| ∼ 1, while for very different dynamics we will obtain a value of C near to zero. Figure 8 shows the synchronization indicator C q 2 i , q 2 j with position variances of (a) the external pair of oscillators i, j = 1, 3 and (b) for i, j = 1, 2 of the open chain with identical couplings and varying the external oscillators frequencies (in the same range as in Fig. 5 ). We see immediately the high resemblance with the R map of Fig. 5 and some interesting differences induced by the shape of the normal modes. Effectively the external pair of oscillators (1, 3) synchronizes (C ∼ 1) in all regions where disparate decay rates (R ≪ 1) are predicted since these two oscillators are linked by the weakly-damped normal mode in these regions (Fig. 8a). As for the internal pair (1, 2), it does depend on the weakly-damped mode in the vicinities of the dashed-dotted hyperbola where synchronization is actually found. On the other hand, near the diagonal the weakly-damped mode approximates to C δ = 1/ √ 2(1, 0, −1) excluding the central oscillator from the induced collective motion. Consistently, the 1, 2 pair (Fig. 8b) does not synchronize for ω 1 ∼ ω 2 (near the diagonal).
We finally point out that, as expected, the synchroniza- tion frequency is that of the weakly-damped mode (Ω 0 ) for the first order moments (position and momenta) and 2Ω 0 for the second order momenta.

V. CONCLUSIONS
Decoherence in an open quantum system can be avoided or reduced by tuning the system parameters in a common environment context. The shape of the interaction Hamiltonian between system and bath can be used in order to engineer the protection of some degrees of freedom from the environmental action. This analysis has been solved in the present work for a system of three coupled harmonic oscillators in contact with a bosonic bath in thermal equilibrium, developing the necessary general relations so as to obtain a NS composed by one or two non-dissipative normal modes.
Different open and close chain configurations have been explored, highlighting the richer variety of NS configurations available when the dissipative system is extended from two to three harmonic oscillators.For a symmetric open chain (equal frequencies of the external pair of oscillators and same coupling) a closed analytical expression for the asymptotic entanglement (logarithmic negativity) has been derived observing the appearance of three different phases depending on temperature and squeezing of the initial state (sudden death of entanglement, a infinite series of sudden death events and revivals and asymptotic robust entanglement). Sudden death of entanglement can be avoided for arbitrarily high bath temperatures by increasing the squeezing in the initial state for both one or two modes NS. Remarkably this critical squeezing in order to avoid sudden death depends logarithmically on temperature.
Other dynamical effects such as the emergence of synchronization of mean values and variances have been ana-lyzed in different situations by simply assuming the relaxation to a thermal state of the normal modes coupled to the environment. Coherent oscillations appear when only a surviving normal mode is present in the dynamics, inducing synchronization in the natural oscillators that depend on it. Interestingly the parameter manifold leading to NSs include several not symmetric configurations: for instance an hyperbolic relation among frequencies can be satisfied for identical couplings in an open chain; in this case both asymptotic entanglement and synchronization are predicted even if all the oscillators natural frequencies are different (a possibility offered by a chain of three oscillators and absent in the case of two). A more general scenario where oscillators networks are considered with dissipation acting globally and locally has been discussed in Ref. [34].
Furthermore an analysis of situations in which the NS conditions are not accomplished at all has been performed. Indeed, important properties can be present although when deviating from NS conditions entanglement does not survive: robust conservation of discord during a long transient dynamics and the emergence of synchronous oscillations are found before thermalization. These effects are interpreted in relation to disparate decay rates for the normal modes. As long as there is a weakly-damped mode surviving among several stronglydamped modes, effects such as robust discord and synchronization arise among the oscillators following this normal mode.
While we have focused along this paper in specific cases of all the possible three oscillator configurations (in which calculations are greatly simplified) the strategy provided here is rather general and applies straightforwardly to other choices of system parameters that produce the decoupling of one or more normal modes. In fact this method can extend this engineering of the normal modes to more complicated systems such as disordered harmonic lattices or complex networks. All the information about bipartite quantum correlations for a Gaussian continuous-variable state is condensed in its covariance matrix defined through the ten second order moments of q (A,B) and p (A,B) (in our case first order moments are initially zero). This bipartite covariance matrix defined for a system of two oscillators A and B, can be written as:

VI. ACKNOWLEDGEMENTS
where α, β, γ are (2 × 2) blocks: α(β) contains the second order moments of oscillator subsystem A (B), and γ contains correlations of both subsystems. The minimum symplectic eigenvalue (of the covariance matrix corresponding to the partially transposed density matrix), necessary to calculate the logarithmic negativity, is given by: with a = 4 det(α), b = 4 det(β) , g = 4 det(γ) and s = 16 det V . Normal modes coupled to the environment will reach in the asymptotic limit a thermal state, given by the Gibbs distribution e − H S T /Tr(e − H S T ). For a normal mode (k) we get: being Ω k the corresponding normal mode frequency and T the reservoir temperature. On the other hand the uncoupled modes evolve freely. This means that the asymptotic covariance matrix can be calculated by expressing all second order moments of natural oscillators in terms of the normal modes and then substituting the asymptotic expressions corresponding to free modes or thermalized ones. The covariance matrix in the asymptotic limit can be separated into three parts corresponding to the contributions of each normal mode. In terms of the blocks we can write: For a non-dissipative normal mode (n) we have: and for a dissipative one (k) we get: While elements in V diss are given by the expressions (A3) those of V no−diss are the ones corresponding to a free evolution of an harmonic oscillator. This analysis gives all the necessary elements in order to calculate asymptotic entanglement for pairs of oscillators in every particular situation in which one or two of the normal modes are un-coupled from the environmental action. Consider the specific case of an open chain (λ 13 = 0) in which we have two equal frequencies (ω 1 = ω 3 ≡ ω) and two equal couplings (λ 12 = λ 23 ≡ λ = 0) (Fig. 1a). In this case we get only one normal mode decoupled from the bath. In order to calculate the expression of the minimum symplectic eigenvalue we have to first calculate the elements of the three normal modes, that are shown here as vector columns: where we have labeled the non-dissipative mode as δ and the other two modes as {±}. The corresponding normal modes frequencies are Ω δ = ω and Ω ± = (ω 2 2 + ω 2 )/2 ± √ ǫ, defining ǫ = ( ) 2 + 2λ 2 , and c ± is nothing but a normalization constant. We can now obtain all the terms appearing in V ± .
The initial condition given in Eq. (24) can be now rewritten in terms of the non-dissipative normal mode as Q 2 δ (0) = 1 2ω e −2r , P 2 δ (0) = ω 2 e 2r and {Q δ , P δ }(0) = 0, and then their free evolution is given by: (e 2r sin 2 (ωt) + e −2r cos 2 (ωt)) where we have already used that Ω δ = ω. By substituting in V no−diss (Eq. A4) we can now obtain the expressions of the determinants and introducing into equation (A2) yields the following expression for the minimum symplectic eigenvalue: which is an oscillatory function with frequency 2ω defined by the following functions: and the dependence on bath temperature and shape of dissipative normal modes is given by: From Eq. (A7), we can obtain the minimum entanglement (given when t = (2n + 1) π 2ω ; n = 1, 2, 3, ...) and the maximum one (given when t = (n + 1) π ω ; n = 1, 2, 3, ..) in order to recover Eq. (25) in section (III A) of the main text with the proper definitions specified there.
On the other hand, if we move to situation represented in Figure 1e by fixing λ =λ 0 (see Eq. 22) we have that the normal modes transform into: being the last one, C CM , the only dissipative mode. The corresponding normal frequencies are respectively: Ω δ = ω, Ω ǫ = 2ω 2 2 − ω 2 and Ω CM = 2ω 2 − ω 2 2 . Naturally we have to restrict ourselves to the regime 2ω 2 3 > ω > ω 2 3 /2 in order for these quantities to be real and positive.
Keeping the same initial condition as in the previous case, we have that nothing changes in the expression of the free evolution of mode δ of Eq. (A6) while the free evolution of mode ǫ is given by: We have assumed the same squeezing parameter r in the central oscillator of the chain (note that in the previous case the initial state of the central oscillator is not relevant and we have not specified it). Following the same procedure as above, we calculate an expression for the minimum symplectic eigenvalue. It is worth noting that in this case we have two contributions to the determinants of the free type V no−diss (Eq. A4) corresponding to the two non dissipative modes, and only a dissipative one V diss (Eq. A5) corresponding to the center of mass mode. The minimum symplectic eigenvalue yields: where we have defined the following quantities in order to simplify the expression. The constant terms: A 0 = cosh(2r) 4(σ Q + σ P ) + J + (Ω 2 ǫ + ω 2 ) B 0 = 64σ P σ Q + 4(ω + ω 2 ) 2 81ωω 2 + 32Ω ǫ ωJ + 3 ωσ P Ω ǫ + Ω ǫ σ Q ω and the oscillating terms: A 1 (t) = 4 cos(2ωt) sinh(2r)(σ Q − σ P ) + + J + cos(2ωt) sinh(2r)(Ω 2 ǫ + ω 2 ) + + J − cos(2Ω ǫ t) cosh(2r)(Ω 2 ǫ − ω 2 ) − − J − cos(2(Ω ǫ − ω)t) sinh(2r) (Ω ǫ + ω) 2 2 − − J − cos(2(Ω ǫ + ω)t) sinh(2r) (Ω ǫ − ω) 2 2 B 1 (t) = cos(2Ω ǫ t) 32Ω ǫ ωJ + 3 where different frequencies coming from the two nondissipative modes are present. We have used J ± = 1 12ω e 2r 2ω2+ω Ω 2 ǫ ± e −2r 2ω+ω2 ωω2 and the two new bath dependent functions are now simply given by the contribution of the center of mass mode's thermal state:

Appendix B: Master Equation
Assuming the general framework provided in section II, we can postulate a simple approach to the dynamical evolution in the weak-coupling limit between system and environment. By using the general Born and Markov approximations as well as initial product state, a Markovian Master Equation (MME) for the reduced density matrix of the open system can be obtained [44] in the normal modes basis. The resulting equation is not of the Lindbland form, thus complete positivity (CP) is not guaranteed [57]. This issue can be solved in two different ways, either considering a rotating-wave-approximation (RWA) in the interaction Hamiltonian (4) x i x j → a i a † j + h.c. or performing a strong-type RWA in the non-Lindbladian master equation by eliminating oscillatory terms of the form e ±i(Ωi±Ωj )t that appear in the interaction picture. The latter is the one we will pursue. The advantages of this method not only reside in obtaining a Master Equation in Lindbland form (thus CP), but also in that dynamical evolution can be solved analytically. However an exhaustive analysis in the case of two harmonic oscillators shows a very well agreement between results using the original non-Lindbladian master equation and the strong RWA here used [33].
The MME for the evolution of the reduced density matrix for a common bath in the strong RWA is then: Here Γ n and D n are constant coefficients (by virtue of the Markov approximation) accounting for the damping and diffusion effects respectively. Note that under this approximation, each normal mode is dissipating separately to the bath, i.e. they have independent decay rates. The bath has been considered to be in a thermal (Gibbs) state at temperature T , and to be composed by a continuum of frequencies characterized by an spectral density J(Ω). For simplicity it has been considered to be Ohmic with a sharp cutoff J(Ω) = 2γ π Ω Θ(Λ − Ω), where Θ(x) is the Heaviside step function, Λ is the largest frequency present in the environment (cutoff frequency) and γ is a constant quantifying the strength of system-environment interaction (thus in the weak-coupling we have always that γ ≪ Ω i ∀i ∈ {1, ..., N }). This assumptions leads to the following definitions of the Master Equation coefficients: Γ n = κ 2 n π 2 J(Ω n ) Ω n = γ κ 2 n (B2) D n = κ 2 n π 2 J(Ω n ) coth Ω n 2T = γ κ 2 n Ω n coth Ω n 2T where we also assume Ω i < Λ ∀i ∈ {1, ..., N }.
As we are interested in classical and quantum correlations of the system oscillators, a description for the evolution of the first and second order moments is necessary. This can be written in a simple form as: where R is a column vector containing the M = (2N + 1)N independent second order moments of the normal modes. The matrix M condenses all the information about their dynamical evolution and N determines the stationary values for long times (whenṘ = 0). The dynamics of R can be solved in terms of the eigenvalues of M: where the i = j eigenvalues determine the evolution of Q 2 i , P 2 i and {Q i , P i } , while the ones with i = j determine that of Q i Q j , P i P j and {Q i , P j } . Note that by virtue of eqs.(B2) and (B3) the decay of the normal modes is entirely governed by the effective couplings mentioned above, thus differences in their magnitude produce disparate temporal scales for the dissipation and diffusion of normal modes.
The stationary state is found to be (R ∞ = M −1 N): being all the other second order moments equal to zero. Note that these expressions for the asymptotic limit recover the thermal state of the system at the bath temperature T given by the Gibbs distribution in Eqs.(A3).