Quantum entanglement at high temperatures? 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 (NESS) 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 T1> T2. For 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 Tc, defined as the temperature above which quantum entanglement vanishes. With the Langevin equations derived we give a full display of how entanglement dynamics in this system depends on T1, T2, the inter-oscillator coupling and the system-bath coupling strengths. For weak oscillator-bath coupling the critical temperature Tc 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.


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, as discussed in e.g., [5,6] B) in a nonequilibrium condition and evolving, possibly to an JHEP11(2015)090 equilibrium state 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. 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 finite temperature quantum entanglement under such conditions. We illustrate these two conditions with two generic models: Whereas Case B) is exemplified by a quantum system made of at least two coupled harmonic oscillators (HO) interacting with one common thermal bath as described above, 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 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 (the opposite situation can more easily lead to 'sudden death', as has been shown for two qubits without direct coupling each interacting with its own bath [7].). 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 timedependent inter-oscillator coupling. Before one can bring these cases under the same roof of nonequilibrium steady state condition 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 1 of NESS only for Case C1 in [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.

JHEP11(2015)090
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 a 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 [8][9][10][11][12]. 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 Ghesquière, Sinayskiy & Petruccione [13].

Comments on prior results on the same problem
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-rotating-wave, or NRW) [14]. 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. condition. Although it works better for strong coupling to the environment the results obtained under these approximations have unphysical behavior at low temperatures. In addition, Ludwig et al. [10] 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., in Point a) above, it would be more informative if situations can be explained explicitly whether entanglement can be generated and sustained at sufficiently low temperature, 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 inter-oscillator 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 a finite 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 wellexplored 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 (PHS) separability criterion [15][16][17] and negativity [18,19] 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,20].
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 JHEP11(2015)090 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 cutoff frequency corresponds to the resolution time limit of the detector. The ultraviolet divergence arises from making the unrealistic assumption that all high frequencies are excited. Generically speaking, for weak system-bath coupling these cutoffdependent terms are subdominant and can be ignored. However in the strong coupling case, they may not be negligibly small. 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 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 [21].
First we summarize the difference in results with a special but commonly used configuration of the shared bath case -zero separation between oscillators [22,23], 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.
The results in the (idealized) zero-separation configuration of the shared bath case can, as explained in [24], at best be interpreted as transients, because in fact 1) both normal modes couple with the bath, and will eventually reach relaxation which renders their latetime dynamics independent of the initial conditions, 2) the field-induced interaction effects JHEP11(2015)090 between oscillators cannot be ignored when they are at a finite separation. This implies that additional complication may arise for late-time entanglement. Thus the main difference between Case C1 and Case B with finite separation arises from the former lacking in this non-Markovian contribution.
1. In both cases, late-time entanglement is independent of time, and insensitive to the initial conditions.
2. The behaviors of late-time entanglement in the shared bath case largely depends on the dominance between the direct inter-oscillator coupling and indirect non-Markovian field-induced interaction, whereas in the private bath case the interoscillator coupling is the sole factor that determines the structure of late-time entanglement.
This paper is organized as follows: In section 2, we briefly discuss the dynamics of the reduced system in the NESS configuration, and introduce the covariance matrix, which constitutes the building blocks of the entanglement criterion/measure. In section 3, we address 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 section 4, we then offer a more intuitive viewpoint to understand how all sorts of interactions can affect entanglement between oscillators. In section 5 we summarize our results in the context of entanglement in nonequilibrium quantum systems. A more thorough exposition of the temperature dependence of the covariance matrix elements at high, zero and low temperature cases are provided in the appendix and the supplemental material.
2 The model and the covariance matrix of the system

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 JHEP11(2015)090 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 where G ii H (t − t ) 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.
In the next section we use an example to illustrate how the covariance matrix elements are calculated, leaving the more complete exposition of them in the appendices.

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.13). 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

JHEP11(2015)090
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.18) 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.15) 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.

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 [17] that in the case of continuous Gaussian variables, the Peres-Horodecki-Simon separability criterion [15,16] 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 3)

JHEP11(2015)090
where ρ is the density matrix of the state of the reduced system. We have assumed R = 0. The column matrix R takes the form R T = (χ 1 , p 1 , χ 2 , p 2 ) in our case, 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 (3.1). A negative value of ζ + thus implies the existence of quantum entanglement.
Although (3.1) 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-Simon separability criterion in terms of the symplectic eigenvalues of the partially transposed covariance matrix V pt . Let η ≷ stand for the symplectic eigenvalues of V pt , the partial transposition of V. Without loss of generality we assume η > is greater than η < . 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 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). It does not affect the identification of the existence of entanglement, but it will mess 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 [18], denoted by N or the logarithmic negativity E N [19]. For the Gaussian states under study they can be respectively defined by 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.

