Entanglement Dynamics of Coupled Quantum Oscillators in Independent NonMarkovian Baths

This work strives to better understand how the entanglement in an open quantum system, here represented by two coupled Brownian oscillators, is affected by a nonMarkovian environment (with memories), here represented by two independent baths each oscillator separately interacts with. We consider two settings, a ‘symmetric’ configuration wherein the parameters of both oscillators and their baths are identical, and an ‘asymmetric’ configuration wherein they are different, in particular, a ‘hybrid’ configuration, where one of the two coupled oscillators interacts with a nonMarkovian bath and the other with a Markovian bath. Upon finding the solutions to the Langevin equations governing the system dynamics and the evolution of the covariance matrix elements entering into its entanglement dynamics, we ask two groups of questions: (Q1) Which time regime does the bath’s nonMarkovianity benefit the system’s entanglement most? The answers we get from detailed numerical studies suggest that (A1) For an initially entangled pair of oscillators, we see that in the intermediate time range, the duration of entanglement is proportional to the memory time, and it lasts a fraction of the relaxation time, but at late times when the dynamics reaches a steady state, the value of the symplectic eigenvalue of the partially transposed covariance matrix barely benefit from the bath nonMarkovianity. For the second group of questions: (Q2) Can the memory of one nonMarkovian bath be passed on to another Markovian bath? And if so, does this memory transfer help to sustain the system’s entanglement dynamics? Our results from numerical studies of the asymmetric hybrid configuration indicate that (A2) A system with a short memory time can acquire improvement when it is coupled to another system with a long memory time, but, at a cost of the latter. The sustainability of the bipartite entanglement is determined by the party which breaks off entanglement most easily.


Introduction
In this work we continue the inquiries two of us conducted in 2015 concerning the quantum thermodynamics and entanglement dynamics of two coupled oscillators, each with its private bath, for possible applications to quantum processes in systems involving two heat baths such as quantum transport and the design of quantum heat engines and other quantum devices. In [1] we show that for bilinear couplings between the oscillators, and between each oscillator and its own bath, described by a scalar field at temperatures T 1 > T 2 , the dynamics of the system can be solved exactly at arbitrary temperatures, even for strong coupling, thanks to the Gaussian structure. In particular, a nonequilibrium steady state (NESS) of this system can indeed exist at late times, its insensitivity to the initial conditions of the system testifies to the uniqueness of this NESS. In [2] we studied how entanglement between the two coupled oscillators depend on the temperatures of the two baths. We find that the valid entanglement criterion in general is not a function of the bath temperature difference, in contrast to thermal transport in the same setting. In fact lowering the temperature of one of the thermal baths does not necessarily help to safeguard the entanglement between the oscillators. Rather, quantum entanglement will disappear if any one of the thermal baths has a temperature higher than a critical temperature T c . The baths in both of these studies are made up of scalar fields which are Ohmic (linear in frequency), and, with a full spectrum, there is no imposed cutoff in its frequency range. With no particular time scale involved they are memoryless. Scalar field baths are Markovian.
In this work we focus on the effects of nonMarkovian (nM) baths on the nonMarkovian dynamics of an open system in the same setup, with constant coupling between the two quantum oscillators, but not limiting our attention to very early or very late times. Rather, we aim at ranges where nonMarkovian dynamics are prominent and examine how the entanglement of the system evolves. In particular, we are interested in how the memory effects associated with the nonMarkovian dynamics affects the system's entanglement throughout the evolution. Our next paper [3] will consider time-dependent coupling between the two oscillators and focus on the conditions which could allow entanglement to survive even at high temperatures, the so-called 'hot entanglement' first discovered by Galve et al. [4]. Since nonMarkovianity is the central theme of both papers we begin with a brief description of nonMarkovianity in open quantum systems.

NonMarkovianity in Open Quantum Systems
'Closed' systems in Nature do not exist, except perhaps the universe as a whole. It is a conceptual idealization. Almost all physical systems possess some environment which they interact with, no matter how weakly. Likewise, treating dynamical processes of a system as Markovian (memoryless), though commonplace, is an approximation which needs justification. In nature, rarely can any system function properly in a sustained way without memory.
Open system approach to the description of Nature has a long history, three centuries at least, in terms of thermodynamics and statistical mechanics, where a bath characterized by a few parameters serves as the environment of a system one is interested in, and clever notions like ensembles are introduced for treating systems in equilibrium conditions or linear response theories for near-equilibrium cases.
Quantum systems add to these the factors of spin-statistics and quantum phase effects like coherence and entanglement. The latter's significance was first put to use in quantum optics in the 50 s/60 s and in quantum information sciences in the 90 s/00 s, both bringing forth application consequences of revolutionary magnitudes.
Placed in this light, one can see the importance of nonMarkovian quantum processes and appreciate both their theoretical and practical values, with regard to the foundations of open quantum systems and applications to the many branches of quantum information sciences and engineering.

