Coupled modes locally interacting with qubits: critical assessment of the rotating wave approximation
CCoupled modes locally interacting with qubits: critical assessment of the rotatingwave approximation
P. C. C´ardenas
Departamento de F´ısica y Matem´aticas, Universidad Aut´onoma de Manizales,Antigua Estaci´on del Ferrocarril, Manizales, Caldas, Colombia
W. S. Teixeira and F. L. Semi˜ao
Centro de Ciˆencias Naturais e Humanas, Universidade Federal do ABC, Santo Andr´e, 09210-170 S˜ao Paulo, Brazil (Dated: June 6, 2017)The interaction of qubits with quantized modes of electromagnetic fields has been largely ad-dressed in the quantum optics literature under the rotating wave approximation (RWA), whererapid oscillating terms in the qubit-mode interaction picture Hamiltonian can be neglected. At thesame time, it is generally accepted that provided the interaction is sufficiently strong or for longtimes, the RWA tends to describe physical phenomena incorrectly. In this work, we extend theinvestigation of the validity of the RWA to a more involved setup where two qubit-mode subsys-tems are brought to interaction through their harmonic coordinates. Our treatment is all analyticthanks to a sequence of carefully chosen unitary transformations which allows us to diagonlize theHamiltonian within and without the RWA. By also considering qubit dephasing, we find that thepurity of the two-qubit state presents non-Markovian features which become more pronounced asthe coupling between the modes gets stronger and the RWA loses its validity. In the same regime,there occurs fast generation of entanglement between the qubits which is also not correctly describedunder the RWA. The setup and results presented here clearly show the limitations of the RWA in ascenario amenable to exact description and free from numerical uncertainties. Consequently, it maybe of interest for the community working with cavity or circuit quantum electrodynamic systems inthe strong coupling regime.
I. INTRODUCTION
Modeling physical phenomena via the coupling of two-level systems (qubits) to quantized harmonic oscillatorshas historically been of great interest in diverse fieldsranging from quantum optics [1–4] and solid-state physics[5, 6] to quantum biology [7–10]. This approach has be-come frequent in modern physics since fully quantum-mechanical descriptions may reveal phenomena not cov-ered by classical or semiclassical approaches. As well-known examples, one has the study of dissipation anddecoherence of a qubit via the spin-boson model [11] orthe presence of collapses and revivals of the atomic pop-ulation inversion in the atom-field interaction [3, 4].For the latter case, the simplest quantum descriptionis given by the exactly solvable Jaynes-Cummings (JC)model [2], which considers a two-level atom weakly cou-pled to a single mode of the electromagnetic field. Usu-ally, the “energy non-conserving’ terms in the atom-fieldinteraction Hamiltonian are neglected through the so-called rotating wave approximation (RWA). However, theuse of the RWA in this problem may not describe dynami-cal properties of the model correctly when the atom-fieldinteraction becomes sufficiently strong [12, 13]. Otherstudies have addressed the limitations of the RWA indiverse configurations [14–17, 19–21]. Only recently ithas been shown that the non-RWA atom-field Hamil-tonian possesses a symmetry rendering the model inte-grable [22]. However, due to the lack of closed form ex-pressions for the eigenstates, one usually has to appeal todifferent effective approaches to treat the problem with- out the RWA. This includes the use perturbation seriesfor path-integrals [12] or particular regimes such as far-from-resonance cases (dispersive limit) [18].Recently, the fabrication of artificial atoms with su-perconducting circuits [23–28] has favored the control ofqubit-oscillator interactions, and thus regimes where theRWA breaks down can be explored experimentally. Incircuit quantum electrodynamics (circuit QED) a qubitor two-level system can be produced, for instance, byusing a thin insulator between two superconducting ma-terials [Josephson junction (JJ)], and controlling eitherthe number of Cooper pairs that tunnel from one sideto the other (charge qubit) or the phase of their wave-functions (phase qubit). Also, by adding one or moreJJ in a superconducting loop [Superconducting Quan-tum Interference Device (SQUID)], a qubit can be pro-duced by controlling the external magnetic flux throughthe SQUID (flux qubit) [23, 24]. On the other hand, aquantum harmonic oscillator naturally represents a sin-gle electromagnetic mode trapped in a transmission line[24, 26].Giving all these developments in the control of sim-ple systems consisting of qubits and bosonic modes, itis natural to search for configurations which allow usto further understand the limits of the usually takenRWA. We explore the validity of the RWA in a setupcomprised of two identical qubit-oscillator systems thatare coupled through their harmonic coordinates. Thissetup is amenable to implementation in superconductingcircuits as discussed in [29], where non-Markovian fea-tures are discussed within the RWA. The present work a r X i v : . [ qu a n t - ph ] J un is organized as follows. We analytically diagonalize thefull non-RWA Hamiltonian in (Sec. III A) and then solvethe master equation which includes dephasing for thequbits (Sec. III B). Predictions contrasting our full treat-ment with the RWA model in [29] are then presented, inparticular for the two-qubit subsystem for which purity(Sec. III C) and entanglement dynamics (Sec. III D) areinvestigated. Finally, in Sec. IV we present our conclu-sions. II. THE MODEL
In this work, we are interested in the setup depictedin Fig.1. It consists of two groups of subsystems ( α and β ), each comprising a qubit and a bosonic mode.These groups will interact through the modes. Insideeach group, the local interaction takes the usual spin-boson form with no transverse field ( (cid:126) = 1) H SB = g ( σ z A + µI A ) (cid:0) a † + a (cid:1) + g ( σ z B + µI B ) (cid:0) b † + b (cid:1) , (1)where σ z A ( B ) are the Pauli matrices for the qubits; a † ( a )and b † ( b ) are the creation (annihilation) operators forthe corresponding bosonic modes, g is a coupling con-stant, and I A ( B ) is the identity operator acting on thecorresponding qubit state space. The µ -terms naturallyappear in some circuit QED architectures when workingout of the so-called degeneracy point [26, 27], so that, forcompleteness, they are included here. More details canbe found in [29].The total energy of the full setup reads H j = H + H SB + H BB j , (2)with the free Hamiltonian H = ω σ z A + σ z B ) + ω (cid:0) a † a + b † b (cid:1) , (3)where ω is the resonance frequency of the qubits and ω isthe angular frequency of the modes. The index j ∈ { , } S α S β gγ Mode a Qubit A gγ Mode b Qubit B λ FIG. 1. (Color online) Two subsystems S α and S β , eachformed by a qubit and a bosonic mode, interact through themodes (coupling strength λ ). Each qubit is also subjected todephasing at rate γ . Inside each subsystem, the qubit andthe mode interact with each other (coupling strength g ). defines the form of the interaction mechanism betweenthe modes which can be either H BB = λ (cid:0) a † + a (cid:1) (cid:0) b † + b (cid:1) (4)or H BB = λ (cid:0) a † b + ab † (cid:1) , (5)with λ being a coupling constant. In the first approach( j = 1), the spatial coordinates of each oscillator are cou-pled in a kind of quadrature-quadrature form. In circuitQED, this can be induced by coupling the two transmis-sion lines to an auxiliary qubit that mediates a geometricsecond-order interaction [30], while in cavity QED thiscan be done by placing a partially reflecting mirror be-tween two optical cavities [31]. On the other hand, thesecond approach ( j = 2) arises from the RWA performedon H BB so that its oscillating terms are neglected inthe interaction picture. For the interaction of a two-levelatom with an electromagnetic mode, which in the RWAgives rise to the JC Hamiltonian, one can only comparethe RWA and non-RWA Hamiltonians through succes-sive approximations or numerics. Here, we will be ableto perform such investigation in a fully analytic mannerby exactly solving˙ ρ j = − i [ H j , ρ j ] + γ σ z A ρ j σ z A + σ z B ρ j σ z B − ρ j ) , (6)for initial states of interest. In Eq.(6), γ is the rate ofpure dephasing caused by independent Markovian bathsacting on the qubits. This is by far the most rele-vant noise when working outside the degeneracy point[32, 33]. Energy relaxation of qubits or the transmis-sion lines (bosonic modes), as well as dephasing on thelatter, can be made negligible compared to dephasing inthe qubits [26, 27]. These experimental facts provide us abackgound to leave aside noise mechanisms besides qubitdephasing as a first approximation. This is especiallyconvenient in our case because our goal is to provide ana-lytical expressions that evidence inadequacy of the RWAin certain regimes. Those neglected noise mechanismswould render the problem unsuitable to analytic treat-ment and can be numerically investigated elsewhere. III. RESULTSA. Diagonalization
In order to analytically solve Eq.(6), we start by di-agonalizing the Hamiltonians H j . This is achieved bythe unitary transformation U j given by U j = P j S j T D j where D j = e δ j ( a † − a + b † − b ) (7)is the displacement operator with δ j = gµ/ (cid:0) ω + 2 − j λ (cid:1) , T = e π ( a † b − ab † ) (8)corresponds to a beam-splitter operation, S j = e − rj +2 ( a − a † ) e − rj − ( b − b † ) (9)is a squeezing operator with r ± = ln (1 ± λ/ω ) and r ± = 0, and P j = e λ j + ( σ zB + σ zA )( a † − a ) e λ j − ( σ zB − σ zA )( b † − b ) (10)is a polaron transformation [34] and λ j ± = ge − r j ± / ( √ j ± ). Notice that the operation S j re-duces to the identity for j = 2. This is so becausethis transformation is responsible for the elimination ofterms with ab and a † b † , not present in Eq. (5). It isimportant to realize that the application of S requires λ < ω/
2, so that the frequencies of the normal modes arereal numbers. Values of λ close to such limit have beenassociated to quantum chaos in nonlinear oscillators [35].With the help of these transformations, we obtain thediagonal Hamiltonian H (cid:48) j = U j H j U † j that reads H (cid:48) j = ω ,j σ z A + σ z B ) + χ j σ z A σ z B + Ω j + a † a + Ω j − b † b, (11)with shifted qubit frequencies ω ,j = ω − gδ j (12)and normal mode frequenciesΩ j ± = ω cosh (cid:0) r j ± (cid:1) ± λe − r j ± . (13)The Hamiltonian (11) is an interesting physical result. Itimplies that, in spite of the model ( j = 1 or j = 2), themodes decouple from the qubits, and the latter interactthrough an Ising-type Hamiltonian with χ j = 2 g (cid:18) e − r j − Ω j − − e − r j + Ω j + (cid:19) . (14)For uncoupled modes ( λ = 0), no effective coupling be-tween the qubits is observed ( χ j = 0). For finite λ , wecan already spot the fundamental differences in the non-RWA ( j = 1) and RWA ( j = 2) descriptions. This canbe seen from the plots in Fig. 2, where we present thedependence of ω ,j , χ j , and Ω j ± on the modes couplingconstant λ . These physical frequencies clearly indicatethat the RWA dismally fails with the increase of λ . B. Dynamics
In the space of Hamiltonian (11), the system densitymatrix is given by ρ (cid:48) j = U j ρ j U † j . A new transforma-tion defined as ρ I j = e iH (cid:48) j t ρ (cid:48) j e − iH (cid:48) j t finally allows one torewrite the master equation (6) in a very compact formas ˙ ρ I j = γ (cid:0) σ z A ρ I j σ z A + σ z B ρ I j σ z B − ρ I j (cid:1) . (15) ω , j / ω χ j / ω Ω j + / ω λ / ω Ω j - / ω FIG. 2. (Color online) Parameters of the diagonal Hamil-tonian (11) as functions of λ (in units of ω ). First row :Effective frequencies of the qubits.
Second row : Effectivecoupling strength of the qubits.
Third and forth rows :Effective frequencies of the modes. In all cases, the solidblue lines represent the complete model ( j = 1), whereas thedashed red lines represent the RWA model ( j = 2). We used ω = ω , g = 0 . ω , and µ = 1. Although the modes do not appear explicitly in Eq. (15),they have not yet been traced out. What happens is thatthey are frozen in this interaction picture and were com-pletely decoupled from the qubits due to the transforma-tion U j = P j S j T D j . Consequently, we can solve Eq. (15)in the qubits subspace and tensor the result with the ini-tial state of the modes. The transformations back to theoriginal picture will then restore the time evolution of thewhole system, entangling modes and qubits. By denot-ing | ψ (cid:105) A ⊗ | φ (cid:105) B as just | ψφ (cid:105) , we then solve Eq. (15) inthe standard basis {| ee (cid:105) , | eg (cid:105) , | ge (cid:105) , | gg (cid:105)} , in which | e ( g ) (cid:105) stands for the excited (ground) state of a single qubit. Inour case, we do not have to move the system state entirelyback to the original picture because we are interested inthe dynamics of the qubits. The polaron operation P † j is the only one required as the other unitary transfor-mations employed to diagonalize Eq. (2) are in fact localin the modes. The density matrix of the qubits is thenobtained as ρ AB j ( t ) = Tr ab (cid:104) P † j e − iH (cid:48) j t ρ I j e iH (cid:48) j t P j (cid:105) , (16)where the partial trace is taken over the modes degreesof freedom.From now on, we will focus on a particular choice of ini-tial states that suitably illustrates the dynamics of statepurity and entanglement for the qubits. On their own,these two quantities carry a lot of information and theiranalysis is then of general importance in quantum in-formation. More important to us, they both depend onthe whole density matrix and not only on its diagonalelements. Quantities such as occupation probabilities ofthe bare states would depend only on diagonal densitymatrix elements. For all these reasons, entropy and en-tanglement are then very good candidates to spot thedifferences between the full model and its RWA version.We consider the qubits to be initially prepared in theeigenstate of the Pauli matrix σ x A ( B ) associated with theeigenvalue +1, i.e., | + (cid:105) A ⊗ | + (cid:105) B ≡ | ++ (cid:105) = 12 ( | ee (cid:105) + | eg (cid:105) + | ge (cid:105) + | gg (cid:105) ) . (17)On the other hand, the modes are initially set in theproduct of coherent states, | α (cid:105) a ⊗ | β (cid:105) b ≡ | α (cid:105) | β (cid:105) = e − | α | | β | ∞ (cid:88) n,m =0 α n β m √ n ! m ! | nm (cid:105) , (18)where α and β are complex amplitudes, and | n (cid:105) ⊗ | m (cid:105) ≡| nm (cid:105) is the two-mode Fock state. Solving Eq. (15) for theinitial state | ψ (0) (cid:105) = | ++ (cid:105) | α (cid:105) | β (cid:105) and using Eq.(16), weobtained the 16 components of ρ AB j ( t ) in the standardbasis as ρ AB,mn j ( t ) = 14 γ mn j ( t )Θ mn,j + ( t )Θ mn,j − ( t ) (19)with m, n ∈ { ee, eg, ge, gg } and m (cid:54) = n . The diago-nal elements m = n are time-independent and given by ρ AB,mm j = 1 /
4. The factors γ mn j ( t ) arise from the non-unitary dynamics followed by each qubit (a dephasing factor), whereas Θ mn,j ± ( t ) are scalar products originatedfrom the partial trace operations. Explicit expressions for γ mn j ( t ) and details about Θ mn,j ± ( t ) can be found in theAppendix. C. Purity of the two-qubit subsystem
In general, quantum information processing requiresthat pure states, like superposition states, remain pureduring time evolution. However, when a quantum systemis interacting with others, its reduced dynamics will, ingeneral, affect the state purity. The same is valid whenthe system is in contact with a bath and, in this case, thelost of purity is typically irreversible. A good measureof how much a state ρ is pure in a d -dimensional statespace is given by a quantity called purity, defined as P =Tr (cid:2) ρ (cid:3) , with 1 /d ≤ P ≤ P = 1 ( P = 1 /d ), the system is in a pure (maximallymixed) state. For our purposes of contrasting RWA andnon-RWA descriptions, the purity is more suitable thanthe full entropy since the former allowed us to get analyticand exact expressions. For the initial conditions given byEqs. (17) and (18), we obtain P j ( t ) = Tr (cid:104) ρ AB j ( t ) (cid:105) = 14 + 2Γ j ( t ) , (20)withΓ j ( t ) = (cid:12)(cid:12) ρ AB,eeeg j ( t ) (cid:12)(cid:12) + (cid:12)(cid:12) ρ AB,eege j ( t ) (cid:12)(cid:12) + (cid:12)(cid:12) ρ AB,eegg j ( t ) (cid:12)(cid:12) + (cid:12)(cid:12) ρ AB,egge j ( t ) (cid:12)(cid:12) + (cid:12)(cid:12) ρ AB,eggg j ( t ) (cid:12)(cid:12) + (cid:12)(cid:12) ρ AB,gegg j ( t ) (cid:12)(cid:12) . (21)Explicitly, the function Γ j ( t ) readsΓ j ( t ) = 4 e − [ f j + ( t )+ f j − ( t )+2 γt ] + e − [ f j + ( t )+ γt ] + e − [ f j − ( t )+ γt ] , (22)with f j ± ( t ) = 16 λ j ± (cid:2) cosh(2 r j ± ) − sinh(2 r j ± ) cos(Ω j ± t ) (cid:3) sin (cid:18) Ω j ± t (cid:19) . (23)Figure 3 compares the dynamics of P j ( t ) for eachmodel under different mode coupling regimes. Evidently,the purity is maximum at t = 0 as the initial state of thequbits is | ++ (cid:105) and not entangled with the modes. Thefirst thing to be noticed is that there is a clear competi-tion between the Markovian dynamics induced by the de-phasing baths and the non-Markovian dynamics inducedby the mode-mode interaction [29]. To be more precise, the oscillations appear as the result of the latter while theenvelop (purity damping) is caused by the baths. Suchfeatures are caused by exponentials of multiples of − γt and f j ± ( t ) in Eq. (22), respectively. Purities for both thecomplete and the RWA models turned out to be indepen-dent on the initial coherent states of the modes. This isso because α and β can be eliminated from the dynam-ics via a time-independent unitary transformation on themodes (basis transformation). λ = ωλ = ωλ = ω j ( t ) ω t FIG. 3. (Color online) Purity P j of the two-qubit subsystemas function of the dimensionless time ω t for different valuesof mode coupling strengths λ . Solid blue lines represent thepurity for the complete model ( j = 1), whereas dashed redlines represent the purity under the RWA model ( j = 2). Weused ω = ω , g = 0 . ω , and γ = 5 × − ω . Let us now closely examine the dependence of thetwo-qubit purity on mode-mode coupling strength λ .From Fig. 3, we can see that, for moderate couplings( λ ≈ . ω ), the predictions of the RWA and non-RWAmodels already disagree considerably. Although theyboth have similar orders of magnitude, the oscillationsare not in phase anymore (compared to small λ ). This isa direct consequence of the deviations in Ω j ± caused es-sentially by the squeezing parameter r j ± . In addition, as λ becomes larger (approaching the limit ω/ j = 1(non-RWA) it is an oscillating function while for j = 2 itis constant and equals one. D. Entanglement of the qubits
We now turn our analysis to the two-qubit entan-glement generation in the studied setup. Some quan-tum information tasks, such as quantum teleportation[37], need the handling of large amounts of entangle-ment to be performed properly, so that it is essentialto determine how entangled a certain system is. In thiswork, we use the concept of entanglement of formationof an arbitrary two-qubit mixed state [38, 39]. First,one defines the so-called concurrence function C ( ρ ) =max { , √ (cid:15) − √ (cid:15) − √ (cid:15) − √ (cid:15) } , where (cid:15) i s are the eigen-values in decreasing order of ρσ y A σ y B ρ ∗ σ y A σ y B , ρ ∗ is thecomplex conjugate of ρ , and σ y the y -Pauli matrix. Then,the entanglement of formation of ρ can be defined as E F ( ρ ) = h (cid:32) (cid:112) − C ( ρ ) (cid:33) , (24)where h ( x ) = − x log x − (1 − x ) log (1 − x ). Both con-currence and entanglement of formation are equal to zero(unity) for a separable (maximally entangled) state.Figure 4 compares the dynamics of E F [ ρ AB j ( t )] for dif-ferent values of qubit dephasing rates γ and mode cou-pling strengths λ . As expected, E F starts from zero as theinitial state is separable. Then E F oscillates and reachesits maximum values for the conditions taken in the firstand third rows ( γ = 0). In this scenario, the state of thetwo-qubit system periodically changes from separable tomaximally entangled states in the absence of dephasing.In the second and forth rows, on the other hand, such os-cillations are damped due to the non-null dephasing rate( γ = 5 × − ω ), and the capability of producing higherpeaks of E F is enhanced in the strong coupling regime.Thus, the non-Markovian aspect of the dynamics acts asan entanglement generator, whereas the Markovian parttends to destroy it over long times.Another feature present in Fig. 4 is that the main fre-quency of oscillation of E F increases with λ , indepen-dently of the model. However, for the full model, thegeneration of entanglement is evidently faster than inthe RWA. This might be interesting for a scenario whereentanglement needs to be preserved in the presence ofstrong dephasing. Fast generation of entanglement hasrecently attracted interest of the community and has al-ready been proposed in other setups [40–42]. This can beobtained, for example, if two non-interacting qubits areweakly coupled to a common Ohmic bath [43, 44].By examining further the short-time behavior of E F (Fig. 5), one can see that it essentially oscillates with thesame fast frequencies of the two-qubit purity (Fig. 3).This is indeed expected since both quantities are indi-rectly related to the local entropies for each qubit sub-system. Differently from what is observed in Fig. 3,the deviations in E F for both models become evidenteven for small values of mode coupling strengths (e.g. λ = 0 . ω ). This is a consequence of the complexityof E F compared to P j . To evaluate the former it is nec-essary density matrix diagonalization and application of λ = ω λ = ωλ = ω ℰ F [ ρ A B j ( t ) ] ω t FIG. 4. (Color online) Long-time behavior of the entanglement of formation E F of the two-qubit subsystem as function of thedimensionless time ω t for different values of mode coupling strengths λ and qubit dephasing rates γ . Solid blue lines represent E F for the complete model ( j = 1), whereas dashed red lines represent E F under the RWA ( j = 2). First and second rows : λ = 0 . ω (left) and λ = 0 . ω (right) with γ = 0 (first row) and γ = 5 × − ω (second row). Third and forth rows : λ = 0 . ω with γ = 0 (third row) and γ = 5 × − ω (forth row). The remaining parameters are ω = ω , g = 0 . ω , and α = β = 2. logarithmic functions while for the latter it is just nec-essary to square it and trace. Also, since what justifiesthe RWA is precisely first order perturbation theory [45],valid for short times, we indeed expect that the strongerthe λ the shorter the time range for which RWA pro-vides a satisfactory answer, and this is clearly seen fromthe plots in Figs. 4 and 5. Therefore, one can concludethat the dynamics of entanglement is more sensible tothe variations of λ than purity is. IV. CONCLUSION
We have provided an illustrative example where theinadequacy of the RWA can be analytically investigated.The system is composed of quantum two-level systemsand harmonic oscillators, both ubiquitous in controlledquantum systems such as trapped ions or circuit QED.In particular, we have focused on a model comprised oftwo qubit-mode subsystems which are brought into in-teraction through their harmonic coordinates. We haveperformed exact diagonalization of the full model and itsRWA version, and analytically solved the master equa-tion for initial states of interest.We then have shown that the modes coupling strength λ plays a fundamental role on the variety of responses ofthe qubits as displayed by state purity and entanglement.The predictions of the complete model and the RWAmodel for short times and λ ≈ . ω agree well for puritybut not for entanglement dynamics. Also, even for shorttimes, as soon as λ ≈ . ω , purity is no longer described correctly by the RWA. At longer times, when two-qubitentanglement is fully generated, the failure of the RWAbecomes more noticeable as it considerably reduces themain frequency of entanglement oscillations. Our studyalso showed that, the stronger the modes are coupled, thelarger the reduction of the degree of purity at short timescale and the faster generation of entanglement. More-over, we verified the competition between Markovian de-phasing on the qubits and non-Markovianity induced bythe modes coupling constant in conformity with Ref. [29]. ACKNOWLEDGMENTS
P.C.C. would like to thank FAPESP for the supportthrough Grant No. 2012/12702-7. W.S.T would like tothank CAPES for the current scholarship. F.L.S. ac-knowledges participation as a member of the BrazilianNational Institute of Science and Technology of QuantumInformation (INCT/IQ). F.L.S. also acknowledges partialsupport from CNPq under Grant No. 307774/2014-7.
Appendix
The time-dependent factors γ mn j ( t ) that appear in thequbits density matrix elements ρ AB,mn j ( t ) in Eq. (19) are λ = ωλ = ωλ = ω ℰ F [ ρ A B j ( t ) ] ω t FIG. 5. (Color online) Short-time behavior of the entangle-ment of formation E F of the two-qubit subsystem as functionof the dimensionless time ω t for different values of modecoupling strengths λ . Solid blue lines represent E F for thecomplete model ( j = 1), whereas dashed red lines represent E F under the RWA ( j = 2). The remaining parameters are ω = ω , g = 0 . ω , γ = 5 × − ω and α = β = 2. explicitly γ eeeg j ( t ) = e − [ i ( ω ,j + χ j )+ γ ] t e i Im [ A j + ( t )+ B j − ( t ) ] ,γ eege j ( t ) = e − [ i ( ω ,j + χ j )+ γ ] t e i Im [ A j + ( t ) − B j + ( t ) ] ,γ eegg j ( t ) = e − iω ,j + γ ) t e i Im (cid:104) λ j + Z ∗ j (cid:16) − e i Ω j + t (cid:17)(cid:105) ,γ egge j ( t ) = e − γt e − i Im (cid:104) λ j − W ∗ j (cid:16) − e i Ω j − t (cid:17)(cid:105) ,γ eggg j ( t ) = e − [ i ( ω ,j − χ j )+ γ ] t e i Im [ A j − ( t ) − B j − ( t ) ] ,γ gegg j ( t ) = e − [ i ( ω ,j − χ j )+ γ ] t e i Im [ A j − ( t )+ B j + ( t ) ] , (A.1)where A j ± ( t ) = λ j + Z ∗ j − λ j + ( Z ∗ j ± λ j + ) e i Ω j + t and B j ± ( t ) = λ j − W ∗ j − λ j − ( W ∗ j ± λ j − ) e i Ω j − t are com-plex functions. Information about the complex ampli-tudes of the initial coherent state of the modes is en-coded in Z j = z j cosh( r j + ) + z ∗ j sinh( r j + ) and W j = w cosh( r j − ) + w ∗ sinh( r j − ), with z j = ( α + β + 2 δ j ) / √ w = ( β − α ) / √
2. Moreover, γ mm j = 1 and, giventhe Hermiticity of the density matrix, the remainingphases are complex conjugates of the ones in Eq. (A.1). Itis interesting to notice that the amplitudes of the coher-ent states appear in the density matrix but, as explainedbefore, do not show up in quantities that are independentof local time-independent transformations.From the partial trace over the bosonic modesin Eq. (16) using the total initial state | ψ (0) (cid:105) = | ++ (cid:105) | α (cid:105) | β (cid:105) , one finds the termsΘ mn,j ± ( t ) = (cid:10) Y m,j ± ( t ) , ξ j ± ( t ) (cid:12)(cid:12) Y n,j ± ( t ) , ξ j ± ( t ) (cid:11) , (A.2)which are scalar products of squeezed coherent stateswith squeezing parameter ξ j ± = − r j ± e − i Ω j ± t and am-plitudes Y ee,j + ( t ) = (cid:0) Z j + 2 λ j + (cid:1) e − i Ω j + t − λ j + ,Y eg,j + ( t ) = Z j e − i Ω j + t = Y ge,j + ( t ) ,Y gg,j + ( t ) = (cid:0) Z j − λ j + (cid:1) e − i Ω j + t + 2 λ j + ,Y ee,j − ( t ) = W j e − i Ω j − t = Y gg,j − ( t ) ,Y eg,j − ( t ) = (cid:0) W j − λ j − (cid:1) e − i Ω j − t + 2 λ j − ,Y ge,j − ( t ) = (cid:0) W j + 2 λ j − (cid:1) e − i Ω j − t − λ j − . (A.3)Notice that in the RWA, i.e., j = 2, Θ mn, ± ( t ) reducesto the overlap of coherent states (cid:10) Y m, ± ( t ) (cid:12)(cid:12) Y n, ± ( t ) (cid:11) . [1] C. Gerry and P. Knight, Introductory Quantum Optics (Cambridge University Press, Cambridge, U.K., 2005).[2] E. T. Jaynes and F. W. Cummings, Proc. IEEE , 89(1963).[3] J. H. Eberly, N. B. Narozhny, and J. J. Sanchez-Mondragon, Phys. Rev. Lett. , 1323 (1980).[4] P. L. Knight and P. M. Radmore, Phys. Rev. A , 676(1982).[5] G. D. Mahan, Many-Particle Physics , 2nd ed. (Plenum, New York, 1990).[6] U. Weiss,
Quantum Dissipative Systems , 3rd ed. (WorldScientifc, Singapore, New Jersey, London, Hong Kong,2008).[7] M. Schr¨oter,
Dissipative Exciton Dynamics in Light-Harvesting Complexes (Springer Spektrum, Heidelberg,Germany, 2015).[8] N. Lambert, Y.-N. Chen, Y.-C. Cheng, C.-M. Li, G.-Y.Chen, and F. Nori, Nat. Phys. , 10 (2013). [9] E. J. O’Reilly and A. Olaya-Castro, Nat. Commun. ,3012 (2014).[10] J. Gilmore and R. McKenzie, J. Phys.: Condens. Matter , 1735 (2005).[11] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A.Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. ,1 (1987).[12] K. Zaheer and M. S. Zubairy, Phys. Rev. A , 1628(1988).[13] M. H. Naderi, J. Phys. A , 055304 (2011).[14] Q.-H. Chen, Y. Yang, T. Liu, and K.-L. Wang, Phys.Rev. A , 052306 (2010).[15] J. Hausinger and M. Grifoni, Phys. Rev. A , 062320(2010).[16] J. Hausinger and M. Grifoni, New J. Phys. , 115015(2008).[17] G. S. Agarwal, Phys. Rev. A , 1778 (1971).[18] D. Zueco, G. M. Reuther, S. Kohler, and P. H¨anggi, Phys.Rev. A , 033846 (2009).[19] J.-Q. Liao, J.-F. Huang, and L. Tian, Phys. Rev. A ,033853 (2016).[20] A. T. Sornborger, A. N. Cleland, and M. R. Geller, Phys.Rev. A , 052315 (2004).[21] S. Agarwal, S. M. H. Rafsanjani, and J. H. Eberly, Phys.Rev. A , 043815 (2012).[22] D. Braak, Phys. Rev. Lett. , 100401 (2011).[23] J. Clarke and F. K. Wilhelm, Nature (London) , 1031(2008).[24] R. J. Schoelkopf and S. M. Girvin, Nature (London) ,664 (2008).[25] A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S.Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J.Schoelkopf, Nature (London) , 162 (2004).[26] A. Blais, R-S. Huang, A. Wallraff, S. M. Girvin, and R.J. Schoelkopf, Phys. Rev. A , 062320 (2004).[27] A. Blais, J. Gambetta, A. Wallraff, D. I. Schuster, S. M.Girvin, M. H. Devoret, and R. J. Schoelkopf, Phys. Rev.A , 032329 (2007).[28] D. I. Schuster, A. A. Houck, J. A. Schreier, A. Wallraff, J. M. Gambetta, A. Blais, L. Frunzio, J. Majer, B. John-son, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf,Nature (London) , 515 (2007).[29] P. C. C´ardenas, M. Paternostro, and, F. L. Semi˜ao, Phys.Rev. A , 022122 (2015).[30] M. Mariantoni, Frank Deppe, A. Marx, R. Gross, F. K.Wilhelm, and E. Solano, Phys. Rev. B , 104508 (2008).[31] M. J. Hartmann, F. G. S. L. Brando, and M. B. Plenio,Nat. Phys. , 849 (2006).[32] D. Vion, A. Aassime, A. Cottet, P. Joyez, H. Pothier,C. Urbina, D. Esteve, and M. Devoret, Science , 886(2002).[33] Y.-D. Wang, Y. Li, F. Xue, C. Bruder, and K. Semba,Phys. Rev. B , 144508 (2009).[34] A. Wurger, Phys. Rev. B , 347 (1998); A. Nazir, Phys.Rev. Lett. , 146404 (2009); C. K. Lee, J. Moix, andJ. Cao, J. Chem. Phys. , 204120 (2012).[35] U. Naether, J. J. Garc´ıa-Ripoll, J. J. Mazo, and D. Zueco,Phys. Rev. Lett. , 074101 (2014).[36] S. M. Barnett, Quantum Information (Oxford UniversityPress, New York, United States, 2009).[37] C. H. Bennett, G. Brassard, C. Cr´epeau, R. Jozsa, A.Peres, and W. K. Wootters, Phys. Rev. Lett. , 1895(1993).[38] W. K. Wootters, Phys. Rev. Lett. , 2245 (1998).[39] R. Horodecki, P. Horodecki, M. Horodecki, and K.Horodecki, Rev. Mod. Phys. , 865 (2009).[40] Z. H. Peng, Y. X. Liu, Y. Nakamura, and J. S. Tsai,Phys. Rev. B , 024537 (2012).[41] A. Posazhennikova, R. Birmuske, M. Bruderer, and W.Belzig, Phys. Rev. A , 042302 (2013).[42] N. Qiu and X.-B. Wang, Phys. Rev. A , 062332 (2013).[43] F. Benatti, R. Floreanini, and U. Marzolino, Europhys.Lett. , 20011 (2009).[44] F. Benatti, R. Floreanini, and U. Marzolino, Phys. Rev.A , 012105 (2010).[45] F. Nicacio and F. L. Semi˜ao, J. Phys. A49