Entanglement Dynamics of Quantum Oscillators Nonlinearly Coupled to Thermal Environments
EEntanglement Dynamics of Quantum OscillatorsNonlinearly Coupled to Thermal Environments
Aurora Voje, Alexander Croy, and Andreas Isacsson ∗ Department of Applied Physics, Chalmers University of Technology, SE-412 96 G¨oteborg, Sweden (Dated: September 21, 2018)We study the asymptotic entanglement of two quantum harmonic oscillators nonlinearly coupledto an environment. Coupling to independent baths and a common bath are investigated. Numericalresults obtained using the Wangsness-Bloch-Redfield method are supplemented by analytical resultsin the rotating wave approximation. The asymptotic negativity as function of temperature, initialsqueezing and coupling strength, is compared to results for systems with linear system-reservoir cou-pling. We find that due to the parity conserving nature of the coupling, the asymptotic entanglementis considerably more robust than for the linearly damped cases. In contrast to linearly damped sys-tems, the asymptotic behavior of entanglement is similar for the two bath configurations in thenonlinearly damped case. This is due to the two-phonon system-bath exchange causing a supressionof information exchange between the oscillators via the bath in the common bath configuration atlow temperatures.
PACS numbers: 03.67.Bg, 03.65.Yz, 85.85.+j
I. INTRODUCTION
Entanglement challenges our comprehension since the1930’s [1, 2], and still remains a highly relevant topic.Problems relating to entanglement creation and manip-ulation are of importance for a broad range of questionsrelated to quantum information science[3], like quan-tum cryptography[4], quantum dense coding[5], quantumcomputation algorithms [6] and quantum state telepor-tation [7–9]. In particular, recent experimental advances[10–13] pave the way for entanglement-based technology.In this paper, asymptotic effects of nonlinear dissi-pation on the entanglement of harmonic oscillators areinvestigated and compared to the widely studied situ-ation of linearly damped (LD) systems. In the latterthe oscillators are linearly coupled to bosonic reservoirs.Such system-reservoir interactions have for instance beeninvestigated within both Markovian [14, 15] and non-Markovian dissipation models [16–19]. For systems ini-tially in squeezed states, high temperature entanglement[20], and the exotic behavior of entanglement sudden dis-appearance and revival (ESDR) [14, 16, 18, 21, 22] havebeen found.Typically dissipation destroys quantum entanglement.However, it is known that by engineering the system-reservoir coupling, entanglement can be generated [23–25]. For example, one possibility to entangle initially sep-arable states is through the introduction of multi-quantadissipation, or nonlinear damping (NLD). Naturally oc-curring NLD has been reported in systems which pos-sess strong intrinsic nonlinearities [26]. Among these arecarbon-based nanomechanical systems like graphene andcarbon nanotubes [25, 27, 28]. Additionally, there havebeen reports of inducing nonlinear dissipation in optome- ∗ Electronic address: [email protected] chanical systems [29, 30], and suggestions for possibleemergence of NLD in solid state quantum devices [31].In Ref. [32] we demonstrated the possibility of entan-glement generation from initially separable states. Here,in the light of previous studies on linearly damped oscilla-tor systems, we investigate how NLD affects the asymp-totic state behavior of initially entangled states. Forthe latter we choose two-mode squeezed vacuum states[13, 14, 33, 34]. These states are entangled, Gaussianstates, which approach the maximally entangled EPR-state by an increase of the squeezing parameter.While NLD is usually accompanied by a conservativenonlinearity, it was found in Ref. [32] that a weak con-servative Duffing nonlinearity did not affect the asymp-totic state behavior when only the lowest lying eigen-states are occupied. Hence, we here limit the study topurely harmonic oscillators, nonlinearly coupled to eitherone common or two individual environments. This ap-proach serves to isolate effect of the nonlinear relaxationbehavior on the entanglement and allows a more trans-parent comparison with the linearly damped systems.Compared to a linearly damped system we find thatthe parity protection inherent to the two-phonon ex-change between the system and the reservoirs, andpresent in the nonlinearly damped systems, changes theasymptotic behavior in several ways. Firstly, the asymp-totic decay of entanglement is considerably slower for twouncoupled oscillators. This is due to the necessity ofsimultaneous excitation processes of the two oscillatorsneeded for thermal dephasing. Secondly, for individualoscillators coupled to a common bath, we do not repro-duce the sharp transition between steady state entangle-ment and disentanglement in the infinite time limit, seenin linearly damped systems. For the linearly dampedsystem, persistent entanglement is connected to the rela-tive oscillator motion degree of freedom being decoupledfrom the bath. For the nonlinearly damped system, nosuch decoupling occurs. Finally, for weakly coupled os- a r X i v : . [ c ond - m a t . m e s - h a ll ] D ec cillators, parity protection in combination with coherentoscillations in the oscillator populations due to the cou-pling leads to disappearance and reappearance of entan-glement reminiscent of ESDR behavior.The organization of this paper is: First, in Sec. II,we present the model Hamiltonians and derive quantummaster equations (QMEs) for two separate harmonic os-cillators coupled to either individual baths or a commonbath. Then, in Sec. III, we present the asymptotic en-tanglement behavior as function of temperature, initialsqueezing and dissipation rate. We compare our resultson nonlinearly damped systems to previous results onlinearly damped systems. In Sec. IV, we comment onsome features of the asymptotic entanglement behaviorof a coupled oscillator system. II. QUANTUM MASTER EQUATIONS FORUNCOUPLED OSCILLATORS
First we consider two different scenarios where a sys-tem of two independent harmonic oscillators with fre-quencies ω are coupled quadratically in position to ei-ther two individual, or one common reservoir of harmonicoscillators. The situation with two weakly coupled oscil-lators is discussed in section IV. For the uncoupled oscil-lators the Hamiltonian is H = H S + H B + H cb , ibSB . Measur-ing length, time and energy in units of (cid:112) (cid:126) / mω , ω − and (cid:126) ω respectively, we have H S = (cid:88) j =1 , (cid:18) p j + 12 ω q j (cid:19) , (1a) H B = (cid:88) j (cid:88) k ω jk b † jk b jk , (1b) H ibSB = (cid:88) j q j (cid:88) k η jk ( b † jk + b jk ) , (1c) H cbSB = (cid:88) j q j (cid:88) k η k ( b † k + b k ) . (1d)Here, p j = i (cid:112) ω / a † j − a j ) and q j = ( a † j + a j ) / (cid:112) (2 ω )denote momentum and oscillation amplitude of oscillator j , respectively, while a j ( a † j ) is the annihilation (creation)operator of the j -th oscillator.The system-bath coupling part of the total Hamilto-nian is denoted by H ibSB for two individual baths, and H cbSB , for the common bath. For the individual bathconfiguration the operator b † jk ( b jk ) creates (destroys) aphonon in state k of reservoir j with the frequency ω jk .The coupling strength of oscillator j to reservoir state k is denoted by η jk . Similarly, for the common bath, theoperator b † k ( b k ) creates (destroys) a phonon in state k ofthe common reservoir with the frequency ω k . The cou-pling strengths of both oscillators to the reservoir state k are denoted by η k .To study the time evolution of the system we numer-ically solve the QMEs for the reduced density matrix ρ in the weak system-reservoir coupling limit. To obtainanalytical results we implement the rotating wave ap-proximation (RWA). Below we summarize the QMEs forboth bath configurations with and without RWA. A. QME for coupling to individual baths
Using the Born-Markov approximation in the interac-tion picture with respect to H S , the general QME for theindividual bath configuration is given by [35] ∂∂t ρ ( t ) = − (cid:80) l,j ∞ (cid:82) dτ (cid:104) S l ( t ) , S j ( t − τ ) ρ ( t ) (cid:105) C lj ( τ ) − (cid:104) S l ( t ) , ρ ( t ) S j ( t − τ ) (cid:105) C jl ( − τ ) . (2)The operators S j ( t ) = e i H S t ( a † j + a j ) e − i H S t and B j ( t ) = (cid:80) k η jk (cid:16) b † jk e i ω jk t + b jk e − i ω jk t (cid:17) allow us to rewrite thecoupling Hamiltonian as H ibSB ( t ) = (cid:80) j =1 , S j ( t ) ⊗ B j ( t )in the interaction picture. Assuming initial thermal equi-librium of the reservoirs, ρ B = ρ B , ⊗ ρ B , , their correla-tion functions C jl ( τ ) = Tr B { B j ( t ) B l ( t − τ ) ρ B } are C jl ( τ ) = δ jl (cid:90) dω π κ j ( ω ) (cid:2) N ( ω ) e i ωτ +( N ( ω ) + 1) e − i ωτ (cid:3) , (3)where N ( ω ) = ( e ω/k B T − − is the Bose-Einstein distri-bution and κ j ( ω ) = 2 π (cid:80) k | η jk | δ ( ω − ω jk ) are the spec-tral densities. The specific form of κ j depends on the mi-croscopic details of the system-reservoir coupling. If κ j issufficiently smooth around the frequencies of interest, theexact frequency dependence is not crucial. To be specific,we use an Ohmic spectral density, κ j ( ω ) = Γ j ω/ (2 ω ),where Γ j is the non-linear dissipation strength of the j -thbath.Further, we define the one-sided Fourier transform ofthe reservoir correlation function12 γ j ( ω ) + i σ j ( ω ) = ∞ (cid:90) dτ e i ωτ C jj ( τ ) . (4)The rates γ j determine the strength of dissipation, while σ j renormalize the system Hamiltonian. For simplicitywe from here on let ω denote the renormalized systemfrequencies and neglect the corresponding small inducedconservative nonlinearity.Using the expression of the bath correlation function(3) one finds γ j (2 ω ) = Γ j [ N (2 ω ) + 1] , (5a) γ j ( − ω ) = Γ j N (2 ω ) . (5b)In the RWA, equation (2) simplifies to˙ ρ = − (cid:88) j =1 , (cid:104) γ j (2 ω ) L [ a † j ] + γ j ( − ω ) L [ a j ] (cid:105) ρ, (6)with L [ X j ] ρ = X j X † j ρ + ρX j X † j − X † j ρX j . (7) B. QME for coupling to a common bath
For the common reservoir configuration the summationin (2) can be omitted and the general form of the commonbath QME is ∂∂t ρ ( t ) = − ∞ (cid:82) dτ (cid:104) S ( t ) , S ( t − τ ) ρ ( t ) (cid:105) C ( τ ) − (cid:104) S ( t ) , ρ ( t ) S ( t − τ ) (cid:105) C ( − τ ) , (8)where the common system and bath operators in (8)are S ( t ) = (cid:80) j =1 , ( a † j e i ω t + a j e − i ω t ) and B ( t ) = (cid:80) k η k (cid:16) b † k e i ω k t + b k e − i ω k t (cid:17) . The correlation reservoirfunction is then given by C ( τ ) = Tr B { B ( t ) B ( t − τ ) ρ B } .In this case the interaction picture RWA QME is˙ ρ = − (cid:88) j =1 , (cid:104) γ j (2 ω ) (cid:0) L [ a † j ] + L [ a † j ] (cid:1) + γ j ( − ω ) (cid:0) L [ a j ] + L [ a j ] (cid:1)(cid:105) ρ (9)with L [ X j ] ρ = X j X † j − ( − j ρ + ρX j X † j − ( − j (10) − X † j − ( − j ρX j . The QME in (9) is similar to (6), but with additionalcross terms by which the two sub-systems are connectedvia the bath.As shown in Ref. [18], coherence and entanglement of asqueezed state is better preserved in a symmetric system.We therefore consider a setup in which the dissipationrates for the two oscillators are set equal in all system-bath configurations, Γ j = Γ . We also define γ (2 ω ) = γ − and γ ( − ω ) = γ . Further, we define basis vectors | n, i (cid:105) = | n (cid:105) ⊗ | i (cid:105) , denoting eigenstates with n quanta inoscillator 1 and i quanta in oscillator 2. III. RESULTS FOR INDEPENDENTOSCILLATORS
In order to compare the asymptotic entanglement ofnonlinearly damped independent oscillators to linearlydamped ones, we solve the QME (2) numerically. We usethe Wangsness-Bloch-Redfield approach in the eigenba-sis of the system Hamiltonian [36–38] in a Hilbert spacetruncated above M = 8 eigenstates for each oscillator.To facilitate the comparison, the system is initializedwith the two-mode squeezed vacuum ρ (0) = ˆ S ( ξ ) | , (cid:105)(cid:104) , | ˆ S † ( ξ ) , (11) S quee z i ng pa r a m e t e r 𝑟 Temperature 𝑘 B 𝑇/𝜔 Temperature 𝑘 B 𝑇/𝜔 S quee z i ng pa r a m e t e r 𝑟 D i s en t ang l e m en t t i m e 𝜏 d i s = Γ 𝑡 d i s FIG. 1: (Color online) Main figure: scaled disentanglementtime τ dis = Γ t dis (colorbar) of nonlinearly damped two-modesqueezed vacuum states as function of temperature T andsqueezing parameter r , for two individual baths with simula-tion time τ sim = 50, damping rate Γ = 10 − ω , and negativ-ity cut-off ε = 10 − . Inset: disentanglement time τ dis (insetcolorbar scale) of linearly damped two-mode squeezed vac-uum states as function of T and r for two individual baths,with simulation parameters as in the main figure. The dashedlines are the contours of τ dis = (cid:2) , , (cid:3) (right to left), as the-oretically predicted in [14]. where the two-mode squeezing operator is ˆ S ( ξ ) = e ξa † a † − ξ ∗ a a and ξ = re i θ . In the Fock basis, using θ = π , equation (11) becomes [39] ρ = 1cosh ( r ) ∞ (cid:88) n,m =0 ( − ( n − m ) [tanh( r )] n + m | n, n (cid:105)(cid:104) m, m | . (12)To quantify the entanglement we use the measure ofnegativity N = ( || ρ T || − /
2, where ρ T denotes thepartial transpose of the bipartite density matrix with re-spect to oscillator one. The negativity corresponds to theabsolute value of the sum of negative eigenvalues of ρ T and vanishes for separable states [40]. A. Asymptotic entanglement for coupling toindividual baths
For the individual bath configuration and a Markovianmodel in the RWA, it was shown in Ref. [14] that for fi-nite temperatures (
T >
0) all linearly damped two-modesqueezed vacuum states disentangle within a finite time,and relax to the ground state. At T = 0 the relaxationto the ground state occurs in the limit of infinite time.These results were further probed with non-Markovian Scaled time 𝜏 = Γ 𝑡
0 50 150 100 10 -25 -20 -15 -10 -5 -0 N e g a t i v i t y 𝒩 𝛾 𝒩 = 𝜕 𝜏 l o g 𝒩 𝑘 B 𝑇/𝜔
0 2 1 3 0 1 2 (a) 𝑇 𝑇 𝑇 𝑇 𝑇 𝑇 𝑇 S l ope 𝛾 𝒩 = | 𝜕 𝜏 l o g 𝒩 | Dissipation rate Γ /𝜔 𝑘 𝐵 𝑇/𝜔 = 0 𝑘 𝐵 𝑇/𝜔 = 3 𝑘 𝐵 𝑇/𝜔 = 2.5 𝑘 𝐵 𝑇/𝜔 = 1.5 𝑘 𝐵 𝑇/𝜔 = 0.5 𝑘 𝐵 𝑇/𝜔 = 2 𝑘 𝐵 𝑇/𝜔 = 1 (b) FIG. 2: (a)
Main figure: Time evolution of the negativity of a nonlinearly damped squeezed two-mode vacuum withsqueezing parameter r = 1 /
20, individual bath configuration, for temperatures 0 ≤ k B T /ω ≤ = 10 − ω [1 / , / , , , T and different Γ overlap when plotted as function of τ = Γ t .Inset: Slope of the negativity, γ N = | ∂ τ log N | , extracted from the second half of the points in the main figure, as functionof T for Γ = 10 − ω [1 / , / , , ,
10] (color code). The dashed line is a slope fit of 2 N (2 ω ). (b) Slope of the negativity γ N = | ∂ τ log N | from the inset in (a) as function of 10 Γ /ω . models in Refs. [16, 17], lending support to the predic-tions of Ref. [14]. Thus, concluding that the Markovianand non-Markovian dynamics coincide for times largerthan reservoir correlation times.For NLD, earlier studies of a single oscillator systemsundergoing NLD show that parity conservation bringsthe system to a final non-classical steady state [23, 25].For a bipartite system with no inter-mode coupling, itfollows that the same parity conservation will, at T = 0,bring the system into a general steady state ρ ( ∞ ) = P | (cid:105)(cid:104) | + P | (cid:105)(cid:104) | + ρ , | (cid:105)(cid:104) | + ρ , | (cid:105)(cid:104) | + (cid:104) ρ , | (cid:105)(cid:104) | + ρ , | (cid:105)(cid:104) | + ρ , | (cid:105)(cid:104) | + ρ , | (cid:105)(cid:104) | + ρ , | (cid:105)(cid:104) | + ρ , | (cid:105)(cid:104) | + H . c . (cid:105) , (13)with matrix elements ρ ni,mj determined by the initialstate. The particular initial state (12) leads to a steadystate of the form (13) where several elements are zero,reducing it to ρ ( ∞ ) = P | (cid:105)(cid:104) | + P | (cid:105)(cid:104) | + ( ρ , | (cid:105)(cid:104) | + H . c) . (14)The element ρ , is important for the asymptotic neg-ativity (entanglement). While initially there are multi-ple off-diagonal elements contributing to the negativity N ( t = 0) = ( e | r | − / N ( ∞ ) = | ρ , ( ∞ ) | T =0 . This can be verified through the characteristic equationfor ρ T of the steady state (14). For a general M × M basis size the characteristic equation is given by( − µ ) M − ( P − µ )( P − µ )( µ − | ρ , | ) = 0 , (15)with only one negative root µ = −| ρ , | .Comparing the nonlinear and linear decays of squeezedstates at T = 0 one finds that the nonlinearly dampedstates remain entangled with a saturating negativity,whereas the linearly damped states asymptotically dis-entangle in the limit of t → ∞ . This can be seen inFig. 1 showing the scaled disentanglement time τ dis (col-orbar) of a nonlinearly damped (main figure) and linearlydamped (inset) two-mode squeezed vacuum as functionof T and r . The results are obtained from numerical sim-ulations. The scaled disentanglement time is defined as τ dis = Γ t dis , where t dis is the time at which N < ε ,and ε is the negativity cut-off. Temperature is mea-sured in units of (cid:126) ω /k B . The main figure and insethave identical simulation parameters, but different disen-tanglement time scales (main and inset colorbars). Thewhite, dashed lines in the inset are the theoretically pre-dicted disentanglement times derived in Ref. [14].The nonlinearly and the linearly damped systems for T > r . The main difference is the disentanglement time scale.The nonlinearly damped states disentangle much slowerthan the linearly damped states. During the chosen evo-lution time all LD states disentangle, while a part of theNLD states remain entangled (white region, main figure).After a longer time evolution all NLD states will eventu-ally disentangle.The time evolution of the negativity while relaxing tothe steady state is shown in Fig. 2 (a). Here we use aconstant value of r = 1 /
20 for several values of T and Γ .The graphs for the same T and different Γ overlap whenexpressed in terms of scaled time units τ = Γ t . Initiallythe negativity has a rapid initial transient during whichthe initial squeezed state reduces to a state of the formgiven by Eq. (14). This is followed by a slow exponentialdecay.To quantify the asymptotic decay, the inset in Fig. 2(a)shows the negativity slope γ N = | ∂ τ log N | , extractedfrom the second half of the data in the main panel ofFig. 2(a), as function of T . The decay rate solely dependson T , which is further corroborated in Fig. 2 (b), showingthe negativity slope γ N from (a) as function of 10 Γ . Asshown in appendix A, in the limit of low temperatures,within the RWA, the slope is given by the expression N ( τ ) = ( | ρ , ( ∞ ) | T =0 ) e − N (2 ω ) τ , shown as the dashedline in the inset to Fig. 2(a).The much slower decay of the negativity in the nonlin-early damped case compared with the linearly dampedsystem can be understood as follows: For finite temper-atures the disentanglement of the nonlinearly dampedsqueezed vacuum states is related to the slow, thermaldephasing of the parity protected matrix element ρ , .However, thermal decoherence of this element requires asimultaneous two-quanta excitation of both oscillators,a process which is less probable compared to linear de-coherence, where neither an individual nor simultaneousexcitation is needed to achieve a de-excitation. B. Asymptotic entanglement for coupling to acommon bath
We now turn to the situation with the two independentoscillators nonlinearly coupled to a common reservoir.For a linearly damped system in this configuration, theasymptotic steady states of an initial two mode squeezedstate (11) can be divided into entangled and separablestates [14, 18, 19]. To which category a state will belong,depends on r and T . This result was obtained in Ref.[14] using Markovian dynamics and RWA and is shownin Fig. 3(inset). An interesting feature is that the systemnever disentangles for T = 0. This entanglement preser-vation can be explained in terms of normal modes, e.g.center of mass and relative coordinates. As only the cen-ter of mass motion is affected by the dissipation to thebath, the relative motion of the oscillators evolves freely.For nonlinear system reservoir coupling, there is no de-coupling of the relative oscillator motion from the bath,and we do not find any finite temperature steady stateentanglement.As seen in the main figure 3, showing the scaled disen-tanglement time (colorbar) of nonlinearly damped two-mode squeezed vacuum states as function of T and r for the common bath configuration, the asymptotic en- Temperature 𝑘 B 𝑇/𝜔 S quee z i ng pa r a m e t e r 𝑟 D i s en t ang l e m en t t i m e 𝜏 d i s = Γ 𝑡 d i s Temperature 𝑘 B 𝑇/𝜔 S quee z i ng pa r a m e t e r 𝑟 FIG. 3: Main figure: disentanglement time τ dis = Γ t dis (colorbar) of nonlinearly damped two-mode squeezed vacuumstates as function of temperature T and squeezing parameter r , common bath, with simulation time τ sim = 50, dampingrate Γ = 10 − ω , and the negativity cut-off ε = 10 − . Inset:Entanglement borderline of the linearly damped two-modesqueezed vacuum states, common bath, in the phase space of r and T . The dashed line is the theoretical prediction in [14]and the dots with errorbars are numerical data. The simula-tion parameters are: τ sim = 50, Γ = 10 − ω and ε = 10 − . tanglement behavior is very similar to the nonlinearlydamped squeezed states in the individual bath configu-ration. Again, the nonlinearly damped states disentangleslower than the linearly damped states in the disentan-gled region in the inset of Fig. 3. Like in the main panelof Fig. 1, not all states have yet disentangled for the cho-sen simulation time (white region), but will do so after alonger evolution.The similarity in the behaviors stem from a supressionof information exhange between the oscillators via thebath in the steady state ρ ( ∞ ) for T ≈
0. In the QME(9) the term L ρ quickly brings the initial state (11) tothe steady state (14), which cannot be further affectedby the term L ρ , as L ρ ( ∞ ) = 0. For higher tempera-tures, the term L ρ will contribute to some informationexchange, but not enough to significantly alter the in-fluence of the L ρ -term. The qualitative evolution of thestate is therefore as for the individual bath configuration.This is supported by the results in Fig. 4, displaying theslow temperature dependent negativity decay for various T and Γ (main figure), and the temperature dependentnegativity decay-rates γ N = | ∂ τ log N | = 2 N (2 ω ) (in-set). The derivation of the exponent can be found inappendix A and is equal to the thermal decay exponentobtained for individual baths. Scaled time 𝜏 = Γ 𝑡
0 50 150 100 10 -25 -20 -15 -10 -5 -0 N e g a t i v i t y 𝒩 𝑘 B 𝑇/𝜔
0 2 1 3 𝛾 𝒩 = 𝜕 𝜏 l o g 𝒩 FIG. 4: (Color online) Time evolution of the negativityof a nonlinearly damped squeezed two-mode vacuum withsqueezing parameter r = 1 /
20, common bath configura-tion, for temperatures 0 ≤ k B T /ω ≤ = 10 − ω [1 / , / , , , γ N extracted from the second half of the points in the mainfigure, as function of T for Γ = 10 − ω [1 / , / , , , N (2 ω ). IV. RESULTS FOR COUPLED OSCILLATORS
Finally, we consider two weakly coupled oscillators.Based on the results of the preceding section, the mainqualitative difference between LD and NLD stems fromthe parity conservation for the individual oscillators.When the oscillators are coupled only the parity of theentire system is conserved. In particular, the element ρ , which we found to give the asymptotic negativity,is no longer protected and can decay via an intermediatetransition to the state ρ , .Since the situation is close to the linear case, we restrictthe discussion to individual baths. The system Hamilto-nian (1a) is adjusted to include an inter-mode couplingterm H S = (cid:88) j =1 , (cid:18) p j + 12 ω j q j (cid:19) + √ ω ω λq q . (16)The corresponding Hamiltonian in the RWA is H S , RWA = (cid:88) j =1 , ω j a † j a j + λ a † j a j − ( − j . (17)By the same procedure as described in Sec. II, the RWA-QME for a symmetric system, ω j = ω , in the weak intermode coupling limit, λ (cid:28) ω , obtains the form [32]˙ ρ = − (cid:88) j =1 , (cid:104) γ j (2 ω ) L [ a † j ] + (18) γ j ( − ω ) L [ a j ] (cid:105) ρ + D ( λ ) ρ . To the lowest order in λ the oscillators are individuallycoupled to their respective reservoirs, and the superop-erator D ( λ ) becomes D ( λ ) ρ = Υ + L [( n − n )] ρ −
12 Υ − (cid:104) ( n − n )( a † a − a † a ) ρ − ( a † a − a † a ) ρ ( n − n )+ ρ ( a † a − a † a ) † ( n − n ) − ( n − n ) ρ ( a † a − a † a ) † (cid:105) . (19)Here Υ ± = γ ( λ ) ± γ ( − λ ) with γ ( λ ) = κ ( λ )[ N ( λ )+ 1] and γ ( − λ ) = κ ( λ ) N ( λ ).The coupling λ plays a dual role of contributing to os-cillator interaction via H S and to decoherence via D .For linearly damped oscillators the decoherence terms in D would only arise if Γ (cid:54) = Γ [41]. There is no suchrestriction however for the nonlinearly damped systemwhere the superoperator D consists of two decoherenceterms, proportional to Υ + and Υ − respectively. The tem-perature dependencies of these terms differ from thoseof γ j . For temperatures exceeding the coupling energy, λ (cid:28) k B T , we have Υ + > Υ − . In the low temperaturelimit, λ (cid:29) k B T , both terms approach equal magnitudesΥ + ≈ Υ − .For weakly coupled oscillators, linearly coupled to in-dividual baths, a phase diagram separating the entan-gled steady states from non-entangled steady states ex-ists (see for instance Ref. [20]). Whether or not the stateremains entangled depends on temperature and strengthof coupling. A more in depth study, using both Marko-vian as well as non-Markovian evolution can be foundin Ref. [17]. As explained in Ref. [17], starting from aninitial two-mode squeezed state, an undamped coupledoscillator system will display coherent oscillations duringits time evolution with corresponding oscillations in neg-ativity. Adding finite linear damping, results in loss ofentanglement in the long time limit and suppression ofcoherent oscillations in the negativity.For the case of nonlinear coupling to individual bathswe find results which are similar to those in Ref. [17], withcoherent oscillations reflected in the decay curve shownin Fig. 5. As can be seen, with increasing temperaturethe oscillations vanish but the initial rapid decay remainsunaltered. The coherent oscillations can be traced backto the time evolution for T = 0 and λ >
0. In this case,after the transient rapid decay, the negativity dynamicsis governed by the evolution of the ρ , matrix element.Since the two-mode squeezed vacuum only has even en-tries, it suffices to analyze the matrix elements with total
0 2 4 6 8 10 12
Scaled time 𝜏 = Γ 𝑡 N e g a t i v i t y 𝒩 𝑘 B 𝑇/𝜔 = 0 𝑘 B 𝑇/𝜔 = 0.1 𝑘 B 𝑇/𝜔 = 0.2 𝑘 B 𝑇/𝜔 = 0.3 𝜌 , 𝜆𝑡 𝜋 𝜋 𝜋 −4 × 10 FIG. 5: Decay of the negativity N as function of scaled time τ = Γ t for an initial state with squeezing parameter r = 0 . λ = Γ = 10 − ω . The cou-pling is reflected in the coherent oscillations superposed onthe exponential decay towards the thermal equilibrium state.The inset shows the time evolution of ρ , of a nonlinearlydamped squeezed two-mode vacuum with squeezing param-eter r = 1 /
20, individual bath configuration, for T = 0 anddamping rate λ = Γ = 10 − ω . amount of two quanta. The evolution of ρ , is influ-enced by the elements ρ , and ρ , which in the RWAapproximation contribute to decoherence and hence theasymptotic negativity decays to zero. The detailed cal-culation is given in Appendix A 3. The resulting timeevolution of ρ , is shown in the inset of Fig. 5.Numerically solving the full QME at T = 0, we seea residual nonzero negativity. This is due to the RWA-Hamiltonian ( H S ∝ H + λ ( a † a + a † a )) not properlyreproducing the correct ground state of the coupled os-cillator system since terms proportional to a ( a † a † ) areneglected. Hence, the numerical results display a smallresidual entanglement. In the inset of Fig. 5 the compar-ison of the results of the numerical simulation and theRWA shows a very good agreement. V. CONCLUSIONS
We have studied the asymptotic behavior of entangle-ment between two harmonic oscillators when they arequadratically coupled to an environment. In particular,we have investigated to what extent phenomena, knownfrom studying the decay of two-mode squeezed states inthe corresponding linearly damped systems, change whendamping is nonlinear. We find that the number par- ity conservation associated with pure nonlinear damp-ing causes significant reduction of the disentanglementrate. Moreover, the equilibrium distribution is differentfrom the standard Bose-distribution. Further, in con-trast to the linearly damped systems, we find no quali-tative difference between oscillators coupled to commonbaths and coupled to individual baths. We attribute thelatter effect to the lack of a conserved quantity (relativeoscillator energy) in the nonlinearly coupled system incompination with a suppressed information exhange be-tween the oscillators at low temperatures. For weaklycoupled oscillators, the number parity is no longer indi-vidually conserved, hence, the system can relax to theground state.The results here are obtained in the Markovian limit.Extending the study to a non-Markovian dynamics couldalter the picture presented here. For instance, it is knownthat for linearly damped oscillator systems studied withnon-Markovian dissipation models [17–19], a more de-tailed and complex picture of the asymptotic behavioremerges. Still, in those studies the overall characteristicsobtained in the Markovian limit, for instance the divi-sion into entangled and separable steady states, remainintact.At present it is not known whether it is possible to real-ize a system where the dominant dissipation mechanismat low excitation levels is purely nonlinear. However, forsome dissipation mechanisms, see for instance Ref. [28],symmetry can dictate that the lowest order coupling tothe environment must be quadratic in the coordinates.Systems with such symmetries are thus strong candidatesfor studying NLD in the quantum regime. Moreover, itwas suggested that engineering of NLD might be feasi-ble [29–31]. Exploiting the reduced disentanglement insystems with NLD is a promising path towards realizingentanglement based technologies.
Acknowledgments
Financial support for this work was provided by theSwedish Research Council VR.
Appendix A1. Negativity exponent - individual baths.
Here it is shown that γ N = 2 N (2 ω ). As argued inSec. III A, N ( τ ) of a nonlinearly damped squeezed stateasymptotically only depends on the density matrix ele-ment ρ , . Equation (9) is used to obtain the respectiveequation of motion, along with equations of motion forelements X = √ ρ , + ρ , ), which influence ρ , .By assuming that the other elements have already de-cayed, one obtains a set of two coupled, first order differ-ential equations˙ ρ , = 2 γ − X − γ ρ , , (A1a)˙ X = 12 γ ρ , − aX, (A1b)where a = ( γ − + 5 γ ) and γ ± is scaled by Γ . Forsolutions of the form of ρ , ∼ Ae r τ + Be r τ on finds, r , = − (16 γ +2) ± γ +1) − γ (21 γ +1)] , (A2)where γ − = γ + 1 was used. For low T , γ = N (2 ω ) <
1, the square root in (A2) can be Taylor ex-panded up to second order to yield r ≈ − N (2 ω ) , r ≈ − N (2 ω ) + 2) , (A3)and ρ , ( τ ) = Ae − N (2 ω ) τ + Be − N (2 ω )+2) τ , (A4)where the first term denotes the slow thermal decay. Thesecond term describes a rapid initial decay of the matrixelement. The amplitudes are given by A = ρ , (0)(11 N (2 ω ) + 2)14 N (2 ω ) + 2 , (A5a) B = ρ , (0) (cid:104) − N (2 ω ) + 214 N (2 ω ) + 2 (cid:105) , (A5b)where ρ , (0) is given by (12). Also, A > r > A (cid:29) B for low T . Hence γ N = 2 N (2 ω ) τ .
2. Negativity exponent - common bath.
Here we verify that the decay of entanglement, as gov-erned by the matrix element ρ , , for the individualbath case is γ N = 2 N (2 ω ). With the same assump-tions as for the individual baths, from equation (9) oneobtains the equations of motion for the matrix elementsresponsible for the negativity in the common bath case˙ ρ , = 2 √ γ − Z − γ ρ , , (A6a)˙ Z = 8 √ γ ρ , − γ − + 3 γ ) Z, (A6b) where Z = ρ , + ρ , + ρ , + ρ , . In this casewe obtain the solution ρ , = Ce − N (2 ω ) τ + De − N (2 ω )+4) τ , (A7)with the amplitudes C = ρ , (0)(15 N (2 ω ) + 4)18 N (2 ω ) + 4 (A8a) D = ρ , (0) (cid:104) − N (2 ω ) + 418 N (2 ω ) + 4 (cid:105) , (A8b)where ρ , ( τ = 0) is given in (12). Also, C > r > C (cid:29) D for low T . Hence the decay rate isagain dominated by γ N = 2 N (2 ω ).
3. Negativity evolution - individual bath, non-zerointer-mode coupling.
For T = 0 and an inter-mode coupling λ > ρ , element.For an initially two-mode squeezed state, the equationsof motion governing the evolution are˙ ρ , = i λ √ Y , (A9a)˙ Y = (i √ λ − √ ρ , − ( γ − + 2Υ)Y , (A9b)where Y = ρ , + ρ , . For a solution of the form ρ , = A + e r + τ + A − e r − τ one finds r ± = −
12 (Γ + 2Υ) ± (cid:112) (Γ + 2Υ) − λ + i2 λ Υ) , (A10)with amplitudes A + = − ρ , (0) r − ( r + − r − ) , (A11a) A − = ρ , (0) (cid:104) r − ( r + − r − ) (cid:105) . (A11b) [1] A. Einstein, B. Podolsky, and N. Rosen, Phys. Rev. ,777 (1935).[2] E. Schr¨odinger, Naturwiss. , 807 (1935).[3] R. Horodecki, P. Horodecki, M. Horodecki, andK. Horodecki, Rev. Mod. Phys. , 865 (2009).[4] A. K. Ekert, Phys. Rev. Lett. , 661 (1991).[5] C. H. Bennett and S. J. Wiesner, Phys. Rev. Lett. ,2881 (1992).[6] P. W. Shor, Phys. Rev. A , R2493 (1995).[7] B. Yurke and D. Stoler, Phys. Rev. A , 2229 (1992). [8] C. H. Bennett, G. Brassard, C. Cr´epeau, R. Jozsa,A. Peres, and W. K. Wootters, Phys. Rev. Lett. , 1895(1993).[9] S. Bose, V. Vedral, and P. L. Knight, Phys. Rev. A ,822 (1998).[10] D. Bouwmeester, J.-W. Pan, K. Mattle, M. Eibl, H. We-infurter, and A. Zeilinger, Nature , 575 (1997).[11] D. Boschi, S. Branca, F. De Martini, L. Hardy, andS. Popescu, Phys. Rev. Lett. , 1121 (1998).[12] A. Furusawa, J. L. Sørensen, S. L. Braunstein, C. A. Fuchs, H. J. Kimble, and E. S. Polzik, Science , 706(1998).[13] W. Pfaff, B. J. Hensen, H. Bernien, S. B. van Dam, M. S.Blok, T. H. Taminiau, M. J. Tiggelman, R. N. Schouten,M. Markham, D. J. Twitchen, et al., Science , 532(2014).[14] J. S. Prauzner-Bechcicki, Journal of Physics A: Mathe-matical and General , L173 (2004).[15] F. Benatti and R. Floreanini, Journal of Physics A:Mathematical and General , 2689 (2006).[16] S. Maniscalco, S. Olivares, and M. G. A. Paris, Phys.Rev. A , 062119 (2007).[17] K.-L. Liu and H.-S. Goan, Phys. Rev. A , 022312(2007).[18] J. P. Paz and A. J. Roncaglia, Phys. Rev. Lett. ,220401 (2008).[19] J. P. Paz and A. J. Roncaglia, Phys. Rev. A , 032102(2009).[20] F. Galve, L. A. Pach´on, and D. Zueco, Phys. Rev. Lett. , 180501 (2010).[21] T. Yu and J. H. Eberly, Phys. Rev. Lett. , 140404(2004).[22] A. Serafini, F. Illuminati, M. G. A. Paris, andS. De Siena, Phys. Rev. A , 022318 (2004).[23] H. D. Simaan and R. Loudon, J. Phys. A , 435 (1978).[24] A. Kronwald, F. Marquardt, and A. A. Clerk, Phys. Rev.A , 063833 (2013).[25] A. Voje, A. Croy, and A. Isacsson, New Journal ofPhysics , 053041 (2013).[26] M. Dykman and M. Krivoglaz, Soviet Scientific Reviews, Section A, Physics Reviews , 265 (1984).[27] A. Eichler, J. Moser, J. Chaste, M. Zdrojek, I. Wilson-Rae, and A. Bachtold, Nat. Nanotechnol. , 339 (2011).[28] A. Croy, D. Midtvedt, A. Isacsson, and J. M. Kinaret,Phys. Rev. B , 235435 (2012).[29] A. Nunnenkamp, K. Børkje, J. G. E. Harris, and S. M.Girvin, Phys. Rev. A , 021806 (2010).[30] W. Leo´nski and A. Miranowicz, J. Opt. B , S37 (2004).[31] M. J. Everitt, T. P. Spiller, G. J. Milburn, R. D. Wilson,and A. M. Zagoskin, Quantum Computing , 1 (2014).[32] A. Voje, A. Isacsson, and A. Croy, Phys. Rev. A ,022309 (2013).[33] R. Simon, Phys. Rev. Lett. , 2726 (2000).[34] S. L. Braunstein and P. van Loock, Rev. Mod. Phys. ,513 (2005).[35] H. P. Breuer and F. Petruccione, The Theory of OpenQuantum Systems (Oxford University Press, USA, 2002).[36] R. K. Wangsness and F. Bloch, Phys. Rev. , 728(1953).[37] F. Bloch, Phys. Rev. , 1206 (1957).[38] A. G. Redfield, Adv. Magn. Reson. , 1 (1965).[39] P. Drummond and Z. Ficek, Quantum Squeezing , Physicsand Astronomy Online Library (Springer, 2004), ISBN9783540659891.[40] G. Vidal and R. F. Werner, Phys. Rev. A , 032314(2002).[41] M. de Ponte, M. de Oliveira, and M. Moussa, Ann. Phys.317