System's NonMarkovian Dynamics Linked to NonMarkovian Baths
Because of the manifold challenges unmatched by, and not encountered in, the more familiar Markovian processes, serious systematic investigations of nonMarkovian processes in open quantum systems (OQS) [5][6][7] after three decades are still in its developmental stage. Notable landmarks were the derivation of a nonMarkovian master equation for quantum Brownian motion (QBM) [8] (and its associated Fokker-Planck-Wigner equation [9] and Langevin equation [10]), mathematical formulations and proposals for the criteria and measures of nonMarkovianity [11,12], and studies of nonMarkovian behavior in various schemes of quantum information processing. This subject is now also enriched by several comprehensive reviews [13][14][15].
Solving the underlying nonMarkovian dynamics enables one to investigate the timeevolution of entanglement in the open system. For specific prior work with similar setups as here with nonMarkovianity emphasis we mention the following: Maniscalco et al. [16,17] study the non-Markovian dynamics of a two-mode bosonic system interacting with two uncorrelated thermal bosonic reservoirs and, from it, the dynamics of entanglement for bipartite Gaussian states. They analyze the effects of short-time system-reservoir correlations on the separability thresholds and show that the relevant parameter is the reservoir spectral density. If the frequencies of the involved modes are within the reservoir spectral density, entanglement persists for a longer time than in a Markovian channel. Liu and Goan [18] investigate the entanglement evolution of two interacting harmonic oscillators under the influence of non-Markovian thermal environments for both cases, namely, each oscillator has its own independent thermal reservoir, or sharing a common reservoir. An and Zhang [19] use the rotating-wave form for the inter-oscillator coupling to investigate the entanglement dynamics of continuous-variable quantum channels in terms of an entangled squeezed state of two cavity fields in a general non-Markovian environment at zero temperature. The influence of environments with different spectral densities, e.g., Ohmic, sub-Ohmic, and super-Ohmic, is numerically studied. Wilson et al. [20] study the nonMarkovian effects on the dynamics of entanglement in an exactly solvable model that involves two independent oscillators, each coupled to its own stochastic noise source. They found that all memory effects enter via the functional form of the energy and hence the time of death and rebirth is controlled by the amount of noisy energy added into each oscillator.
We also mention in passing the other commonly studied set-up, which is to place the two coupled oscillators in a shared bath. It is more involved due to the addition of environment-induced nonMarkovian interaction between the two oscillators. For entanglement dynamics under constant inter-oscillator coupling, see e.g., [21,22], where a comparison table with a handful oft-cited works can be found. We plan to address hot entanglement under parametric coupling in a common bath setting in a future work. Instead of appealing to formal narratives or abstract proofs we would like to seek explicit solutions of the reduced system dynamics to get some direct feels of the physics. To do this we need to cover time ranges beyond the initial transient stage where most of the papers addressing the nonMarkovian criteria and measures are focused upon. Behavior of the nonMarkovian open quantum system in the mid-time-range is needed to assess how long and in what degree an initially entangled state can sustain its entanglement against environmental degradation. The initial transient stage often contains artefacts due to somewhat contrived initial conditions such as the jolt, resulting from the assumption of a product initial state. The late time behavior is often easier to obtain from asymptotic analytical solutions, but then in most cases without any external drive, the initial entanglement is likely to have died off.
In this and a follow-up paper [3] we assume the same two coupled quantum oscillators and two thermal baths set-up as described in the two 2015 papers mentioned earlier but place our focus on how a nonMarkovian (nM) environment may alter the entanglement dynamics of the coupled oscillator system. We will consider the simplest nonMarkovianity configuration, where the memory time of the bath is realized by the cutoff scale in the spectral density. It quantifies the time scale over which the system's evolution depends on its past history. The state of the system at earlier times is remembered by the bath, and that information is continuously fed back to the system at later moments which affects the system's subsequent evolution.

NonMarkovian Baths with Different Spectral Density Functions
To identify what imparts the nonMarkovian (nM) properties of a bath note that two factors enter in the spectral density function J(ω, Λ). We can go by what is described in [8]: "Different environments are classified according to the functional form of the spectral density I(ω). On physical grounds, one expects the spectral density to go to zero for very high frequencies. Let us introduce a certain cutoff frequency Λ (a property of the environment) such that I(ω) → 0 for ω Λ. The environment is classified as Ohmic if in the physical range of frequencies (ω < Λ) the spectral density is such that I(ω) ∝ ω, as supraohmic if I(ω) ∝ ω n if n > 1, or as subohmic if n < 1".
Let us write the spectral density J(ω) of the bath in the frequency domain ω in the generic form (cf., e.g., Equation (I.3) in [8]) J(ω) = ω n P Λ (ω) involving two factors. First, the power of ω in it distinguishes between Ohmic and non-Ohmic cases, as quoted above. An environment modeled by a massless scalar field in four dimensional spacetime is an example of Ohmic bath. The broader class of non-Ohmic spectrum is subtler in the sense that it may induce instability in the oscillator's dynamics. An example is radiation reaction in moving charges in an electromagnetic field, with the appearance of higher derivative terms and runaway solutions. The advantage of taking a nonMarkovian approach, treating the electromagnetic field as a supra-Ohmic bath, in mitigating these 'pathologies' is discussed in a recent paper [23].

Memory in NonMarkovian Baths Linked to Finite Cut-off Frequency
The second factor, the cutoff scale Λ, is what we focus on in this paper. In this context, the cutoff scale serves dual purposes: straightforwardly it suppress the high frequency contribution of the bath modes, so we require that in the spectral density that the cutoff function P Λ (ω) should fall off to zero sufficiently fast, when ω is greater than the cutoff scale Λ. For simplicity, we assume that Λ is the only scale in P Λ (ω). Then another requirement that P Λ (ω) → 1 as Λ → ∞ will give the spectral density of a Markovian bath in that limit.
When the cutoff scale in the spectral density of the bath takes on a finite value, the Langevin equation of motion, which dictates the evolution of the reduced system, will contain a nonlocal integral expression. The nonloal expression determines how the system's evolution depends on its past history, and the cutoff scale measures the time scale this dependence lasts. The inverse of the cutoff scale gives the memory time and this is how it is linked to the nonMarkovianity of the bath.
Several types of cut-off functions have been analyzed in [23], namely, the Lorentzian and an exponentially decaying spectral density. The undesirable features of imposing a hard-cutoff as commonly practiced is also addressed. In this paper we suppose that P Λ (ω) has a double Lorentzian form, P Λ (ω) = Λ 4 /(Λ 2 + ω 2 ) 2 . In our next paper for parametrically-driven coupled systems we shall study nonMarkovian baths of Lorentzian form, previously also investigated by Estrada and Pachon [24].

Symmetric versus Asymmetric Set-Ups
In this work we assume that the two private baths attached to each oscillator is Ohmic and has a double-Lorentzian spectrum. We consider the qualitative behavior difference between two categories of parameter choices which we refer to as 'symmetric' versus 'asymmetric' settings or configurations. By a symmetric setting, we mean that both oscillators and the attached private baths have the same physical parameters, while in the asymmetric case, any of the parameters can be different. Here we focus on the asymmetric setting caused by different cutoff scales in the bath's spectral density. For the special case of asymmetric setting when one private bath is nonMarkovian and the other is Markovian we call the reduced system a hybrid system.

Issues of Interest
Here, in the two coupled oscillators with independent (uncorrelated and private) baths setting, we pose the following questions:

1.
How does the nonMarkovianity in the bath affect the system's entanglement? 2.
Which time regime does bath nonMarkovianity benefit the system's entanglement most?

