Exact dynamics of quantum correlations of two qubits coupled to bosonic baths

The counter rotating-wave term (CRT) effects from the system-bath coherence on the dynamics of quantum correlation of two qubits in two independent baths and a common bath are systematically investigated. The hierarchy approach is extended to solve the relevant spin boson models with the Lorentzian spectrum, the exact dynamics for the quantum entanglement and quantum discord (QD) are evaluated, and the comparisons with previous ones under the rotating-wave approximation are performed. For the two independent baths, beyond the weak system-bath coupling, the CRT essentially changes the evolution of both entanglement and QD. With the increase of the coupling, the revival of the entanglement is suppressed dramatically and finally disappears, and the QD becomes smaller monotonically. For the common bath, the entanglement is also suppressed by the CRT, but the QD shows quite different behaviors, if initiated from the correlated Bell states. In the non-Markovian regime, the QD is almost not influenced by the CRT and generally finite in the long time evolution at any coupling, while in the Markovian regime, it is significantly enhanced with the strong coupling.

ence, which is of fundamental importance in the quantum information processing and measurement [25,26]. Specially, Yu and Eberly [27] found that the dynamical behavior of the entanglement is significantly different from the single qubit in a pair of qubits separately coupling with Markovian baths, which can be modeled from the well-known spin boson models [28]. It shows finite time disentanglement called as entanglement sudden death (ESD) [27]. However, this was investigated under quantum master equation with the Born and Markov approximation, which assumes that system-bath interaction is weak and the relaxation of the bath is much faster than that of the system. Then further works have been done in the non-Markovian regime in the framework of the rotating-wave approximation (RWA) to reveal novel characters of the quantum correlations compared to the Markovian effect [29][30][31][32][33][34], which is only valid under the condition of weak system-bath interaction. Besides, entanglement dynamics beyond the weak coupling spin boson models has also been exploited based on the optimal polaron transformation [35], where the higher order terms are neglected. It is pointed out that the corresponding accuracy in the wide parameter zone should be carefully analyzed [36]. Therefore, It is highly desirable to solve the dynamics without performing any approximation.
Beyond both the Born-Markov approximation and RWA, Tanimura et al. [37] have developed an efficient hierarchy approach, which has been later extensively applied to some chemical and biophysical systems [38,39]. Recently, this method has also been extended to study the dynamics of entanglement for two qubits coupled to a common bath [40,41]. However, the dynamics of QD, one recently mostly studied quantum correlations has not been investigated beyond the Born-Markov approximation and the RWA. Moreover, the exact study of the dy-namics of both entanglement and QD for two qubits coupled to their own baths are not available in the literature, to the best of our knowledge. This case is designed for two remote qubits each interacting locally with its own environment. Such a physical condition is also relevant to the quantum information science and quantum computation.
In the present paper, to give a comprehensive picture of the two qubits coupled to their own bath and common bath, we will extend the hierarchy equation to these two kinds of baths to study the dynamics of the pairwise entanglement and the QD in both Markovian and non-Markovian regimes. The effect of the counter-rotating term on these dynamics in a wide coupling regime will be systematically studied. This paper is organized as follows. In Sec. II we briefly review the definition of the QD. In Sec. III we exactly solve the reduced qubits density matrix of two independent spin boson models, and investigate its quantum correlation compared with that in RWA. In Sec. IV the quantum correlation of two separate qubits interacting with a common bath is analyzed. Finally, a summary is given in Sec. V.