JHEP11(2015)090
Comparing the negativity (3.5) with the Peres-Horodecki-Simon criterion (3.4), 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 [25][26][27] 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 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 [27]. For symmetric states, the predictions from both measures coincide.

Late time behavior of the entanglement measures
At late time when the system reaches the steady state, the covariance matrix takes the form 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 section 3.1, 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 JHEP11(2015)090 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 (3.9) as, (3.11) 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 JHEP11(2015)090 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 [28], where alternatively det V can be written as det 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-Horodecki-Simon criterion can be used to identify 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.

JHEP11(2015)090
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. In the next few sections we will compare the predictions of the separability criterion with those of the negativity η < . We will see that as far as the identification of entanglement is concerned, both still yield consistent prognosis at finite temperature, but the separability criterion fails to show reliable quantification of entanglement at higher 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 aymptotically once they are coupled. On the other hand, when we consider contributions due to the finite JHEP11(2015)090 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 (3.18) with ζ + = 0, yielding (3.20) 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 (3.21), 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 figure 3, we let the cutoff JHEP11(2015)090 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 figure 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 effects of temperature and show the inadequacy of the separability criterion at finite temperatures.

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 JHEP11(2015)090 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 figure 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., [29]). 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 JHEP11(2015)090 Figure 5. The symplectic eigenvalue η < is plotted with respect to the temperatures of two private baths. The black curve η < = 0 divides the separable state (η < > 0, pink shade) from the entangled state (η < < 0, green shade). Essentially the curve traces along the region βω = O(1). This result can be easily mapped to the logarithmic negativity by E N (ρ) = max 0, − ln 2η < . The oscillating frequency ω and the cutoff Λ are chosen to be 5 and 10000, respectively. The damping constant γ is 0.1, and the inter-oscillator coupling σ is 20. 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 (3.18). 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 . 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 temperature is higher for stronger σ but weaker γ. The oscillating frequency ω and the cutoff Λ are chosen to be 5 and 10000, respectively where we have used the L'Hôpital's rule to evaluate the limit of such an expression lim x→0 x ln 1 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. 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 (3.21), 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 figure 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 ζ + = ζ  . 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 (3.26) 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 (3.26) is that the critical temperature is lowered when the oscillator-bath interaction gets stronger and vice versa.

JHEP11(2015)090
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) 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 σ (0) 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 (B.5). Solving η < = 1/2 leads to . (3.30) 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. In addition, we have learned that at finite temperature the separability criterion fails to provide a quantifiable description of entanglement. Thus, in the next section, we will inspect the high temperature entanglement based on the negativity exclusively.

High temperature βω 1
We now turn our attention to the high temperature regime and ask if entanglement at high temperatures is at all possible. From the plot of the symplectic eigenvalue η < against the bath temperatures β 1 and β 2 in figure 5 we see that the surface η < forms a very flat basin which is symmetric with 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 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 figure 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 (3.31) for β c is not JHEP11(2015)090 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.
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 2 √ 3/ √ 3ω 2 + 4σ. 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 figure 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 oscillator-JHEP11(2015)090 bath 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ë 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