3.
Can the memory of one nonMarkoviann bath be passed on to another Markovian bath? 4.
Does this memory transfer help to sustain the system's entanglement dynamics?
The answers to these questions regarding the baths' nonMarkovianity and the system's entanglement dynamics can be found in the last section Section 6. Before drawing conclusions on these issues we need to proceed in three stages: First, provide some background and explain the characteristics of the non-Markovian dynamics for the configuration under consideration. Then we commence our investigation with these questions in mind for the two different configurations, which we call symmetric and asymmetric. In performing numerical investigation we choose parameter ranges also with these issues in mind. Finally, after sufficiently useful results are obtained for the system's entanglement dynamics, we examine their behaviors, extract their meanings and attempt to answer the above list of questions.
The paper is organized as follows: We first summarize the formulation of the coupled nonMarkovian system in Section 2. Then in Section 3 we identify the characteristics of the nonMarkovian dynamics by means of numerical method. It is most easily done in the symmetric configuration setting, where we focus on the effects of the memory time scale on the evolution of the system from solutions to the equation of motion and the evolution of the covariance matrix. This clarifies the role of memory effects in nonMarkovian dynamics. We then look into the entanglement dynamics from the initial moment to the stage when the coupled system has relaxed to a steady state, and pin down the time regime when the entanglement in the system benefits from the bath nonMarkovianity most. During the course of investigation, we notice some similarity in the behavior between a nonMarkovian system strongly coupled to the thermal bath and a Markovian system weakly coupled to the bath. Something similar to this has been used by some authors as a way to 'simplify' the often hard-to-grasp nonMarkovian behaviors. In Section 4, we examine and compare these two cases, emphasizing their fundamental differences. In Section 5 we address an interesting issue in a coupled nonMarkovian system. By using an asymmetric configuration setting of a hybrid system we ask, can the memory effect be transferred between the coupled system and if so, what is the consequence to the system's entanglement? We conclude with a discussion of these issues in the last section, Section 6.

NonMarkovian Dynamics
Suppose we have a pair of coupled harmonic oscillators, each of which has its own private bath. Let χ i (t) be the canonical position of the ith oscillator, whose private bath φ i has an initial temperature β −1 i . We further assume that two oscillators have same physical frequency ω P . The coupling strength between the oscillators is denoted by σ. In addition, the oscillators have the same mass m, coupling strength e with each individual's bath, but each bath may not have the same cutoff scale in its spectral density.
The action describing such an interacting system is given by where φ i are two private baths with the spatial coordinate x i , and the ith oscillator is located at the origin z i . The private bath is initially prepared in a thermal states, but in the beginning, both baths are independent and uncorrelated. Then the equations of motion in the compact matrix form is [1,2,25] with Ξ(t) = {χ 1 (t), χ 2 (t)} T and Φ(t) = {φ 1 (t), φ 2 (t)} T , where T denotes the matrix transpose, and On the righthand side of (2), the operator Φ describes a noise force due to quantum and thermal fluctuations of the bath. It satisfies the Gaussian statistics The expectation value is taken with respect to the initial thermal state of the bath. The to formulate the fluctuation-dissipation relation of the ith bath [25], so they are also respectively called the noise kernel and the dissipation kernel. We use the convention to define the Fourier transformation of the function f (τ). We introduce the kernel function Γ (Φ) (τ) to isolate the contributions to renormalizations (When we carry out integration by parts on the integral expression in (2) with the identification (3), it yields a correction to Ω 2 B in (2). When the cutoff scale is infinite as in the field theory, we conventionally call this procedure renormalization (of the frequency).) of or the corrections to the parameters due to nonMarkovianity of the private baths. Integration by parts leads tö with γ = e 2 /(8πm), and Its formal solution can be formally found by the Laplace transformation of (8), and the Laplace transform of the solution is given bỹ Explicitly,D 2 is given bỹ where Thus due to the nonvanishing inter-oscillator coupling σ,D 2 will intermingle dynamics of both oscillators. This is better seen from the example Explicitly the displacement of oscillator 1 is also imparted by the thermal noise from the private bath of oscillator 2, indirectly via oscillator 2 through the inter-oscillator coupling. We may identify the normal modes by diagonalizingD 2 ; however the transformation matrix to the normal modes tends to be rather complicated in the asymmetric case, and in the time domain such a transformation is not local in time, so we will not pursue this approach, except for the symmetric setting.