II. QUANTUM DISCORD
The QD is one main route to fully measure the nonclassical correlation, which can be extended from the classical information theory [5,[14][15][16]. It is interpreted as the difference of two quantum mutual correlations of subsystems A and B, before and after local measurement operated on subsystem B [14,15]. The total quantum correlation of two subsystems is determined by their joint density matrix ρ as I(ρ) = S(ρ A ) + S(ρ B ) − S(ρ), where ρ A(B) = Tr B(A) {ρ}. The von Neumann entropy S(ρ a ) = −Tr{ρ a ln ρ a }. The other quantum mutual information J (ρ), which is equivalent to I(ρ) in classical information frame, is derived by locally measuring subsystem B with a complete set of orthonormal projectors {Π B k }, where k is one outcome state of B. After this measurement, the joint density is reduced to conditional counterpart as where the corresponding probability is p k = Tr{(I⊗Π B k )ρ}. Then the mutual information based on specific performance is shown as Q = S(ρ A ) − k p k S(ρ k ). Hence, we should maximize it to capture all classical correlation by obtaining J (ρ) = max {Π B k } {Q}. Finally, the QD is defined as The quantum nature of the system can be explicitly observed from the QD. It is only zero in pure classical correlation condition, and has the same values as the entanglement in pure state. While it remains finite even in separable mixed states, where the corresponding entanglement may completely disappear [14].

III. TWO INDEPENDENT SPIN-BOSON MODELS
The Hamiltonian of two independent qubits coupled two independent bosonic baths A and B is given by H = where a = A, B, σ a β (β = x, y, z) is the Pauli operator with Zeeman energy ω a of spin β, b † a,k (b a,k ) creates (annihilates) one boson with frequency ω a k in bath a, and g a k is the coupling strength between the system and the bath.
The spectra density of the bosonic bath is selected as the following Lorentz distribution which can describe an bosonic field inside an imperfect cavity mode ω a 0 with the system-bath coupling strength λ a . The corresponding correlation function at zero temperature is shown as (A13) Then by applying the Feynman-Vernon influence functional procedure in Appendix, the hierarchy equation is derived by (A16) The commutation relations A × B = AB−BA and A o B = AB + BA. The unit vectors e 1 = (1, 0), e 2 = (0, 1), and n = (n 1 , n 2 ), m = (m 1 , m 2 ). Q a (a = A, B) in the system-bath interaction represents σ a x . At the initial state, only ρ (0,0),(0,0) (0) is not zero. For the simplest condition in the following calculations, the two subsystems are identical as ω a = ω 0 , λ a = λ, γ a = γ.
For comparison, we thereby briefly review the main results derived in the RWA, where system-bath interaction becomes H R sb = a,k (g a k σ a + b a,k + g a * k σ a − b † a,k ). When the initial system reduced density matrix in the basis then we have where with R = λγ/2 − γ 2 /4. It is known that the bath relaxation time is τ b ≈1/γ, and the system-bath correlation time is τ r ≈1/λ. From the time dependence of the function P t , we immediately know that the Markovian regime corresponds to γ/2 > λ where P t decays exponentially, while the non-Markovian regime to γ/2 < λ where P t shows oscillating behavior describing the coherent process between the system and the bath. This non-Markovian condition has been realized in cavity quantum electrodynamics having Rydberg atoms inside the Fabry-Perot cavities with γ/λ≈0.1 [42]. We should point out here that the Markovian and non-Markovian regimes are only defined in the framework of the RWA, which can not be well defined in the exact treatment without the RWA. Throughout the present paper, we still use such notations. Doing so is only for comparisons and analysis.
As usual, the initial Bell states with anti-correlated spins and correlated spins will be considered in the present paper which both obey X form in Eq. (6) during the evolution. So it is straightforward to derive the concurrence, a pairwise entanglement [29] C Φ (t) = max{0, 2α 1 − α 2 P t }, The analytical expression for QD in this case is quite complicated and also only consists of α and P t , which is implied in Refs. [33,43,44].
It is shown analytically above that the dynamics of both entanglement and QD under the RWA can be distinguished in the non-Markovian and Markovian regimes, characterized by the ratio γ/λ. In each regime, the essential feature only depends on the weights α 2 of the initial states, and independent of the system-bath coupling strength λ. In other words, in the framework of the RWA, all essential properties for the physical observable are independent of the system-bath coupling strength, which is not always true with the increasing coupling realized in many recent experiments [45]. Hence, it is necessary to go beyond the RWA.

A. Exact dynamics and comparisons
The dynamics of the pairwise concurrence of the two qubits in their own baths without the RWA can be obtained exactly by the hierarchy approach outlined above. Initiated from the anti-correlated Bell state (9) with α = 1/ √ 2, the time evolution of the concurrence from weak, intermediate, to strong coupling cases are presented in the upper panel of Fig. 1. The concurrence under the RWA has been investigated exactly in Ref. [29] both in the Markovian and non-Markovian regimes. The relevant results for the same parameters are also reproduced in the upper panel of Fig. 1 for comparison.
In the weak system-bath coupling (λ = 0.02), the rotating-wave term dominates the entanglement evolution as indicated in Fig. 1(a), and the results including the CRT show negligible deviation from those in RWA. In the non-Markovian regime, the concurrence exhibits periodic oscillations with amplitude damping. The corresponding time of the zero points is t n = [nπ − arctan(R/2γ)]/R (n = 1, 2, · · · ). Whether it only vanishes at these discrete moments is very subtle, the periods of the time for zero entanglement is too short to be visible in this case . While in the Markovian regime, the concurrence decays exponentially and vanish asymptotically. The decay rate is larger than the non-Markovian one. So at such a weak coupling, we have not found any evident effect from the CRT, so the previous RWA description is really valid.
What happens if we increase the system-bath coupling to the intermediate regime, say λ = 0.5. It is interesting to observe in Fig. 1(b) for λ = 0.5 that the dynamic behavior is evidently different from that in the RWA. From Eq. (11), one can see that if initiated from anticorrelated Bell state, in the RWA the ESD never happens, which was first found in Ref. [29] in the same system. But without the RWA, in both the Markovian (f = 5) and the non-Markovian (f = 0.1) regimes, the concurrence exhibits ESD obviously. It follows that the RWA can not describe the essential feature of the evolution of the entanglement beyond the weak coupling regime and is therefore broken down under strong coupling.
Importantly, this physical condition can be realized in the recent strong coupling experiments [45]. In the non-Markovian limit (γ/λ≪1), the spectral function of baths is reduced to the single mode case by J(ω)≈ λγ 2 δ(ω − ω 0 ). Hence, the system can be approximately simplified to two independent Rabi models. The corresponding effective system-bath coupling constant is given by g = λγ 2 , so g/ω 0 = 0.11 for f = 0.1, which can be practically realized in the recent experiments of the superconducting qubits coupled to LC resonators, where the coupling constant g/ω 0 already reached to ten percentages [45].
In the Markovian regime, the entanglement vanishes quickly and permanently beyond the weak coupling as shown in Fig. 1 (b) and 1(c), mainly due to the generation of extra phonons with increasing coupling, which disturb the pairwise entanglement. No revival of the entanglement appears in this case. Interestingly, the evolution behavior is almost unchanged as the increasing coupling strength for λ ≥ 0.5, which follows that λ = 0.5 is strong enough in this case. However, in the non-Markovian regime, the entanglement dynamics show essential different behavior with the system-bath coupling. The revival of the entanglement after the ESD occurs at λ = 0.5, but this revival phenomena never happens in the strong coupling regime, as shown in Fig. 1(c) for λ = 2. We propose that the more extra phonons activated by the CRT would suppress the feedback from the bath to the qubit [48] even in the non-Markovian regime, resulting in the permanent disentanglement. The effective systembath coupling strength for f = 0.1 in Fig. 1(c) can be estimated as g/ω 0 = 0.45, which is hopefully realized in the near future [46,47].
Then we investigate the dynamics of the QD for two qubits in two independent baths to explore the influence of the CRT. The QD in the RWA has been discussed by Wang et al. [33] and Fanchini et al. [43]. In the lower panel of Fig. 1(d), we display the evolution of the QD for the same parameters as the upper panel by both the hierarchy approach and RWA. In the weak coupling regime, it behaves similar to the concurrence in Fig. 1(a) for both the Markovian and non-Markovian evolutions, where the CRT can be ignored. Beyond the weak coupling regime, as shown in Fig. 1(e) and 1(f), the CRT drives the QD to deviate from those in RWA obviously. Specifically for the non-Markovian case, the QD decays dramatically and shows weak revival signal, which mainly attributes to the excitation of the phonons in baths. For the Markovian case, the QD decreases faster than that in RWA and decay to zero asymptotically.
To see the comprehensive effect of the system-bath interaction, we plot the dynamics of the quantum correlation as a function of both t and λ in Fig. 2 for initial anti-correlated Bell state. Both the entanglement and the discord are suppressed monotonically with the increasing interaction. The revival phenomena for entanglement is not present completely as the coupling strength exceeds a critical value, e.g. λ c is about 1.5 for the parameters in Fig. 2 (a). The region of ESD becomes wider with the increase of the system-bath coupling. But Fig. 2 (b) reveals that the QD never vanishes suddenly, consistent with the previous observation concluded from an arbi- trary Markovian evolution [51]. Therefore, such a limitation is removed in the present observation in the exact evolution which definitely includes the non-Markovian effect. No sudden death of the QD is universal in this sense, in sharp contrast with the entanglement. More over, in our exact solution for the full model without the RWA, we actually have not found the real zero QD at any time, despite arbitrary small. This nonzero nature for the QD maybe intrinsic at least for the spin-boson model. Previous observation in the RWA that the QD can vanish at some discreet times [33,43] may be an artificial result due to the absence of the CRT.
The essential feature for the dynamics of the quantum entanglement and discord under the RWA is not altered qualitatively with the system-bath coupling in the same non-Markovian or Markovian regime. However, it is not that case in the real exact solution without the RWA. For the weak system-bath interaction, we have shown above that the RWA can basically give the right description of the quantum correlation of two separate spin-boson models for both Markovian and non-Markovian regimes. While in the strong coupling regime, the CRT should be necessarily included for correctly depicting the phonons generation from the baths. And the quantum correlation shows significantly different behavior from that in RWA, especially in non-Markovian case. Hence, we focus on the dynamics of the quantum correlation in the strong coupling and non-Markovian regime in the following with various mostly used initial states.

B. Non-Markovian dynamics at strong coupling
First, we examine the dependence of the dynamics of the quantum correlation on the value of α for the initial Bell state in Eq. (9), which characterize weights of the two superposed anti-correlated spin states. Since the quantum correlation is similar in the regime α 2 ∈(0, 1/2) to that in (1/2, 1), the first zone is selected to study the correlation evolution in Fig. 3(a)-3(c). In RWA, both the concurrence and the QD show damping oscillations with stronger amplitude as α 2 increases. The other evident feature is that they only vanish at some discrete moments. When the CRT is included, the concurrence exhibits ESD after finite time evolution in the whole regime. The starting time of ESD is much earlier than first vanishing moment in the RWA. The QD first decreases almost to zero and then revives with amplitude much too weaker than that in the RWA, mainly attributed to the suppress from the phonons. Fig. 3(d)-3(f) presents the results for the another initial Bell state in Eq. (10). It is known that the discord dynamics for α 2 ∈(0, 1) is similar to the concurrence for α 2 ∈(1/2, 1). Hence, we also select these initial states for comparison as in Ref. [43]. The exact dynamics shows that the ESD occurs very quickly and irrelevant with the value of α. The QD first decreases almost to zero, and then revives weakly. The revival amplitude becomes stronger with the value of α. For the same α, the revival amplitude is much too weaker than that in the RWA, due to the CRT effect, as demonstrated in the inset.
Finally, we assume the extended Werner like state (ρ AB (0) = r|ξ ξ| + 1−r 4 I) as the initial state and see what happen to the evolution of these quantum correlations. The main results are list in Fig. 4. As one classical kind of the mixed states, it is important in quantum information processing [52][53][54]. As the initial state with |ξ = |Φ( 1 √ 2 ) , the concurrence in RWA shows complete death in r∈(0, 1/2), whereas the discord in the regime r∈(0, 1) behaves similar with the concurrence in r∈(1/2, 1) [33]. Hence, we study the quantum correlation in the regime r∈(1/2, 1) to explore the effect of the CRT in Fig. 4(a)-4(c). The entanglement vanishes suddenly and permanently for all values of r, RWA ones show damping oscillations for large r (e.g. r = 0.7, 0.9). The QD is also strongly suppressed first, and then revives weakly. The revival amplitude increase slightly with r. It is again shown that the QD is not so fragile as the entanglement. If |Ψ( 1 √ 2 ) is replaced by |ξ in the initial extended Werner like state, no essential different behavior can be observed except that the revival amplitude becomes much too weaker, as demonstrated in Fig. 4 For the three kinds of initial states studied above, it is generally find that the larger initial quantum correlations is, the stronger they evolves at the same moments.

IV. TWO QUBITS COUPLED TO ONE COMMON BOSONIC BATH
The Hamiltonian of two qubits coupled to the common bath H = H s + H b + H sb reads where the notations are the same as those in Eq. (2). If the spectrum of the common bath is also Lorentz distribution J(ω) = 1 2π λγ 2 (ω−ω0) 2 +γ 2 , and the correlation function of the bath should be also C(t) = λγ 2 e −(γ+iω0)t . By choosing the proper auxiliary influence functional in Appendix, the hierarchy equation is given by (A17) [41] where n = (n 1 , n 2 ), the unit vector e 1 = (1, 0) and e 2 = (0, 1). The reduced density matrix for qubits is ρ AB (t) = ρ (0,0) (t), by which we can derive the quantum correlation numerically.
In RWA, the system-bath interaction becomes H R sb = (σ A If we represent the Hamiltonian in the basis {|0 = |00 , | + = 1 2 (|10 + |01 ), | − = 1 √ 2 (|10 − |01 ), |2 = |11 }, | − is found to be isolate from others. Hence we can rewrite the Hamiltonian in three level form [55] For the single particle excitation or double excitation case, the motion equation of the system density matrix can be derived by the called pseudomode mater equation in the interacting picture [43,[55][56][57][58][59] ∂ρ where a † (a) is the creator (annihilator) of the bosonic pseudomode, and ρ I (t) = e iH0t ρ(t)e −iH0t . By integrating the bosonic part of the density matrix, the system reduced density is given by ρ AB (t) = Tr a {ρ(t)}. For single excitation, an alternatively exact solution has also been done in Ref. [31]. If the initial state is chosen by X form in Eq. (6), the concurrence is derived by [59] C(t) = 2 max{0, |ρ 23 For the common bath, the entanglement evolution initiated from the anti-correlated Bell state without the RWA has been studied by Ma et al. [41]. On the other hand, in the previous studies of the dynamics of QD under the RWA, the correlated Bell state is usually selected as the initial state [33,43,55,59]. Therefore in the present exact study without the RWA, we will focus on the initial correlated Bell state.
First, we exhibit the dynamics of the pairwise concurrence of the two qubits coupled to a common bath with and without the RWA initiated from the Bell state in Eq. (9) with α = 1/ √ 3 in the upper panel of Fig. 5 for various coupling strength in both Markovian and non-Markovian regimes. Similarly, for the weak coupling (λ = 0.02), the CRTs take little effect for both dynamics in different evolution regimes, showing the same results as in Ref. [43]. In the intermediate coupling regime (λ = 0.5), the ESD occurs more quickly than that in the RWA in the Markovian regime, and the evolution shows slight deviations from the RWA ones with damping oscillations in the non-Markovian regime. When the  FIG. 5: Dynamics of quantum correlation of two qubits coupled to the common bath with initial correlated Bell state: strong coupling regime is reached, the CRT takes remarkable effects in both evolution regimes, where concurrence vanishes completely and permanently without the RWA. Hence, the CRT also suppress the entanglement dramatically as the system-bath coupling for the common bath, similar to those for the independent baths. Then we turn to the relevant QD, as shown in the lower panel of Fig. 5. It is interesting to note that in the non-Markovian regime, the CRTs take very limited effect with the increasing coupling. Even at the strong coupling, only slightly deviations are observed. More interestingly, the QD takes finite value at any time, and take a moderate finite value asymptotically in the long time limit, for arbitrary coupling strength. In the Markovian regime, the QD tends to a remarkable value in the intermediate and strong coupling regime. Actually, even at weak coupling like λ = 0.05, the asymptotical value of QD takes a considerable non-zero value.
From the numerical iteration of the hierarchy equation in Eq. (14) with Bell correlated initial state under Markov regime (e.g. f = 10 in Fig. 5), we find that the form of the two qubits reduced density matrix under the same basis as in Eq. (6) always to be X type as In the weak system-bath coupling (say λ = 0.02), a(t) and b(t) decrease gradually to 0 after the finite time evolution, whereas d(t) rises to reach 1. Hence both concurrence and the discord decay to 0 exhibit monotonically as indicated in Fig. 5 that in the RWA. When the coupling strength becomes stronger, e.g. λ = 0.5, 2.0, the phonons are excited considerably due to the strong coupling [48]. According to the CRT k g * , the qubits are also excited accompanying the extra excitations of the phonons, which stabilize a(t) and b(t) in long time evolution. It is interesting to observe that they tend to 1/3 and 1/6 in the strong coupling limit at t→∞. The corresponding reduced density matrix is then described by ρ(t→∞) = 1 3 (|00 00| + |11 11| + | + + |). which is not altered by different α 2 . Under this separable mixed state, the disentanglement occurs naturally from Eq. (17), and the QD is stabilized at 0.23.
To see the overall effect of the system-bath interaction on the dynamics of the quantum correlation beyond the weak coupling regime, we draw the corresponding 3D plots for λ ≥ 0.1 in Fig. 6 in the non-Markovian regime. For the entanglement, one can find that there are two regimes. One is damping oscillation at λ < 1. The other is the regime where the entanglement revives following the ESD, and the revival will be suppressed and finally disappear with the further increase of the coupling regime. Very interestingly, after weak oscillations, the QD is generally robust in the long time evolution in arbitrary coupling regime, which is very useful as the quantum information resource.

V. CONCLUSION
In summary, without performing any approximation, we exactly investigate quantum correlations in two typical kinds of spin-boson models by applying the hierar-chy approach. The corresponding RWA results are also reproduced for comparison. The Markovian and non-Markovian effects of two kinds baths in strong systembath interaction on the dynamics of the quantum correlation are revealed in detail, which demonstrates that the CRT should be necessarily included beyond the weak coupling regime to capture the essential features of quantum correlations correctly.
For two independent baths, the CRT suppresses both the entanglement and the QD dramatically, due to the activations of phonons in baths as the system-bath coupling increases. If initiated from the anti-correlated Bell state, The ESD is driven by extra phonons activated with the increasing system-bath coupling, which is absent in the RWA study. It follows that previous picture under the RWA is essentially modified with the increasing coupling. The QD is found to be always higher than zero, despite extremely small, also different from that in the RWA. The dynamics with other initial conditions is also discussed at the strong coupling and non-Markovian regime. The general trend is that the larger initial quantum correlations, the stronger they evolve.
For the model of two qubits coupled to the common bath, the suppression of the entanglement by the CRT is also observed. It is interesting to find that the nonclassical correlation reserved by the QD can be enhanced as the coupling strength, in contrast with the entanglement. It follows that while the entanglement becomes fragile with the increasing coupling to the surrounding environment, the QD is however robust in the long time evolution. In recently experiments, the strong system-bath coupling has been reached in circuit QED systems [45], where g/ω 0 ≈0.1 for multi modes and the RWA description is broken down. The present study for the evolution of the quantum correlations beyond the RWA may motivate the corresponding experimental studies based on these strong coupling systems.

VI. ACKNOWLEDGEMENT
It is known that ρ I s (t) = Tr b {ρ I (t)}. Hence by integrating the bath degree of freedom, the reduced system density is Considering the propagator in functional form and the initial factorized state ρ(0) = ρ s (0)⊗ρ b (ρ b = exp(−βH b )/Z is the equilibrium state), the reduced density matrix is reexpressed by Then the functional field is introduced by with corresponding bath correlation function Here we should note that the variable Q[α] is a real number. Then we generalize the results to the multimode case, where the system-bath interaction now is H sb = aQ aFa , andF a = k c a kx a k . The influence functional becomes The correlation function in independent baths is In the schrodinger picture, the motion equation of the system density is ∂t .
The key step here is to connect time differentiation of the density matrix with that of the influence functional.
The proper way is employing the hierarchy equation by introducing the auxiliary influence functional While such auxiliary choice is closely dependent on the specific form of the bath correlation function. The hierarchy equation used in the present paper is described in the following. For two independent baths with Lorentz type spectral function J a (ω) = 1 2π λaγ 2 a (ω−ω a 0 ) 2 +γ 2 a , the correlation function simply is according to the following formula we introduce the auxiliary influence functional of each subsystem as where Q a = σ a x , (a = A, B) and index k = (k 1 , k 2 ). Then the whole influence functional becomes Finally, the evolution of the density matrix reads ∂ ∂t ρ n, m (t) = −(iH × S + n · ν A + m · ν B )ρ n, m (t) (A16) + k=1,2 (−1) k Q × A ρ n+ e k , m (t) where unit vector e 1 = (1, 0) and e 2 = (0, 1) ν a = (ν a + , ν a − ), with ν a + = γ a + iω a 0 and ν a − = γ a − iω a 0 . A × B = AB − BA and A o B = AB + BA.
For common bosonic bath with the similar Lorentz distribution, if we choose the auxiliary influence functional (A17)