JHEP11(2015)090
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. (C.1) and (C.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 . Alternatively let us take a look at the formal expression of the symplectic eigenvalue η < . In the case β 1 = β 2 , from (3.15), we see due to the fact that there is no cross-correlation between two normal modes. 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: entanglement in nonequilibrium quantum systems
To obtain a fuller perspective of 'hot entanglement' and entanglement in nonequilibrium quantum systems it is useful to compare the results obtained here for a quantum system (of two coupled oscillators) in nonequilibrium steady state NESS (each oscillator has its own bath -Case C1 ) with the same system in a shared thermal bath (Case B ) studied in [21], which we did in section 1.4, and that which is kept in thermal equilibrium at all times (Case A).
The launchpad for this investigation began with a study of the effects of direct coupling and spatial separation on the late-time entanglement between two oscillators in a JHEP11(2015)090 shared zero-temperature bath [24], extending the earlier work of [29] for early time entanglement between two oscillators without direct coupling. The nonequilibrium dynamics of this system in a shared thermal bath Case B can be found in [21], where the parameter ranges for entanglement to survive at a finite temperature are given. It is shown there that generically when two coupled oscillators separated at a fixed distance evolve under the influence of a shared thermal bath, their dynamics is usually non-Markovian due to the field-induced interactions between the two oscillators (or qubits, e.g., [30,31]). The asymptotic correlations / entanglement between the oscillators depend greatly on the relative strengths between their direct coupling and the indirect field-induced interaction which is non-Markovian in nature. When direct coupling dominates over induced interaction, the late-time entanglement tends to survive better under 1) stronger inter-oscillator coupling, 2) weaker oscillator-bath interaction and at 3) a longer distance between them. On the other hand, in the regime dominated by the field-induced interactions, these factors enter in opposite ways for asymptotic entanglement. In the case of weak oscillator-bath coupling, direct coupling is more essential in maintaining late-time entanglement at higher temperature because the field-induced interaction is heavily suppressed. Thus the critical temperature is still bounded by the inverse of the center-of-mass mode frequency, but tends to be lower for shorter separation due to adverse effects from the field-induced interaction. Late-time entanglement can sustain over much larger separation, on account of direct coupling, than the scale of the order of the inverse cutoff frequency inherited from the thermal bath.
These results are dissimilar to, but not in conflict with, the results for Case A, since in Case A direct coupling is limited to the neighboring sites. Hence the two-site entanglement is expected to fall off rapidly beyond the nearest neighbors. In addition, as we implicitly pointed out earlier that in the regime that direct coupling is negligible and the field-induced interaction governs, late-time entanglement deteriorates swiftly with increasing separation. This is not unexpected since it is known [32,33] that in the weak oscillator-bath coupling and zero direct inter-oscillator coupling limit, both configurations will yield similar results. For stronger oscillator-bath interaction, the field-induced interaction can sustain over a very long history in the evolution of both oscillators. Deviation in results between Case A and Case B will become more apparent. Nonetheless a strong oscillator-bath interaction can likely induce dynamical instability in the oscillators, a case worthy of closer analysis later. Now in the case (Case C1 ) studied here for two bilinearly coupled oscillators each interacting with its own bath, the qualitative features of entanglement dynamics for quantum systems in nonequilibrium steady state (NESS) can be summarized as follows: 1. Quantum entanglement in systems of this setup is harder to sustain at finite temperatures. Thermal fluctuations from the baths disrupt 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 entangle-JHEP11(2015)090 ment 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.
In a future communication we will report on our findings of hot entanglement in Case C2 situations.

JHEP11(2015)090
In this appendix we will outline the approaches we used to evaluate the elements of the covariance matrix at late times for three situations: (1) high temperature limit, (2) zero temperature case and (3) low temperature regime in section A-B. In section C, we summarize the temperature dependence of the covariance matrix elements, and discuss their general features at finite temperature. We leave the detailed evaluations of the late-time expressions of the covariance matrix in the supplemental material.

A 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.18) 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., [10]).
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.15) 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 the supplemental material. Nonetheless for the following discussions we will bring forward the results

JHEP11(2015)090
for V 22 and V 13 here. When both private baths have the same temperature β −1 , we have from (C.1), (C.3) and (C.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. (A.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

JHEP11(2015)090
when β 1 = β = β 2 . Including this contribution we see the total energy for the two-oscillator system in the high-temperature limit becomes 13) 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 (A.13) with (A.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 [34] that both results can be equivalent only in the limit of vanishingly weak oscillator-bath coupling.

B 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.

B.1 V 11 at zero temperature
The vacuum contribution of G kk H (κ) is so we need the following two integrals to evaluate the vacuum contribution of V 11 , 3) The sum of J 1 and J 2 can be expressed as where the dimensionless function f (z) is defined by

JHEP11(2015)090
It is clear that Ω ± are the resonance frequencies of the two normal modes. Therefore from (2.15), 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 the supplemental material.
B.2 V 11 at low temperature βω 1 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.15), 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 (B.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 (B.7) should be proportional to n −2 β −2 . Their overall contributions will introduce an addition factor ∞ n=1 1 n 2 β 2 = π 2 6 1 β 2 , (B.11) to (B.10). Therefore after taking this into account, we obtain the low temperature correction to V 11 as follow: We leave the derivations of the zero and the low temperature expressions for V 13 , V 14 , V 22 , V 24 in the supplemental material.
C Temperature dependence of the covariance matrix Because elements of the covariance matrix at finite temperature may prove useful for more general purposes, we collect their expressions for both high and low temperatures for the system at late times when it reaches a NESS, the existence of which for this setup is shown in our earlier paper [1].
Here we summarize the temperature dependence of the elements of the covariance matrix.
θ(β j Λ − 1) ln β j Λ + m 2 8ω 2 γ 2 + σ 2 4ω 2 γ 2 + σ 2 1 β 1 + σ 2 4ω 2 γ 2 + σ 2 1 β 2 , β 1,2 ω 1 . (C.6) 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.18). 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 < Λ.
In figure 10 we show the comparisons of 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 vacuum contribution, fits nicely with the numerical JHEP11(2015)090 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.
6. V 24 = 1 2 {p 1 (∞), p 2 (∞)} : In this case, since the leading contribution of the high-temperature approximation vanishes, we have to include the next-order term.
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 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.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.