Symmetric Setting
We first start with the symmetric setting. In this case, both oscillator are assumed to have the same physical frequencies ω 1,P = ω 2,P = ω P and damping constants γ 1 = γ 2 = γ. Likewise, both private baths will have the same power spectral densities and the cutoff scales . This setting is particularly convenient to investigate the effects of the bath's nonMarkovianity, manifested in this case by the cutoff scale or the memory time.
In the symmetric setting, the normal modes are nothing but the center-of-mass and the relative superpositions of the original two modes, Then we havë where ω 2 ± = ω 2 ± σ represents the oscillating frequencies of the normal modes, and φ ± is given by Here we remind that both private baths, associated with each oscillator, are initially uncorrelated and do not have direct coupling. Thus each mode acts as an Ohmic, nonMarkovian oscillator with different oscillating frequencies. Since for the equation of motion of the form (15) takes the generic form we will first use this equation to discuss the effects of the bath nonMarkovianity on the system's dynamics.
The general solution of χ(t) is given by the inverse Laplace transform of in whichΓ We will choose the double Lorentzian form for the bath spectral density. The two fundamental solutions d 1 (t) and d 2 (t) are two particularly useful homogeneous solutions to the equation of motion (17). They satisfy the initial conditions d 1 (0) = 1,ḋ 1 (0) = 0, d 2 (0) = 0 andḋ 2 (0) = 1, and are used to construct the complete solution to (17). For the double Lorentzian bath spectral density, the Laplace transformd 2 (z) has four poles: two negative real poles and two complex poles, whose real parts are real. Since for a given set of γ, ω P , and Λ, the two real poles are more negative than the real parts of the complex poles, so at late times, the damping behavior of d 2 (t) is controlled by the real part of the complex poles. Thus, we may construct an effective damping constant, which is approximately given by when γ/Λ 1, γ/ω P 1 in the long memory time limit. It implies that these two modes will have two different effective damping constants due to the difference in the normal mode frequencies, even though they have the same bath configuration. Figure 1 shows how the effective damping constant depends the oscillating frequency and the cutoff scale. We can immediately identify two features: (1) for a given cutoff scale, the damping is weaker, that is, smaller effective damping constant, when the frequency is larger. It is clearly seen in the second row of Figure 1, and the difference is more significant for smaller cutoff scales. This feature is related to the second feature: (2) for a given oscillating frequency, the effective damping constant γ EFF decreases if we lower the cutoff scale Λ, and such change is more dramatic when Λ is smaller than the oscillating frequency ω P . The weakening of the effective damping constant can be understood from two aspects. Since the inverse cutoff scale can be interpreted as the memory time of the bath, which moderates the duration the system, with which the bath interacts, depends on its past history. From this viewpoint, when the memory time is longer than the typical period 2π/ω P of the oscillatory motion, or when the cutoff scale Λ is smaller than the oscillating frequency ω P , the system dynamics from the previous cycle is still in good coherence with the current one. Thus the motion is progressively superposed to compete against decaying due to dissipation. This effect is expected to be more prominent when the memory is much longer, because the contributions from more earlier cycles will coherently join in [23]. Alternatively, the weakening effect can also be understood from resonance absorption [26]. The equation of motion (17) more or less describes a driven, damped harmonic oscillator, although this may not be obvious when the nonlocal integral expression is present. However, we may start from the limit Λ → ∞, where the integral expression gives 2γχ(t) and Equation (17) reduces to the standard form. It is known that such a system has peaked power spectrum, centered at around ω ∼ ω P with a width of the order O(γ). Now heuristically suppose we only allow bath modes whose frequencies are lower than the location of the resonance peak, ω P . These modes will not be very efficient in exchanging energy between the oscillator and the bath, so neither can they effectively drive the oscillator, nor can drain the energy out of the oscillator. Therefore it leads to relatively weak damping to the oscillator's motion. By these arguments, we expect that for a fixed Λ, in particular when Λ < ω P , a larger value of ω P means that the bath modes which participate in the interactions are further away from the resonance peak, thus rendering less effective energy exchange. Or, when the memory time becomes even longer than the system's period, thus allowing more coherent cycles from the past history of the oscillator's motion. These make the effective damping weaker. Figure 2 shows the cutoff scale dependence of the effective damping constant by numerical method. We see that when Λ 3 ω P , the effective damping constant γ EFF monotonically decreases with Λ, and is smaller than its Markovian counterpart γ = 0.3 ω P . On the other hand, in this case when Λ is greater than 20 ω P , the nonMarkovian effect is barely present. This raises an interesting observation. Here we choose γ = 0.3 ω P , which typically falls in the strong oscillator-bath coupling regime. By introducing nonMarkovianity via engineering the bath's spectrum, we seem to effectively render the reduced dynamics of the oscillator to behave like a weak coupling case. Or, by tuning the bath spectrum, we may be able to vary the effective coupling between the oscillator and the bath. The Λ dependence of the effective damping constant γ EFF of the oscillator, coupled to a thermal bath with the double Lorentzian spectral density. Here we choose ω P = 1 ω P , and γ = 0.3. In this case when Λ < 3.57 ω P , the effective damping constant is already weaker than γ.
Based on the above arguments and the idea of effective damping constant, we may be led to ask an important strategic question: Can we find a weak-coupling Markovian effective system to mimic the strong-coupling nonMarkovian system? It can greatly facilitate the theoretical study of the nonMarkovian system, because the latter is computationally challenging and demanding. In contrast, the former is analytically and exactly solvable. The answer depends. We will come back to this in Section 4.
To compute the entanglement measure (Relevant material can be found in Appendix A) via negativity, we start from the covariance matrix elements of the normal modes, since each mode is "decoupled" (Strictly speaking this is not entirely correct because φ ± are not independent). For the discussion of sustainability of entanglement, we suppose two oscillators are initially prepared in a two-mode squeezed vacuum state (see Appendix B for the discussion of the two-mode squeezed state). When the squeeze parameter η is not zero, the initial state is already entangled, and the degree of entanglement grows with increasing η. With respect to the normal modes, the nonzero covariance matrix elements of the coupled oscillator in this initial state are These two conditions ensure the covariance matrix with respect to the normal mode remains blockwise diagonal at all times Thus we only need to compute a smaller set of covariance matrix elements, Here d (±) i are the counterparts of the fundamental solutions, discussed in (18), of the normal modes.
Then we can restore the covariance matrix elements of the canonical variables of two modes by suitable superpositions of the covariance matrix elements of the normal modes. This step is essential because entanglement is partition dependent [27]. The entanglement measure, negativity, of a bipartite system depends on how we partition the bipartite system. Its values vary if we use the normal modes, instead of the canonical variables of the bipartite system, to construct the measure. This can be traced to the fact that the partial transposition does not belong to the symplectic transformation (see Appendix A for a brief discussion).
Before we proceed, we take a look at the time evolution of the selected covariance matrix element χ 2 + (t) for different choices of the bath temperatures and the bath cutoff scales. In Figure 3, different rows correspond to different initial bath temperatures β −1 , while different columns are associated with different bath memory times Λ −1 . We immediately see that χ 2 + (t) relaxes more slowly for longer memory times, consistent with the time evolution of the solution shown in Figure 1. We observe that each covariance matrix element can be decomposed into two components: (1) the intrinsic (active) component depends on the initial condition of the oscillator, and is damped by dissipation to the bath, and (2) the induced (passive) component is generated by the quantum and thermal fluctuations from the bath, independent of the oscillator's initial condition. The latter is the only component that will survive at late times, so the late-time value of the covariance matrix elements will not depend on the initial conditions of the system. Both components obey different statistics and are not correlated for the linear system. The memory endowed by the bath has different effects in both components. In the intrinsic component, the nonMarkovian effect leads to a weaker damping rate, so we expect that the information on the initial conditions of the system will linger for an extended period of time. On the other hand, in the passive component, the coherent superposition of the oscillator's motion over previous cycles, due to the memory effect, will compete with the accumulative intervention, resulting from the thermal fluctuations of the bath.
In each column, when the intrinsic component dies down, the value of the χ 2 + (t) dispersion will not decay to zero; otherwise it will violate uncertainty principle. In fact it will more of less come a constant due to balance between the thermal fluctuations and damping, and this constant value grows with the bath temperature because at higher temperatures, the bath fluctuations are stronger and becomes more random, less correlated, but damping in this case is independent of the bath temperature.
A longer memory time implies a longer effective relaxation time scale, and thus a longer duration for the effects of thermal fluctuations to accumulate. This probably explains the larger values of the dispersions χ 2 + (t) at late times when the cutoff scale Λ is smaller in each row of Figure 3. However, this difference is rather insignificant at higher bath temperature, signifying the ineffectiveness of nonMarkovian effects at higher temperatures.
We will transform the covariance matrix σ ± (t) back to that with respect to the modes χ 1,2 by where U 2 is a global symplectic matrix, that is, U 2 ∈ Sp(4, R). The symplectic eigenvalues of the partially transposed covariance matrix σ PT is then given by with ∆(σ PT ), det σ PT spelled out below in terms of covariance matrix elements of the normal modes as Thus, the entanglement of the system at any given moment t is determined by the value λ PT − (t) < 1/2. The time evolution of the symplectic eigenvalue λ PT − is shown in Figure 4. Since it contains products of various elements of the covariance matrix, we cannot simply decompose λ PT − into the intrinsic and the induced components as we did for each covariance matrix elements. However, we still expect that the initial entanglement between the oscillators will be gradually lost because contributions related to the oscillators' initial conditions will be exponentially small, and that the late-time behavior of the symplectic eigenvalue λ PT − is predominantly controlled by the bath. Thus the existence of the entanglement at late times will be independent of the initial conditions of the oscillators, but determined by the configuration of the bath such as the memory time, temperature and the oscillator-bath coupling strength.  Since we choose an entangled initial state, we see that the memory effect is much more significant during the relaxation stage. The duration τ ENT for both oscillators to remain entangled is correlated with the length of memory time. In Figure 4 the coupled oscillators become separable when λ PT − ≥ 1/2. Comparing the plots in each row, we observe that the entanglement interval during the relaxation stage increases with the memory time τ MEM . In other words, the entanglement is more robust against the thermal fluctuations and can be sustained at higher bath temperatures, owing to the nonMarkovian effect due to the bath's memory. This point is particularly clearly seen when we compare the Λ = 5 ω P , β = 10 ω −1 P and Λ = 1 ω P , β = 1 ω −1 P cases. In the latter, the bath temperature has 10 times higher than the former, but the entanglement dies out at approximately the same time around t = 2 ω −1 P , as shown in Figure 4. It is expected to be more dramatic for even smaller cutoff scale, as is inferred from the behavior of oscillator dynamics, shown in Figure 1. The result that τ ENT ∼ O(1) ω −1 P may not seem significant as it appears. However, recall that the oscillators are strongly coupled to the bath because γ = 0.5 ω P . In Figure 5, we will see that τ ENT for the corresponding Markovian bath (green curve) is merely τ ENT 0.197 ω −1 P . Furthermore, let us put it in a comparative perspective. For example, τ ENT 5.86 ω −1 P for the Λ = 1 ω P , β = 1 ω P case. The effective damping constants for the two normal modes are approximately given by γ  The time evolution of symplectic eigenvalue λ PT − (t) of the partially transposed covariance matrix σ PT (t). The system is entangled when λ PT − (t) < 1/2. The notation (Λ i , Λ j ) denote the cutoff scales of the private baths attached to the coupled oscillators. For example, the asymmetric setting (orange curve) discussed in this section corresponds to the curve with the tag (Λ 1 , Λ 2 ). Two other symmetric settings (Λ 1 , Λ 1 ) (blue curve) and (Λ 2 , Λ 2 ) (green curve) are used to contrast the effect of memory transfer on entanglement dynamics. The relevant parameters are m = 1, ω P = 1, γ = 0.5, σ = 0.2, β = 10, Λ 1 = 1 and Λ 2 = 20, and they are expressed in the unit of ω P .
In Figure 4, we also see that although the oscillators become disentangled after τ ENT , determined by ω P , γ, Λ and β, the entanglement can resurrect at a later time, when the memory time is sufficiently long and the bath temperature is low enough, as shown in the case Λ = 1, β = 10. However, this late-time entanglement does not benefit from the initial squeezing because its effect has been dissipated away. The plots show that although the nonMarkovian effect may improve entanglement at late times, the benefit is marginal. To understand this better, let us focus on the behavior of the covariance matrix elements and the symplectic eigenvalue η PT − after the oscillators are relaxed to the steady state. The late-time result in this case is particularly simple, since σ χ + p + (∞) and σ χ − p − (∞) vanish. We are left with only four elements σ χ ± χ ± (∞) and σ p ± p ± (∞), and they are given by Figure 6 shows that the late-time value of σ χ + χ + (t) decreases with the larger cutoff scale Λ or the shorter memory time. It is consistent with earlier observations, and it can be clearly seen that the variation of σ χ + χ + (t) is less significant, implying that the nonMarkovian effect is much weaker at higher bath temperature. This is also demonstrated in Figure 7, where the curves corresponding to different cutoff scales essentially converge when the bath temperature T B is of the order ω P . Note that the late-time values of the covariance matrix elements are mainly governed by the bath. When the bath has a higher initial temperature β −1 , it imparts stronger thermal fluctuations to the system. Thus the values of σ χ + χ + (∞) quickly increase with the bath temperature. In Figure 8, we show the dependence of the symplectic eigenvalue λ PT − on the bath temperature and the bath cutoff scale at late times, after the dynamics of the oscillators are fully relaxed. Recall that a large damping constant in the strong oscillator-bath regime is chosen, so the entanglement between the oscillators is typically difficult to maintain at late times. We see a trend, though not very significant, that the longer memory time τ MEM = Λ −1 or the lower bath temperature T B = 1/β is prone to keep the oscillators' entanglement at late times. In the left panel, with a longer memory time, the entanglement still exists at higher bath temperatures although this critical temperature does not change much with Λ and is still of the order O(ω −1 P ). On the other hand, the right panel shows that with a lower bath temperature, the entanglement may still survive for shorter memory times. These examples illustrate the coherent superposition induced by the nonMarkovian memory effect is capable to counter, at least marginally, the debilitating effects accumulated over the whole course of evolution due to the thermal fluctuations.

