Partly invariant steady state of two interacting open quantum systems
aa r X i v : . [ qu a n t - ph ] D ec Partly invariant steady state of two interacting open quantum systems
J´ozsef Zsolt Bern´ad ∗ and Juan Mauricio Torres Institut f¨ur Angewandte Physik, Technische Universit¨at Darmstadt, D-64289, Germany (Dated: July 17, 2018)We investigate two interacting open quantum systems whose time evolutions are governed byMarkovian master equations. We show a class of coupled systems whose interaction leaves invariantthe steady state of one of the systems, i.e., only one of the reduced steady states is sensitive to theinteractions. A detailed proof with the help of the Trotter product formula is presented. We applythis general statement to a few models, one of which is the optomechanical coupling model wherean optical cavity is coupled to a small mechanical oscillator.
PACS numbers: 03.65.Yz
I. INTRODUCTION
Quantum Markovian master equations are the simplestbut also the most used cases of study for the dynam-ics of open quantum systems. The concept of open sys-tems plays an important role since perfect isolation inthe quantum realm and access to all degrees of freedomis not possible. The main source of inspiration was thetheory of lasers where models describing the interactionof a central subsystem weakly coupled to an uncontrol-lable environment were studied. In the last two decadesthere has been great interest in the study of compos-ite quantum systems in order to create entanglement be-tween them [1] with the ultimate aim to realize quantuminformation procedures [2]. However, there is always anenvironment causing decoherence which is also the bot-tleneck to realize highly efficient quantum informationprotocols. Therefore the theory of open quantum sys-tems was frequently applied to various models [3].There is a subclass of models which can be decomposedinto many interacting parts which in turn are coupled toindependent environments. The best example is the caseof two-level atoms interacting with a single mode of theradiation field and coupled to independent reservoirs, likethe model of the many-atom laser [4] or the one-atommaser [5]. The time evolution of the composite system isderived from a microscopic model or is assembled on theground of phenomenological considerations from equa-tions which describe independently the parts. The lattercase usually assumes a low coupling strength between thetwo parties which allows a separate description of dissi-pative effects. In both cases the possible treatments ofenvironments consist in Markovian and non-Markovianapproaches [6]. Obtaining the full solution to the dy-namical equation is a very complicated task even with aMarkovian approach [5, 7] and sometimes one has to besatisfied by identifying the steady state [8] or its symme-try properties [9]. In general, the role of the interactionin a composed system was studied typically for two cou-pled harmonic oscillators [10, 11] and for two interacting ∗ Electronic address: [email protected] identical atoms [12].In this paper, we focus on the case when two indepen-dent quantum systems A and B are subject to Markovianmaster equations and their interaction is studied throughan equation which is obtained by adding the two masterequations to a Hamiltonian interaction. For physical sys-tems, this consideration usually implies that the interac-tion strength between A and B has to be the smallestparameter in the system. We concentrate on a specialclass of interaction Hamiltonians that cause dephasingon system A and lead to a generator (Liouvillian) of thecomplete system that commutes with the Liouvillian ofsystem A ’s master equation. The motivating example isthe optomechanical coupling between a single-mode elec-tromagnetic field and a small mechanical oscillator in thepresence of independent decoherence mechanisms [13]. Itwill be demonstrated that the considered class of inter-action Hamiltonians leads to an invariant steady state ofsystem A , i.e., tracing out system B in the steady state ofthe composite system one always obtains the same steadystate of system A which is independent of the interactionHamiltonian. In contrast, the steady state of system B may depend on the interaction. In the case of finite di-mensional systems, theorems of linear algebra can be em-ployed to prove the statement. However, when dealingwith infinite dimensional systems a different and moresophisticated method is required. Therefore, the proofof the invariant steady state is based on the propertiesof quantum dynamical semigroups where we employ theTrotter product formula and identify the growth boundsof the master equations. While the proof of our statementis rather general we will explicitly show it in some simplemodels. The steady state of system B does not have anyparticular property and its solution is model dependent.Therefore in our examples we study the steady state ofsystem B by analyzing its average excitation. Our workpresents a general approach towards composite systemsgoverned by independent Markovian master equations.This paper is organized as follows. In Sec. II we discussthe motivation of our analysis. In Sec. III we show thegeneral statement about the steady state. In Sec. IV, weapply these findings to three different systems: two inter-acting spins, a spin interacting with a harmonic oscillatorand two coupled harmonic oscillators. A discussion aboutour findings is summarized up in Sec. V. In Appendix Awe briefly discuss the Liouvillians of the Markovian mas-ter equations involving some technical problems regard-ing infinite dimensional Hilbert spaces. In Appendix Bwe analyse the form of the equilibrium state of a dampedharmonic oscillator in the number state representation. II. MOTIVATION
There are many physical situations in which the dy-namics of the system can be decomposed into two par-ties subject to a Hamiltonian interaction. The individualsub-systems might include a relaxation mechanism thatunder the Markovian assumption leads to a master equa-tion for the density matrix ˆ ̺ of the composite systemthat can be written as d ˆ ̺dt = L ˆ ̺ = L A ˆ ̺ + L B ˆ ̺ + i [ˆ ̺, ˆ H I ] . (1)The superoperators L A and L B are generators of the in-dividual systems and ˆ H I is the interaction Hamiltonian.One system that allows this type of separation and hasgained a lot of recent interest is the model of optomechan-ical coupling, where the radiation field and a mechanicaloscillator are coupled via radiation pressure. In the sin-gle mode scenario, numerical simulations show that thesteady state of the radiation field remains independentof the coupling strength and of the parameters of themechanical oscillator. This interesting effect can be gen-eralized to two systems where one of them is composedof spins and harmonic oscillators and the other systemstays arbitrary and the interaction between them holdsthe same symmetries as the optomechanical coupling. Inthe following we motivate the problem by showing theoptomechanical model. We consider the simplest Hamil-tonian describing the interaction between a single modeof the radiation field and a vibrational mode of a me-chanical oscillator. This interaction is derived for a cavitywith a moving mirror [14, 15] and the quantization of theradiation field is considered with time-varying boundaryconditions. The single mode assumption for both quan-tized field and mirror motion was successfully applied todescribe most of the experiments to date and leads to thefollowing interaction Hamiltonian ( ~ = 1)ˆ H I = g ˆ a † ˆ a (ˆ b † + ˆ b ) . (2)Thereby, ˆ a (ˆ a † ) is the annihilation (creation) operator ofthe single mode radiation field with frequency ω and ˆ b (ˆ b † ) is the annihilation (creation) operator of the mov-ing mirror’s mode with frequency ν . The parameter g describes the optomechanical coupling strength. This in-teraction Hamiltonian commutes with the free Hamilto-nian of the radiation field ω ˆ a † ˆ a but it does not commutewith the free Hamiltonian of the mirror ν ˆ b † ˆ b . Account-ing to an irreversible damping mechanism the Liouville operators of the individual systems can be cast in thefollowing form L A ˆ ̺ = iω [ˆ ̺, ˆ a † ˆ a ] + κ (¯ n + 1) D [ˆ a ]ˆ ̺ + κ ¯ n D [ˆ a † ]ˆ ̺, (3) L B ˆ ̺ = iν [ˆ ̺, ˆ b † ˆ b ] + γ ( ¯ m + 1) D [ˆ b ]ˆ ̺ + γ ¯ m D [ˆ b † ]ˆ ̺, which are written in terms of the Lindbladian dissipator D [ˆ a ]ˆ ̺ = ˆ a ˆ ̺ ˆ a † − (ˆ a † ˆ a ˆ ̺ + ˆ ̺ ˆ a † ˆ a ) . (4)The constant κ is the loss rate of the radiation field and γ is the damping rate of the mirror’s mode. Both areassumed to be coupled to a thermal reservoir with anaverage number of photons ¯ n and phonons ¯ m .An interesting feature can be distinguished fairly easilyin this model by evaluating the mean occupation photonand phonon numbers at the steady state of the dynamics.This calculation is explained in detail in Sec. IV and herewe merely show the result, namely h ˆ a † ˆ a i st = ¯ n, (5) h ˆ b † ˆ b i st = ¯ m + 4¯ n g γ + 4 ν + 4¯ n (¯ n + 1)(2 κ + γ ) g γ [(2 κ + γ ) + 4 ν ] . The most striking effect is that the mean photonic occu-pation number remains unchanged, i.e. it is insensitiveto the coupling to the mechanical oscillator. In contrastthe mean phonon number grows with increasing photonnumber ¯ n and coupling strength g . In the special casewhen the photonic bath is at zero temperature (¯ n = 0),the occupation numbers of the cavity and the movingmirror remain unchanged. It is not a complicated taskto show that the composite steady state is insensitive tothe interaction in this situation. This is actually the ap-proximate situation of current experimental scenarios, asfor optical frequencies ¯ n can be regarded to be zero. Inaddition, the phononic contributions due to the couplingin the second line of Eq. (5) are negligible in the low cou-pling regime for low values of ¯ n . Nevertheless, with thenew generations of optomechanical systems reaching thestrong coupling regime and for microwave frequencies,this effect can become important. It is worth mentioningthat if the coupling is too strong, a more detailed modelof the loss mechanisms has to be considered such as inRef. [16] where the eigenstates of the Hamiltonian in Eq.(2) where used to derive the corresponding master equa-tion. It can be shown that even in this case the effectpersists: the mean photon number remains unchanged.We will comment more on this at the end of Sec. III.In an attempt to explain the effect one can realize thate − i ˆ H I t (cid:0) L A ˆ ̺ (cid:1) e i ˆ H I t = L A (cid:16) e − i ˆ H I t ˆ ̺ e i ˆ H I t (cid:17) . (6)It turns out that by identifying this special characteris-tic one is able to generalize this behaviour. There aretwo main purposes of this work. One is to demonstratethat for ¯ n = 0 the steady state of the optical mode re-mains unchanged in contrast to the changing mechanicalmode. Second, we will extend the result to a broaderclass of systems that share the same symmetry relations.We will show that these arguments result in an invariantsteady state for the system with symmetries and this isindependent of the details of the other system. III. STEADY STATE ANALYSIS
In this section we generalize the statement presentedin Sec. II and prove it by using the properties of dy-namical semigroups. First, let us consider two quan-tum systems A and B , where each of these systemsundergo a Markovian non-unitary dynamics, describedby the Gorini-Kossakowski-Sudarshan-Lindblad genera-tor [17–19]. The non-unitary evolution is a trace of theinteraction of each system with an independent environ-ment. These environments are considered to be such thatthe evolution of systems A and B has a Markovian irre-versible behaviour. The dynamical equation of the jointsystem d ˆ ̺/dt = L ˆ ̺ is governed by the Liouville operator L ˆ ̺ = L A ˆ ̺ + L B ˆ ̺ − i [ ˆ V A ˆ V B , ˆ ̺ ] . (7)We assume that system A is composed of N P harmonicoscillators and N S spins and that its dynamics is de-scribed by the Lindblad master equation L A ˆ ρ A = − i [ ˆ H A , ˆ ρ A ] + ˜ L A ˆ ρ A , ˜ L A = N P X p =1 L ( p )HO + N S X s =1 L ( s )S . (8)The operator ˆ ρ A denotes the reduced density matrix ofsystem A alone and the individual Liouvillian for eachoscillator has the form L ( p )HO = κ p (¯ n p + 1) D [ˆ a p ] + κ p ¯ n p D [ˆ a † p ] , (9)which is written in terms of the dissipator defined in Eq.(4). ˆ a p (ˆ a † p ) is the annihilation (creation) operator of the p th harmonic oscillator. Similarly we have for each spin L ( s )S = γ s (1 − ¯ m s ) D [ˆ σ ( s ) − ] + γ ¯ m s D [ˆ σ ( s )+ ] , (10)with the ladder operators ˆ σ ( s ) ± = (ˆ σ ( s ) x ± i ˆ σ ( s ) y ) / s thspin. All the constants γ p and κ s have dimension of fre-quency and are assumed to be non-zero. The dimen-sionless parameters { ¯ n p ≥ } p =1 ,..,N p are proportional tothe temperature of the individual thermal reservoir towhich each oscillator is coupled and 0 ≤ ¯ m s ≤ s ∈ { , ..., N s } .It is a challenge in the description of a harmonic os-cillator that we have to work with unbounded operators,which cannot be defined everywhere. This issue is inher-ited by the Liouvillian in (9) which is defined on elementsof a Banach space. In order to discuss the properties ofthis Liouvillian we require tools from the theory of dy-namical semigroups. These tools are supported by Ap-pendix A where we show the spectrum of the Liouvillians and we discuss the definition of the adjoint. Examplesrelated to the domain of the harmonic oscillator’s Liou-villian are also presented.We assume for system A the following form of the in-teraction Hamiltonianˆ V A = N P X p =1 ˆ a † p ˆ a p + N S X s =1 ˆ σ + s ˆ σ − s . (11)The operator ˆ V A can be regarded as the excitation op-erator of system A and therefore it is a dephasing typeinteraction. Our last assumption is that the Hamiltonianof system A commutes with the interaction Hamiltonian,i.e. [ ˆ H A , ˆ V A ] = 0 . (12)Examples of ˆ H A are the Jaynes-Cummings, Ising, andTavis-Cummings Hamiltonians. If these conditions aremet, we can state that the steady state of the reducedsystem A is independent of the interaction with system B . This can be restated as follows: if ˆ ̺ st is the steadystate of the complete system and solution of the equation L ˆ ̺ st = 0, then L A Tr B { ˆ ̺ st } = 0 , (13)and therefore ˆ ρ A st = Tr A { ˆ ̺ st } with L A ˆ ρ A st = 0. This state-ment does not hold in general for system B as L B Tr A { ˆ ̺ st } 6 = 0 . (14)To prove this statement let us start by rewriting theLiouvillian L in terms of superoperators in the followingform L = L B + ˜ L A − i ˜ C A − i A A C B − i C A A B (15)where C i ˆ ̺ = [ V i , ˆ ̺ ] , A i ˆ ̺ = { V i , ˆ ̺ } , i ∈ { A, B } , ˜ C A ˆ ̺ = [ H A , ˆ ̺ ] . (16)One can note that the superoperator C A commutes withthe Liouvillian, [ C A , L ] = 0 . (17)The commutation with the first and last three terms inEq. (15) is evident from Eqs. (12) and (16), and be-cause superoperators acting separately on systems A and B commute. The commutation with ˜ L A can be shownusing the commutation relation between the interactionHamiltonian ˆ V A and the Lindblad operators.In order to exploit the commutation property between L and C A , it is convenient to work in the eigenbasis ofthe superoperator C A , which is defined in (16). As it isformed by the operator ˆ V A it is convenient to identifyfirst the eigenvalue equation of ˆ V A , i.e.ˆ V A | n, j i = n | n, j i , n, j ∈ N . (18)Due to (12) we haveˆ V A ˆ H A | n, j i = n ˆ H A | n, j i (19)which means that ˆ H A does not couple vectors with dif-ferent n .The kets | n, j i represent states where the number ofphotons plus the number of excited spins sum to n . Thesestates have a degeneracy d n = d n ( N P , N S ) that is enu-merated by the index j = 1 , , . . . d n and depends onthe number of spins and oscillators in system A . For in-stance, in the Jaynes-Cummings model these states are | n, i = | n i ⊗ | e i and | n, i = | n + 1 i ⊗ | g i . Using thestates in (18) one can find that the eigenvectors of thesuperoperator C A are given byˆ σ l,ν = (cid:12)(cid:12)(cid:12) n + | l | + l , j E D n + | l |− l , k (cid:12)(cid:12)(cid:12) , ν = { n, j, k } , (20)with l ∈ Z and where we have introduced the collectiveindex ν . These eigenvectors are rank 1 operators and itcan be verified that they solve the eigenvalue equation C A ˆ σ l,ν = [ ˆ V A , ˆ σ l,ν ] = l ˆ σ l,ν (21)with the eigenvalue l ∈ Z . As these eigenvectors aredegenerate, any superposition of σ l,ν with fixed l is alsoan eigenvector of C A . Furthermore, they are orthogonalwith respect to the Hilbert-Schmidt inner product(ˆ σ l,ν | ˆ σ l ′ ,ν ′ ) = Tr A n ˆ σ † l,ν ˆ σ l ′ ,ν ′ o = δ l,l ′ δ ν,ν ′ , (22)where we have introduced for convenience a ket-bra no-tation | ˆ σ l,ν ) and (ˆ σ l,ν | to write the inner product. Let B ( H ) be the set of all bounded linear operators on aHilbert space H . Given the fact that the states | n, j i are complete in the Hilbert space H A it follows that theeigenvectors ˆ σ l,ν form a complete orthonormal basis forthe Hilbert space of Hilbert-Schmidt operators B ( H A ) := (cid:8) ˆ x ∈ B ( H A ) : Tr (cid:8) ˆ x † ˆ x (cid:9) < ∞ (cid:9) (23)which contains the set of all density matrices of system A . Therefore, in the previously defined ket-bra notationone can write the completeness relation as X l P l = I , with P l = X ν | ˆ σ l,ν )(ˆ σ l,ν | . (24)Thereby we have introduced the projectors P l thatproject onto the subspace of fixed value of the eigen-value l of C A . Using the completeness relation one canseparate the Liouville operator in different subspaces as L = P l P l L where we have P l L = LP l . (25) P l L acts only onto the subspace containing eigenelementsˆ σ l,ν of C A for fixed value of l . In particular this superop-erator can be written in terms of the projectors as P l C A = C A P l = l P l . (26) In search of the steady state structure restricted to sys-tem A we consider the reduced dynamical equation andwe note that for the subspace labeled by l = 0 one has P Tr B {L ˆ ̺ } = (cid:16) ˜ L A − i ˜ C A (cid:17) P Tr B { ˆ ̺ } . (27)We have used Eqs. (15) and (26), and the relationsTr B (cid:8) L B ˆ ̺ (cid:9) = Tr B (cid:8) C B ˆ ̺ (cid:9) = 0. This shows that if ˆ ̺ st is the steady state of L and ˆ ρ A st is the steady state of L A then P ˆ ρ A st = P Tr B { ˆ ̺ st } . (28)This still does not imply the statement in Eq. (13) aswe have to show that the equality also holds for the sub-spaces l = 0.For the remaining part of the proof we use a moretechnical approach based on known results of quantumdynamical semigroups [20]. In this sense, let us recallthat the Liouville operator L in Eq. (7) is the genera-tor of the quantum dynamical semigroup e L t . These twooperators act on density matrices ˆ ̺ which are trace-classoperators in the Hilbert space H = H A ⊗ H B . The Ba-nach space of all self-adjoint trace-class linear operatorscontaining the set of density matrices can be defined as[20] L ( H ) := { ˆ ̺ ∈ B ( H ) : ˆ ̺ † = ˆ ̺, k ˆ ̺ k < ∞} . (29)with the norm || ˆ ̺ || = Tr np ˆ ̺ † ˆ ̺ o . (30)It was shown by Kossakowski (see Definition 1 and The-orem 1 in [20]) that the Liouville operator L generatesa contraction semigroup on L ( H ). This means that thenorm of the evolution operator is bounded as k e L t k op ≤ e L t ˆ ̺ ∈ L ( H ) is bounded as k e L t ˆ ̺ k ≤ k e L t k op k ˆ ̺ k ≤ k ̺ k , ∀ ̺ ∈ L ( H ) (31)where we introduced the operator norm kOk op = sup {kO ˆ x k : ˆ x ∈ L ( H ) , k ˆ x k = 1 } . (32)Our goal now is to find a smaller bound than the boundgiven in Eq. (31) for the elements ˆ ̺ ( t ) = e L t ˆ ̺ but pro-jected to the subspaces with fixed value | l | . For this pur-pose we consider the projectors Q l = P l + P − l , l = 0 , (33)which act as the identity operator for system B . Bytaking fixed values of | l | one ensures that the self-adjointproperty is preserved, i.e. ( Q l ˆ ρ ) † = Q l ˆ ρ whenever ˆ ρ † = ˆ ρ .This does not hold when projecting with P l =0 . Therefore Q l ˆ ̺ ∈ L ( H ) and the operator norm of each projector k Q l k op = 1. The projector Q l projects all elements ˆ x ∈ L ( H ) onto a closed linear subspace of L ( H ) which is alsoa Banach space with norm k . k : L l ( H ) := {Q l ˆ x : ˆ x ∈ L ( H ) } . (34)In the next step we take advantage of the fact thatthe generators L and ˜ L A do not mix spaces of differentvalue of l , i.e. Q l L = LQ l . This allows us to define theoperators restricted to L l ( H ) L l := Q l L = LQ l , ˜ L Al := Q l ˜ L A = ˜ L A Q l . (35)Now we can evaluate the norm of the projected elementsas kQ l ˆ ̺ ( t ) k ≤ k e L l t k l op kQ l ˆ ̺ (0) k . (36)where we also used that Q l = Q l . The operator norm k . k l op is induced by the Banach space L l ( H ). In order tofind a better bound we take advantage of our knowledgeof the generator ˜ L A and separate the dynamics with theuse of the Trotter product formula (see Corollary 5 . e L l t = lim N →∞ (cid:16) e ˜ L Al t/N e ( L l − ˜ L Al ) t/N (cid:17) N . (37)As discussed in Appendix A the Liouvillian ˜ L A hasonly a point spectrum . Therefore in the Banach space L l ( H ) the dynamical semigroup of system A has the op-erator norm k exp( ˜ L Al t ) k l op ≤ e η l t (see Definition 2 . .
12 in Chap. IV of Ref. [21]) where η l ∈ R − isthe largest eigenvalue of ˜ L Al in the Banach space L l ( H )( l = 0). This follows from the fact that there is only onezero eigenvalue of ˜ L A with corresponding eigenvector ˆ ρ A st that has no projection on the subspaces with l = 0, i.e. Q l ˆ ρ A st = 0 and P ˆ ρ A st = ˆ ρ A st . The rest of the eigenvalues of˜ L A are negative. This allows us to find a bound indepen-dent of N for the evolution operator that we expandedusing the Trotter formula in Eq. (37). Using the submul-tiplicative property of the operator norm k . k l op and thefact that ( L l − ˜ L Al ) also generates a strongly continuousone-parameter contracting semigroup, we have k e ˜ L Al t/N e ( L l − ˜ L Al ) t/N k l op ≤ k e ˜ L Al t/N k l op k e ( L l − ˜ L Al ) t/N k l op ≤ k e ˜ L Al t/N k l op ≤ e η l t/N (38)which yields k lim N →∞ (cid:16) e ˜ L Al t/N e ( L l − ˜ L Al ) t/N (cid:17) N k l op ≤ e η l t . (39)This upper bound is independent of N which implies thefollowing bound for the norm ||Q l ˆ ̺ ( t ) || ≤ e η l t || ̺ (0) || ⇒ lim t →∞ Q l ˆ ̺ ( t ) = 0 . (40)Thus the steady state has projection 0 on the subspaces l = 0. This fact does not depend on the interaction Hamiltonian ˆ V A ˆ V B . This is a stronger result than theoriginal requirement for our statement, becauseTr B {Q l ˆ ̺ ( t → ∞ ) } = 0 ⇒ Tr B { ˆ ̺ ( t → ∞ ) } = ˆ ρ A st , (41)but we found that Q l ˆ ̺ ( t → ∞ ) = 0.The above approach can also be used for system A alone in the case when we know the spectrum of the Li-ouville operator, but due to the internal interactions itis hard to solve the dynamics or to find the steady state.As an example someone can think about the Jaynes-Cummings model where the two-level system and theharmonic oscillator have an independent damping mech-anism. Despite the fact that finding the analytical formof the steady state is a challenging task, in an appropri-ate representation we have to deal only with the blockdiagonal part of the density matrix.We end this section with a comment on the choice ofseparate damping mechanisms. This is a valid descrip-tion for weakly interacting systems, where the environ-ments are insensitive to the dressed eigenstates of coupledsystems A and B . For strong coupling, the physical de-scription is not valid anymore and a more detailed lossymodel has to be derived. For instance, in Ref. [16] thiswas investigated for the optomechanical model and themaster equation was derived using the dressed state ofthe Hamiltonian in Eq. (2). The result is a master equa-tion that differs from Eq. (1) and (3) by an additional dis-sipative term proportional to D [ˆ a † ˆ a ] and where the me-chanical dissipators are replaced by D [ˆ b ] → D [ˆ b − g/ν ˆ a † ˆ a ]and D [ˆ b † ] → D [ˆ b † − g/ν ˆ a † ˆ a ]. With these corrections, thecondition in Eq. (17) is still fulfilled and therefore thesteady of system A , i.e. the optical mode, remains invari-ant as well. In general, our method applies to any twodissipative systems A and B with a separable Hamilto-nian interaction ˆ V A ˆ V B and where: the assumptions (12)and (17) are satisfied, and one can identify a generator˜ L A acting solely on system A whose zeroth eigenvalue be-longs only to the eigenspace l = 0 of C A (see Eq. (21)).There are techniques to infer about the spectrum of a fi-nite dimensional operator. This implies that our resultscan be applied to finite dimensional systems includingcorrections due to the strong coupling limit. However,in the case of infinite dimensional Hilbert spaces withnon-zero temperature baths, the choice is restricted be-cause only the spectral decomposition of the optical mas-ter equation (damped harmonic oscillator) is known. IV. EXAMPLESA. A case study of two interacting spins
The situation can be easily understood for the sim-plest case of two interacting spin systems ( A and B ) withindependent Markovian dynamics. We consider an in-teraction Hamiltonian which satisfies the symmetries ofthe optomechanical system. Due to the simple structureof the Hilbert spaces the steady state of the whole sys-tem has an analytical formula. The Hilbert spaces are H A = H B = C with orthonormal basis {| i A , | i A } forsystem A and {| i B , | i B } for system B . Each systemalone is governed by the master equation d ˆ ρ i /dt = L i ˆ ρ i with i ∈ { A, B } expressed in terms of the the Liouvillians L i ˆ ρ i = iω [ˆ ρ i , ˆ σ iz ] + γ i (1 − s i ) D [ˆ σ i − ]ˆ ρ i + γ i s i D [ˆ σ i + ]ˆ ρ i (42)where ω i is the frequency difference between the energylevels | i i and | i i . γ i has the dimension of frequencyand 0 ≤ s i ≤
1. The spin operators are given byˆ σ iz = | i i h | i − | i i h | i , ˆ σ i − = (ˆ σ i + ) † = | i i h | i . (43)The steady state of each separate system has the sameform ˆ ρ i st = (1 − s i ) | i i h | i + s i | i i h | i . (44)By considering the interaction Hamiltonianˆ H I = Ωˆ σ Az (cid:0) ˆ σ B + + ˆ σ B − (cid:1) , (45)one can realize that the interaction Hamiltonian com-mutes with the free Hamiltonian of system A , i.e.[ ˆ H I , ˆ σ Az ] = 0 . (46)The steady state ˆ ̺ st of the whole system can be analyti-cally evaluated due to the small dimension of the Hilbertspace. We do not show the explicit expressions here butsimply emphasize that it has the formˆ ̺ st = ρ , , , ρ , , , ρ , , , ρ , , , ρ , , , ρ , , , ρ , , , ρ , , , , (47)where ρ i,j,k,l = h i | A h j | B ˆ ̺ st | k i A | l i B and i, j, k, l ∈ { , } .The reduced steady state of system A is given byTr B { ˆ ̺ st } = (cid:18) − s A s A (cid:19) = ˆ ρ A st , (48)whereas for system B it readsTr A { ˆ ̺ st } = (cid:18) B B B ∗ − B (cid:19) = ˆ ρ B st . (49)The steady state of system B is given by the real pa-rameter B and the complex parameter B both of thembeing functions of ω , Ω, γ i and s i ( i ∈ { A, B } ). In thecase of strong interaction between systems A and B thefollowing relations holdlim Ω →∞ B = 1 / , lim Ω →∞ B = 0 , (50)resulting in a maximally mixed state for Tr A { ˆ ̺ st } .In Fig. 1 we found that the population inversion h ˆ σ Bz i st = Tr { ˆ σ Bz ˆ ̺ st } of system B has a monotone decreas-ing behaviour as a function of the interaction strength Ω, W (cid:144) Γ A X Σ` z \ st BA FIG. 1: The population inversion of system A and B in thesteady state as a function of Ω /γ A . We set γ A /γ B = 1, s B = 0 . s A = 0 .
8. The solid line shows that in system A the average excitation difference between the energy levels | i A and | i A remains independent of interaction strength Ω,while the dashed and the dotted line related to system B shows a decrease. The parameters are ω/γ A = 1 (dottedcurve) and ω/γ A = 10 (dashed curve). and furthermore its dependence on the frequency differ-ence between the energy levels ω governs the rate of thedecrease. The decrease is slowed down by the increase ofthe frequency difference ω . System A returns to the samesteady state independently of the interaction strength Ω.The steady state of system B depends on the interactionstrength which can outplay the relaxation mechanism ofsystem B . This preliminary result already hints thatthis model can only explain physical situations whereasΩ ≪ ω . B. A spin system interacting with a harmonicoscillator
In this section we consider an interaction Hamiltonianbetween a finite and an infinite dimensional system. TheHilbert spaces are H A = C and H B = Γ s ( C ) (sym-metric Fock space) with orthonormal basis {| i A , | i A } for system A and {| n i B } n ∈ N for system B . System A isgoverned by the master equation d ˆ ρ A dt = L A ˆ ρ A (51)= iω A [ˆ ρ A , ˆ σ z ] + γ A (1 − s ) D [ˆ σ − ]ˆ ρ A + γ A s D [ˆ σ + ]ˆ ρ A , where we use the spin operators ˆ σ z = | i A h | A − | i A h | A and ˆ σ − = ˆ σ † + = | i A h | A . The parameter ω A stands forthe frequency difference between the levels | i A , | i A , γ A has the dimension of frequency and 0 ≤ s ≤
1. System B is governed by the master equation d ˆ ρ B dt = L B ˆ ρ B (52)= iω B [ˆ ρ B , ˆ b † ˆ b ] + γ B (1 + ¯ n ) D [ˆ b ]ˆ ρ B + γ B ¯ n D [ˆ b † ]ˆ ρ B , where ˆ b (ˆ b † ) is the annihilation (creation) operator ofthe harmonic oscillator with frequency ω B . ¯ n ≥ γ A hasthe dimension of frequency.The steady states of the two systems are (see AppendixB) ˆ ρ A st = (1 − s ) | i A h | A + s | i A h | A , ˆ ρ B st = 1¯ n + 1 ∞ X n =0 (cid:18) ¯ n ¯ n + 1 (cid:19) n | n i B h n | B . (53)We assume that the two systems interact through theHamiltonian ˆ H I = Ωˆ σ z (cid:16) ˆ b † + ˆ b (cid:17) , (54)and the evolution of the whole system is given by themaster equation d ˆ ̺dt = − i [ ˆ H I , ˆ ̺ ] + L A ˆ ̺ + L B ˆ ̺. (55)A similar type of optical master equation has been sug-gested to describe the properties of defect centers in crys-talline structures [22]. The eigenvectors of ˆ V A = ˆ σ z are | i A and | i A and therefore we can construct the eigen-vectors of the superoperator C A ˆ ρ = [ˆ σ z , ˆ ρ ]:ˆ σ l =1 ,ν =0 = | i A h | A , ˆ σ l = − ,ν =0 = | i A h | A , ˆ σ l =0 ,ν =0 = | i A h | A , ˆ σ l =0 ,ν =1 = | i A h | A . (56)In the case when l = ± L A ˆ ρ = L A ˆ ρ + iω A [ˆ σ z , ˆ ρ ] has only one eigenvalue η l = ± = − γ A / ̺ = Z (cid:18) P , ( α, α ∗ ) P , ( α, α ∗ ) P , ( α, α ∗ ) P , ( α, α ∗ ) (cid:19) | α i B h α | B d α, (57)where d α = d Re { α } d Im { α } and P i,j ( α, α ∗ )( i, j ∈{ , } ) are quasi-distributions of Glauber-Sudarshan type[23, 24]. We are going to focus only on the off-diagonalelement P , ( α, α ∗ ) and using the following identitiesˆ b † (cid:18) e | α | | α i B (cid:19) = ∂∂α (cid:18) e | α | | α i B (cid:19) , ˆ b | α i B = α | α i B , (58)the equation of motion (55) for the representation (57)yields dP , dt = (cid:18) iω A − γ A (cid:19) P , + D [ P , ] − i Ω (cid:18) ∂∂α ∗ + ∂∂α − α ∗ − α (cid:19) P , , (59) where D [ P , ] = − iω B (cid:18) ∂∂α ∗ α ∗ − ∂∂α α (cid:19) P , (60)+ γ B h (cid:18) ∂∂α ∗ α ∗ + ∂∂α α (cid:19) + 2¯ n ∂ ∂α∂α ∗ i P , . In order to obtain the above differential equations anintegration by parts was taken in the variables α and α ∗ .The structure of Eq. (59) suggests a Gaussian ansatz ofthe form P , = e − a ( t )+ b ( t ) α + c ( t ) α ∗ − d ( t ) | α | . (61)Substituting Eq. (61) into the Eq. (59), and a compari-son of the terms proportional to 1, α , α ∗ and | α | leadsto the following set of differential equations:˙ a = − iω A − i Ω ( c + b ) − γ B ¯ n ( bc − d ) − γ B + γ A / , ˙ b = iω B b − i Ω ( d + 2) + γ B / b − γ B ¯ n bd, ˙ c = − iω B c − i Ω ( d + 2) + γ B / c − γ B ¯ n cd, ˙ d = γ B d − γ B ¯ n d . (62)The solution of these equations is very complicated andis subject to the initial conditions. However, we are in-terested in the limit t → ∞ and b ( t ), c ( t ) and d ( t ) tendto a finite limit, but the normalization of Eq. (61) ischaracterized by a ( t ) ≈ γ A t/
2. Therefore in the longtime limit P , ( t → ∞ ) = 0 and the decay to this limit isupper bounded by e − γ A t/ . The same arguments applyto P , ( t ). It is immediate that the steady state of thespin system remains invariantTr B { ˆ ̺ ( t → ∞ ) } = (1 − s ) | i A h | A + s | i A h | A . (63)In order to gain some insight into the evolution of thediagonal elements, we study the average excitation of theharmonic oscillator h ˆ b † ˆ b i = Tr { ˆ b † ˆ b ˆ ̺ ( t ) } and the popu-lation inversion of the spin system h ˆ σ z i = Tr { ˆ σ z ˆ ̺ ( t ) } .Making use of Eq. (55) we can derive the following equa-tions of motion for h ˆ b † ˆ b i and for h ˆ σ z i ddt h ˆ b † ˆ b i = −
2Ω Im[ h ˆ σ z ˆ b i ] − γ B h ˆ b † ˆ b i + γ B ¯ nddt h ˆ σ z ˆ b i = (cid:18) − iω B − γ B − γ A (cid:19) h ˆ σ z ˆ b i − i Ω+ γ A (2 s − h ˆ b i ,ddt h ˆ b i = (cid:18) − iω B − γ B (cid:19) h ˆ b i − i Ω h ˆ σ z i ,ddt h ˆ σ z i = − γ A h ˆ σ z i + γ A (2 s − , (64)where Im[ z ] denotes the imaginary part of the complexnumber z .The obtained equations form a closed set of differentialequations and they can be solved for a set of given initialconditions. However, we are interested in the behaviourof these quantities in the steady state ˆ ̺ st . Setting theleft hand side of Eq. (64) to zero we find the followingsolutions to the population inversion h ˆ σ z i st = Tr { ˆ σ z ˆ ̺ st } = 2 s − , (65)and to the average excitation of the harmonic oscillator h ˆ b † ˆ b i st = ¯ n + (2 s − ( γ B ) + 4 ω B +4(1 − s ) s γ A + γ B γ B (2 γ A + γ B ) + 4 ω B . (66)The average excitation of the harmonic oscillator h ˆ b † ˆ b i st shows an increasing dependency on the interactionstrength Ω. The most peculiar scenario occurs when s = ¯ n = 0, in which case the non-interacting systemsrelax to their ground state ˆ ρ i st = | i i h | i ( i ∈ { A, B } ).While this excitation extracting mechanism dominatesboth systems, due to the interaction Hamiltonian in Eq.(54) the average excitation of the harmonic oscillator forΩ > h ˆ b † ˆ b i st = 4Ω ( γ B ) + 4 ω B . (67) W (cid:144) Γ A X b ` Ö b ` \ st FIG. 2: The average excitation of the harmonic oscillator inthe steady state as a function of Ω /γ A . We set γ B /γ A = 1,¯ n = 0 and s = 1 /
2. The solid line with ω B /γ A = 1 and thedashed line with ω B /γ A = 5 shows an increase in the averageexcitation, but this increase is mediated by the frequency ofthe harmonic oscillator ω B . Fig. 2 shows that the average excitation of the har-monic oscillator starts at the value ¯ n and increases asa function of Ω. This increase depends on the relax-ation parameters of the spin system and of the har-monic oscillator. The found behaviour suggests thatthe model can only describe physical scenarios whenever ω B ≫ Ω p γ A /γ B .We also investigated numerically the decaying mecha-nism of the whole system. We considered the initial con-dition to be | i A + | i A √ | i B . Thereby choosing ¯ n = 0 and Γ A t ÈÈ Ρ ` A H t L - Ρ ` A s t ÈÈ FIG. 3: The norm distance between the evolved state ˆ ρ A ( t ) =Tr B { ˆ ̺ ( t ) } and the steady state of system A as a function of γ A t . We set γ B /γ A = 1, ¯ n = 0, s = 1 / ω A /γ A = ω B /γ A = 10. The solid line with Ω /γ A = 5 and the dashedline with Ω /γ A = 0 shows the difference between coupled anduncoupled system. s = 1 / A h | ˆ ̺ ( t ) | i A and A h | ˆ ̺ ( t ) | i A .We simulated the harmonic oscillator with ten basis vec-tors which are enough because the initial condition isthe steady state for system B and the chosen couplingstrength Ω is not big enough to drive the harmonic os-cillator towards higher excitations than 10 and the con-tributions from these states are negligible. Fig. 3 showsthat the decay e − γ A t/ of the uncoupled system A , i.e.,Ω = 0, is the upper bound for the off-diagonal elementswhich is in agreement with our general findings in Sec.III. C. The optomechanical system
In this section we study more in detail the model whichmotivated our work, see Sec. II. The system is composedof two harmonic oscillators, i.e., two infinite systems, anda direct approach to find out the steady state of the jointsystem requires the solutions of an infinite number ofrecurrence relations. Therefore our result in Sec. III isuseful to identify some properties of the steady state. TheHilbert spaces are H A = Γ s ( C ) and H B = Γ s ( C ) withorthonormal basis {| n i A } n ∈ N for system A (optical mode)and {| n i B } n ∈ N for system B (mechanical mode). Theuncoupled systems are governed by the master equation d ˆ ρ A dt = L A ˆ ρ A , d ˆ ρ B dt = L B ˆ ρ B , (68)with the Liouville operators defined in Eq. (3). We re-call that ˆ a (ˆ a † ) and ˆ b (ˆ b † )are the annihilation (creation)operators of the quantum harmonic oscillators A and B with frequency ω and frequency ν respectively. The coef-ficients with dimensions of frequency characterizing thedissipative dynamics are κ and γ for the optical andmechanical mode respectively. The temperature of thephotonic reservoir is proportional to ¯ n and to ¯ m in thephononic case. The steady states of the separate systemsare (see Appendix B)ˆ ρ A st = 1¯ n + 1 ∞ X n =0 (cid:18) ¯ n ¯ n + 1 (cid:19) n | n i A h n | A , ˆ ρ B st = 1¯ m + 1 ∞ X n =0 (cid:18) ¯ m ¯ m + 1 (cid:19) n | n i B h n | B . (69)The two systems interact with each other through theoptomechanical interaction ˆ H I = g ˆ a † ˆ a (cid:16) ˆ b † + ˆ b (cid:17) and theevolution of the whole system is given by the masterequation (1).Connecting with the formalism in Sec. III we notethat in this case ˆ V A = ˆ a † ˆ a and for all n ∈ N | n i A is aneigenvector with eigenvalue n . The eigenvectors of thesuperoperator C A ˆ ρ = [ˆ a † ˆ a, ˆ ρ ] areˆ σ l,n = ( | n + l i A h n | A , l ≥ | n i A h n + l | A , l < n ∈ N . (70)In the subspaces defined by the eigenvectors l = 0 thesmallest eigenvalue is η l = − κ | l | / B { ˆ ̺ st } = ˆ ρ A st . (71)In this particular example it is hard to work out an al-ternative proof for the above relation as we did in sec-tion IV B. In general Tr A { ˆ ̺ st } 6 = ˆ ρ B st and using the sameapproach as in the previous section we are going to an-alyze the average excitations h ˆ a † ˆ a i = Tr { ˆ a † ˆ a ˆ ̺ ( t ) } and h ˆ b † ˆ b i = Tr { ˆ b † ˆ b ˆ ̺ ( t ) } . The time evolution of the densitymatrix ˆ ̺ ( t ) is governed by Eq. (2) and the followingequations of motion can be derived for h ˆ a † ˆ a i and h ˆ b † ˆ b i ddt h ˆ b † ˆ b i = − g Im[ h ˆ a † ˆ a ˆ b i ] − γ h ˆ b † ˆ b i + γ ¯ m,ddt h ˆ a † ˆ a ˆ b i = κ ¯ n h ˆ b i − γ + 2( κ + iν )2 h ˆ a † ˆ a ˆ b i − ig h (cid:0) ˆ a † ˆ a (cid:1) i ddt h ˆ b i = − iν h ˆ b i − ig h ˆ a † ˆ a i − γ h ˆ b i ,ddt h (cid:0) ˆ a † ˆ a (cid:1) i = κ ¯ n − κ h (cid:0) ˆ a † ˆ a (cid:1) i + κ (4¯ n + 1) h ˆ a † ˆ a i ,ddt h ˆ a † ˆ a i = − κ h ˆ a † ˆ a i + κ ¯ n. (72)The steady state solutions of Eq. (72) are h ˆ a † ˆ a i st = ¯ n, (73) h ˆ b † ˆ b i st = ¯ m + 4¯ n g γ + 4 ν + 4¯ n (¯ n + 1)(2 κ + γ ) g γ [(2 κ + γ ) + 4 ν ] . The above relation shows that the average excitation ofsystem B is increasing as a function of g . The only case when it remains the same is for ¯ n = 0, i.e., system A is relaxed to the ground state. So, whenever system A is relaxed to a finite temperature Gibbs state the aver-age excitation number of system B is increased abovethe average excitation given by its own damping mech-anism. These results resemble the same effect found inSec. IV B. The parameter region where this model worksis constrained by the relation ν ≫ g ¯ n p κ/γ . This is usu-ally fulfilled, because in experimental scenarios ¯ n ≈ Κ t ÈÈ Ρ ` A H t L - Ρ ` A s t ÈÈ FIG. 4: The norm distance between the evolved state ˆ ρ A ( t ) =Tr B { ˆ ̺ ( t ) } and the steady state of system A as a function of κt . The dashed curve corresponds to the free evolution of A ,i.e. g = 0. The full line with g = 0 . κ shows the differencebetween the coupled and the uncoupled system. The initialstate for both systems A and B is a coherent state | α i , with α = 0 .
15. The parameters are: ω = 10 κ , ν = 1 . κ , γ = 0 . κ ,¯ n = 0 .
015 and ¯ m = 0 . Fig. 4 shows the relaxation of the optical mode to thesteady state. We numerically evaluate the trace norm k ˆ ρ A ( t ) − ˆ ρ A st k to measure the distance of the evolved stateto the steady state in Eq. 69. The dashed line shows thefree scenario where the optical mode is solely dampedfrom an initial coherent state | α i with α = 0 .
15. Thefull line shows the reduced density matrix of the opticalfield after the optomechanical interaction. Both modesare taken initially in a coherent state with α = 0 .
15. Onecan note that the trace norm vanishes for both situations,meaning that the optical system always relaxes to thesteady state ˆ ρ A st . Even though the optical steady state isthe same, the relaxation mechanism changes due to theinteraction with the mechanical mode. V. CONCLUDING REMARKS
We have discussed the steady state of two interactingquantum systems each of them having an independentnon-unitary Markovian evolution. The interaction be-tween these two systems was considered in a way thatit commutes with one of the system Hamiltonians andthat the unitary operator generated by it leaves the Li-0ouvillian of the non-unitary dynamics invariant. It isdemonstrated that the steady state of the whole systemhas an interesting property, the trace over the other sys-tem’s degrees of freedom results in a steady state for thecommuting system which is also the steady state of thissystem without the interaction. In plain words, the com-muting system relaxes always to the same steady state,independently of the interaction strength.The method applied here is based on the knowledge ofthe spectrum of the commuting system’s Liouville oper-ator. The eigenvalue equations of these generators areknown only in a few cases. Therefore we restricted ourstatement to the case where the system with the com-muting properties is composed of spins and harmonic os-cillators with damping mechanisms defined in Eqs. (9),and (10), though the presented method is general enoughto deal with other systems as well. In the case when thissystem is finite (more than two levels) the solutions tothe eigenvalue equations can be done either with analyt-ical formulas or with numerical simulations. However,we chose the model of a two-level (spin) system withdecay processes due to its frequent application in var-ious physical scenarios. In the case of infinite systemsour knowledge is restricted to the harmonic oscillator.The decoherence mechanism of a harmonic oscillator canbe described by an optical master equation whose eigen-value equations are known [5]. This is the reason whythis model is our second choice. These damping mod-els with their Liouvillian generators are considered inthe decomposition of the commuting system, althoughwe allow internal interactions between the systems. Forexample: in the case of a spin and a harmonic oscilla-tor we have the two independent damping mechanismsand a Jaynes-Cummings type interaction between them;for two spins and a harmonic oscillator we have three in-dependent damping mechanisms and a Tavis-Cummingstype interaction between them; for a spin and two har-monic oscillators we have again three independent damp-ing mechanisms and a dipole coupling in the Lamb-Dickeregime, an interaction used to describe trapped ions withcenter-of-mass motion inside an optical cavity. In ouranalysis we do not have to specify the other system. Wecomment here that in the case of a Caldeira-Leggett mas-ter equation [26] the solution to the eigenvalue equationis still unknown. The proof uses theorems from the the-ory of one-parameter dynamical semigroups. The reasonto apply these results is the Liouville operator of the har-monic oscillator which is an unbounded operator actingon the Banach spaces of trace-class operators. We usethe fact that dynamical semigroups are contractive semi-groups on the subspace of trace-class operators whoseelements are self-adjoint operators. We slice up this sub-space into further closed subspaces and derive in eachof them an inequality for the norm of the evolved den-sity matrix with the help of the Trotter product formula.These inequalities allow us to identify which elements ofthe density matrix tend to zero in the limit of infinitetime. The projections which splice up the subspace are defined by those vectors which are eigenvectors of theinteraction’s superoperator.We considered three examples of this scenario: a spin-spin, a spin-harmonic oscillator, and a two harmonic os-cillator system. In the first example every detail aboutthe time evolution and the steady state can be workedout analytically due to the small dimension of the Hilbertspace. For the other two systems the situation is morecomplicated and in order to understand the evolution ofthe system we derived closed differential equations forthe averaged excitations (harmonic oscillator) and thepopulation inversion (spin). We focused on the steadystate solutions which show that the interaction increasesthe average excitation number only for the system with-out the commuting properties. In the case of the op-tomechanical model, it turns out that this effect is unde-tectable in current experimental realizations due to theweak coupling between the optomechanical resonator andthe single mode of the radiation field and the low valueof thermal photons in optical frequencies.Finally let us make some comments on the limitationsand possible extensions of our approach. The results ofthis work are based on a class of coupled systems withindependent Markovian master equations, models whichcan explain certain physical scenarios like the actual op-tomechanical experiments. However, it is known thatstrong interaction between two coupled dissipating sub-systems leads to serious defects, whenever independentlyderived non-unitary evolutions are added to the unitaryevolution [27–35]. If these defects are present in the so-lutions, then one must start the theoretical aproach ona microscopic level and derive a new dissipative model.The structure of these master equations depends on theinteraction between the two subsystems. For example, inRef. [27] the two subsystems are harmonic oscillators andtheir creation and annihilation operators are mixed up inthe generator of the master equation, which leads to anopen problem of determining the spectrum of this Liou-villian. Our statements are based on the requirement ofindependent master equations in order to keep the non-commuting system arbitrary, however when the details ofthis system are known and they fulfill the assumptions ofSec. III then the steady state of the commuting systemremains invariant, see for example the strongly coupleddissipative optomechanical model [16].As a final word we think that the presented analysismay offer interesting perspectives for actual and futuremodels with Markovian master equations.
Acknowledgement
This work is supported by the BMBF project Q.com.We are grateful to F. Sokoli, N. Trautmann, H. Frydrych,and G. Alber for stimulating discussions.1
Appendix A: Eigensystem of the Liouville operator
In this appendix we present the eigensystem of adamped two-level system and a damped harmonic os-cillator. Due to the fact that the Hilbert space of a har-monic oscillator is infinite dimensional we are going tosummarize some basic facts about the adjoint of a Ba-nach space operator, a notion needed for the definitionof the eigenvectors.
1. Eigensystem of L S In the main text we considered the following masterequation for a spin system L S ˆ ρ = γ (1 − ¯ m ) (cid:0) ˆ σ − ˆ ρ ˆ σ + − { ˆ σ + ˆ σ − , ˆ ρ } (cid:1) + γ ¯ m (cid:0) ˆ σ + ˆ ρ ˆ σ − − { ˆ σ − ˆ σ + , ˆ ρ } (cid:1) . (A1)The corresponding Liouville operator is non-Hermitian,but it can still be diagonalized leading to the dual eigen-value equation L S ˆ ρ i = λ i ˆ ρ i , L † S ˇ ρ i = λ ∗ i ˇ ρ i . (A2)The elements ˆ ρ i are the four right eigenvectors given byˆ ρ = ˆ I + m − ˆ σ z , ˆ ρ z = ˆ σ z , ˆ ρ ± = ˆ σ ± , (A3)with the corresponding eigenvalues λ = 0 , , λ z = − γ, λ ± = γ/ . (A4)The adjoint of the Liouville operator is L † S ˆ ρ = γ (1 − ¯ m ) (cid:0) ˆ σ + ˆ ρ ˆ σ − − { ˆ σ − ˆ σ + , ˆ ρ } (cid:1) + γ ¯ m (cid:0) ˆ σ + ˆ ρ ˆ σ − − { ˆ σ − ˆ σ + , ˆ ρ } (cid:1) , (A5)and its eigenvectors are the following four elementsˇ ρ = ˆ I , ˇ ρ z = ˆ σ z − m − ˆ I , ˇ ρ ± = ˆ σ ± . (A6)The eigenvectors of L S and L † S belong to the space of2 × M ( C ) and form a biorthogonal systemas it can be seen that they fulfill the relation h ˇ ρ i , ˆ ρ j i = Tr n ˇ ρ † i ˆ ρ j o = δ i,j (A7)with the Kronecker delta for all i, j ∈ { , z, + , −} .
2. Eigensystem of L HO The master equation describing a Markovian dampedharmonic oscillator can be written in the following form L HO ˆ ρ = γ (¯ n + 1) (cid:0) ˆ a ˆ ρ ˆ a † − { ˆ a † ˆ a, ˆ ρ } (cid:1) + γ ¯ n (cid:0) ˆ a † ˆ ρ ˆ a − { ˆ a ˆ a † , ˆ ρ } (cid:1) . (A8) The eigenvalue equations of L HO and L † HO read L HO ˆ ρ n,k = λ n,k ˆ ρ n,k , L † HO ˇ ρ n,k = λ ∗ n,k ˇ ρ n,k . (A9)The first equation is solved by the following right eigen-vectors (see Ref. [5])ˆ ρ n,k = ˆ a † | k | + k ( − n (¯ n +1) k +1 : L ( | k | ) n (cid:16) ˆ a † ˆ a ¯ n +1 (cid:17) e − ˆ a † ˆ a ¯ n +1 : ˆ a † | k |− k (A10)where L ( k ) n is a generalized Laguerre polynomial and: f (ˆ a † ˆ a ) : is the normal ordering of f , a function of thenumber operator ˆ a † ˆ a . The eigenvalues are λ n,k = − γ (cid:18) n + | k | (cid:19) , n ∈ N , k ∈ Z . (A11) L † HO is a representation of the adjoint of L HO with thehelp of the trace as a linear functional and can be foundto be L † HO ˆ ρ = γ (¯ n + 1) (cid:0) ˆ a † ˆ ρ ˆ a − { ˆ a † ˆ a, ˆ ρ } (cid:1) + γ ¯ n (cid:0) ˆ a ˆ ρ ˆ a † − { ˆ a ˆ a † , ˆ ρ } (cid:1) , (A12)and its eigenvectors are given byˇ ρ n,k = ( − ¯ n ) n n !(¯ n +1) n ( n + k )! ˆ a | k | + k : L ( | k | ) n (cid:16) ˆ a † ˆ a ¯ n (cid:17) : ˆ a | k |− k . (A13)It can be verified that they are orthogonal with respectto the inner product h ˇ ρ n,k , ˆ ρ n ′ ,k ′ i = Tr n ˇ ρ † n,k ˆ ρ n ′ ,k ′ o = δ n,n ′ δ k,k ′ , (A14)and this is consistent with hL † HO ˇ ρ n,k , ˆ ρ n,k i = h ˇ ρ n,k , L HO ˆ ρ n,k i = λ n,k . (A15)
3. Comments on the eigensystem of L H0 Banach space structure.-
In the case of the dampedharmonic oscillator one has to deal with an infinite di-mensional Hilbert space. The Liouville operator is a mapon the set of trace class operators which do not form aHilbert space, but do form a Banach space. Geomet-rical problems are inherent in the concept of a Banachspace and therefore a slightly different notion of adjointoperators is needed.Let
E, F be Banach spaces and if T : E → F is abounded linear operator then we define its adjoint by T ′ : F ′ → E ′ , T ′ ( ϕ ) → ϕ ◦ T, ϕ ∈ F ′ (A16)where T ′ is mapping bounded linear functionals on F tobounded linear functionals on E and E ′ ( F ′ ) is the dualof E ( F ). The set of trace class operators on a Hilbertspace H is defined as B ( H ) := { ˆ x ∈ B ( H ) : Tr n √ ˆ x † ˆ x o < ∞} . (A17)2If L : B ( H ) → B ( H ) is a bounded linear operator thenthe adjoint L ′ is defined on the dual B ( H ) ′ . Due to thefact that the set of trace class operator form a two-sidedideal in B ( H ) (see Theorem VI.19 in Ref. [36]) we havean isometric isomorphism B ( H ) ′ ∼ = B ( H ) . (A18)Thus, we may translate the adjoint L ′ to an operator L † : B ( H ) → B ( H ) satisfyingTr { ( L † ˆ y ) † ˆ ρ } = Tr { ˆ y † ( L ˆ ρ ) } (A19)for all ˆ y ∈ B ( H ) and ˆ ρ ∈ B ( H ).We can already conclude that the eigenvectors of L ′ HO can be understood as the linear functional ϕ for a fixedˇ ρ n,k as ϕ [ˇ ρ n,k ](ˆ ρ ) = Tr { ˆ ̺ † n,k ˆ ρ } , ∀ ˆ ρ ∈ Dom( L HO ) . (A20) Domain.-
When dealing with infinite dimensional sys-tems, the domain of the operators is a non-trivial con-cept that has to be handled carefully as we are goingto elucidate in the following lines. The Hilbert spaceof the harmonic oscillator is the symmetric Fock space H = Γ s ( C ). Every element can be expressed in terms ofthe orthonormal basis {| n i} n ∈ N as x = ∞ X n =0 x n | n i ∈ H , ∞ X n =0 | x n | < ∞ . (A21)The creation and the annihilation operators (ˆ a † , ˆ a ) havethe domainDom(ˆ a † ) = Dom(ˆ a ) = n x : ∞ X n =1 n | x n | < ∞ o (A22)which is dense in H . The number operator ˆ a † ˆ a has adifferent and more restricted domain given byDom(ˆ a † ˆ a ) = n x : ∞ X n =1 n | x n | < ∞ o . (A23)In order to exemplify this difference, let us define thefollowing three vectors in H : x = ∞ X n =0 n | n i , x = ∞ X n =0 n / | n i ,x = ∞ X n =0 n | n i . (A24)These vectors have the following properties x / ∈ Dom(ˆ a ) , x / ∈ Dom(ˆ a † ˆ a ) ,x ∈ Dom(ˆ a ) , x / ∈ Dom(ˆ a † ˆ a ) ,x ∈ Dom(ˆ a ) , x ∈ Dom(ˆ a † ˆ a ) . (A25)This is because even when x is a normalized vector,ˆ ax / ∈ H and ˆ a † ˆ ax / ∈ H are not normalized anymore and therefore they are not in H . In the case of x , itcan be noted that ˆ ax ∈ H but ˆ a † ˆ ax / ∈ H . With thisobservation we can conclude that Dom(ˆ a † ˆ a ) ⊂ Dom(ˆ a ) ⊂H .In the case of the Liouville operator and for the sake ofsimplicity, let us consider the case of zero a temperaturebath (¯ n = 0). In this case the Liouvillian reads L ˆ ρ = ˆ a ˆ ρ ˆ a † − (cid:0) ˆ a † ˆ a ˆ ρ + ˆ ρ ˆ a † ˆ a (cid:1) . (A26)In order to show that not all elements of B ( H ) are in thedomain of L , we consider the following density matrixˆ ρ = 6 π ∞ X n,m =1 nm | n ih m | , (A27)The normalization condition is fulfilled as P ∞ n =1 n − = π /
6. Applying L one getsˆ ρ ⋆ = L ˆ ρ = 6 π ∞ X n =0 √ n + 1 ( | n ih | + | ih n | ) ++ 6 π ∞ X n,m =1 p ( n + 1)( m + 1) − n + m nm ! | n ih m | . Now, let us remember that the normed space of Hilbert-Schmidt operators is a subset of all bounded linear oper-ators defined on a Hilbert space H and at the same timeit contains the set trace class operators, i.e. B ( H ) ⊆ B ( H ) ⊆ B ( H ) . (A28)In order to determine if ˆ ρ ⋆ is a Hilbert-Schmidt operatorwe evaluateTr { ˆ ρ † ⋆ ˆ ρ ⋆ } = ∞ X n,m =1 p ( n + 1)( m + 1) − n + m nm ! + 2 ∞ X n =2 n + 1 (A29)which is clearly diverging. Therefore ˆ ρ ⋆ / ∈ B ( H ) im-plying that ˆ ρ ⋆ / ∈ B ( H ). So L in Eq. (A26) cannot bedefined on the whole space B ( H ), but only on its domainDom( L ). The properties of Dom( L HO ) are not studiedhere, but one can at least state thatDom( L HO ) ⊂ B ( H ) , which implies that the bounded linear functionals onDom( L HO ) regarding to the norm of B ( H ) are more nu-merous than the bounded linear functionals on B ( H ). Linear independence and completeness.-
The eigen-vectors ˆ ρ n,k due to the normal ordered form : L ( k ) n ( z ˆ a † ˆ a )e − z ˆ a † ˆ a : with z ∈ (0 ,
1] are trace class oper-ators. Nonetheless, (ˆ a † ˆ a ) m ˆ ρ n,k is also a trace class op-erator for any finite value of m . All the eigenvalues of3 L HO are distinct which implies that the eigenvectors ˆ ρ n,k are linearly independent. We have a countable infiniteset of eigenvectors whose linear independence propertyis equivalent to show that all possible finite selection ofeigenvectors are linearly independent. In order to showthis, let { ˆ ρ i (1) , ˆ ρ i (2) , ..., ˆ ρ i ( k ) } be a finite set of eigenvec-tors with corresponding eigenvalues { λ i (1) , λ i (2) , ..., λ i ( k ) } and let c , c , ..., c k be scalars such that c ˆ ρ i (1) + c ˆ ρ i (2) + ... + c k ˆ ρ i ( k ) = 0 . (A30) L HO has non-degenerate eigenvalues and applying N ( L HO ) = ( L HO − λ i (2) )( L HO − λ i (3) ) ... ( L HO − λ i ( k ) )to the above equation results that all terms except thefirst are annihilated, resulting c N ( λ i (1) )ˆ ρ i (1) = 0 , N ( λ Ai (1) ) = 0 , (A31)and since ˆ ρ i (1) = 0 it follows that c = 0. Similarly, c = ... = c k = 0, which proves the linear independenceof { ˆ ρ i (1) , ˆ ρ i (2) , ..., ˆ ρ i ( k ) } . This proof can be applied to allpossible finite sets of eigenvectors. The eigenvectors arenot orthogonal in regard to the Hilbert-Schmidt innerproduct. The set of eigenvectors has maximal possiblecardinality, so { ˆ ρ n,k } n ∈ N ,k ∈ Z is a basis. Spectrum.-
Finite rank operators are dense in the setof trace class operators (see the Corollary after TheoremVI.21 in Ref. [36]). The eigenvectors ˆ ρ n,k are finite rankoperators and they form a basis which implies that L HO is a closed densely defined operator. The point spec-trum of L † HO coincides with the point spectrum of L HO and this means that L HO has no residual spectrum (seeProposition 1 .
12 in Chap. IV of Ref. [21]). Now weknow the point spectrum and the residual spectrum, butwe have not yet discussed all possible approximate eigen-values. λ is an approximate eigenvalue if there exists asequence { ˆ x n } n ∈ N ⊂ Dom( L HO ) such that k ˆ x n k = 1and lim n →∞ kL HO ˆ x n − λ ˆ x n k = 0. The latter equationcan be expanded in regard to the basis { ˆ ρ n,k } n ∈ N ,k ∈ Z andyields lim n →∞ k X n,k ( λ n,k − λ ) Tr { ˆ ̺ n,k ˆ x n } ˆ ρ n,k k = 0 , which is fulfilled only when λ n,k = λ , meaning that L HO has only a point spectrum. Appendix B: Steady state of a damped harmonicoscillator
Let us consider the model of a damped harmonic os-cillator given by the master equation d ˆ ρdt = − γ (cid:0) ˆ a † ˆ a ˆ ρ − a ˆ ρ ˆ a † + ˆ ρ ˆ a † ˆ a (cid:1) − γ (cid:0) ˆ a ˆ a † ˆ ρ − a † ˆ ρ ˆ a + ˆ ρ ˆ a ˆ a † (cid:1) (B1)where ˆ a (ˆ a † ) is the annihilation (creation) operator ofthe harmonic oscillator. γ and γ have the dimension of a frequency. To ensure that Eq. (B1) describes theevolution of a density matrix we require the condition γ < γ .To investigate in more detail the steady state of thedamping, we consider the density matrix of the field innumber state representation | n i ( n ∈ N )ˆ ρ st = ∞ X n,m =0 ρ n,m | n ih m | . (B2)We set the left-hand side of Eq. (B1) to zero, and takingmatrix elements of the density matrix we get0 = − γ m + n ) ρ n,m + γ p ( n + 1)( m + 1) ρ n +1 ,m +1 − γ m + 1 + n + 1) ρ n,m + γ √ nmρ n − ,m − . (B3)The matrix elements are coupled via a second order re-currence relation to their nearest neighbours in each ofthe diagonals of the density matrix.First, we consider the main diagonal ( m = n ) and Eq.(B3) simplifies to ρ n +1 ,n +1 = γ n + γ ( n + 1) γ ( n + 1) ρ n,n − γ nγ ( n + 1) ρ n − ,n − . (B4)The solution can be determined by induction, thus yield-ing ρ n,n = ( γ /γ ) n ρ , , n ≥ , (B5)The normalization condition Tr { ˆ ρ } = 1 results in ρ , ∞ X n =0 (cid:18) γ γ (cid:19) n = 1 . (B6)With the geometric series ∞ X n =0 (cid:18) γ γ (cid:19) n = 11 − γ /γ , (B7)one finds that ρ , = ( γ − γ ) /γ .Now, we consider only the upper diagonals m > n due to the self-adjoint property of the density matrix ˆ ρ s .Defining l = m − n ≥ γ = ǫγ with ǫ ∈ [0 ,
1) Eq.(B3) takes the form ρ n,n + l = r ǫl,n ρ ,l ,r ǫl,n +1 = ( f nl + ǫg nl ) r ǫl,n − ǫh nl r ǫl,n − (B8)where we introduced the following abbreviations f nl = ( n + l/ / p ( n + 1)( n + l + 1) ,g nl = ( n + 1 + l/ / p ( n + 1)( n + l + 1) ,h nl = p n ( n + l ) / p ( n + 1)( n + l + 1) , (B9)4which are all positive and fulfill the following inequalities f nl − h nl > , g nl > , ∀ n ≥ , l ≥ . (B10)For the first few terms it can be shown that they havesolutions of the form r ǫl,n = n X i =0 a n,li ǫ i , (B11)with coefficients obeying the properties a n,li > , a n − ,li < a n,li +1 , ∀ n ≥ i ≥ . (B12)We apply a proof by induction to show that r ǫl,n +1 takesthe form of Eq. (B11) with the conditions of Eq. (B12).First we use Eq. (B8) which results in r ǫl,n +1 = n X i =0 f nl a n,li ǫ i + n X i =0 g nl a n,li ǫ i +1 − n − X i =0 h nl a n − ,li ǫ i +1 , (B13)and translates to equations for the coefficients: a n +1 ,l = f nl a n,l , (B14) a n +1 ,li = f nl a n,li + g nl a n,li − − h nl a n − ,li − , i ≥ . Now, applying the conditions of Eqs. (B10), and (B12)we find that a n +1 ,li > , a n,li < a n +1 ,li +1 (B15)and it follows immediately that every r ǫl,n can be cast inthe form of Eq. (B11).As all coefficients are positive, the sum is always largerthan the zeroth order coefficient r l,n = n − Y j =0 f nj = Γ (cid:0) n + l (cid:1) Γ (cid:0) l (cid:1) s l ! n !( n + l )! , (B16)which we have expressed in terms of the gamma functionΓ( z ), and multiplied by ρ ,l solves the first order recur-rence relation obtained from the Eq. (B3) in the casewhen ǫ = 0 which implies γ = 0.These considerations allow us to set a lower bound tothe following sum: | ∞ X n =0 ρ n,n + l | ≥ | ρ ,l | ∞ X n =0 r l,n ≥ | ρ ,l | ∞ X n =0 r ,n , (B17)as the relation r l,n ≥ r ,n can be verified by noting thateach term in the product of Eq. (B16) has the property f nl ≥ f n which can be checked from Eq. (B9). Stirling’sformula in an inequality form [37, 38] √ πn (cid:16) ne (cid:17) n < n ! < √ πn (cid:16) ne (cid:17) n e n (B18)and the relation Γ( n + 1 /
2) = √ π (2 n )! / (4 n n !) yield r ,n > e − n p πn ( n + 1) > √ πe / ( n + 1) . Thus, we have found that all the solutions to Eq. (B3)for l = 0 fulfill the inequality | ∞ X n =0 ρ n,n + l | > | ρ ,l | ∞ X n =1 √ π e / ( n + 1) ! . (B19)The right-hand side of the above inequality is divergentfor | ρ ,l | > H with an orthonormal basis ( e i ). The set of Hilbert-Schmidt operators B ( H ) with the Hilbert-Schmidt innerproduct h ˆ X, ˆ Y i = X i h ˆ Xe i , ˆ Y e i i = Tr { ˆ X † ˆ Y } , ∀ ˆ X, ˆ Y ∈ B ( H )(B20)is a Hilbert space and it is also a two-sided ideal in B ( H )(see Theorem VI.22 in Ref. [36]): || ˆ A ˆ X || , || ˆ X ˆ A || < ∞ , ˆ X ∈ B ( H ) , ˆ A ∈ B ( H ) . (B21)In our case the Hilbert space H = Γ s ( C ) and any den-sity matrix ˆ ρ ∈ B ( H ) ( || ˆ ρ || ≤
1) and ˆ ρ / ∈ B ( H )( || ˆ ρ / || = 1).Now let us consider the right shift operatorˆ S R | n i = | n + 1 i , n ≥ , ˆ S R ∈ B ( H ) . (B22)Any l th power of the right shift operator ˆ S lR is also abounded operator. Then | Tr( ˆ S lR ˆ ρ s ) | = | ∞ X n =0 ρ n,n + l | (B23)and ˆ S lR ˆ ρ s ∈ B ( H ). The Cauchy-Schwarz-Bunyakovskiinequality for the Hilbert-Schmidt inner product gives: | Tr( ˆ X † Y ) | = |h ˆ X, ˆ Y i| ≤ || X || || Y || , ∀ ˆ X, ˆ Y ∈ B ( H ) . (B24)We observe that for the Hilbert-Schmidt operators ˆ X =ˆ ρ / and ˆ Y = ˆ S lR ˆ ρ / | Tr( ˆ S lR ˆ ρ s ) | ≤ || ˆ ρ / || || ˆ S lR ˆ ρ / || = 1 , (B25)where we used the relation (cid:16) ˆ S lR (cid:17) † ˆ S lR = ˆ I , the identityoperator on H . Therefore, Eq. (B19) obviously leads to1 ≥ | ∞ X n =0 ρ n,n + l | > | ρ ,l | + 1 √ π e / | ρ ,l | ∞ X n =1 n + 1 . (B26)which is true if and only if ρ ,l = 0. If ρ ,l = 0, thenall the upper diagonal elements of ˆ ρ st are zero. It isimmediate from the self-adjoint property of ˆ ρ st that alllower diagonal elements are also equal to zero. Finally,we can conclude that ˆ ρ st has only non-zero elements inthe diagonal which are given in Eq. (B5).5 [1] R. Horodecki, P. Horodecki, M. Horodecki, and K.Horodecki, Rev. Mod. Phys. , 865 (2009).[2] M. A. Nielsen and I. L. Chuang, Quantum Computationand Quantum Information (Cambridge University Press,Cambridge, UK, 2000).[3] L. Aolita, F. de Melo, and L. Davidovich, Rep. Prog.Phys. , 042001 (2015).[4] W. Weidlich, H. Risken and H. Haken, Z. Phys. , 396(1967).[5] H.-J. Briegel and B.-G. Englert, Phys. Rev. A , 3311(1993).[6] H.-P. Breuer and F. Petruccione, The theory of openquantum systems (Oxford University Press, Oxford2002).[7] J. M. Torres, Phys. Rev. A , 052133 (2014).[8] T. Prosen, Phys. Rev. Lett. , 030603 (2014).[9] V. V. Albert and L. Jiang, Phys. Rev. A , 022118(2014).[10] A. Chimonidou and E. C. G. Sudarshan, Phys. Rev. A , 032121 (2008).[11] I. Dorofeyev, Can. J. Phys. , 537 (2013).[12] P. W. Milonni and P. L. Knight, Phys. Rev. A , 1096(1974).[13] S. Mancini, V. Giovannetti, D. Vitali, and P. Tombesi,Phys. Rev. Lett. , 120401 (2002).[14] C. K. Law, Phys. Rev. A , 433 (1994).[15] C. K. Law, Phys. Rev. A , 2537 (1995).[16] D. Hu, S.-Y. Huang, J.-Q. Liao, L. Tian, and H.-S. Goan,Phys. Rev. A , 013812 (2015).[17] A. Kossakowski, Rep. Math. Phys. , 247 (1972).[18] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J.Math. Phys. , 821 (1976).[19] G. Lindblad, Commun. Math. Phys. , 119 (1976). [20] R. S. Ingarden and A. Kossakowski, Ann. Phys. , 451-485 (1975).[21] K.-J. Engel and R. Nagel, One-Parameter Semigroups forLinear Evolution Equations (Springer New York, 2000).[22] R. Betzholz, J. M. Torres and M. Bienert, Phys. Rev. A , 063818 (2014).[23] R. J. Glauber, Phys. Rev. Lett. , 84 (1963).[24] E. C. G. Sudarshan, Phys. Rev. Lett. , 277 (1963).[25] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt,Rev. Mod. Phys. , 1391 (2014).[26] A. O. Caldeira and A. J. Leggett, Ann. Phys. (N.Y.) ,374 (1983).[27] D. F. Walls, Z. Phys.