Canonical versus noncanonical equilibration dynamics of open quantum systems
Chun-Jie Yang, Jun-Hong An, Hong-Gang Luo, Yading Li, C. H. Oh
aa r X i v : . [ qu a n t - ph ] A ug Canonical versus noncanonical equilibration dynamics of open quantum systems
Chun-Jie Yang,
1, 2
Jun-Hong An,
1, 3, ∗ Hong-Gang Luo,
1, 4
Yading Li, † and C. H. Oh ‡ Center for Interdisciplinary Studies, School of Physical Science and Technology, Lanzhou University, Lanzhou 730000, China Qian Xuesen Laboratory of Space Technology, China Academy of Space Technologyin China Aerospace Science and Technology Corporation, Beijing 100094, China Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543, Singapore Beijing Computational Science Research Center, Beijing 100084, China
In statistical mechanics, any quantum system in equilibrium with its weakly coupled reservoir isdescribed by a canonical state at the same temperature as the reservoir. Here, by studying the equi-libration dynamics of a harmonic oscillator interacting with a reservoir, we evaluate microscopicallythe condition under which the equilibration to a canonical state is valid. It is revealed that the non-Markovian effect and the availability of a stationary state of the total system play a profound role inthe equilibration. In the Markovian limit, the conventional canonical state can be recovered. In thenon-Markovian regime, when the stationary state is absent, the system equilibrates to a generalizedcanonical state at an effective temperature; whenever the stationary state is present, the equilibriumstate of the system cannot be described by any canonical state anymore. Our finding of the physicalcondition on such noncanonical equilibration might have significant impact on statistical physics. Aphysical scheme based on circuit QED is proposed to test our results.
PACS numbers: 03.65.Yz, 05.30.Rt, 05.30.Ch
I. INTRODUCTION
Recently, how quantum systems, open [1–8] or closed[9–12], thermalize microscopically from a certain initialstate to a canonical state has attracted much attention.In statistical mechanics, based on some basic postulates,such as the equal a priori probability postulate or theassumption of ergodicity, it can be proven kinematicallythat any quantum system in equilibrium with its inter-acting reservoir is described by a canonical state at thesame temperature as the reservoir [13]. However, thesebasic postulates rely on the subjective lack of knowledgeof the system. It was revealed that they could be aban-doned by examining the entanglement induced by the in-teraction between the system and the reservoir, by whichthe canonical state can be achieved without referring toensembles or time averages [1].Another way to relax these assumptions in statisticalmechanics from the perspective of the microscopic theoryis to study the dynamics of open systems [14–19]. It is ex-pected that the equilibrium behavior of any (many-body)system in statistical mechanics should be the steady-statelimit of the nonequilibrium dynamics. In this way, thestudy of the equilibration dynamics of an open systemcan give a bridge between an arbitrary initial state of thesystem and its statistical equilibrium state. It has beenproven that [14, 17], in the limit of the vanishing couplingand the memoryless correlation function of the reservoir,the steady state of the interested quantum system is acanonical state. However, in real physical situations, e.g.,photonic crystal [20] and semiconductor phonon [21] en- ∗ [email protected] † [email protected] ‡ [email protected] vironments, low-dimensional solid-state system [22], andeven purely optical systems [23, 24], these conditions aregenerally invalid and the strong memory effect of thereservoir may even cause the anomalous decoherence be-havior, i.e., the halting of decoherence with the increaseof the coupling between the system and the reservoir[8, 25, 26]. Thus it is desirable to know how and whena nonequilibrium dynamics results in a relaxation of thesystem to the equilibrium states universally with respectto widely differing initial conditions and compatible tostatistical mechanics.In this work, we are interested in two questions. (1)What state does the system dynamically evolve to? (2)What is the physical condition under which the canonicalstate in statistical mechanics is valid? To answer thesequestions, we study the exact equilibration dynamics of aharmonic oscillator interacting with a thermal reservoir.A generalized canonical state with an effective tempera-ture, which reduces to the conventional one in the vanish-ing limit of the coupling strength, and its validity condi-tion are obtained microscopically. The mechanism of itsbreakdown is due to the formation of a stationary stateof the whole system [28, 29]. Our results reveal the com-plete breakdown of the equilibrium statistical mechanicswhen the stationary state is present. The possible testifi-cation of our prediction in a coupled cavity array systemin circuit QED is also proposed.Our paper is organized as follows. In Sec. II, weshow the model of a harmonic oscillator interacting witha bosonic reservoir at finite temperature and its ex-act master equation derived by Feynman-Vernon’s influ-ence functional method. In Sec. III, the canonical andnoncanonical thermalization dynamics of the harmonic-oscillator system is revealed. In Sec. IV, we proposean experimentally realizable scheme to test our result onnoncanonical thermalization. Finally, a summary is givenin Sec. V. II. MODEL AND EXACT DYNAMICS
We study a harmonic oscillator coupled to a reservoir.The Hamiltonian of the total system is ˆ H = ˆ H s + ˆ H r + ˆ H i with (¯ h = 1)ˆ H s = ω ˆ a † ˆ a, ˆ H r = X k ω k ˆ b † k ˆ b k , ˆ H i = X k g k (ˆ a ˆ b † k + h.c.) , (1)where ˆ a and ˆ b k (ˆ a † and ˆ b † k ) are, respectively, the annihila-tion (creation) operators of the oscillator with frequency ω and the k th reservoir mode with frequency ω k , and g k is the coupling strength. g k can be characterized by thespectral density J ( ω ) ≡ X k | g k | δ ( ω − ω k ) = ηω ( ω/ω c ) s − e − ω/ω c , (2)where ω c is a cutoff frequency, and η is a dimensionlesscoupling constant. The reservoir is classified as Ohmic if s = 1, sub-Ohmic if 0 < s <
1, and super-Ohmic if s > a ˆ b k and ˆ a † ˆ b † k have been discarded. Note that this reduc-tion is more than a simple approximation to the quan-tum Brownian motion [2, 3] since our model describesan interaction conserving the total number of excitationsˆ N tot ≡ ˆ a † ˆ a + P k ˆ b † k ˆ b k , a symmetry that can be imposedby nature right from the start to describe relevant ex-periments. Explicitly, our system is highly pertinent toquantum-optical setting where the system oscillator candescribe the quantized fields in cavity [33] or in circuit[34] QED, mechanical oscillators in optomechanics [35],and atomic ensemble in the large- N limit [36].Setting ρ tot (0) = ρ (0) ⊗ e − β ˆ H r Tr r e − β ˆ H r with β = ( k B T ) − ,we can derive the exact master equation by Feynman-Vernon’s influence functional method [1, 2, 4, 39]˙ ρ ( t ) = − i Ω( t )[ˆ a † ˆ a, ρ ( t )]+[Γ( t ) + Γ β ( t )2 ][2ˆ aρ ( t )ˆ a † − ˆ a † ˆ aρ ( t ) − ρ ( t )ˆ a † ˆ a ]+ Γ β ( t )2 [2ˆ a † ρ ( t )ˆ a − ˆ a ˆ a † ρ ( t ) − ρ ( t )ˆ a ˆ a † ] , (3)where Ω( t ) = − Im[ ˙ u ( t ) /u ( t )] is the renormalized fre-quency, Γ( t ) = − Re[ ˙ u ( t ) /u ( t )] and Γ β ( t ) = ˙ v ( t ) +2 v ( t )Γ( t ) are the dissipation and noise coefficients. u ( t )and v ( t ) are determined by˙ u ( t ) + iω u ( t ) + Z t dt µ ( t − t ) u ( t ) = 0 , (4) v ( t ) = Z t dt Z t dt u ∗ ( t ) ν ( t − t ) u ( t ) , (5)under u (0) = 1. Physically, u ( t ), which is tempera-ture independent and induced purely by the quantum vacuum effect, plays a role of dissipation, while v ( t ),which is dependent on the temperature and induced bythermal effect, plays a role of fluctuation in the dy-namics. Equation (S28) can serve as a nonequilibriumversion of the fluctuation-dissipation relation, which en-sures the positivity of ρ ( t ) during the evolution. Thekernel functions µ ( x ) = R ∞ dωJ ( ω ) e − iωx and ν ( x ) = R ∞ dωJ ( ω )¯ n ( ω ) e − iωx with ¯ n ( ω ) = e βω − (see Supple-mental Material [40]). Equation (3) keeps the same Lind-blad form as the corresponding Markovian master equa-tion, but the coefficients are time dependent. All thenon-Markovian effect from the reservoir has been incor-porated into these coefficients self-consistently. It canrecover the Markovian one under the conditions that thesystem-reservoir coupling is very weak and the charac-teristic time scale of the reservoir correlation function ismuch smaller than the typical time scale of the system.Then the second-order perturbative solutions of Eqs. (4)and (S28) read u ( t ) = e − ( κ + iω ′ ) t and ˙ v ( t ) = 2 κ ¯ n ( ω )with κ = πJ ( ω ) and ω ′ = ω + P R J ( ω ) ω − ω dω . Thusthe coefficients reduce to Γ( t ) = κ , Ω( t ) = ω ′ , andΓ β ( t ) = 2 κ ¯ n ( ω ), which equal the ones in the Markovianmaster equation (see Supplemental Material [40]). III. THERMALIZATION DYNAMICS
With the exact master equation (3) at hand, the ther-malization dynamics of the system can be studied readily.According to the detailed balance condition, the steadystate for positive Γ( t ) and Γ β ( t ) can be constructed as(see Supplemental Material [40]) ρ ( ∞ ) = ∞ X n =0 [Γ β ( ∞ ) / (2Γ( ∞ ))] n [1 + Γ β ( ∞ ) / (2Γ( ∞ ))] n +1 | n ih n | , (6)which defines a unique canonical state for any initial state ρ con ≡ ρ ( ∞ ) = e − β eff ˆ H s / Tr s e − β eff ˆ H s , ( β eff = 1 /k B T eff ) , (7)with an effective temperature T eff = ω k B { ln[1+ ∞ )Γ β ( ∞ ) ] } − . T eff reduces exactly to the reservoir temperature T in theMarkovian limit. Thus the canonical ensemble assump-tion in statistical mechanics, i.e., a system in contactwith a reservoir is described by the canonical state withthe same temperature as the reservoir, can be dynam-ically confirmed only under the Markovian approxima-tion. In the non-Markovian case, the equilibrated tem-perature T eff shows a dramatic deviation from T .To get a qualitative understanding on T eff , we solve u ( t ) by Laplace transform and get u ( t ) = L − [ s + iω +˜ µ ( s ) ]with ˜ µ ( s ) = R ∞ J ( ω ) s + iω . Using Cauchy residue theorem,the inverse Laplace transform can be calculated by find-ing the poles of the integrand, i.e., y ( E ) ≡ ω − Z ∞ J ( ω ) ω − E dω = E, ( E = is ) . (8) H c L Η T e ff (cid:144) T FIG. 1. (Color online) Γ( t ) (a) and Γ β ( t ) (b) in different η .(c) T eff evaluated when t = 100 /ω . ω c = 8 . ω , β = 0 . /ω ,and s = 1 have been used. Note that the roots of Eq. (8) match well with the energyspectrum of Eq. (1) in the single-excitation subspace (seeSupplemental Material [40]). It is understandable basedon the fact that the vacuum-induced dissipation u ( t ) isdominated by the single-excitation process in Eq. (1).Since y ( E ) is a monotonically decreasing function in theregion E ∈ ( −∞ ,
0) for Eq. (2), Eq. (8) has only one neg-ative root if y (0) <
0. No further discrete root exists inthe region E ∈ (0 , + ∞ ) because it would make y ( E ) di-vergent. The obtained u ( t ) contributed from this discretenegative root, i.e. e − iEt , corresponds to a stationary-state solution to Eq. (4) [28, 29], which has a vanishingΓ( t ) and means a halt of decoherence [8, 25, 26]. By thiscriterion, the stationary state for Eq. (2) is formed when ω − ηω c γ ( s ) ≤ , (9)where γ ( s ) is Euler’s Γ function. The vanishing Γ( t )would result in the divergence of T eff . Hence, we arguethat the canonical equilibration would be invalid whenthe stationary state is formed. In this case, the dynam-ics will keep the superposition amplitude of the formedstationary state unchanged, which causes a notable con-tribution of the stationary state in the equilibrium state.It is the reason for the breakdown of the canonical equi-libration.To verify the validity of Eq. (7) and evaluate the con-sequence of the formed stationary state on the equilibra-tion of the system, we present numerical simulations onthe equilibration dynamics. For concreteness, we choosethe Ohmic spectral density, i.e., s = 1 in Eq. (2), to dothe numerics. Then we can obtain from the condition(9) that the stationary state is formed when ω ≤ ηω c .Figures 1 and 2 show the obtained Γ( t ), Γ β ( t ), and T eff indifferent parameter regimes. We can see that when thestationary state is absent, i.e., η < .
125 in Fig. 1 and ω c < ω in Fig. 2, both Γ( t ) and Γ β ( t ) tend to posi-tive constant values asymptotically. The positivity of the H c L Ω c (cid:144) Ω T e ff (cid:144) T FIG. 2. (Color online) Γ( t ) (a) and Γ β ( t ) (b) in different ω c .(c) T eff evaluated when t = 100 /ω . η = 0 . β = 0 . /ω , and s = 1 have been used. two coefficients guarantees the overall thermalization ofthe system during the dynamics. Consequently, the sys-tem approaches a canonical state governed by Eq. (7),as shown in Figs. 1(c) and 2(c), irrespective of in whatstate the system initially resides. Here it is qualitativelyconsistent with the one under the Born-Markovian ap-proximation. However, whenever the stationary state isformed, i.e., η ≥ .
125 in Fig. 1 and ω c ≥ ω in Fig. 2,Γ( t ) and Γ β ( t ) possess transient negative values and ap-proach zero asymptotically, which confirms our expecta-tion based on the stationary-state analysis. The negativecoefficients reveal a nonzero non-Markovianity [9], whichin turn gives an insight to the dynamical consequence ofthe formed stationary state (see Supplemental Material[40]). The vanishing of the two coefficients causes thedissipationless dynamics [25, 42], which is totally differ-ent from the Markovian case, and the ill definition of thecanonical state (7). In Figs. 1(c) and 2(c), we really seethat T eff changes to be divergent when the parameterstend to the critical point forming the stationary state.Another interesting observation is that even when theequilibration to a canonical state is valid in the absenceof the stationary state, T eff still shows a dramaticallyquantitative difference to T . They match well only when η is vanishingly small [see Fig. 1(c)], and when ω c iscomparable with ω [see Fig. 2(c)].We next investigate an actual example by employingan explicit initial state of the system harmonic oscilla-tor as ρ (0) = e −| α | | α ih α | . In the framework of theinfluence-functional method, we can readily calculate itstime evolution state and the expectation value of thebosonic number operator N ( t ) = | α u ( t ) | + v ( t ). Fig-ures 3(a) and 4(a) show the evolution of N ( t ) in the fullparameter regime. We can see clearly that N ( t ) behavesdivergently at the critical point forming the stationarystate, i.e., η = 0 .
125 in Fig. 3(a) and ω c = 10 ω inFig. 4(a). We plot in Figs. 3(b) and 4(b) the ratio N ( t )to N con ≡ Tr s [ˆ a † ˆ aρ con ] in the parameter regime where =0.01/ =0.05/ =0.1/ l n N () FIG. 3. (Color online) N ( t ) (a) and ˜ n ( t ) ≡ N ( t ) /N con (b) indifferent η . | α | = 1 . N ( ∞ ) and the corresponding Markovian valuesshown by dashed, dotted, and dot-dashed lines in different T . =0.01/ =0.05/ =0.1/ l n N () c / FIG. 4. (Color online) N ( t ) (a) and ˜ n ( t ) ≡ N ( t ) /N con (b) indifferent ω c . | α | = 1 . N ( ∞ ) and the corresponding Markovian valuesshown by dashed, dotted, and dot-dashed lines in different T . the stationary state is absent. We can see that the ra-tio in both of the results exclusively tends to 1 in thelong-time limit. It verifies the validity of Eq. (7) to de-scribe the equilibrium state of the system. We show inFigs. 3(c) and 4(c) the long-time values of bosonic num-ber N in different initial temperatures of the reservoir.We can see that, irrespective of the initial temperature, N ( ∞ ) diverges at the same critical point. It means thatthe condition for the formation of the stationary state isindependent of the initial temperature of the reservoir,which confirms our analytical criterion (9).To verify the uniqueness of the equilibrium state whenthe stationary state is absent, we plot in Fig. 5 the timeevolution of N ( t ) in different initial condition α . Figure5(a) indicates that, when the stationary state is absent, N ( t ) tends to the unique steady value N con . However,when the stationary state is formed, the steady value of FIG. 5. (Color online) N ( t ) in different α when η = 0 .
05 (a)and 0 . C − N L − N C c C c C L C c C c C N L N C s L s S (a) FIG. 6. (Color online) (a) Scheme on a TLR (formed by C s and L s ) with ω = 1 / (2 √ L s C s ) interacting with an ar-ray of TLRs (formed by C i and L i ) with identical ω C =1 / (2 √ L i C i ) coupled via C c in strength ξ = 2 Z C c ω C ( Z being the impedance). Γ( t ) (b) and Γ β ( t ) (c) in different ω . ξ = 0 . ω C , g = 0 . ω C , β = 0 . /ω C , and N = 200 havebeen used. N ( t ) is not unique and dependent on the initial condition α , which also proves the breakdown of the thermaliza-tion to a canonical state in this situation. This corre-sponds to a situation which cannot be described by theequilibrium statistical mechanics. The similar result canalso be achieved in different ω c regimes with definite η .The initial-state dependence of Fig. 5(b) also reveals theincorrectness of the result of Ref. [17], which is initial-state independent, in describing our system. IV. PHYSICAL REALIZATION
A promising physical system to test our predictionis an array of coupled cavities, which can now be re-alized experimentally in an optical-ring resonator sys-tem [43], in a microdisk cavity system [44], in a pho-tonic crystal system [45–47], and synthesized in opticalwaveguides [48, 49]. In Fig. 6, we depict a schemerealized in a circuit QED system [50–52]. Here, thesystem is realized by the transmission line resonator(TLR), and the reservoir is realized by an array of ca-pacitance coupled TLRs, which in turn couples to thesystem via inductance. The reservoir is governed byˆ H r = ω C P N/ j = − N/ ˆ b † j ˆ b j + ( ξ P N/ j = − N/ ˆ b † j +1 ˆ b j + h.c.) andthe interaction is ˆ H i = g ˆ a † ˆ b + h.c.. ˆ H r is recast intoˆ H r = P k ǫ k ˆ b † k ˆ b k by a Fourier transform ˆ b j = P k ˆ b k e ikjx with ǫ k = ω C + 2 ξ cos kx and x being the spatial sep-aration of the two neighbor resonators. One can noticethat the dispersion relation of the reservoirs shows fi-nite band width 4 ξ centered at ω C , which can inducea strong non-Markovian effect even in the weak- andintermediate-coupling regimes. The stationary state isformed if ω falls in the gap of the spectrum of the reser-voir, i.e., ω < ω C − ξ [26]. We really see in Figs. 6(b)and (c) that both Γ( t ) and Γ β ( t ) calculated based on anexperimentally accessible parameter regime [50, 51] ap-proach zero in this situation, where the equilibration toa canonical state is expected to be invalid. V. CONCLUSION
In summary, we have studied the non-Markovian equi-libration dynamics of a harmonic oscillator coupled to areservoir. A generalized canonical state is constructed microscopically. It recovers the conventional canonicalstate in the Markovian approximation, but gets to be in-valid whenever the stationary state is formed. Our workgives a bridge between the ensemble theory of the equi-librium statistical mechanics and the non-equilibrium dy-namics, thus builds a unified framework to treat the is-sues in equilibrium and nonequilibrium statistical me-chanics. In particular, the explicit criterion for the for-mation of the stationary state may supply a guideline toexperiments on exploring the noncanonical equilibrationof open system via engineering the specific structure ofthe reservoir [53, 54].
ACKNOWLEDGEMENTS
This work is supported by the Fundamental ResearchFunds for the Central Universities, by the SpecializedResearch Fund for the Doctoral Program of Higher Edu-cation, by the Program for NCET, by the NSF of China(Grant No. 11175072, No. 11174115, No. 11325417, andNo. 61271145), and by the National Research Founda-tion and Ministry of Education, Singapore (Grant No.WBS: R-710-000-008-271).
Note added : Recently, we became aware of two relatedworks [55, 56]. [1] S. Popescu, A. J. Short, and A. Winter, Nat. Phys. ,754 (2006).[2] N. Linden, S. Popescu, A. J. Short, and A. Winter, Phys.Rev. E. , 061103 (2009).[3] H. Tasaki, Phys. Rev. Lett. , 1373 (1998).[4] S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zangh`ı,Phys. Rev. Lett. , 050403 (2006).[5] P. Reimann, New J. Phys. , 055027(2010).[6] O. Lychkovskiy, Phys. Rev. E , 011123 (2010).[7] C. K. Lee, J. Cao, and J. Gong, Phys. Rev. E. , 021109(2012).[8] S. Genway, A. F. Ho, and D. K. K. Lee, Phys. Rev. Lett. , 130408 (2013).[9] M. Rigol, V. Dunjko, and M. Olshanii, Nature (London) , 854 (2008).[10] A. V. Ponomarev, S. Denisov, and P. H¨anggi, Phys. Rev.Lett. , 010405 (2011).[11] M. C. Ba˜nuls, J. I. Cirac, and M. B. Hastings, Phys. Rev.Lett. , 050405 (2011).[12] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat-tore, Rev. Mod. Phys. 83, 863 (2011).[13] L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, London, 1958).[14] E. Geva, E. Rosenman, and D. Tannor, J. Chem. Phys. , 1380 (2000).[15] M. Esposito and P. Gaspard, Phys. Rev. E. , 066113(2003).[16] T. Mori and S. Miyashita, J. Phys. Soc. Jpn. , 124005(2008).[17] Y. Suba¸si, C. H. Fleming, J. M. Taylor, and B. L. Hu,Phys. Rev. E , 061132 (2012). [18] D. Pagel, A. Alvermann, and H. Fehske, Phys. Rev. E , 012127 (2013).[19] A. Rancon and J. Bonart, Eur. Phys. Lett. , 50010(2013).[20] U. Hoeppe, C. Wolff, J. K¨uchenmeister, J. Niegemann,M. Drescher, H. Benner, and K. Busch, Phys. Rev. Lett. , 043603 (2012).[21] H. Tahara, Y. Ogawa, and F. Minami, Phys. Rev. Lett. , 037402 (2011).[22] C. Galland, A. H¨ogele, H. E. T¨ureci, and A. Imamoˇglu,Phys. Rev. Lett. , 067402 (2008).[23] B.-H. Liu, L. Li, Y.-F. Huang, C.-F. Li, G.-C. Guo, E.-M. Laine, H.-P. Breuer, and J. Piilo, Nat. Phys. , 931(2011).[24] K. H. Madsen, S. Ates, T. Lund-Hansen, A. L¨offler, S.Reitzenstein, A. Forchel, and P. Lodahl, Phys. Rev. Lett. , 233601 (2011).[25] W.-M. Zhang, P.-Y. Lo, H.-N. Xiong, Matisse Wei-YuanTu, and F. Nori, Phys. Rev. Lett. , 170402 (2012).[26] Y.-Q. L¨u, J.-H. An, X.-M. Chen, H.-G. Luo, and C. H.Oh, Phys. Rev. A , 012129 (2013).[27] H.-B. Liu, J.-H. An, C. Chen, Q.-J. Tong, H.-G. Luo,and C. H. Oh, Phys. Rev. A , 052139 (2013).[28] M. Miyamoto, Phys. Rev. A 72, 063405 (2005).[29] F. Dreisow, A. Szameit, M. Heinrich, T. Pertsch, S.Nolte, A. T¨unnermann, and S. Longhi, Phys. Rev. Lett. , 143602 (2008).[30] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A.Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. ,1 (1987).[31] A. Caldera and A. Leggett, Physica A , 587 (1983). [32] B. L. Hu, J. P. Paz, and Y. Zhang, Phys. Rev. D. ,2843 (1992).[33] X. Zhou, I. Dotsenko, B. Peaudecerf, T. Rybarczyk, C.Sayrin, S. Gleyzes, J. M. Raimond, M. Brune, and S.Haroche, Phys. Rev. Lett. , 243602 (2012).[34] J. M. Fink, M. G¨oppl, M. Baur, R. Bianchetti, P. J.Leek, A. Blais, and A. Wallraff, Nature (London) ,315 (2008).[35] T. J. Kippenberg and K. J. Vahala, Opt. Express ,17172 (2007).[36] K. Hammerer, M. Aspelmeyer, E. S. Polzik, and P. Zoller,Phys. Rev. Lett. , 020501 (2009).[37] R. P. Feynman and F. L. Vernon, Ann. Phys. (NY) ,118 (1963).[38] J.-H. An and W.-M. Zhang, Phys. Rev. A , 042127(2007).[39] H.-N. Xiong, W.-M. Zhang, X. Wang, and M.-H. Wu,Phys. Rev. A , 012105 (2010).[40] See the Supplemental Material for a detailed derivation ofthe exact master equation, the reduction to its Markovianlimit, the energy spectrum of the total system in thesingle-excitation subspace, the derivation of Eq. (6), andrelationship between the formed stationary state and thenon-Markovianity.[41] H.-P. Breuer, E.-M. Laine, and J. Piilo, Phys. Rev. Lett. , 210401 (2009).[42] Q.-J. Tong, J.-H. An, H.-G. Luo, and C. H. Oh, Phys.Rev. A , 052330 (2010).[43] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor,Nat. Photon. , 1001 (2013).[44] P. E. Barclay, K. Srinivasan, O. Painter, B. Lev, and H.Mabuchi, Appl. Phys. Lett. , 131108 (2006).[45] A. Majumdar, A. Rundquist, M. Bajcsy, V. D. Dasika,S. R. Bank, and J. Vuˇckovi´c, Phys. Rev. B , 195312(2012).[46] K. Hennessy, A. Badolato, M. Winger, D. Gerace, M.Atature, S. Gulde, S. F¨alt, E. L. Hu, and A. Imamoˇglu,Nature (London) , 896 (2007).[47] M. Notomi, E. Kuramochi, and T. Tanabe, Nat. Photon. , 741 (2008).[48] M. Verbin, O. Zilberberg, Y. E. Kraus, Y. Lahini, andY. Silberberg, Phys. Rev. Lett. , 076403 (2013).[49] M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer,D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Sza-meit, Nature (London) 496, 196 (2013).[50] D. L. Underwood, W. E. Shanks, J. Koch, and A. A.Houck, Phys. Rev. A , 023837 (2012).[51] A. A. Houck, H. E. Tureci, and J. Koch, Nat. Phys. 8,292 (2012).[52] D. J. Egger and F. K. Wilhelm, Phys. Rev. Lett. ,163601 (2013).[53] J. T. Barreiro, M. M¨uller, P. Schindler, D. Nigg, T. Monz,M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R.Blatt, Nature (London) , 486 (2011).[54] K. W. Murch, U. Vool, D. Zhou, S. J. Weber, S. M.Girvin, and I. Siddiqi, Phys. Rev. Lett. , 183602(2012).[55] C. Y. Cai, L.-P. Yang, and C. P. Sun, Phys. Rev. A ,012128 (2014).[56] H.-N. Xiong, P.-Y. Lo, W.-M. Zhang, F. Nori, and D. H.Feng, arXiv:1311.1282. SUPPLEMENTAL MATERIAL
S-1. THE EXACT DECOHERENCE DYNAMICS
In this section, we give the detailed derivation of ex-act master equation using Feynman-Vernon’s influence-functional theory [1–4] in the coherent-state representa-tion [5].
A. Influence functional method
The Hamiltonian of the total system reads ˆ H = ˆ H s +ˆ H r + ˆ H i with (¯ h = 1)ˆ H s = ω ˆ a † ˆ a, ˆ H r = X k ω k ˆ b † k ˆ b k , ˆ H i = X k g k (ˆ a ˆ b † k + h.c.) , (S1)The total density matrix obeys i ˙ ρ tot ( t ) = [ ˆ H, ρ tot ( t )], wehave the formal solution as ρ tot ( t ) = e − i ˆ Ht ρ tot (0) e i ˆ Ht . (S2)In the coherent-state representation, Eq. (S2) is ex-panded as h ¯ α f , ¯z f | ρ tot ( t ) | α ′ f , z f i = Z dµ ( z i ) dµ ( α i ) dµ ( z ′ i ) dµ ( α ′ i ) ×h ¯ α f , ¯z f ; t | α i , z i ; 0 ih ¯ α i , ¯z i | ρ tot (0) | α ′ i , z ′ i i×h ¯ α ′ i , ¯z ′ i ; 0 | α ′ f , z f ; t i , (S3)We assume that the initial density matrix factors into asystem part and an environment part ρ tot (0) = ρ (0) ⊗ ρ r (0). The dynamics of the system can be obtained byintegrating over the environmental variables, ρ (¯ α f , α ′ f ; t ) = Z dµ ( α i ) dµ ( α ′ i ) J (¯ α f , α ′ f ; t | ¯ α ′ i , α i ; 0) × ρ (¯ α i , α ′ i ; 0) , (S4)where ρ (¯ α f , α ′ f ; t ) ≡ R dµ ( z f ) h ¯ α f , ¯z f | ρ tot ( t ) | α ′ f , z f i .The effective propagation functional of the density ma-trix is defined as J (¯ α f , α ′ f ; t | ¯ α ′ i , α i ; 0) = R dµ ( z f ) dµ ( z i ) dµ ( z ′ i ) × h ¯ α f , ¯z f ; t | α i , z i ; 0 i ρ r ( ¯z i , z ′ i ; 0) h ¯ α ′ i , ¯z ′ i ; 0 | α ′ f , z f ; t i . (S5)Eq. (S5) contains two propagators of the total sys-tem: the forward and backward propagators, whichcan be expressed as path integrals. To evaluate theforward propagator operator e − i ˆ Ht between the initial( | α i , z i i ) and the final ( h ¯ α f , ¯z f | ) coherent states, one gen-erally divides the time interval t f − t i into N subinter-vals. Then by inserting the resolution of identity, i.e. R dµ ( α ) dµ ( z ) | α, z ih ¯ α, ¯z | = 1 with the integration mea-sures dµ ( α ) = e − ¯ αα d ¯ αdα πi and dµ ( z ) = Q k e − ¯ z k z k d ¯ z k dz k πi ,( N −
1) times between each subintervals and taking the limit of large N , the path integral representation of thepropagator can be obtained h ¯ α f , ¯z f ; t | α i , z i ; 0 i = Z D z D α exp( iS s [¯ α, α ]+ iS i [ ¯z , z , ¯ α, α ] + iS r [ ¯z , z ]) , (S6)where iS s [¯ α, α ] = ¯ α ( t ) α ( t ) + ¯ α (0) α (0)2 − Z t dτ [ ¯ α ( τ ) ˙ α ( τ ) − ˙¯ α ( τ ) α ( τ )2 + iH s (¯ α, α )] , (S7) iS r [ ¯z , z ] = X k { ¯ z k ( t ) z k ( t ) + ¯ z k (0) z k (0)2 − Z t dτ [ ¯ z k ( τ ) ˙ z k ( τ ) − ˙¯ z k ( τ ) z k ( τ )2 + iH r ( ¯z , z )] , (S8) iS i [ ¯z , z , ¯ α, α ] = − i Z t dτ H i (¯ α, α, ¯z , z ) . (S9)are the (complex) actions corresponding to the system,reservoir, and interaction Hamiltonian ˆ H s , ˆ H r , and ˆ H i ,respectively [6]. All the functional integrations are evalu-ated over paths ¯z ( τ ), z ( τ ), ¯ α ( τ ), and α ( τ ) with endpoints ¯z ( t ) = ¯z f , z (0) = z i , ¯ α ( t ) = α f , and α (0) = α i . Substi-tuting Eq. (S6) and a similar expression for the backwardpropagator into Eq. (S5), we obtain J (¯ α f , α ′ f ; t | ¯ α ′ i , α i ; 0) = Z D αD α ′ exp { i ( S s [¯ α, α ] − S ∗ s [¯ α ′ , α ′ ]) }F [¯ α, α, ¯ α ′ , α ′ ] , (S10)where F [¯ α, α, ¯ α ′ , α ′ ] = Z dµ ( z f ) dµ ( z i ) dµ ( z ′ i ) D z D z ′ × exp { i ( S r [ ¯z , z ] + S i [ ¯z , z , ¯ α, α ] − S ∗ r [ ¯z ′ , z ′ ] − S ∗ i [ ¯z ′ , z ′ , ¯ α ′ , α ′ ]) } ρ r ( ¯z i , z ′ i ; 0) (S11)is the influence functional containing all the environmen-tal effects on the system. B. Evaluation of the influence functional
The path integral of the environmental part in Eq.(S11) can be evaluated by the saddle point method un-der the boundary conditions z k (0) = z ki , ¯ z k ( t ) = ¯ z kf .We have the equations of motion for the stationary pathas ˙ z k ( τ ) + iω k z k ( τ ) = − ig ∗ k α ( τ ) , (S12)˙¯ z k ( τ ) − iω k ¯ z k ( τ ) = ig k ¯ α ( τ ) , (S13)where the paths of ¯ α , α are taken as external sources.Substituting the formal solutions of Eqs. (S12) and (S13) z k ( τ ) = z ki e − iω k τ − ig ∗ k R τ dt ′ e − iω k ( τ − t ′ ) α ( t ′ )(S14)¯ z k ( τ ) = ¯ z kf e − iω k ( t − τ ) − ig k R tτ dt ′ e iω k ( τ − t ′ ) ¯ α ( t ′ )(S15)into Eq. (S6), one can obtain the result of forward prop-agator. It is noted that the prefactor under the contribu-tion of the stationary path in the coherent-state pathintegral is equal to one and the saddle point methodto the evaluation of the environmental part here is ex-act. The similar expression for the backward propagator h ¯ α ′ i , ¯z ′ i ; 0 | α ′ f , z f ; t i can also be determined. Assuming the initial state of the reservoir is ρ r (0) = e − β ˆ H r Tr r e − β ˆ H r , we have ρ r ( ¯z i , z ′ i ; 0) = Y k (1 − e − β ¯ hω k ) exp( − | z ′ ik | + | ¯ z ik | e − β ¯ hω k ¯ z ik z ′ ik ) . (S16)Substituting this initial state, the forward and backwardpropagators into Eq. (S11), and using the Gaussian inte-gral identity R d zπ e − γ ¯ zz + δz f (¯ z ) = γ f ( δγ ) repeatedly, theinfluence functional can be evaluated as F [ α, α, α ′ , α ′ ] = exp { Z t dτ Z t dτ ′ [ ν ( τ − τ ′ )(¯ α ′ ( τ ) − ¯ α ( τ ))( α ( τ ′ ) − α ′ ( τ ′ )) + µ ( τ − τ ′ )¯ α ′ ( τ ) α ( τ ′ )] − Z t dτ Z τ dτ ′ [ µ ( τ − τ ′ )¯ α ( τ ) α ( τ ′ ) + µ ∗ ( τ − τ ′ )¯ α ′ ( τ ′ ) α ′ ( τ )] } , (S17)where ν ( τ − τ ′ ) = P k | g k | ¯ n ( ω k ) e − iω k ( τ − τ ′ ) and µ ( τ − τ ′ ) = P k | g k | e − iω k ( τ − τ ′ ) . C. Evaluation of the effective propagationfunctional
With Eq. (S17), Eq. (S10) can be recast into J (¯ α f , α ′ f ; t | ¯ α ′ i , α i ; 0) = Z Dµ ( α ) Dµ ( α ′ ) exp {
12 [¯ α ( t ) α ( t ) + ¯ α (0) α (0) + ¯ α ′ ( t ) α ′ ( t ) + ¯ α ′ (0) α ′ (0)] − A [¯ α, α ; ¯ α ′ , α ′ ] } A [¯ α, α ; ¯ α ′ , α ′ ] = Z t dτ { ¯ α ( τ ) ˙ α ( τ ) − ˙¯ α ( τ ) α ( τ )2 + ˙¯ α ′ ( τ ) α ′ ( τ ) − ¯ α ′ ( τ ) ˙ α ′ ( τ )2 + iω [¯ α ( τ ) α ( τ ) − ¯ α ′ ( τ ) α ′ ( τ )] − Z t dτ ′ [ ν ( τ − τ ′ )(¯ α ′ ( τ ) − ¯ α ( τ ))( α ( τ ′ ) − α ′ ( τ ′ )) + µ ( τ − τ ′ )¯ α ′ ( τ ) α ( τ ′ )]+ Z τ dτ ′ [ µ ( τ − τ ′ )¯ α ( τ ) α ( τ ′ ) + µ ∗ ( τ − τ ′ )¯ α ′ ( τ ′ ) α ′ ( τ )] } (S18)According to the main idea of saddle point method to evaluate the path integral, we can perform the path integralby expanding the effective action A around the saddle point paths or the classical paths A [¯ α c + δ ¯ α, α c + δα ; ¯ α ′ c + δ ¯ α ′ c , α ′ c + δα ′ ] = A [¯ α c , α c ; ¯ α ′ c , α ′ c ] + δA + δ A. (S19)By requesting δA = 0, we obtain the equations of motion satisfied by the variables in the classical paths˙ α c ( τ ) + iω α c ( τ ) + Z t dτ ′ ν ( τ − τ ′ )( α c ( τ ′ ) − α ′ c ( τ ′ )) + Z τ dτ ′ µ ( τ − τ ′ ) α c ( τ ′ ) = 0 , (S20)˙¯ α ′ c ( τ ) − iω ¯ α ′ c ( τ ) + Z t dτ ′ ν ∗ ( τ − τ ′ )(¯ α ′ c ( τ ′ ) − ¯ α c ( τ ′ )) + Z τ dτ ′ µ ∗ ( τ − τ ′ )¯ α ′ c ( τ ′ ) = 0 , (S21)˙¯ α c ( τ ) − iω ¯ α c ( τ ) + Z t dτ ′ ν ∗ ( τ − τ ′ )(¯ α ′ c ( τ ′ ) − ¯ α c ( τ ′ )) + Z t dτ ′ µ ∗ ( τ − τ ′ )¯ α ′ c ( τ ′ ) − Z tτ dτ ′ µ ∗ ( τ − τ ′ )¯ α c ( τ ′ ) = 0 , (S22)˙ α ′ c ( τ ) + iω α ′ ( τ ) + Z t dτ ′ ν ( τ − τ ′ )( α c ( τ ′ ) − α ′ c ( τ ′ )) + Z t dτ ′ µ ( τ − τ ′ ) α c ( τ ′ ) − Z tτ dτ ′ µ ( τ − τ ′ ) α c ( τ ′ ) = 0 . (S23)We also find that the second-order term δ A is inde-pendent on the variables in the classical paths, which means that it contributes only a time-dependent con-stant to the path integral. Our final task is to evaluate A [¯ α c , α c ; ¯ α ′ c , α ′ c ] using Eqs. (S20, S21, S22, S23). It isinteresting to find that A [¯ α c , α c ; ¯ α ′ c , α ′ c ] = 0 . (S24)Accordingly, we finally have J (¯ α f , α ′ f ; t | ¯ α ′ i , α i ; 0) = M ( t ) exp {
12 [¯ α f α ( t ) + ¯ α (0) α i + ¯ α ′ ( t ) α ′ f + ¯ α ′ i α ′ (0)] } , (S25)where M ( t ) contributed from the second-order term canbe determined by normalization.
1. The form of α ( t ) and α ′ (0) Defining χ ( τ ) = α ( τ ) − α ′ ( τ ) ≡ u ( τ ) χ ( t ) with u ( t ) =1, we can obtain from Eqs. (S20) and (S23)˙ u ( τ ) + iω u ( τ ) − Z tτ dτ ′ µ ( τ − τ ′ ) u ( τ ′ ) = 0 . (S26)We can see from Eq. (S20) that the initial condition of α ( τ ) can be fixed by α ( τ ) = u ( τ ) α i − v ( τ )[ α ( t ) − α ′ f ] with u (0) = 1 and v (0) = 0, and˙ u ( τ ) + iω u ( τ ) + Z τ µ ( τ − τ ′ ) u ( τ ′ ) dτ ′ = 0 , (S27)˙ v ( τ ) + iω v ( τ ) − Z t ν ( τ − τ ′ ) u ( τ ′ ) dτ ′ + Z τ µ ( τ − τ ′ ) v ( τ ′ ) dτ ′ = 0 . (S28)Comparing Eq. (S27) with Eq. (S26), we have u ( τ ) = u ∗ ( t − τ ). Thus we have α ′ ( τ ) = u ( τ ) α i − v ( τ )[ α ( t ) − α ′ f ] − u ( τ )[ α ( t ) − α ′ f ], which induces α ′ (0) = α i − u ∗ ( t )[ α ( t ) − α ′ f ] . (S29)Similarly, α ( t ) = u ( t ) α i − v ( t )[ α ( t ) − α ′ f ], which induces α ( t ) = u ( t ) α i + v ( t ) α ′ f v ( t ) . (S30)Substituing Eq. (S30) into Eq. (S29), we get α ′ (0) = [1 − | u ( t ) | v ( t ) ] α i + u ∗ ( t )1 + v ( t ) α ′ f . (S31)
2. The form of ¯ α (0) and ¯ α ′ ( t ) Defining ς ( τ ) = ¯ α ′ ( τ ) − ¯ α ( τ ) ≡ u ′ ( τ ) ς ( t ) with u ′ ( t ) =1, we can obtain from Eqs. (S21) and (S22)˙ u ′ ( τ ) − iω u ′ ( τ ) − Z tτ dτ ′ µ ∗ ( τ − τ ′ ) u ′ ( τ ′ ) = 0 . (S32) One can see that u ′ ( τ ) = u ∗ ( τ ) = u ( t − τ ). We can seefrom Eq. (S21) that the initial condition of ¯ α ′ ( τ ) can befixed by ¯ α ′ ( τ ) = u ′ ( τ )¯ α ′ i − v ′ ( τ )[¯ α ′ ( t ) − ¯ α f ] with u ′ (0) = 1and v ′ (0) = 0, and˙ u ′ ( τ ) − iω u ′ ( τ ) + Z τ µ ∗ ( τ − τ ′ ) u ′ ( τ ′ ) dτ ′ = 0 , (S33)˙ v ′ ( τ ) − iω v ′ ( τ ) − Z t ν ∗ ( τ − τ ′ ) v ′ ( τ ′ ) dτ ′ + Z τ µ ∗ ( τ − τ ′ ) v ′ ( τ ′ ) dτ ′ = 0 . (S34)Comparing Eqs. (S33) and (S34) with Eqs. (S27) and(S28), we have u ′ ( τ ) = u ∗ ( τ ), v ′ ( τ ) = v ∗ ( τ ). So ¯ α ′ ( τ ) = u ∗ ( τ )¯ α ′ i − v ∗ ( τ )[¯ α ′ ( t ) − ¯ α f ], which induces¯ α ′ ( t ) = u ∗ ( t )¯ α ′ i + v ∗ ( t )¯ α f v ∗ ( t ) . (S35)Similarly, ¯ α ( τ ) = u ∗ ( τ )¯ α ′ i − [ v ∗ ( t ) + u ( t − τ )][¯ α ′ ( t ) − ¯ α f ],which induces¯ α (0) = [1 − | u ( t ) | v ∗ ( t ) ]¯ α ′ i + u ( t )1 + v ∗ ( t ) ¯ α f . (S36)Substituting Eqs. (S30, S31, S35, S36) into Eq. (S25),we get J (¯ α f , α ′ f ; t | ¯ α ′ i , α i ; 0) = M ( t ) exp[ J ( t )¯ α f α i − J ( t ) ¯ α f α ′ f − J ( t ) ¯ α ′ i α i + J ∗ ( t )¯ α ′ i α ′ f ] , (S37)where M ( t ) = 11 + v ( t ) , J ( t ) = u ( t )1 + v ( t ) , (S38) J ( t ) = − v ( t )1 + v ( t ) , J ( t ) = | u ( t ) | v ( t ) − . (S39) D. Exact master equation
Eq. (S4) with the final form of the effective propoga-tion functional given by Eq. (S37) is the exact solution ofthe decoherence dynamics of the system. On one hand,the evolution of any initial state can be calculated by exe-cuting the integration in Eq. (S4). On the other hand, anexact master equation can be obtained by making timederivative to Eq. (S4). In the following, we derive themaster equation. First, we have the following functionalidentities from Eq. (S37) α i J = 1 J ( t ) [ ∂∂ ¯ α f + J ( t ) α ′ f ] J , (S40)¯ α ′ i J = 1 J ∗ ( t ) [ ∂∂α ′ f + J ( t ) ¯ α f ] J , (S41)which are useful to eliminating the dependence of thetime derivative of the effective propagation functional onthe initial variable α i and ¯ α ′ i . Then making the timederivative to Eq. (S37) and the substitution of Eqs. (S40,S41) and remembering Eq. (S4), we obtain the evolutionequation0˙ ρ (¯ α f , α ′ f ; t ) = n − Γ β ( t ) − (cid:2) Γ β ( t ) + Γ ( t ) + i Ω ( t ) (cid:3) ¯ α f ∂∂ ¯ α f + Γ β ( t ) ¯ α f α ′ f + [Γ β ( t ) + 2Γ ( t )] ∂ ∂ ¯ α f ∂α ′ f − (cid:2) Γ β ( t ) + Γ ( t ) − i Ω ( t ) (cid:3) α ′ f ∂∂α ′ f o ρ (¯ α f , α ′ f ; t ) (S42)where Γ ( t ) + i Ω ( t ) = − ˙ u ( t ) u ( t ) , Γ β ( t ) = ˙ v ( t ) + 2 v ( t )Γ( t ) . (S43)To obtain the operator equation of the reduced density matrix, it will be useful to introduce the following functionaldifferential relations. It can be proven the following identities for arbitrary system states in the coherent-staterepresentation¯ α f ∂ρ (¯ α f , α ′ f ; t ) ∂ ¯ α f ↔ ˆ a † ˆ aρ ( t ) , ¯ α f α ′ f ρ (¯ α f , α ′ f ; t ) ↔ ˆ a † ρ ( t )ˆ a, ∂ ρ (¯ α f , α ′ f ; t ) ∂ ¯ α f ∂α ′ f ↔ ˆ aρ ( t )ˆ a † , α ′ f ∂ρ (¯ α f , α ′ f ; t ) ∂α ′ f ↔ ρ ( t )ˆ a † ˆ a, (S44)with which we arrive at our final operator form of the non-Markovian master equation˙ ρ ( t ) = − i Ω ( t ) [ˆ a † ˆ a, ρ ( t )]+Γ ( t ) (cid:2) aρ ( t )ˆ a † − ˆ a † ˆ aρ ( t ) − ρ ( t )ˆ a † ˆ a (cid:3) +Γ β ( t ) (cid:2) ˆ a † ρ ( t )ˆ a + ˆ aρ ( t )ˆ a † − ˆ a † ˆ aρ ( t ) − ρ ( t )ˆ a ˆ a † (cid:3) . (S45) E. The evolved state of the coherent state as theinitial state
Consider explicitly that the system is initially in a co-herent state e −| α | / | α i , we have ρ (¯ α i , α ′ i ; 0) = exp( − | α | + ¯ α i α + ¯ α α ′ i ) . (S46)Substituting (S46) into Eq. (S37) and performing theGaussian integration, we can get ρ (¯ α f , α ′ f ; t ) = M ( t ) exp[ − (1 + J ( t )) | α | + J ( t )¯ α f α − J ( t )¯ α f α ′ f + J ∗ ( t )¯ α α ′ f ] . (S47)Remembering that the original density matrix is ex-pressed as ρ ( t ) = R dµ ( α f ) dµ ( α ′ f ) ρ (¯ α f , α ′ f ; t ) | α f ih ¯ α ′ f | ,we can readily calculate the expectation value of thebosonic number operator as N ( t ) ≡ Tr s [ˆ a † ˆ aρ ( t )] = | α u ( t ) | + v ( t ) . (S48) S-2. THE PERTURBATIVE SOLUTIONS OF u ( t ) AND v ( t ) In this section, we give the proof of the recovery of ex-act master equation to the conventional Markovian mas-ter equation under the second-order perturbation to thecoupling strength between the system and the reservoir.The solution of Eq. (S27) can be expanded as the orderof the coupling strength labeled by λ , u ( t ) = u (0) ( t ) + λ u (2) ( t ) + O ( | g k | ) , (S49)substituting which into Eq. (S27) we can get u (0) ( t ) = e − iω t and˙ u (2) ( t ) = − iω u (2) ( t ) − Z t µ ( t − t ) u (0) ( t ) dt . If the characteristic time scale of the environmentalcorrelation µ ( t − t ) is much smaller than the typicaltime scale t of the system, then we can extend the up-per limit of the integration into + ∞ . Using the iden-tity R + ∞ dt e − i ( ω − ω ) t = πδ ( ω − ω ) + i P ω − ω , we get u (2) ( t ) = − ( κ + i ∆ ω ) te − iω t , where κ = πJ ( ω ) and∆ ω = P R J ( ω ) ω − ω dω . Going back to Eq. (S49), we get u ( t ) = e − ( κ + iω ′ ) t , ω ′ = ω + ∆ ω, (S50)where λ has been set back to one.The second-order approximate solution of Eq. (S28)reads v ( t ) = Z t dt Z t dt u ∗ (0) ( t ) ν ( t − t ) u (0) ( t )= Z t dt Z t dt e − i ( ω − ω )( t − t ) ν ( t − t ) . (S51)Since we are only concerned about its time derivative, so˙ v ( t ) = Z ∞ dω Z t dt [ e − i ( ω − ω )( t − t ) + e i ( ω − ω )( t − t ) ] × J ( ω )¯ n ( ω ) ≃ κ ¯ n ( ω ) , (S52)where we have also used the extension of the upper limitof the integration to + ∞ . Substituting Eqs. (S50, S52)into Eq. (S43), we obtain the second-order perturbativevalues of of the parameters in the master equation asΓ( t ) = κ, Ω( t ) = ω ′ β ( t ) = 2 κ ¯ n ( ω ) , (S53)where the second term of Γ β ( t ) has been neglected be-cause it is the high-order term. These values are justthe parameters appearing in the master equation derivedunder the conventional Markovian approximation.1 S-3. ENERGY SPECTRUM INSINGLE-EXCITATION SUBSPACE
Here we want to further emphasize that the roots ofEq. (8) in the main text give the energy spectrum of thetotal system in the single-excitation subspace. There-fore, we call the state represented by the discrete root E as the stationary state of the total system. Further-more, since the single-excitation process in our Hamilto-nian is solely contributed by the vacuum of the environ-ment, this gives further evidence on calling u ( t ) as thevacuum-induced dissipation function. The eigensolutionin the single-excitation subspace can be obtained in thefollowing. The total excitation number of our model, i.e.,ˆ N tot = ˆ a † ˆ a + P k ˆ b † k ˆ b is conserved, which means that theHilbert space is divided into independent subspaces withdefinite h ˆ N i . The eigenstate of the total system in thesubspace with h ˆ N tot i = 1 can be expanded as | ϕ i = c | , { k }i + X k d k | , k i , (S54)substituting which into the eigenequation, we get c ω + X k g k d k = E c , (S55) c g k + d k ω k = E d k . (S56)From Eq. (S56), we have d k = g k c E − ω k , substituting whichinto Eq. (S55), we obtain ω + X k | g k | E − ω k = E . (S57)It matches well with Eq. (8) in the main text in thecontinuous limit of the environmental frequency. Theformation of the stationary state is in a same mechanismas the spin-boson model case in Ref. [7, 8].The formation of this stationary state would cause thetransient negative of the decay rate, which reveals thatthe system has a strong non-Markovian effect. Such ef-fect can be quantified by the non-Markovianity [9]. Thisquantitative characterization in turn can supplies someinsight to the profound impact of the formed station-ary state on the dynamics of the system. The non-Markovianity is defined as N = max ρ (0) ,ρ (0) Z σ> σ ( t, ρ , ( t )) dt, (S58)where σ ( t, ρ , ( t )) = ddt D ( ρ ( t ) , ρ ( t )) is the change rateof the trace distance D ( ρ ( t ) , ρ ( t )) between two densitymatrix ρ ( t ) and ρ ( t ). The trace distance is definedas D ( ρ , ρ ) = Tr | ρ − ρ | with | ρ | ≡ p ρ † ρ . The in-tegration is performed over all the time duration with σ ( t, ρ , ( t )) > ρ , (0). It has been proven inRef. [9] that for zero-temperature environment and the FIG. S1. The behavior of the non-Markovianity N with thechange of η in the case of the Ohmic spectral density. ω c =8 . ω and β = ∞ have been used. The stationary state ispresent when η > . /ω , which guarantees to include allthe time region with σ ( t, ρ , ( t )) > one-excitation-number involved at most in the dynamicsof the whole system, N takes the maximal value when ρ (0) = | ih | and ρ (0) = | ih | .In Fig. S1, we plot the non-Markovianity calculated inthe zero-temperature case, where the effect of the formedstationary state can be manifested clearly. We can seethat when the stationary state is absent for the param-eter regime η < . N equals to zero persistently.Whenever the stationary state is formed, i.e., η > . N takes non-zero value. It is understandable based onthe fact the presence of the stationary state would causethe decay coefficients negative. The negative decay im-plies that the energy or information flows back from theenvironment to the system, which results in the increaseof the trace distance between the two initial states. S-4. THE DERIVATION OF EQ. (6)
Eq. (6) is valid in the parameter regime without thestationary state, where both Γ( t ) > β ( t ) > N = ˆ n r − ˆ n l and { ˆ K = ˆ n r − ˆ n l , ˆ K + = ˆ a r † ˆ a l , ˆ K − = ˆ a r ˆ a l † } ,where ˆ • r and ˆ • l operates, respectively, on the right andleft vector, we can transform Eq. (3) into a Schr¨odinger-equation-like form˙ ρ ( t ) = L ( t ) ρ ( t ) , (S59) L ( t ) = − i Ω( t )( ˆ N + 1) + A ( t ) ˆ K + + B ( t ) ˆ K − − [ A ( t ) + B ( t )] ˆ K − A ( t ) − B ( t )2 , (S60)where A ( t ) = 2 × Γ β ( t )2 and B ( t ) = 2 × [Γ( t ) + Γ β ( t )2 ]. The2actions of these operators on the state areˆ N | m ih n | = ( m − n − | n ih n | , (S61)ˆ K | m ih n | = m + n + 12 | m ih n | , (S62)ˆ K + | m ih n | = p ( m + 1)( n + 1) | m + 1 ih n + 1 | , (S63)ˆ K − | m ih n | = √ mn | m − ih n − | . (S64)ˆ N forms a u(1) algebra and { ˆ K = ˆ n r − ˆ n l , ˆ K + =ˆ a r † ˆ a l , ˆ K − = ˆ a r ˆ a l † } form a su(1,1) algebra,[ ˆ K , ˆ K ± ] = ± ˆ K ± , [ ˆ K − , ˆ K + ] = 2 ˆ K , (S65)[ ˆ N , ˆ K l ] = 0 , ( l = ± , . (S66)With the dynamical symmetry structure at hands andby introducing a time-dependent similarity transforma-tion ˆ U = e α + ( t ) ˆ K + e α − ( t ) ˆ K − with α ± ( t ) satisfying dα + ( t ) dt = B ( t )[ α + ( t ) − A ( t ) B ( t ) ][ α + ( t ) − , (S67) dα − ( t ) dt = B ( t )[1 − α + ( t ) α − ( t )]+[ A ( t ) + B ( t )] α − ( t ) , (S68)under the initial condition α ± (0) = 0, the time-dependent master equation (S59) can be diagonalizedinto ˙¯ ρ ( t ) = ˆ¯ L ( t )¯ ρ ( t ) , (S69)ˆ¯ L ( t ) = ˆ U − ( t ) ˆ L ( t ) ˆ U ( t ) − ˆ U − d ˆ Udt = − i Ω( t )( ˆ N + 1) + γ ( t ) ˆ K − A ( t ) − B ( t )2 (S70)with γ ( t ) = 2 B ( t ) α + ( t ) − [ A ( t ) + B ( t )] and ¯ ρ ( t ) =ˆ U − ρ ( t ). Then the solution of ¯ ρ ( t ) as well as ρ ( t ) canbe obtained as ρ ( t ) = ˆ U e R t ¯ L ( τ ) dτ ρ (0)= e R t γ ( τ ) − [ A ( τ ) − B ( τ )]2 dτ e α + ( t ) ˆ K + e α − ( t ) ˆ K − X n,k ρ n,k × e R t [ ( n + k ) γ ( τ )2 − i ( n − k )Ω( τ )] dτ | n ih k | , (S71)where the initial state ρ (0) = P n,k ρ n,k | n ih k | has beenused. The asymptotic form of the obtained solution (S71)can be analyzed by the following method. We study firstthe asymptotic forms of α ± ( t ). From Eq. (S67), we have dα + ( t ) dt > , if 0 < α + ( t ) < A ( t ) B ( t ) dα + ( t ) dt < , if A ( t ) B ( t ) < α + ( t ) < A ( t ) and B ( t ) satisfy A ( t ) > B ( t ) > α + ( t ) evolves asymptotically from its initialvalue α + (0) = 0 to the unique stable value α + ( ∞ ) = A ( ∞ ) B ( ∞ ) , (S73)under the conditions A ( t ) = 2 × Γ β ( t )2 > B ( t ) =2 × [Γ( t ) + Γ β ( t )2 ] >
0, which is always satisfied when thestationary state is absent [see Fig. 1(a,b), Fig. 2(a,b),and Fig. 6(a,b)]. To check the asymptotic form of α − ( t ),we define y ( t ) = α − ( t ) e R t γ ( τ ) dτ , whose time derivativereads dy ( t ) dt = [ dα − ( t ) dt + α − ( t ) γ ( t )] e R t γ ( τ ) dτ = B ( t ) e R t γ ( τ ) dτ , (S74)where Eq. (S68) has been used. Because of γ ( ∞ ) = 2 B ( ∞ ) α + ( ∞ ) − [ A ( ∞ )+ B ( ∞ )] = − ∞ ) < , (S75)we thus havelim t →∞ dy ( t ) dt = 0 ⇒ y ( ∞ ) = Const. , (S76)which in turn implies α − ( ∞ ) = ∞ . (S77)With the asymptotic forms (S73, S76, S77) at hands,we are ready to evaluate the asymptotic form of Eq.(S71). Eq. (S71) can be expanded as ρ ( t ) = e R t γ ( τ ) − [ A ( τ ) − B ( τ )]2 dτ e α + ( t ) ˆ K + X n,k ρ n,k e R t [ ( n + k ) γ ( τ )2 − i ( n − k )Ω( τ )] dτ k X m =0 α m − ( t ) m ! s n ! k !( n − m )!( k − m )! | n − m ih k − m | . (S78)3Setting m = n + k − p and remembering the form of y ( t ), we have ρ ( t ) = e R t γ ( τ ) − [ A ( τ ) − B ( τ )]2 dτ e α + ( t ) ˆ K + X n,k ρ n,k e − i R t ( n − k )Ω( τ ) dτ y ( t ) n + k × n + k X p = n − k α − p − ( t )( n + k − p )! s n ! k !( p + n − k )!( p − n − k )! | p + n − k ih p − n − k | Due to the asymptotic form (S77), only the term p = 0,i.e., n = k , in the summation over p survives in the long-time limit. Thuslim t →∞ ρ ( t ) = e α + ( ∞ ) ˆ K + X n ρ n,n (Const.) n | ih | = Ce α + ( ∞ ) ˆ K + | ih | , (S79)where the asymptotic forms (S76) and (S75) have beenused and C can be determined by normalization. The final form of Eq. (S79) takes the form ρ ( ∞ ) = C ∞ X n =0 [ A ( ∞ ) B ( ∞ ) ] n | n ih n | = C ∞ X n =0 [ Γ β ( ∞ ) / (2Γ( ∞ ))1 + Γ β ( ∞ ) / (2Γ( ∞ )) ] n | n ih n | . (S80)Remembering both Γ( ∞ ) and Γ β ( ∞ ) are positive val-ues when the stationary state is absent, we can readilycalculate the normalization factor C = β ( ∞ ) / (2Γ( ∞ )) .In summary, in obtaining Eq. (6), we are working inthe parameter regime where the stationary state is absentsuch that Γ( t ) > β ( t ) > t ) and Γ β ( t ) may transiently takes negative values andasymptotically tends to zero. In this regime, Eq. (6)does not work. It implies the breakdown of the canonicalequilibration in our system. [1] R. P. Feynman and F. L. Vernon, Ann. Phys. (N. Y.) ,118 (1963).[2] A. Caldera and A. Leggett, Physica A , 587 (1983).[3] B. L. Hu, J. P. Paz, and Y. Zhang Phys. Rev. D. ,2843 (1992).[4] J.-H. An and W. M. Zhang, Phys. Rev. A , 042127(2007).[5] W. M. Zhang, D. H. Feng, and R. Gilmore, Rev. Mod.Phys. , 867 (1990).[6] J.-H. An, Y. Yeo, and C. H. Oh, Ann. Phys. , 1737 (2009).[7] Q.-J. Tong, J.-H. An, H.-G. Luo, and C. H. Oh, Phys.Rev. A , 052330 (2010).[8] H.-B. Liu, J.-H. An, C. Chen, Q.-J. Tong, H.-G. Luo,and C. H. Oh, Phys. Rev. A , 052139 (2013).[9] H.-P. Breuer, E.-M. Laine, and J. Piilo, Phys. Rev. Lett. , 210401 (2009).[10] S. J. Wang, M. C. Nemes, A. N. Salgueiro, and H. A.Weidenm¨uller, Phys. Rev. A66