Strongly Coupled NonMarkovian Dynamics vs. Weakly Coupled Markovian Dynamics
In the previous section, we note that the dynamics of a coupled system strongly coupled to a nonMarkovian bath seems to behave in some aspects like that of a system weakly coupled to a Markovian bath, and raised the question whether we can use weakly coupled Markovian linear open systems to approximate strongly coupled nonMarkovian linear open systems? This issue is of significance because one may argue that even though many open system processes in the real world involve memories and are thus fundamentally nonMarkovian, we may in practice describe them by using a simpler effectively Markovian formulation. How valid are such prescriptions of convenience? We want to take a closer look at this issue with the help of a concrete example.
For the normal modes of the coupled oscillators described by (15), we suppose there correspond effective equations of motion of the normal modes for the weakly coupled Markovian oscillators For example, when ω P = 1, γ = 0.5 and Λ = 1, the effective damping constants of the normal modes are γ    Figure 9 shows the time evolution of d (+) 1 (t), the covariance matrix element σ χ + χ + (t), the uncertainty relation I + (t) for the normal mode +, and the symplectic eigenvalue η PT − . The blue curve corresponds to the strongly coupled nonMarkovian system (15), while the orange curve represents the effective, weakly coupled Markovian system (31). The fundamental solution in both cases look quite similar. They decay approximately with the same rate, but their phases vary with time. This minor difference starts showing its repercussion effects in the covariance matrix element σ χ + χ + (t), which accounts for the uncertainty of the χ + mode, or the degree of coherence. Both curves notably disagree before relaxation. The situation deteriorates even more for I + , the Robertson-Schrödinger relation of the mode +. Note that both modes are decoupled in the sense of (23). The plot shows that the nonMarkovian description of the same system seems to have a better control of coherence, compared with the effective Markovian description. Finally when we examine the symplectic eigenvalue λ PT − , they give quite a distinct prediction. The original nonMarkovian description shows that the entanglement between the oscillators can be maintained up to τ ENT 5.863 ω −1 P , but the effective Markovian formalism predicts roughly 0.6 τ ENT . Moreover, at late times, the original nonMarkovian description shows the presence of residual entanglement in the Λ = 1 ω P and β = 10 ω −1 P case, which is in the strong oscillator-bath coupling, low temperature regime. Since the disparity is more than marginal, the Markovian approximation, instead, predicts that the state of the system is separable. 1 (t), the covariance matrix element σ χ + χ + (t), the uncertainty relation I + (t) for the normal mode +, and the symplectic eigenvalue η PT − . The blue curve corresponding to the nonMarkovian system (15), while the orange curve is associated with the effective Markovian system (31). For d (+) 1 (t), the green curve describes the envelope of the damped oscillatory evolution. We choose m = 1, ω P = 1, γ = 0.5, σ = 0.2, Λ = 1, and β = 10.
The results in Figure 9, structured as a hierarchy from the evolutionary phase of the canonical variables, to its dispersion, then to the uncertainty of the constituent of the system, and finally into the integrated correlation among the system, indicate that the approximated Markovian description fails to precisely grasp the phase information embedded in nonMarkovian dynamics of the reduced system of the coupled oscillators. The approximated description tends to give poor predictions of the quantities that involving phases, or coherence. Thus the effective, weak coupling, Markovian description cannot be a sufficiently accurate substitute of the original strong coupling, nonMarkovian, linear open systems, even though the former has extraordinary convenience in computations.
Next we turn to the final issue discussed in this paper: whether the memory effect can be transferred from one party to the other in a bipartite system?

