Quantum Entanglement at High Temperatures? II. Bosonic Systems in Nonequilibrium Steady State

This is the second of a series of three papers examining how viable it is for entanglement to be sustained at high temperatures for quantum systems in thermal equilibrium (Case A), in nonequilibrium (Case B) and in nonequilibrium steady state conditions (Case C). The system we analyze here consists of two coupled quantum harmonic oscillators each interacting with its own bath described by a scalar field, set at temperatures $T_1>T_2$. For \textit{constant bilinear inter-oscillator coupling} studied here (Case C1) owing to the Gaussian nature, the problem can be solved exactly at arbitrary temperatures even for strong coupling. 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 NESS setting [1]. Thus lowering the temperature of one of the thermal baths does not necessarily help to safeguard the entanglement between the oscillators. Indeed, quantum entanglement will disappear if any one of the thermal baths has a temperature higher than the critical temperature $T_c$. With the Langevin equations derived we give a full display of how entanglement dynamics in this system depends on $T_{1}$, $T_{2}$ , the inter-oscillator coupling and the system-bath coupling strengths. For weak oscillator-bath coupling the critical temperature $T_c$ is about the order of the inverse oscillator frequency, but for strong oscillator-bath coupling it will depend on the bath cutoff frequency. We conclude that in most realistic circumstances, for bosonic systems in NESS with constant bilinear coupling, `hot entanglement' is largely a fiction. In Paper III we will examine the case (C2) of \textit{time-dependent driven coupling } which contains the parametric pumping type described in [2] wherein entanglement was first shown to sustain at high temperatures.


Contents
A Expressions for V 13 (t), V 12 (t), V 14 (t), V 22 (t), V 32 (t), V 24 (t) 40 A.1 V 13 (t) = χ 1 (t), χ 2 (t) /2: 40 A.2 V 12 (t) = χ 1 (t), p 1 (t) /2: 1 Introduction Recently Galve et al [2] (see also [3]) pointed out the possibility of keeping quantum entanglement alive in a system at high temperatures by driving the system of two oscillators with a time-dependent interaction term. This has generated a great deal of interest in understanding the underlying issues and the basic mechanisms of obtaining the so-called 'hot entanglement' [4]). The word 'hot' conveys three layers of meaning in three different contexts, referring to quantum systems A) kept in thermal equilibrium at all times, B) in a nonequilibrium condition and evolving and C) in a nonequilibrium steady state at late times. Thus before making sweeping statements one needs to discern and analyze systems under at least these three separate situations for the behavior of quantum entanglement therein.
We have analyzed Case B) described above in some detail in our first paper, obtaining the parameter ranges for entanglement to survive at a finite temperature and comparing with the results for Case A obtained earlier in e.g., [5,6]. Our results indicate that, generically, when two coupled oscillators separated at a fixed distance evolve under the influence of a shared thermal bath, their dynamics is usually highly non-Markovian. The asymptotic correlation /entanglement between the oscillators tends to survive better under 1) stronger inter-oscillator coupling, 2) weaker oscillator-bath interaction and at 3) a shorter distance between them. In the case of weak oscillator-bath coupling, the critical temperature is still bounded by the inverse oscillator's natural frequency, but tends to be lower than that the critical temperature in Case C, due to the finite separation between the oscillators. The largest separation before the entanglement drops significantly is of the order of the inverse cutoff frequency inherited in the thermal bath, and the distance will decrease with higher bath temperature. In this limit, the results are similar to Case A. This is not unexpected since it is known [26,27] that in the weak coupling limit, both configurations will yield similar results; furthermore, the non-Markovian mutual interaction between the oscillators is minimal in the weak oscillator-bath coupling regime. For stronger oscillator-bath interaction, the mutual interaction can sustain over a very long history in the evolution of both oscillators. Deviation in results between Case A and Case B will emerge. Nonetheless a strong oscillator-bath interaction can likely induce dynamical instability in the oscillators, a case worthy of closer analysis later.
In this paper we analyze condition C) where the system can maintain a nonequilibrium steady state (NESS) at late times. Since NESS is a distinctly generic state, playing an important role for nonequilibrium systems as fundamental as the equilibrium state in quantum statistical mechanics, it is important to clarify the behavior of high temperature quantum entanglement under such conditions. We illustrate these two conditions with two generic models: Case B) is exemplified by a quantum system made of at least two harmonic oscillators (HO) interacting with one common thermal bath; Case C) is exemplified by a quantum system composed of two coupled harmonic oscillators each interacting with its own (private) thermal bath. We wish to inquire about how entanglement initially present between the two quantum oscillators evolves in time, and calculate at what temperature (approaching from below) it begins to die out.
To identify the root cause of quantum entanglement existing at high temperatures, if it does at all, one needs to identify the determining factors. Coupling in the system is certainly an important factor. Intuitively the stronger the coupling in the system, the weaker the coupling of the system to the baths, the better preserved the entanglement will be. If the coupling can be tuned to "cruise alongside" how entanglement evolves in time, to even amplify it along the way, the better the chance of keeping the entanglement alive. To see these effects more clearly we further divide the nonequilibrium steady state cases into two subcases, C1 and C2. Case C1 is for time-independent inter-oscillator coupling, and Case C2 for time-dependent inter-oscillator coupling. Before one can bring these cases under the same roof of nonequilibrium steady state condiiton one needs to prove or demonstrate that indeed a steady state exists at late times in these setups. We have so far shown the existence of NESS only for Case C1 in [1]. 1 .
Before we treat the Case C1 scenario in full which is the main goal of this paper, we first give a brief description of a Case C2 model to mark the differences so the results of our work can be placed in perspective.
As a model for Case C2 the system is made up of two quantum oscillators interacting with each other via a time-dependent (sinusoidal) coupling. Unlike Case C1 where the temperatures of the two baths are different, here they could be the same. In fact the temperature of the thermal bath and how strong the oscillators are coupled to the baths are not important. The nonequilibrium condition is provided by the external driving agent. Driving leads to production of entanglement even at very high temperatures. For instance, even with a weak environmental coupling, a strong driving amplitude still provides a higher critical temperature.
The physics for these two cases albeit both in NESS is also very different. As explained in [2] , it is the squeezing of the system provided by the external agent and the parametric amplification (pumping) which can offset the thermalization /equilibration process naturally expected for the systems interacting with a bath and dominant at high temperatures.
Parametric driving is what sustains the entanglement in the system. We will study this case in our sequel paper.

Time-independent bilinear inter-oscillator coupling
In the case of a chain of quantum harmonic oscillators coupled bilinearly with each other and with the baths the dynamics of the total system admits a complete solution, by virtue of its Gaussian nature, for all temperatures and for strong coupling within the system and with the baths. This model has been studied by many authors [10-12, 14, 15]. In our recent work [1] functional methods are used to provide an explicit demonstration of the existence of a nonequilibrium steady state. Here we apply the results obtained therein to a study of quantum entanglement in NESS, with the aim of quantifying the claims made in the literature [4], alluding to the possibility of entanglement survival at high temperatures for systems in NESS. Note the present setup of bilinear coupling is different from that of Galve et al [2] where the interaction between the two oscillators is via parametric pumping. For this setup a recent paper closest to our intent is that by Ghesquire, Sinayskiy & Petruccione [7].

Comments on Claims by Other Authors
We make a brief summary of what GSP have done and what claims they made below. For the same model as mentioned above, namely, two bilinearly coupled quantum harmonic oscillators each interacting with its own bath GSP derived a perturbative 'pre-Lindblad' master equation without invoking the rotating wave approximation (non-rotatingwave, or NRW) [8]. They consider two situations: For the study of entanglement they consider the high temperature regime in their Eq. (3) valid for both strong and weak interaction strength with the baths. For the consideration of entropy dynamics related to equilibration issues they take the weak system-bath coupling limit and arrived at their Eq. (4). We will only be concerned with the entanglement issue here. GSP made the following claims: a) Entanglement persists for longer times at lower temperatures. b) In the weak system-bath coupling limit, the late time steady state developed is independent of the initial conditions. c) For the equilibrium case, there exists a critical temperature which is consistent with the result of [5] in the limit.
We limit to two comments regarding their method and claims here. The major differences will become clear in our results with quantitative representation via graphs found in later sections.
1) Regarding the method and approximations: A perturbative 'pre-Lindblad' master equation, even without the RWA, does not in general satisfy the complete positivity condition. Although it works better for strong coupling to the environment the results obtained under these approximations have unphysical behavior at low temperatures. For example, Ludwig et al [12] pointed out the effect from the environment cutoff has to be handled with care.
2) The claim statements are too general -they may not hold for specific conditions. They need be qualified more carefully by specifying the range of (in)validity of the approximations introduced. E.g., Point a) above is sort of expected, but does it also imply that entanglement can be generated and be sustained if the temperature of both baths are sufficiently low, even though the system state is initially separable? Point b) regarding the existence of a NESS -it has been demonstrated for arbitrary strength in bilinear interoscillator coupling and for arbitrary temperatures of the two baths [1]. Point c) There is a distinction between i) a system of two coupled oscillators each with its private bath under NESS studied here, setting the two baths to be at the same temperature (presumably what their 'equilibrium' condition entails) and ii) the system in one common thermal bath (what we call Case A). The situation is a lot more complex -see discussions in the last section of this paper.
The above questions and a broader set of issues will be addressed in a fuller treatment of this generic (bilinear coupling) NESS model in the sections below.

Our Method and Key Findings
The model we use in this work to describe entanglement dynamics at high temperature, namely, two coupled oscillators each interacting with its private bath at different temperatures, has been treated in full in our earlier paper [1], where one can find more technical details of the whole framework. Entanglement in a harmonic chain is also a well-explored subject . The Gaussian nature of this model allows us to obtain exact solutions for arbitrary coupling strengths and temperatures . The central quantity to calculate is the covariance matrix at finite temperature and at late times, where it has been shown that the system approaches NESS. The Peres-Horodecki-Simon entanglement criterion [16][17][18] can be calculated without approximation. This approach has been shown to be totally equivalent to that of directly deriving the reduced density operator of the system [1,19]. A short way to report on our findings is that quantum entanglement will disappear when the bath temperatures become higher than a critical temperature (T c = 1/β c ). Also not surprisingly, asymptotic entanglement is easier to sustain for stronger inter-oscillator coupling and weaker oscillator-bath coupling. The true gain of this investigation is a full display via the Langevin equations we derived of the dependence of entanglement dynamics on the three parameters in this model, temperatures (T 1 , T 2 ) of the baths, the intra-system (inter-oscillator) coupling σ and the system-bath coupling strengths γ. Their interplay is presented in the plots, where the critical temperature dependence on different coupling strengths can be easily seen. For the special case when both baths have the same temperature, we show that the critical temperature, above which the system becomes separable, satisfies β c ω ∼ 2 1 + 4σ/3ω 2 −1 for weak oscillator-bath coupling, ω being the oscillator natural frequency. It is consistent with the general expectation that β c ω ∼ O(1) in the vanishing inter-oscillator coupling σ limit. In the opposite limit, when the oscillator-bath coupling is strong, correction terms with bath cutoff frequency dependence will show up. This is a noteworthy point in a lesser-explored regime, namely, one needs to be mindful of the choice of the environment cut-off frequency in the treatment of strong system-environment coupling.
A cautionary remark is in place here about entanglement measures used for quantum systems at finite temperature: Although the Peres-Horodecki-Simon (PHS) criterion is totally valid to identify the existence of entanglement in a quantum system, it does not serve as a quantifiable measure. We find from explicit calculations that at finite temperature it does not necessarily vary monotonically with the parameters in our system, namely, the temperature or coupling constants. One should exercise caution in using the PHS criterion for a physical understanding of thermal entanglement. In contrast negativity is a valid measure to quantify the dependence of quantum entanglement on these physical parameters.

Differences from the common bath case
To highlight the qualitative features in the behavior of the separability criterion it is useful to contrast the private bath case (Case C1) studied here and the shared bath case (Case B) studied in Paper I. A more detailed description can be found in the last section: 1. The initial Gaussian conditions will be irrelevant in the private bath case, but remain significant in the shared bath case, so the state of entanglement is sensitive to the initial conditions.
2. At late times the entanglement measure for the private bath case is time-independent, but it continues to oscillate in time.
3. The inter-oscillator coupling (σ > 0) plays a more important role in the private bath case than in the shared bath case.
4. In the private bath case, entanglement is easier to survive for stronger inter-oscillator coupling and weaker oscillator-bath coupling, but in the shared bath case, both factors seem to be overshadowed by the intrinsic quantum dynamics of the system which depends on the initial conditions of the oscillators.
This paper is organized as follows: In Sec. 2, we briefly discuss the dynamics of the reduced system in the NESS configuration, and introduce the separability/entanglement criterion. In particular we pay attention to the covariance matrix, which constitutes the building blocks of the separability criterion. In Sec. 3 and 4, we highlight the calculations of the covariance matrix elements at high, zero and low temperature cases. We further examine the temperature dependence of the covariance matrix elements and the validity of the relevant approximations in Section 5. In Sec. 6 we investigate the separability criterion at different temperature regimes in detail and point out its non-monotonic behavior. Because of this we adopt instead negativity as a valid measure of entanglement for quantitative analysis of quantum systems at finite temperature. We derive some relations between the critical temperature and various coupling constants. In Sec. 7, we then offer a more intuitive viewpoint to understand how all sorts of interactions can affect entanglement between oscillators. Finally we summarize our results and compare them with the case of the shared bath in Sec. 8.

The Model
Consider two coupled harmonic oscillators of equal mass m and (bare) natural frequency ω b coupled to each other with strength σ, each of which interacting with its own thermal bath with coupling constant e. We refer to the two oscillators together as the system, and the two baths together as the environment. This setup is a prototype used often for the investigation of nonequilibrium steady state (NESS), the existence of which is shown in a recent paper [1] (see also the references therein). In the Langevin equation approach the two oscillators' amplitude χ 1 , χ 2 satisfy the following equations of motion: where γ is the damping constant related to e by γ = e 2 /(8πm), and ω is the renormalized frequency (wherein the correction from the interaction with the environment has been considered before), and ξ 1 , ξ 2 are the stochastic forces acting on Oscillators 1, 2 (O 1,2 ) respectively. Note they are not specified by hand but determined self-consistently. An overdot denotes taking the time derivative of a variable. The initial state of the oscillator is described by a Gaussian wavepacket and both oscillators are prepared in the same initial configuration. The two private baths (B 1,2 ) are modeled by massless scalar fields at different temperatures β −1 i . In the matrix notation, these two Langevin equations are condensed into one, namely, where The solutions to this equation are given by, where χ χ χ(0),χ χ χ(0) represent the initial configuration of the oscillators. The fundamental solution matrices D 1 , D 2 are a special set of homogeneous solutions to the Langevin equation (2.3), In particular, the Fourier transformation of is equal to θ(τ ) D 2 (τ ), that is, The function θ(τ ) is the unit-step function. Unless mentioned otherwise, we will not distinguish θ(τ ) D 2 (τ ) from D 2 (τ ) for all practical purposes, and denote −κ 2 I + Ω Ω Ω 2 − i 2κ I −1 by D 2 (κ).
The force term ξ i (t) is a stochastic c-number with the statistical properties , is the Hadamard function of the bath scalar field, associated with the ith oscillator [1]. This stochastic force in essence represents the quantum fluctuations of the private bath at a finite temperature.

Entanglement Measures
For continuous-variable systems, the entanglement measure based on the density matrix is not conveniently calculable because the density matrix in this case is infinite-dimensional. However, it has been shown [18] that in the case of continuous Gaussian variables, the Peres-Horodecki separability criterion [16,17] can be reformulated in terms of the covariance matrix of the bipartite system, Here the matrices A, B, C are the block matrices in the covariance matrix V, while the covariance matrix V itself is defined by the canonical variables of the two subsystems where ρ is the density matrix of the state we are interested in. We have assumed R = 0. The column matrix R takes the form R T = (χ 1 , p 1 , χ 2 , p 2 ), and p i is the canonical momentum conjugate to χ i associated with the subsystem i. The angular brackets denote taking the quantum expectation value. In our case, once we have the covariance matrix for the coupled harmonic oscillators in the NESS configuration, we may construct ζ + according to (2.10). A negative value of ζ + thus implies the existence of quantum entanglement.
Although (2.10) constitutes only the second moments of the canonical variables, it offers a complete description of the Gaussian system since for a Gaussian system, all higher moments can be expressed in terms of the second moments. Oftentimes it is instructive to write the Peres-Horodecki separability criterion in terms of the symplectic eigenvalues of the partially transposed covariance matrix. Let η ≷ stand for the symplectic eigenvalues of V pt , the partial transposition of V. Without loss of generality we assume η > is greater than η < . In fact they can be found by solving the eigenvalues of the matrix i Ω Ω Ω · V pt , with Ω Ω Ω = 2 k=1 J. The resulting eigenvalues will appear in pairs by the form ±η > , ±η < , so the symplectic eigenvalues of V pt are given by the absolute value of the eigenvalues of i Ω Ω Ω · V pt . When we write V pt into the Williamson's form, the separability criterion V pt When η < < 1/2, entanglement occurs. Notice that η > is assumed to be larger than η < , so η > is always greater than 1/2. We observe that although a violation of the Peres-Horodecki-Simon separability criterion signals the existence of entanglement, it is not a good measure for a quantitative description of entanglement, in that the criterion includes a unwelcome factor (η > − 1/2), which does not affect the identification of the existence of entanglement, it messes up the correct evaluation of entanglement. This can be understood if we examine the behavior of the symplectic eigenvalues η ≷ about η < ∼ 1/2. For definiteness, we assume that the symplectic eigenvalues are similar monotonic functions of the parameters of the entangled system. We can easily see that if η > changes too fast in the vicinity of η < = 1/2, the product η 2 > − 1/4 η 2 < − 1/4 will not be monotonic there. As is perhaps better known, a simple calculable measure of entanglement which also provides quantifiable accuracy is negativity [20], denoted by N or its logarithm (strictly speaking logarithmic negativity is not merely the logarithm of negativity, although it is related to) [21], the logarithmic negativity E N . For the Gaussian states under study they can be respectively defined by 14) in terms of the symplectic eigenvalue η < of the partially transposed covariance matrix. When η < < 1/2, the Gaussian state ρ is entangled and both measures take nonzero values between 0 + to +∞. In addition, the logarithmic negativity has a convenient feature of being additive. Comparing the negativity (2.14) with the Peres-Horodecki-Simon criterion (2.13), we observe that they are all based on the smaller symplectic eigenvalue η < of the partially transposed covariance matrix V pt , so they will give the same prediction on the existence of entanglement. However, the separability criterion carries an additional undesired factor (η > − 1/2), which may inadvertently scale (η < − 1/2). Thus the separability criterion is not suitable for quantifying entanglement.
Finally, we remark on a subtlety of the entanglement measure. It has been pointed out [22][23][24] that different measures may give different ordering of density operators with respect to the amount of entanglement. To be more specific, given two density matrix ρ 1 and ρ 2 , we can have E 1 (ρ 1 ) ≤ E 1 (ρ 2 ) for one entanglement measure, while E 2 (ρ 1 ) ≥ E 2 (ρ 2 ) for another. In particular, negativity and Gaussian entanglement of formation, the latter forming an upper bound to the true entanglement of formation, have been found to be inequivalent for asymmetric Gaussian states [24]. For symmetric states, the predictions from both measures coincide.
The next few sections will be dedicated to the calculation of elements of the covariance matrix.

Elements of the Covariance Matrix
We use (2.5) to find the elements of the covariance matrix V. Assume that the initial state of each oscillator is depicted by a Gaussian wave packet of the same shape, at rest initially at the bottom of the harmonic potential associated with each oscillator, such that with p i = mχ i . Thus these two oscillators are initially in a separable state. From the solutions (2.5) one can identify the role of the interaction, either between the oscillators or between the oscillator and its private bath, in creating or sustaining the quantum entanglement in the system. To calculate the elements of the covariance matrix V one can show, for example, that When the dynamics of the system evolves into relaxation as t → ∞, the first two terms on the righthand side will be exponentially small if the coupling constant between the oscillator and the bath is not vanishing. Thus at late time χ i (t), χ j (t) /2 simplifies to where we have used the fact that D 2 (τ ) = 0 if τ < 0. Since the Fourier transform of D 2 (s) is defined by we use the property D 2 (−κ) = D * 2 (κ) to arrive at (2.18).
At this point, let us look at a more specific example: the element V 11 (t) = χ 1 (t), χ 1 (t) /2 = χ 2 1 (t) . At late time it takes on the value V 11 , and with γ = e 2 /(8πm). The frequencies ω 2 ± = ω 2 ± σ are the oscillating frequencies of the normal modes, which can be constructed from the superpositions of (2.1) and (2.2). The Fourier transformation of the Hadamard function takes the form (2.23) The term κ/4π represents the vacuum zero-point contribution. The off-diagonal terms of G H are zero because both private baths are not correlated. From the late-time value V 11 of the amplitude uncertainty of O 1 , we observe the following distinct features: (1) it approaches a constant independent of time, (2) its integral expression (2.20) takes a form similar to the Landauer formula, where | D 11 2 (κ)| 2 plays a role of the transmission coefficient, and (3) it depends on both thermal baths even though O 1 does not have a direct contact with B 2 . The last property would not be unexpected because the coupling between the oscillators will bring in correlations between O 1 and B 2 , and vice versa, between O 2 and B 1 . In fact, these features hold quite generally for the all elements of the covariance matrix in the current NESS configuration.
The definition of the covariance matrix V and the expressions for its elements, and their corresponding values at late times are derived in Appendix A-D.
In the next sections we will explicitly evaluate the elements of the covariance matrix for three situations: (1) high temperature limit, (2) zero temperature case and (3) low temperature regime.

The Covariance Matrix at High Temperatures
We consider the high temperature limit βω 1 of the elements of the covariance matrix. In this limit, the Hadamard function of the bath (2.23) is approximately given by We see that its vacuum contribution is relatively negligible, and can be neglected for most cases. However, extra discretion is advised for the evaluation of the momentum uncertainty where the vacuum contribution of the bath can be significant when the coupling between the oscillator and the bath is sufficiently strong. Thus the result can depend on the cutoff scale of the environment field) (see, e.g., [12]).
Here we merely highlight the calculation for the element V 11 at late time. To obtain the high temperature limit of V 11 (∞), that is, V 11 , essentially we evaluate the following two integrals In terms of I 1 and I 2 , we see from (2.20) that the high temperature limit of the element V 11 at late time is given by Here we would like to point out that when the mutual interaction σ is large, in particular when σ → ω, the fluctuations of the oscillator grow significantly. This will be traced back to the small values of ω − . We will come back to this feature in due course. Derivations of the high temperature forms of V 13 , V 14 , V 22 , V 24 are given in Appendix B. Nonetheless for the following discussions we will bring forward the results for V 22 and V 13 here. When both private baths have the same temperature β −1 , we have from (5.1), (5.3) and (5.5) in the weak oscillator-bath coupling limit. This implies that the average harmonic potential energy of Oscillator 1 (O 1 ) will be It is a bit off from the value 1/2β one would expect from the equipartition theorem for a free harmonic oscillator in the high temperature limit. This difference will disappear when the mutual coupling σ between the two oscillators are turned off.
Eq. (3.6) on the other hand tells us the corresponding average value of the kinetic energy in the high temperature limit, is the same as the value obtained from the classical equipartition theorem. We observe that in the high temperature limit the mean kinetic energy is not equal to the mean harmonic potential energy in general, and the sum of the kinetic energy and the harmonic potential energy is not equal to kT : Let us compare this with the average total energy of a free harmonic oscillator in a closed system, in the high temperature limit and E n = n + 1 2 ω. The deviation can be accounted for by the fact that some portion of the total energy of both oscillators is stored in the mutual interaction between O 1 and O 2 . Accordingly the missing piece should come from the expectation value of mσ χ 1 χ 2 . Its contribution to the mechanical energy is when β 1 = β = β 2 . Including this contribution we see the total energy for the two-oscillator system in the high-temperature limit becomes which is that obtained by the classical equipartition theorem for two coupled linear oscillators. This also serves as a consistency check of our calculation. Finally we comment on two issues. First, weak oscillator-bath coupling enables us to ignore the cutoff-dependent effect from the bath. This may not be true in the strong coupling case. Second, despite the resemblance of (3.13) with (3.11), they are quite different in the physical context. The former is considered in the context of open systems while the latter under the assumption of a closed system. It has been shown [13] that both results can be equivalent only in the limit of vanishingly weak oscillator-bath coupling.

The Covariance Matrix at Zero and Low Temperatures
Here we evaluate the vacuum contributions and the low temperature correction of the covariance matrix elements. Due to the zero-point fluctuations of all bath modes, the vacuum contributions of some covariance matrix elements can be divergent. Suitable cutoffs need be introduced to regularize them, with due consideration of the particulars of the bath the system interacts with.
Let us examine, for example, V 11 = V 11 (∞) and work out its zero and low temperature expressions.
so we need the following two integrals to evaluate the vacuum contribution of V 11 , The sum of J 1 and J 2 can be expressed as where the dimensionless function f (z) is defined by It is clear that Ω ± are the resonance frequencies of the two normal modes. Therefore from (2.20), the zero-temperature (vacuum) contribution of V 11 is given by We observe that the vacuum contribution can be clearly separated into decoupled components of two normal modes, with oscillating frequency ω ± respectively. This is another general feature of this system. The zero-temperature expressions for V 13 , V 14 , V 22 , V 24 are given in Appendix C.
The low temperature corrections to the covariance matrix elements basically result from the corresponding correction in the Hadamard function, because the fundamental solution matrix D 1,2 does not depend on temperature. This is a consequence of the fact that the retarded Green's function of the scalar field, which accounts for dissipation in the Langevin equation, is state-independent. As is seen from (2.20), we need the following two integrals to evaluate the low temperature correction of V 11 , We then have the low temperature correction to V 11 given by However this is merely the contribution from the first term in the summation of all finite temperature corrections in (4.7). Since the remaining terms (n > 1) will have algebraically comparable contributions, we have to take them into consideration. We note that the leading term in K 1,2 has a temperature dependence β −2 in the low temperature limit. Thus we expect that the leading contribution for the remainder of the finite temperature corrections in (4.7) should be proportional to n −2 β −2 . Their overall contributions will introduce an addition factor to (4.10). Therefore after taking this into account, we obtain the low temperature correction to V 11 as follow: The low temperature expression of V 11 is then given by the sum of (4.6) and (4.12), We leave the derivations of the zero and the low temperature expressions for Mathematically speaking, the inclusion of the unit-step function θ(βΛ − 1) is to account for the vacuum contribution of the bath modes in the case βΛ > 1, because when βκ > 1, the Hadamard function G kk H (κ) takes the low-temperature form as shown in (2.23). On the other hand when βΛ < 1, the high-temperature approximation of G kk H (κ) is entirely valid up to the cutoff scale, the cutoff-dependent term being subdominant. However on physical grounds, since the cutoff scale by construction is the highest energy scale compatible with the model, the thermal excitation energy thus must be smaller than the cutoff scale. It then implies that even in the high temperature limit, we still have ω β −1 < Λ. Here we show the high/low temperature approximations of V 22 with a numerical calculation. In particular we explicitly highlight the role of the vacuum contribution, that is, the cutoff dependent terms, even in the high temperature approximation for strong oscillator-bath interaction. The pink curve in the plot on the right shows that if the vacuum contribution of the bath is not taken into account, the analytical high-temperature approximation will be way off from the numerical result (the purple curve) in the region βω ∼ O(1). On the other hand, the red curve, which includes the vaccuum contribution, fits nicely with the numerical result. The parameters are chosen to be γ = 0.2, σ = 18, ω = 5, and Λ = 1000. The plot on the left is drawn for weak oscillator-bath interaction γ = 0.2, i.e. γ/ω 1. The cutoff-dependence is seen as dispensable.
In this case, since the leading contribution of the high-temperature approximation vanishes, we have to include the next-order term.
it is similar to V 11 except that we swap β 1 and β 2 , Some comments are in place here: Both oscillator are initially prepared in a state of nonoverlapping Gaussian wavepackets with the same width ς. As they come into interaction with their own private baths, the evolution of each individual oscillator is then driven by its bath and the other oscillator it is directly coupled with. Due to the dissipative self-force on the oscillator arising from its interaction with its own bath, the intrinsic information of the initial state is dispersed away exponentially fast as the system evolves in time. In the end when t → ∞, the values of the dynamical variables of the oscillator are determined by its private bath and by the other oscillator. We want to bring up this point because even when the oscillator-bath coupling constant γ approaches zero, not all of the asymptotic values of the covariant matrix elements are zero. In this limit their values are independent of the parameter ς characterizing the initial state, so they are not related to the intrinsic evolution that begins with the initial configuration. Instead they are the induced components as a consequence of the interaction between the oscillator and the bath. In other words, the results of the covariance matrix in the vanishing γ limit should be understood by the limiting procedures of taking t → ∞ first and then taking γ → 0. This is a good point to discuss in more details in what is meant by the high/low temperature approximations. We only cover the generic situation and discard some extreme cases such as ω − , Ω − , σ → 0, so we assume that σ 1 2 , ω ± and Ω ± are about the same order of magnitude as the oscillating frequency ω. The cutoff frequency is assumed to be much larger than ω, i.e., Λ ω. The magnitude of the parameters γ and β 1,2 are not restricted as long as Ω − remains well-defined. We use V 11 and V 22 as illustrating examples, 1. V 11 : as far as the order of magnitude is concerned, we see vacumm: Roughly speaking, the high temperature limit refers to the case βω 1; on the other hand a consistent low temperature correction requires γ ω 1 (βω) 2 1 , which can be weaker than the naive low temperature limit βω 1, especially in the weak coupling limit γ/ω 1. It implies that in the weak oscillator-bath coupling limit, the low temperature correction has a much wider range of validity. In the strong coupling regime γ ω, the low temperature correction is remains fully valid for the βω 1.
Here, additional subtlety arises due to the presence of the cutoff frequency Λ. The importance of the cutoff-dependent term relies on how the factor γ ω ln Λ ω is compared with unity. In the weak coupling limit, the cutoff dependent term is negligible, so we can safely ignore it unless the cutoff frequency is extremely high, such as Λ O(ω e ω γ ) .
In the strong coupling regime γ ω, we see that the cutoff-dependent term still has a comparable magnitude relative to the high temperature approximation in the interval of the high-to-low temperature transition βω O(1). This interval has a special significance because, as we will see later, this is the region where thermal entanglement may disappear in the nonequilibrium steady state configuration.
Thus at this point, generically speaking, the high-temperature limit refers to βω 1 while the low-temperature limit refers to βω 1. For weak oscillator-bath coupling, the low temperature correction has a wider range of validity than is implied by βω 1 due to the additional factor γ/ω in the corresponding expression. In addition, the cutoff is completely negligible in normal circumstances. By contrast, in the strong coupling regime, the cutoff-dependent contribution enters in determining the critical temperature of thermal entanglement.
6 Entanglement of System in Nonequilibrium Steady State

Late Time Behavior of the Covariance Matrix
At late time when the system reaches the steady state, the covariance matrix takes the form 3) The determinants of the matrices A, B are related to the generalized uncertainty relation for each single subsystem, which also takes into account the correlation between canonical variables. The matrix C contains the cross-correlation among canonical variables between two subsystems. As is briefly discussed in Sec. 2.2, the knowledge of the covariance matrix enables us to use the Peres-Horodecki-Simon separability criterion to determine the quantum entanglement. In fact the separability criterion can be combined with the generalized uncertainty relation to form an unified statement The expression containing the − sign represents the uncertainty relation while that with the + sign represents the separability criterion. We immediately see that in the current case det A and det B are always positive definite by construction. In addition, the expression Tr A · J · C · J · B · J · C T · J , once written explicitly in term of the covariance matrix elements, is found to be always positive. Thus when we rewrite (6.4) as, we immediately recognize that ζ ± actually contains two positive but competing components. This makes it difficult to determine the sign of ζ ± . However, we can use the following argument: Suppose that the uncertainty relation ζ − ≥ 0 always holds. Since the condition ζ − < 0 implies that det C must be negative. Therefore the appearance of negative values of det C may help to signify the existence of entanglement. The sign of det C is less clear, depending on how negative V 13 V 24 can be allowed. Although this is not a sufficient condition, it highlights the role of cross-correlations in entanglement. Before we proceed to evaluate ζ + , we observe that among the elements of the covariance matrix, two of them, V 22 and V 44 , have dependence on the cutoff frequency Λ, which is the highest energy scale that is consistent with the theory. Thus we expect that ζ ± , and in particular, the separability criterion, will depend on the cutoff scale. Since the cutoffdependent term always has the form γ ln Λ, where γ is the system-environment coupling constant, it implies that this cutoff dependence will be suppressed in the weak coupling limit. However, when the system interacts strongly with the environment in the sense that γ/ω is close to unity, the contribution from the factor γ ln Λ can be significant, and can make the separability criterion ambiguous.
Likewise in terms of the matrices A, B, and C we can construct the symplectic eigenvalues η ≷ of the partial transpose V pt of the covariance matrix V [25], where alternatively det V can be written as det A det B + (det C) 2 − Tr{A · J · C · J · B · J · C T · J}. This enables us to calculate the quantitative entanglement measures like negativity or logarithmic negativity for the Gaussian state.
In the sections that follow, we will refer to the special case when both thermal reservoirs have the same temperature. In this case, the Gaussian state becomes symmetric, so (logarithmic) negativity will give an unambiguous ordering of density matrices, in comparison with other quantitative entanglement measures. Since we have A = B, and the matrix C becomes diagonal, the symplectic eigenvalues η ≷ takes a particularly neat form We readily see that are associated with the dynamics of the normal modes of the joint system. This elicits a transparent connection between the entanglement behavior and the underlying dynamics.

Entanglement Behavior
As stated earlier the Peres-Simon-Horodecki criterion can be used to identity the existence of entanglement of the Gaussian states, but it may be inadequate to provide a quantitative description of entanglement, in particular, for a quantum system at finite temperature. We will show later that it does not vary monotonically with temperature and coupling constants. This discrepancy comes from an additional factor in the criterion. It has no effect on the identification of entanglement but it will give an unwarranted bias on the values, rendering it inappropriate for quantifying entanglement. While the separability criterion can be used for the system under study at zero temperature, we need a different measure to quantify finite-temperature entanglement, namely, negativity.  Figure 2. The separability criterion ζ + is plotted against the mutual coupling strength σ between the oscillators at zero temperature. Larger values of the damping constant γ will move the curve upwards and make the entanglement between the two oscillators harder to sustain. The oscillating frequency ω and the cutoff frequency Λ are chosen to be 5 and 10000, respectively

zero temperature
We first examine the separability criterion ζ + at zero temperature. The whole expression for ζ + can be exactly found but it is tremendously large. Here we present the leading terms in the weak oscillator-bath coupling limit, i.e., γ < ω ± is the smallest parameter at hand, Note that it is not sufficient to expand ζ + to first order in γ because they all depend on (ω + − ω − ). This factor will make the first-order expansion of ζ + vanish when σ → 0 no matter what value γ has. In fact ζ + has a finite value when γ = 0, so we have to include terms which are at least of second order in γ.
In addition, as far as the leading contribution is concerned, we see that (ω + − ω − ) 2 is always positive, so ζ + is negative for all nonzero mutual coupling strength σ between the two oscillators. This implies that the oscillators will become entangled aympotically once they are coupled. On the other hand, when we consider contributions due to the finite value of the damping constant γ, we find that, in particular in the limit σ → 0, we have ω + → ω − and The separability criterion ζ + is positive for σ = 0 when γ = 0. With increasing σ, the value of ζ + gradually falls below zero at some critical value of σ c . Thus we see that the curve of the separability criterion will move upwards with larger values of the damping constant γ, that is, with stronger interaction between the oscillator and its private bath. Furthermore, it also indicates that these oscillators are not always entangled, and they can be separable for some choices of γ and σ. For a specific value of γ, the mutual coupling strength σ must be greater than the critical value to render both oscillators entangled. In other words, the bonding between two oscillators has to be strong enough to overcome the incoherent disturbance from their respective baths in order to maintain their entanglement. The larger the values of σ the more stable the mutual entanglement is. Therefore we see that the oscillator-bath interaction and the coupling between oscillators play competing roles in sustaining the entanglement. We now derive a relation between the critical values of different couplings. For the case of a small damping constant γ, the critical value σ c can be obtained by solving (6.13) with ζ + = 0, yielding Inverting it leads to Therefore in the weak coupling regime where γ ω ln Λ ω is small but not vanishing, we find that if the two oscillators are initially prepared in a disentangled state, they can become entangled for sufficiently strong direct mutual coupling between the oscillators. Otherwise, they remain asymptotically in a separable state when the mutual coupling is weak. Finally, we add some remarks on the cutoff dependence of the separability criterion. From (6.16), we see that the dependence on the cutoff always occurs as long as γ = 0. This implies  some discretion is needed in the treatment of the cutoff scale. If we ignore this contribution, then one will encounter the following unphysical situation: If the oscillators are initially in a separable state , prepared as Gaussian wavepackets, their final state is always, at least marginally, entangled even though the mutual interaction is turned off. In contrast, if the cutoff contribution is accounted for, then the final state of the oscillators will not be entangled unless their mutual interaction is strong enough. Secondly, the cutoff scale always appears in the form ln Λ, so the separability criterion is not very sensitive to the choice of the cutoff scale unless it takes some extreme values. In Fig 3, we let the cutoff scale Λ go up to a very high value relative to ω. We see that the separability criterion ζ + become always positive above a critical value Λ c , and Λ c is highly sensitive to the choice of γ. Comparing the two plots in Fig 3, we see a mere change in γ causes a dramatic shift in the value of Λ c . Generally speaking, only for very weak oscillator-bath coupling will the cutoff-dependent terms play a subdominant role in the separability criterion.
So far we have presented the general features in how the separability criterion depends on the interactions. We now investigate the role of temperature in the criterion.

low temperature βω 1
Generally speaking, with increased temperature, thermal fluctuations will become increasingly important in affecting the dynamics of the oscillators from their respective baths. Quantum coherence is expected to deteriorate accordingly. We expect similar degradation may occur in entanglement. Thus it is reasonable to conjecture that once the temperatures of the baths are raised above a certain critical value, the degradation can be so severe that the oscillators become completely separable. However, the situation is more complicated for the present setup because two independent thermal baths are involved. It turns out that lowering the temperature of one of the thermal baths does not necessarily guarantee entanglement between the oscillators. Thus the concept of a universal critical temperature is less well-defined in multiple bath situations.
Here, we discuss the functional dependence of the separability criterion ζ + on temperatures. To begin with, let us suppose that it takes on a generic form ζ + = ζ(β 1 , β 2 ). When a steady state is reached, the separability criterion should be invariant under the exchange of β 1 and β 2 because the configuration of the total system is designed to be symmetric when we swap one oscillator and its private bath with the other oscillator and its private bath. This implies that ζ + (β 1 , β 2 ) = ζ + (β 2 , β 1 ). However, it is unlikely that the temperature dependence of the separability criterion can be reduced to a function of |β 1 − β 2 | solely. If ζ + were a function of |β 1 − β 2 |, it would imply that the separability criterion could be independent of temperature for the case β 1 = β 2 where it would further suggest that both oscillators should be either separable or entangled for all temperatures. We have shown that at least they can not always be separable because in the zero temperature case, we found that both oscillators can be entangled for certain choices of parameters. On the other hand, it is hard to believe that both oscillators remain entangled even at very high temperature. Thus we rule out the possibility that the separability criterion may depend on |β 1 − β 2 |. The same features are also shared by the symplectic eigenvalue η < , as can be seen in Fig 5, but there are two differences: η < is monotonic with respect to the parameters of the joint system and it does not rise up as steeply as the separability criterion in the high temperature regime. The latter is related to the extra factor (η 2 > − 1/4) in the criterion. Furthermore we observe that even for β 1 = β 2 where the reduced system is described @by?@ asymmetric Gaussian states, the symplectic eigenvalue η < , thus negativity, still gives a consistent and physical picture with respective to the ordering of the density matrix in terms of the relevant parameters in question.
Thus, to define more precisely a critical temperature β c , we will look at the special case of β 1 = β 2 . In this case both private reservoirs have the same temperature, yet they are totally uncorrelated. This setup is still distinct from the case that two oscillators share a common bath. In the shared bath case, the oscillators can influence each other indirectly through their interaction with the same bath, whereby non-Markovian effects enter in their dynamics, with dependence on their spatial separation (see, e.g., [9]) . This type of effects are absent in the private bath configuration; nonetheless, other than the direct influence from its own bath, each oscillator can still experience, by means of direct coupling between the two oscillators, the action of the other bath associated with the other oscillator. Therefore the equal-temperature private baths and the single common bath configurations are not the same, but, as we shall see, there are some similar features. Moreover, in this special case the two-mode Gaussian state becomes symmetric so the negativity can give an unambiguous comparison of entanglement between states.
In the low temperature limit, we find the finite temperature correction to the separability criterion is given by where ζ (0) + is the zero-temperature result in (6.13). It is interesting to note that the correction may change sign as the inter-oscillator coupling σ varies from zero to its upper limit. The upper limit of σ is constrained by the condition Ω − = ω 2 − − γ 2 = 0, so σ max ∼ O(ω 2 ). When σ = 0, the term linear in γ vanishes due to ω + = ω − there, but the term quadratic in γ is positive. Hence the correction starts off with a positive value. On the other hand in the limit ω − → γ (i.e., σ → σ max ), we find the finite temperature correction gradually becomes negative where we have used the L'Hôpital's rule to evaluate the limit of such an expression This feature reveals the non-monotonicity of the separability criterion at finite temperature. We stress that this errant behavior does not affect us from reading off the critical values of the parameters. 5.× 10 -6 Figure 6. The separability criterion ζ + is plotted against the inter-oscillator coupling σ and the oscillator-bath coupling γ at low temperature. In each plot, we show the ζ + curve for three different bath temperatures. We see that the critical temperpareue is higher for stronger σ but weaker γ.
The oscillating frequency ω and the cutoff Λ are chosen to be 5 and 10000, respectively Next we look into the effect of finite temperature correction on the critical value of σ c where ζ + transits from a positive to a negative value. In the zero temperature case, we have found this critical value in (6.16), now denoted by σ which is always positive. It means that this correction shifts the curve ζ + upwards by about the order of magnitude O(γ 2 ). It thus implies that the critical values of σ will increase because in general ζ + decreases with σ, as seen in Fig 6. In addition, a higher bath temperature results in a larger correction, and in turn causes σ c to be even greater. Therefore thermal fluctuations from the baths make entanglement harder to maintain. The higher the bath temperature, the more severely the entanglement will deteriorate. This is consistent with our expectation. However, we may be concerned with a possible loophole related to what we found earlier that the finite temperature correction to ζ + may change sign with increasing σ. If it occurred before σ c , we may encounter the opposite conclusion that the lower bath temperature will instead do more harm to the quantum entanglement in the system. We will argue that this is not the case. Let ζ (β) + be the low-temperature correction, so that ζ + = ζ From (6.17) we can derive a relation among the critical values of γ, σ and β for the small γ cases. Similar to (6.15), we can show in the low temperature limit that  Figure 7. The trend of η < with respect to the inter-oscillator coupling σ and the inverse temperature β when both private baths have the same temperature β −1 . It can be translated into the logarithmic negativity by E N (ρ) = max 0, − ln 2η < . We also draw a reference line η < = 1/2, the region below which represents the entangled final state of the joint system. In addition, all these curves are monotonic with respect to the parameters in discussion. The oscillating frequency ω and the cutoff Λ are chosen to be 5 and 10000, respectively. The damping constant γ is 0.5.
Note that the expression in the second pair of square brackets is negative in the low temperature case. We have argued earlier that at low temperature we don't need a strong inter-oscillator coupling to safeguard quantum entanglement, so the curve ζ + can vanish for the small values of σ. Furthermore if Λ ω, we have ln Λ ω (ln Λ ω − 1). Thus we may safely conclude that for small σ and Λ ω, and that the second pair of square brackets in (6.21) is negative in the low temperature case. Alternatively we may roughly see this based on the arguments that for the expansion to be valid we need σ/ω 2 < 1 so the second term should be smaller in magnitude than the first term in the second pair of square brackets. The physical implication of (6.21) is that the critical temperature is lowered when the oscillator-bath interaction gets stronger and vice versa.
Presently we have seen that the low-temperature correction can change sign for sufficiently large mutual coupling; however, this does not affect its usefulness to identify of the existence of entanglement. This unwelcome feature only makes murky the quantitative description of entanglement based on the separability criterion, and it can be traced back to the fact that the separability criterion contains not only η < but also η > , whose existence distorts the information about entanglement, delivered by η < . By comparison negativity is freed from this nonintuitive, unphysical behavior. Let us analyze the finite temperature correction of η < . In the same configuration, it takes a much simpler form with ω 2 ± = ω 2 ±σ. We immediately see that it is always positive and monotonically increasing for all permissible values of σ. Moreover, the finite temperature correction of η < is a monotonically decreasing function of β. This, unlike the separability criterion, give a plausible and physically intuitive description of the extent the state is entangled. Furthermore, since the analytical expression of η < is much simpler than that of the separability criterion, it simplifies the analysis on the critical parameters. Expand out σ c = σ (0) c + γ σ (1) c + · · · and plug this expression back into the symplectic eigenvalue η < = η (0) c satisfies η (0) < = 1/2 at zero temperature. We find which is positive-definite for all permissible ranges of the coupling constant σ and temperature β −1 . Note that the dependence on the cutoff scale is hidden in the expression of σ c . Finally we calculate the critical temperature via the criterion η < = 1/2 in the low temperature regime. Expanding η < with respect to large β gives where f (z) is defined in (4.5). Solving η < = 1/2 leads to .

(6.25)
The inverse critical temperature β c grows with increasing γ but falls off with increasing σ. Therefore, we can see that in the low temperature regime the critical temperature β c at which η < = 1/2 is higher for stronger inter-oscillator interaction, and for weaker oscillatorbath coupling γ. This is totally in line with our intuition that the temperature and the oscillator-bath coupling γ will corroborate to disrupt the quantum coherence between the oscillators and make them harder to remain entangled, while the inter-oscillator coupling will tend to enhance the coordination of both oscillators so their entanglement become more robust. respect to β 1 and β 2 . The surface η < will mildly rise up in the vicinity of β 1,2 ω = O(1) when we approach from the low temperature end. Next we see it sharply climbs up, crossing the dividing curve η < = 1/2 in the region β 1,2 ω = O(1), and enters the high temperature regime. Thus we can make a first observation that, roughly speaking, the dividing curve of η < = 1/2 follows β 1 ω = O(1) and then turns to β 2 ω = O(1). Secondly it implies that separability is determined by the temperature of the warmer bath, instead of the temperature difference, as was mentioned in the previous section. Thirdly, since from earlier discussion we know entanglement tends to survive at higher temperature if the mutual coupling between oscillators is stronger, we use the high temperature approximation to find the critical temperature in the strong σ regime. As shown in Fig 8, we compare the numerical calculation of η < with its low and high temperature approximations, and see that in the large σ case the high-temperature approximation yields a very consistent behavior of η < in the vicinity of βω ∼ O(1), in comparison with the numerical results.
In the high temperature approximation, the symplectic eigenvalue η < is given by The cutoff-dependent factor in those higher order expressions is less important in the weak γ limit because it always appears with the small parameter γ/ω. The critical temperature occurs at η < = 1 2 . Directly solving a transcendental equation like (6.26) for β c is not possible. Nonetheless since ln βΛ always pairs up with γ, we can use the iteration method to derive β c . If we first ignore terms of the order O(γ) and higher, we find β c is given by  Figure 9. We plot the separability criterion ζ + and the symplectic eigenvalue η < together. They crisscross at the critical temperature, and therefor give the same information about the existence of entanglement. However, separability criterion falls off and rises up with increasing β. Note that we shift the values of ζ < downward by 1/2, That is, what we plot in fact is η < − 1/2. The oscillating frequency ω, the cutoff Λ, the damping constant γ and the inter-oscillator coupling σ are chosen to be 5, 10000, 0.2 and 24 respectively.
Substituting it back to seek a correction of order O(γ) we obtain It is indeed consistent with the statement that β c ω = O(1), and it rules out the possibility of the existence of entanglement in the regime βω 1. Again it reveals the fact that with larger inter-oscillator coupling σ we see a higher critical temperature; on the other hand, stronger oscillator-bath interaction γ will cause the critical temperature to decrease.
The same results can be found if we investigate the high temperature approximation of the separability criterion. This is no surprise since we have previously discussed that the separability criterion is perfectly valid for identification of entanglement except for a quantitative measure of entanglement. For example, as shown in Fig. 9, the separability criterion ζ + and the symplectic eigenvalue η < give the same prediction about the location of the critical temperature, but the separability criterion is not a monotonic function of the temperature, which makes it inappropriate as an entanglement measure.
With temperature measured in ratio to the oscillator's natural frequency βω we can conclude that there is no high temperature entanglement in Case C1, namely, between two oscillators each interacting with its own bath.

Intuitive Understanding of Entanglement Behavior
So far we have taken quite some labor to show that asymptotic entanglement between oscillators are easier to sustain for stronger inter-oscillator coupling but weaker oscillatorbath interaction.
Here we would like to offer a physically more transparent illustration as to the competing roles between these two kinds of interactions in terms of normal modes of the oscillator. The Langevin equations (2.1) and (2.2) can be easily decoupled by forming two normal modes and the corresponding equations of motion arë 2) Since we are interested in the late-time dynamics, we will not write down the homogeneous solutions to the Langevin equations. Following the earlier discussions we find that the inhomogeneous solutions are given by The frequency Ω ± is the resonance frequency of the normal modes χ ± . Hence the stronger inter-oscillator coupling σ implies smaller values of ω − but larger values of ω + , and in turn smaller Ω − and larger Ω + . Since the amplitude of the normal modes χ ± is related to the ratio γ/Ω ± , stronger inter-oscillator interaction will induce a larger amplitude of the mode χ − , which will grow with increasing values of σ, meanwhile it causes the mode χ + to oscillate subdominantly and its amplitude decreases with σ. This is intuitively plausible since, e.g., for a very soft spring, or for a particle in a very shallow harmonic potential, a small disturbance could easily induce a large displacement in its motion. Thus in these circumstances it tends to have a large position uncertainty. Furthermore, the consequence from the normal-mode dynamics hints at the fact that when we form the displacements of two oscillators by superposing the normal modes the mode χ + can be overshadowed by χ − . The original displacements χ 1 , χ 2 are more or less determined solely by the mode χ − , in particular in the strong mutual coupling limit Ω − → 0. Furthermore, in this limit, χ 1 and χ 2 will be out of phase by π. Likewise, following similar arguments and taking care of contributions from the resonance, we see that in the case of the conjugate momentum p 1 , p 2 of the two oscillators, the contribution of p + can dominate over that of p − for strong mutual coupling between the oscillators. Furthermore, the phase difference between χ 1 and χ 2 is reflected by the fact that in this limit we should have V 11 ∼ −V 13 . It is particularly easy to see this for the special case β 1 = β 2 . The formal late-time expressions of V 11 and V 13 in this case are In addition Eqs. (5.1) and (5.3) also explicitly demonstrate the same relation. Similarly p + is the dominant mode in the conjugate momenta p 1 , p 2 , so we may expect V 22 ∼ +V 24 and this is clear from for the case β 1 = β 2 . Now let us take a look at the formal expression of the symplectic eigenvalue η < . From (6.10), we see due to the fact that there is no cross-correlation between two normal modes in the case β 1 = β 2 . It is clearly seen that η < is composed of subdominant contributions only, which all have smaller uncertainty with larger σ. Moreover they decrease with increasing values of the inter-oscillator coupling σ. This makes possible that the symplectic eigenvalue η < can fall off with strong inter-oscillator coupling, thus that entanglement can be sustained at higher temperature.
In summary we show that by analyzing the behaviors of the normal mode frequencies with respect to various couplings and parameters of interest, we may get an intuitive understanding of the general features of the entanglement between two oscillators in relation to these parameters.

Summary of Results for Entanglement in Systems in NESS
Having shown the quantitative details in the last section we now provide a summary of the qualitative features of entanglement dynamics in the case (Case C1) studied here for quantum systems in NESS. For two bilinearly coupled oscillators each interacting with its own bath, we find: 1. Quantum entanglement in systems of this setup is harder to sustain at finite temperatures. Thermal fluctuations from the baths disrrupt the coherence between the oscillators.
2. Both the separability criterion and the negativity are perfectly good indicators for the existence of entanglement. However, the former is not necessarily a monotonic function of the parameters in the configuration, so it does not qualify as an entanglement measure. It cannot give a consistent, quantitative comparison between different entangled configurations.
3. The entanglement criterion ζ + or the symplectic eigenvalue η < in general is not a function of the bath temperature difference, in contrast to thermal transport in the same setting [1].
• Lowering the temperature of one of the thermal baths does not necessarily help to keep the entanglement between the oscillators.
• The notion of a critical temperature, where ζ + = 0 or η < = 1/2, is better defined when two private baths have the same temperature.
4. Entanglement between the two oscillators is reduced for stronger oscillator-bath interaction, but enhanced for larger inter-oscillator coupling. They play competing roles as far as their influence on entanglement is concerned.
• strong inter-oscillator coupling better links the dynamics of the two oscillators, and thus improves the coherence between them.
• uncorrelated environment fluctuations corrupts the correlations between the oscillators; stronger oscillator-bath interaction will compound this effect.
• For strong oscillator-bath coupling the critical temperature depends on the damping constant γ and the environment cutoff frequency Λ.
• The effect of environment cutoff cannot be ignored in the low temperature and the strong oscillator-bath coupling cases.
6. Asymptotic quantum entanglement disappears in the high temperature regime βω 1. There is no hot entanglement in systems (with bilinear constant coupling) under NESS conditions.

Comparison: System in a Private Bath vs in a Common Bath
It is useful to make a comparison of the case studied here (Case C1) with what we have studied in Paper I (Case B), namely, a system of two bilinearly coupled oscillators interacting with one common bath.
Case B: common bath. In Paper I we have studied the case of two coupled oscillators at a finite spatial separation, both interacting with a common thermal field bath, which is a finite temperature generalization of the work [9,30]. For comparison with Case C1 studied here we only need to consider the limiting case when the two oscillators are placed next to each other 2 , thus ignoring the spatial variation of entanglement. The action for this setup takes the form Since the two oscillators share the same bath, we can decompose the two harmonic oscillator variables into the fast mode (or center of mass) and the slow mode (or relative coordinate) variables, χ + = 1 2 (χ 1 + χ 2 ), χ − = χ 1 − χ 2 whence the action becomes where ω 2 ± = ω 2 ± σ. We see that the slow mode χ − is decoupled from the bath, while the fast mode χ + now interacts with an effective bath, described by the same scalar field but with the reduced amplitude, φ/ √ 2 and with an effective coupling strength enhanced to √ 2 e. The Langevin equations for the fast and slow variables arë Note there is a subtle issue in this pair of Langevin equations. Although we still write the oscillating frequencies of the fast and the slow mode as ω ± , now they have quite different physical contents. We observe that the fast mode couples with the bath, so its oscillating frequency ω + will acquire a correction due to its interaction with the bath. This correction is absent for the slow mode. Nonetheless because the correction depends on the environment cutoff and it is of the order δω 2 /ω 2 ∼ O(γΛ/ω 2 ), we expect that in the weak oscillator-bath coupling limit, this correction is moderate. In fact it has been shown in [9] that the oscillator-bath coupling constant should be reasonably small or else in may induce instability due to the Coulomb-like interaction at short distances. The stochastic force ξ in this case still possesses the statistical properties For the last four elements, the term caused by the environment vanishes in the limit t → ∞.
The parameter λ is used as a marker for the intrinsic components, and can be set to unity with no consequence. This is in contrast to the components induced by the environment which have e 2 dependence. One feature in the common bath case noticeably different from the private bath case is that the elements of the covariance matrix still contain the information about the initial conditions even though the system has evolved to late time. This is a consequence of the fact that one of the normal modes is completely decoupled from the bath such that part of the initial information is retained in the system and is not lost into the environment. On the contrary, for the private bath case, both the slow and the fast modes are coupled to the environment, and it causes the dispersion of the initial information into the environment. We also note that the components induced by the environment are typically smaller by an order O(γ) than the components intrinsic to the quantum evolution of the oscillators.
Moreover, we have shown that in the private bath case, stronger inter-oscillator coupling renders the oscillating frequency of the slow mode smaller than that of the fast mode. It implies that when the interaction between the oscillators are comparable with the original oscillating frequency ω the slow mode will dominate over the fast mode. From (7.7), we see that the late-time correlation between χ 1 and χ 2 is prone to be negative, meaning that χ 1 tends to be anti-correlated with χ 2 . This is not the case for the common bath case. When both oscillators share a common bath, only the fast mode couples with the bath. The coupling between the oscillators plays a minor role because W + is always of the order ω. At late time t γ −1 , we see that both χ 1 and χ 2 are more or less led by the fast mode, apart from the intrinsic quantum evolution of the system inherited in the slow mode. Hence the bath tends to drive two neighboring oscillators into correlation, meaning that the correlation between χ 1 and χ 2 induced by the shared bath tends to be positive.
At this point, some discretion is advised. First, we observe from (8.14) that the correlation caused by the intrinsic quantum evolution of the system carries a negative sign, and will counteract with the correlation induced by the environment, so the total correlation between χ 1 and χ 2 may not always be positive definite at late time. This will also introduce additional complexity in the entanglement for the shared bath case. Secondly, unlike the private bath case where the elements of the covariance matrix approach a time-independent constant when the NESS is reached, the corresponding elements in the shared bath case remain oscillatory in time reflective of the intrinsic quantum dynamics of the system. We end with a summary of the qualitative behavior of entanglement in a system of two coupled oscillators comparing between the two cases: in one case this system interacts with a common bath, and in the other case, with their own private baths, but kept at the same temperature.
1. From the structure of the normal modes, we find that • private bath: both degrees of freedom are coupled to the bath, so they behave like a pair of damped driven oscillators with different oscillating frequency.
• common bath: only the fast mode is coupled to the bath, the slow mode is totally decoupled from the bath and it acts like a free oscillator.
2. If we separate the elements of the covariance matrix into the intrinsic and the induced components; the former depends on the initial conditions of the oscillators but the latter is entirely driven by the environment, independent of the initial conditions of the oscillators. We see that • for the mode coupled to the bath, the intrinsic component will be exponentially small at late time, and only the induced component survives.
• for the free mode, the intrinsic component oscillates all the time and there is no induced component.
3. This implies that the initial Gaussian conditions will be washed out for the mode coupled to the bath, but they will survive at late times for the uncoupled mode.
• private bath: the initial Gaussian conditions will be irrelevant to the asymptotic entanglement, • common bath: they remain significant, so the final state of entanglement depends on the choice of the initial conditions. 4. At late times the entanglement measure for the private bath case is time-independent, but for the common bath it continues oscillating in time.
5. The amplitude of the driven mode is related to the mode frequency. The smaller the frequency is, the larger the driven amplitude will be.
• private bath: slow mode will have larger driven amplitude than that of the fast mode, so the dynamics of the original canonical variables, which are the superposition of these two modes, will be dominated by the slow mode, especially when the mutual interaction is strong.
• common bath: since there is one driven mode and it is the fast mode, the driven amplitude does not change too much as the mutual coupling strength varies. However, the asymptotic dynamics is determined by the relative magnitude between the slow mode (intrinsic component only) and the fast mode (induced component only).
if the fast mode dominates, then the asymptotic elements of the covariance matrix will be more or less constant in time with small ripples.
if the slow mode dominates, then they will oscillate in time.
6. The inter-oscillator coupling (σ > 0) plays a more important role in the private bath case, but a minor role in the shared bath case.
7. In the private bath case, entanglement is easier to survive for stronger inter-oscillator and weak oscillator-bath coupling, but in the shared bath case, both factors can be overshadowed by the intrinsic components, which are sensitive to the initial conditions of the oscillators.
8. The asymptotic entanglement criterion in the common bath case can thus be broken into three components: one involving the fast mode only, one with slow mode only and the cross term.
• if the fast-mode part is subdominant, then the resulting entanglement criterion will oscillate with time, and that can cause sudden death [28] and revival [29][30][31] (SDR).
• if the fast-mode part is dominant, then there is no SDR phase.

From (A.1) we need the integral
to evaluate V 13 at high temperature. Thus we have the high temperature limit of V 13 given by 3) It means that in this configuration, χ 1 and χ 2 anti-correlated and this anti-correlation grows with the mutual coupling strength σ.

C.1 V 13
We first evaluate the integral We thus obtain V 13 The elements V 14 vanishes because the contributions from both thermal baths cancel.
2. γ → 0: the leading contribution of the momentum uncertainty at late time in this case is lim

C.4 V 24
We need the integral Thus we have the vacuum component of V 24 given by We first evaluate the integral This implies that the low temperature correction to V 13 is given by For V 14 we need the integral Here the finite temperature correction behaves like β −4 , so we will acquire a factor  Before evaluating V 22 , we first evaluate the following two integrals We have the low temperature correction to V 22 given by The corresponding finite temperature correction to the kinetic energy of Oscillator 1 is (D.9)

D.4 V 24
Here we need the integral