Asymmetric Setting
Suppose in the system of two coupled oscillators, one of which (oscillator 1) has a finite memory time due to small cutoff scale of its private bath, but the other (oscillator 2) has a negligible memory time due to the large cutoff scale. From the previous discussions, we learn that when they stand alone, oscillator 1 has a much smaller effective damping constant, resulting from the memory effect, in comparison with oscillator 2, so they will relax with different paces. Now when we couple them together, how does the coupled system evolve? Is the evolution dominated by the memoryless system, or by the nonMarkovian system? Or, will the memory effects in oscillator 1 be transferred to oscillator 2, so that it is shared among them, causing the coupled system to evolve in a cooperative way? Further, how does this affect the entanglement dynamics.
Since both private baths have different cutoff scales, the reduced system of coupled oscillators has an asymmetric configuration. We will start from (8), (10) and (11) with for i = 1, 2. Since working in the normal modes does not reduce computation hurdles, we will directly compute the covariance matrix elements of the canonical variables of the coupled system. For example, σ χ 1 χ 1 (t) will take the form H,0 (s − s ) H,0 (s − s ) .
As before, we suppose that the initial state of the coupled oscillator is a two-mode squeezed vacuum state. We plot the time evolution of the covariance matrix element σ χ i χ j (t) in Figure 10. In the right plot, the blue curve corresponds to σ χ 1 χ 1 (t) when the private baths of two coupled oscillators have the same cutoff scale Λ 1 = 1 ω P , while the green curve represents σ χ 1 χ 1 (t) when the private baths have the cutoff scale Λ 2 = 20 ω P . They can be respectively compared to σ χ 1 χ 1 (t) and σ χ 2 χ 2 (t) in the left plot. For the memoryless case in the right plot (green curve), the covariance matrix element σ χ 1 χ 1 (t) decays very fast, and its time evolution is almost fully relaxed when t ∼ 4 ω −1 P , which is about twice the relaxation time scale γ −1 . On the other hand, in the finite memory case, the blue curve falls off rather slowly. Its oscillatory behavior remains visible when t = 30 ω −1 P . The left plot of Figure 10 describes the above asymmetric setting where the oscillator 1's private bath has the cutoff scale Λ 1 , but the oscillator 2's bath has Λ 2 . We immediately see that for such a hybrid system, σ χ 2 χ 2 (t) (green curve in the left plot) does not decay as fast as the green curve in the right plot. Nonetheless, σ χ 1 χ 1 (t) of oscillator 1 in the hybrid system falls off faster than its counterpart (blue curve) in the right plot. These results imply transfer of the memory effect. Oscillator 2, which is essentially memoryless, benefits from such a transfer because its covariance matrix element shows a prolonged relaxation time (to roughly t ∼ 5 ω −1 P ), but this transfer shortens the memory time of oscillator 1. Thus in the hybrid system, motion of both oscillators in fact takes a cooperative way, and is neither dominated by the memoryless component nor governed by the one with a finite memory. Will their entanglement dynamics shows a similar feature? Figure 5 gives the time evolution of the symplectic eigenvalue λ PT − (t) of the partially transposed covariance matrix σ PT (t). The asymmetric setting is described by the orange curve, where two private baths have different cutoff scales Λ 1 = 1 ω P , Λ 2 = 20 ω P such that one of the oscillator is Markovian and basically memoryless. The blue curve describes the symmetric setting when both private baths have the same cutoff scale Λ 1 . That is, both oscillators are nonMarkovian and have equal memory time. In contrast to the previous case, we have coupled (almost) Markovian oscillators (Λ 2 ) in the green curve. Since the damping constant γ is chosen to be 0.5, we find λ PT − (t) associated with the memoryless case settles down to an equilibrium value 0.929 very quickly, compared to the other two cases. Its value shoots up past 1/2 when τ ENT = t 0.197 ω −1 P , even though the bath temperature T B = 0.1 ω P , is rather low. In this case, once the state becomes separable, it remains disentangled. For the other extreme when both oscillators have the same, long memory time (blue curve), the initial entanglement is sustained up to τ ENT 5.863 ω −1 P , and at late times the curve falls below the 1/2 level, whence the entanglement between the oscillators is revived. Finally, in the hybrid case, the duration the system remains entangled happens to fall between two cases of the symmetric setting. We find τ ENT 1.168 ω −1 P , much shorter than the blue curve case but slightly better than the green curve case. At late times, we observe that the curve gradually dips down to 0.571, still above the 1/2 level, so the system is not entangled at late time. This sloping-down behavior, also seen in the blue curve, may imply the memory effect in one of the private baths is still in effect, because the corresponding bath induces a longer effective relaxation time. We observe that for the hybrid system during the relaxation regime, the duration the quantum entanglement is sustained over falls much shorter than the case when both oscillators have the same long memory time, and improves slightly compared to the Markovian, memoryless case. A similar phenomenon has been observed in [2]. There, two coupled oscillators in the two uncorrelated private bath setting are attached to their own Markovian private baths, which have different bath temperatures. In that case, we find that if the temperature in one of the private baths is raised beyond the critical temperature of the order ω −1 P , then the entanglement between the oscillator is lost even though the other bath is kept at temperature much lower than the critical value. Similarly, here since one of the oscillator is memoryless, the coupled system formed purely by this oscillator is supposed to have a short entanglement time τ ENT (the green curve). Thus even it is coupled to a nonMarkovian oscillator, as in the hybrid case, the entanglement duration still fails to gain any improvement. That is, the entanglement duration of a bipartite system is predominantly controlled by the component that least favors the sustainability of system entanglement.

Discussion
In this paper we present a detailed analysis of the effects of nonMarkovian dynamics on the entanglement dynamics between two coupled oscillators, each of which has its own private bath. In particular, we are interested in how the memory effects in the baths affect the full course of entanglement evolution-from the initial time to the moment the reduced system relaxes to the equilibrium state. We break down the analysis in parts to examine the imprints the nonMarkovian baths leave on the solutions of Langevin equation, time evolution of the covariance matrix, and on the negativity, serving as our entanglement measure. Based on the results we obtained we now can answer the questions posed in the Introduction in the following: 1.
How does the nonMarkovianity in the bath affect the system's entanglement? 2.
Which time regime does the bath's nonMarkovianity benefit the system's entanglement most? 3.
Can the memory of one nonMarkovian bath be passed on to another Markovian bath? 4.
Does this memory transfer help to sustain the system's entanglement dynamics?

How Does the NonMarkovianity in the Bath Affect the System's Entanglement?
This is most easily seen when both oscillator-private bath pairs have the same setting. For the homogeneous solution to the Langevin equations, Figure 1 shows that the solution damps at a slower rate, determined by the cutoff scale Λ in the bath's spectral density than what would be given by the damping constant γ. In this paper it is chosen in the strong coupling regime. For the range of chosen parameters, the actual damping rate can be well approximated by the effective damping constant γ EFF in (21). Roughly speaking, a smaller cutoff scale leads to a lower value of the effective damping constant. Thus the original system that is strongly coupled to the nonMarkovian thermal bath seems to behave like a system weakly coupled to a Markovian bath. The phenomenon of the weakened damping may be understood by comparing the memory time and the dynamical time scale of motion. When the memory time is much longer than the period of the oscillator, the contributions from many previous cycles coherently add up to the current state of motion, and thus compensate for the attenuation caused by damping. Alternatively, a cutoff scale, especially a scale smaller than the resonance frequency of the oscillator driven by the quantum fluctuations of the bath, implies that the oscillator does not lose its energy back to the bath as effectively. This characteristic of a weakened damping is likewise passed on to the time evolution of the covariance matrix elements and the entanglement measure, as shown in Figures 3 and 4. These two plots, as well as, Figure 6 and 8 show that the nonMarkovian effect becomes inconsequential at higher bath temperature.

Which Time Regime Does the nNonMarkovianity in the Bath Benefit the System's Entanglement Most?
Coherent superposition of the history of the system's dynamics effectively reduces the damping, and prolongs the relaxation time scale. However, as the system approaches relaxation, its dynamics is gradually taken over by its private thermal bath. The role of the initial conditions is exponentially suppressed in this case. The prolonged relaxation thus lead to the extended intervention of the thermal fluctuations of the bath. This accumulation seems to almost cancel out the benefit from the memory effect. This can be seen from Figure 3. These observations are reflected in the time evolution of entanglement in Figure 4. Suppose that the oscillators are initially entangled, we see that in the intermediate time range, the duration of entanglement is proportional to the memory time, and it lasts about 1/3 of the relaxation time, but at late times when the dynamics reaches the steady state, the value of the symplectic eigenvalue λ PT − of the partially transposed covariance matrix barely benefit from the bath nonMarkovianity, shown in Figure 8. Thus bath nonMarkovianity has a limited capacity to improve the late-time entanglement.

Can the Memory of One NonMarkovian Bath Be Passed on to Another Markovian Bath?
Note this question is not about whether the memory can be faithfully mapped from one party to the other. Rather, we ask whether a system with a shorter memory time can be inducted to a longer memory time when it is coupled to another system that has a long memory time. As far as the dynamics is concerned, the answer seems to be yes. In Figure 10, where oscillator 2, which is (almost) Markovian, is coupled to oscillator 1, which possess a long memory time. The results in the right panel of Figure 10 serves as a contrast, showing the time evolution of the position uncertainty when both oscillators are oscillator 1 (blue curve) or both oscillators are oscillator 2 (green curve). We see that the green curve falls to a constant very fast within the time scale γ −1 . The green curve on the left panel of Figure 10 on the other hand shows the counterpart when oscillator 2 is coupled to oscillator 1. It does not decay as fast as the former, indicating that in this case the effective damping is weaker, or in other words, the "effective" memory time becomes longer. In this sense, the memory of oscillator 2 indeed improves. But, there is a catch. The effective memory time of oscillator 1 deteriorates when it is coupled to oscillator 2, compared to the case when it is coupled to an identical copy (comparing the blue curves in both panels of Figure 10). Therefore, a system with a short memory time can acquire improvement when it is coupled to another system with a long memory time, but, at the cost of the latter.

Does This Memory Transfer Help to Sustain the System's Entanglement Dynamics?
When systems of different memory times are coupled together, from the conclusion brought forth above, we note that the two oscillators still have quite distinct time scales in their dynamics, even though one has improved a bit and the other degraded a bit. The consequence, as shown in Figure 5, is that the entanglement time is improved, compared to the case when both oscillators are oscillator 2, but is shortened, compared to the other extreme. Overall, our results show the sustainability of the bipartite entanglement is determined by the party which breaks off entanglement most easily.
One more question arises when we delineate the role of bath nonMarkovianity. Earlier, we note that the solution to the equation of motion for a system strongly coupled to a nonMarkovian system behaves similar to that for a system weakly coupled to a Markovian bath in the sense that they can have the same decaying rate with minor difference in the phases (shown in top left plot in Figure 9). We may ask whether we can effectively use the latter to approximate the former because the numerical computation of a Markovian system is much less demanding, and the theoretical analysis is simpler? Or pushing it to the conceptual extreme, is it possible that the system we thought to be Markovian is actually a mere mirage of a nonMarkovian system with stronger coupling? In either case, it is then important to examine the similarity and disparity of both descriptions. The additional plots in Figure 9 show that the minor phase difference in the solution d 1 (t) results in broader dissimilarity in the dispersion, uncertainty relation and then bipartite entanglement. This reflects subtleties of the phase carried in a nonMarkovian system. We assert that the effective Markovian description is not sufficiently accurate to capture the needed phase information, so it cannot faithfully mimic a nonMarkovian system.  Together we arrive at an expression of the symplectic eigenvalues of σ in terms of two invariants ∆(σ) and det σ for the covariance matrix σ. Equation (A7) implies that the uncertainty principle becomes The symplectic eigenvalues λ ± are positive because σ is a positive-definite, symmetric real matrix. For the two-mode Gaussian system, we may write the covariance matrix into a canonical form [29] via suitable local symplectic transformation Sp(2, R) ⊗ Sp(2, R), that is, local changes of the canonical variables. This allows us to readily find the Sp(2, R) ⊗ Sp(2, R) invariants [30]. Then we have I 1 = a 2 , For a bipartite system AB, if the state ρ is separable, then the partial transpose of ρ is still non-negative, and vice versa, a bona fide density matrix [31]. We have a corresponding uncertainty relation, The covariance matrix σ PT is the partial transpose of the original covariance matrix via σ PT = Λ · σ · Λ. When the transposition is carried out with respect to B, the matrix Λ has the form Λ = diag(+1, +1, +1, −1), and then σ PT has the form with det C = − det C. Note that Λ is not a symplectic matrix, that is, Λ · Ω · Λ = Ω. We can likewise write down the corresponding Sp(2, R) ⊗ Sp(2, R) invariants, denoted by an extra prime, for the partially transposed covariance matrix. Since Λ 2 = I, we find that they are related to the invariants associated with σ by and the corresponding uncertainty relation is given by Since we arrive at this expression by assuming the state ρ is separable, Equation (A19) serves as the necessary criterion of separability. Following the earlier arguments, we find it implies The criterion (A19) happens to be sufficient for the two-mode Gaussian state. Thus when λ PT − < 1/2, the state is not separable, that is, entangled. This allows us to introduce a quantifiable entanglement measure [32,33], logarithmic negativity, E N by [34] E N = max{0, − ln(2λ PT − )} .
The Gaussian state is entangled when E N > 0. As far as the existence of entanglement is concerned, logarithmic negativity provides the same amount of information as λ PT − does. However, it has the special feature that it indiscriminately gives zero when the state is separable or when λ PT − ≥ 1/2, and thus clips off quite an amount of dynamical information hidden in λ PT − such as effective squeezing and the effective temperature discussed in the main text. Thus we will use λ PT − to gain fuller information about the system dynamics.

Appendix B. Two-Mode Squeezed State
The two-mode squeezed thermal state ρ (β) TMSQ is defined by The operator S 2 is the two-mode squeeze operator S 2 = exp ζ * a 1 a 2 − ζ a † 1 a † 2 with ζ = η e iθ , η ≥ 0 and 0 ≤ θ < 2π. Its action on the annihilation operator, say a 1 of mode 1, gives S † 2 a 1 S 2 = cosh η a 1 − e iθ sinh η a † 2 . (A24) Note that the two-mode squeeze operator is symmetric in a 1 and a 2 , so a similar result for a 2 can be obtained with 1 ↔ 2. The annihilation operators, a 1 and a 2 , of mode 1 and mode 2 satisfy the canonical commutation relation [a j , a † k ] = δ jk with i, j = 1, 2. The nontrivial moments of a i in the two-mode squeezed thermal state are (β) TMSQ = n 1 + 1 cosh 2 η +n 2 sinh 2 η , a † 1 a 1 (β) TMSQ =n 1 cosh 2 η + n 2 + 1 sinh 2 η , a 2 a † 2 (β) TMSQ = n 2 + 1 cosh 2 η +n 1 sinh 2 η , a † 2 a 2 (β) TMSQ =n 2 cosh 2 η + n 1 + 1 sinh 2 η , a 1 a 2 (β) TMSQ = − e +iθ 2 n 1 +n 2 + 1 sinh 2η , a † 1 a † 2 (β) TMSQ = − e −iθ 2 n 1 +n 2 + 1 sinh 2η , wheren i is the average particle number density in the thermal state ρ β . Now given the canonical pair (χ j , p k ) with [χ j , p k ] = i δ jk , which are expanded by a i , a † i via we find the corresponding covariance matrix elements given by This expression is of particular significance. It tells that in the two mode squeezed thermal state, squeezing and thermal fluctuations play competing roles in sustaining quantum entanglement. The effect of thermal fluctuations are reflected in the mean occupation number of the mode, which will grows proportional to the bath temperature in the high temperature limit because stronger thermal fluctuations tend to excite the mode to higher levels. Furthermore, Equation (A45) indicates that even at high temperature regime, the two modes are still able to maintain entanglement when the squeezing is sufficiently large, compared with the contribution from the thermal fluctuations. That is, since the entanglement occurs when λ PT − < 1/2, we find that it is possible only when the squeezing parameter is at least as large as e −2η n + 1 2 < 1 2 , ⇒ η > 1 2 ln 2n + 1 = 1 2 ln coth βω 2 . (A46) In the high temperature limit, the lower bound of the squeeze parameter η to generate entanglement is where T B = β −1 is the bath's initial temperature.