Third quantization: a general method to solve master equations for quadratic open Fermi systems

Third quan tization: a general metho d to solv e master equation s for quadratic op en F ermi systems T oma ˇ z Prosen Department of ph ys ic s, FMF, University of Ljubljana, Jadransk a 19 , SI-10 00 Ljubljana, Slov enia Abstract. The Lindblad master equation for an arbitr a ry quadr atic sys tem of n fermions is s o lved explicitly in terms o f diagonalizatio n of a 4 n × 4 n matrix, provided that all Lindblad bath o p erators are line ar in the fermionic v ariables. The metho d is applied to the explicit construction of no n-equilibrium steady states and the ca lculation of asymptotic relax ation rates in the fa r from equilibrium pr oblem of heat a nd spin transp ort in a near est neighbor Heisenberg X Y spin 1 / 2 chain in a transverse magnetic field. P A CS num ber s: 02.30.Ik, 03.65.Yz, 0 5.30.Fk, 75.10.Pq A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 2 1. In tro duction Understanding time ev olution of an op en quantum system of many in teracting particles is of primary imp ortance in fundamen tal problems of quan tum phy sics, suc h as decoherence [1, 2] and closely related quan tum measuremen t problem [3, 4], quantum computation [5, 6], or the problem of computation of non-e quilibrium ste ady states (NESS) in q uantum statistical mec hanics [7 , 8, 9]. Ev en though application of the metho ds of Hamiltonian dynamical sy stems and ergodic theory to quan tum systems out of equilibrium give s many promising results [10, 11, 12], the field of op en qu antum systems is still lack ing non- trivial explicitly solv able mo dels, as compared to studies of closed (isolated) quantum systems where we know a la rge b o dy of the so-called c om p letely inte gr able systems [1 3 , 14]. Examples of explicitly solv able mo dels of master equations for op en quan tum systems are limited to quite restricted mo dels of a single particle, single spin o r ha rmonic oscillators (see e.g. [15, 16, 17]). In t his pap er w e sho w that the generator of the master equation of a general quadratic system of n in teracting fermions whic h are coupled to a g eneral set of Mark o vian baths, sp ecified in terms of Lindblad op erato rs which are linear in the fermionic v aria bles - the so called quan tum Liouville sup er-op erator (or Liouvillean) - can b e explicitly diagonalized in terms of 2 n norm al master mo des , i.e. an ticomm uting sup er-op erators w hic h act on the F o ck space of densit y op erators. This can b e understo o d a s a complex (non-canonical) v ersion of the Bo goliub ov transformation [18 ] lifted on the op erator space, and has v ery p o w erful consequences: (i) The NESS of the master equation can b e understo o d as the ‘ground state’ normal mo de of the L io uvillean, whereas the long time r elaxatio n rate is give n b y the eigenmo de closest to the real axis. (ii) The co v ariance matrix of NESS can b e computed explicitly in terms of the eigen v ectors of 4 n × 4 n an tisymmetric complex matrix. It can b e used to completely express ph ysical observ ables in NESS, such as particle/spin densities, curren ts, etc. W e demonstrate the p o w er of this nov el method b y applying it to the problem of heat and spin tr a nsp ort far f rom equilibrium in nearest neigh b or Heisen b erg XY spin 1 / 2 c hains sub ject to a t r a nsv erse magnetic field. As a result w e repro duce b al listic tr ans p ort in the i nte gr able spatially homogeneous case (see e.g. [19, 20, 21, 22, 23, 24, 12] for related recen t s tudies of quantum thermal conductivit y in one dimension), and predic t ide a l ly insulating b eha viour (at all temp eratures) in a disordered case of spatially random in teractions/field. Apart from obtaining n umerical results which go by fa r b eyond what w as so far accessible b y direct n umerical solution of the many-particle Lindblad equation, either directly or b y means of quan tum tra jectories [15], w e also obtain t wo notable analytical results in the spatially homo g eneous (non-disordered) case: (i) W e compute the sp ectral ga p of the Liouvillean i.e . the rate of of relaxatio n to t he NESS and show that it scales with t he inv erse cub e of the chain length. (ii) W e construct evanesc ent normal master mo des o f the Liouvillean, for long c hains, b y whic h w e explain quan titativ ely the exp o nen tial fallo ff of energy densit y or temp erature profiles near the bath con tacts. A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 3 The paper is organized as follo ws. In sec tion 2 w e shall outline a g eneral metho d for the diagonalizatio n of the Liouvillean sup er-op erator for finite quadratic op en F ermi systems and an explicit construction o f NESS. In section 3 we illustrate the metho d b y working out a simple ex ample of a single fermion or a tw o lev el quantum system in a bath. In se ction 4 w e demonstrate the usefulness of the new metho d b y applying it to quantum transp ort in XY spin c hains. In section 5 we discuss p ossible alternativ e applications and generalizations of the metho d a nd reach some conclusions. 2. General metho d o f solution for the Lindblad equation The general master equation g ov erning time ev olution of the densit y matrix ρ ( t ) of an op en quantum system, preserving trace and p ositivit y of ρ , can b e written in the Lindblad form [25, 17] as (we set ~ = 1) d ρ d t = ˆ L ρ := − i[ H, ρ ] + X µ  2 L µ ρL † µ − { L † µ L µ , ρ }  (1) where H is a Hermitian op erator (Hamiltonian), [ x, y ] := xy − y x , { x, y } := x y + y x , and L µ are arbitrary o p erators represen ting couplings to differen t baths (at p ossibly differen t v alues of t hermo dynamic p o ten tials). W e are no w going to describ e a general metho d of explicit solution of (1) for a quadr atic system o f n fermions (or spins 1 / 2 ) with line ar bath op erators H = 2 n X j,k =1 w j H j k w k = w · H w (2) L µ = 2 n X j =1 l µ,j w j = l µ · w (3) where w j , j = 1 , 2 , . . . , 2 n , a re abstract Hermitian Ma jorana op erators satisfying the an ti-commutation relations { w j , w k } = 2 δ j,k j, k = 1 , 2 , . . . , 2 n (4) and generate a Clifford alg ebra. Th us, 2 n × 2 n matrix H can be c hosen to be an tisymmetric H T = − H . Throughout this pap er x = ( x 1 , x 2 , . . . ) T will designate a v ector (column) of appropriate scalar v alued or op erator v alued sym b ols x k . Tw o notable examples, to whic h our formalism is immediately applicable, are: (i) canonical fermions c m , m = 1 , 2 , . . . , n , w 2 m − 1 = c m + c † m w 2 m = i( c m − c † m ) (5) or (ii) spins 1 / 2 with canonical Pauli op erat o rs ~ σ m , m = 1 , 2 , . . . , n , w 2 m − 1 = σ x m Y m ′ 0 . Then, the rate of expo nen tial relaxation to NESS is giv en by the sp e ctr al gap ∆ of the Liouvillean whic h equals ∆ = 2 R e β 2 n . A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 10 Theorem 3: Assume that the rapidit y sp ectrum do es not con tain 0, i.e. β 2 n 6 = 0. Then the exp ectation v alue of an y quadratic obse rv able w j w k in a (unique) NESS can b e explicitly computed as h w j w k i NESS = δ j,k + h 1 | ˆ c j ˆ c k | NESS i = (46) = δ j,k + 1 2 2 n X m =1  v 2 m, 2 j − 1 v 2 m − 1 , 2 k − 1 − v 2 m, 2 j v 2 m − 1 , 2 k − i v 2 m, 2 j v 2 m − 1 , 2 k − 1 − i v 2 m, 2 j − 1 v 2 m − 1 , 2 k  (47) The statements of theorems 1 and 2 simply follow from exact and explicit sp ectral decomp osition (42,43,44). The pro of of theorem 3 is also straigh tfo rw ard: The first expression (46) follow s from the definition of the annihilation maps (9) and the exp licit represen tation of the densit y op erator of NESS, ρ NESS , in the cano nical basis P α . The second, ve ry useful equalit y (47) is then obtained b y expressing ˆ c j thru (24) in terms of NMM maps (35) and using the annihilat ion relations (39). The quadratic correlator of theorem 3 co v ers man y ph ysically in teresting observ ables suc h as densities or curren ts. Ho w ev er if one needs exp ectation v alues of more general observ ables, e.g. an expectation v alue of a hig h order monomial h P α i NESS = h 1 | ˆ c α 1 1 ˆ c α 2 2 · · · ˆ c α 2 n 2 n | NESS i , then one ma y use a Wick the or em a nd rewrite it as a sum of pro ducts of pair-wise con tractions (46). 3. T rivial example: A single fermion in a bath In order to illustrate the metho d and demonstrate conv enience of the results deriv ed in the previous section we first w o rk out a simple example of a single fermion n = 1 (or equiv alen tly , an arbitra ry qubit, a t w o-leve l quan tum system), in a thermal bath. W e tak e the most general single fermion Ha milto nian H = − i hw 1 w 2 + const = 2 hc † c + const ′ and the follo wing bath op erators (see e.g. [16, 29]) L 1 = 1 2 p Γ 1 ( w 1 − i w 2 ) = p Γ 1 c L 2 = 1 2 p Γ 2 ( w 1 + i w 2 ) = p Γ 2 c † (48) where the ratio of coupling constan ts determine the bat h temp erature T , Γ 2 / Γ 1 = exp( − 2 h/T ). Leav ing out the details of a straigh tforward calculation, simply follo wing the steps of the previous section, we arriv e at the follow ing shap e matrix of t he Liouvillean (26) A = − h R + B Γ + , Γ − A 0 = Γ + (49) where R :=     0 0 1 0 0 0 0 1 − 1 0 0 0 0 − 1 0 0     B Γ + , Γ − :=     0 i 2 Γ + − i 2 Γ − 1 2 Γ − − i 2 Γ + 0 1 2 Γ − i 2 Γ − i 2 Γ − − 1 2 Γ − 0 i 2 Γ + − 1 2 Γ − − i 2 Γ − − i 2 Γ + 0     (50) A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 11 and Γ ± := Γ 2 ± Γ 1 . F urther, w e compute NMM rapidities β 1 , 2 = 1 2 Γ + ± i h and the eigen v ector matrix V =     Γ − Γ + − 1 i Γ − Γ + + i − i Γ − Γ + + i Γ − Γ + + 1 − 1 4 − i 4 − i 4 1 4 Γ − Γ + + 1 i Γ − Γ + − i i Γ − Γ + + i − Γ − Γ + + 1 1 4 i 4 − i 4 1 4     (51) Then, using theorem 3 we compute t he exp ectation v alue of o ccupation n umber h c † c i = 1 2 − i 2 h w 1 w 2 i = Γ 2 / (Γ 1 + Γ 2 ) whic h is what w e exp ect in canonical equilibrium. 4. Non-trivial example: transp ort in quan tum spin cha ins Here w e w ork out a ph ysically more in teresting example, namely w e construct NESS for the magnetic and heat transp or t of a Heisen b erg XY spin 1 / 2 c hain, with arbitrary – either homogeneous or p ositionally dependent (e.g. disordered) – nearest neigh b our in teraction H = n − 1 X m =1  J x m σ x m σ x m +1 + J y m σ y m σ y m +1  + n X m =1 h m σ z m (52) whic h is coupled t o t wo thermal/magnetic baths at the ends o f the c hain, generated by t w o pairs of canonical Lindblad op erat o rs [29] (similar to (48)) L 1 = 1 2 q Γ L 1 σ − 1 L 3 = 1 2 q Γ R 1 σ − n L 2 = 1 2 q Γ L 2 σ + 1 L 4 = 1 2 q Γ R 2 σ + n (53) where σ ± m = σ x m ± i σ y m and Γ L , R 1 , 2 are positive coupling constan ts related to bath temp eratures/magnetizations, for example if spins w ere non-interacting t he bath temp eratures T L , R w ould b e give n with Γ L , R 2 / Γ L , R 1 = exp( − 2 h 1 , n /T L , R ). Applying t he inv erse of Jordan- Wig ner transformation (6), σ x m = ( − i) m − 1 Q 2 m − 1 j =1 w j , σ y m = ( − i) m − 1 ( Q 2 m − 2 j =1 w j ) w 2 m , we rewrite (5 2,53) in terms of Ma jorana fermions H = − i n − 1 X m =1 ( J x m w 2 m w 2 m +1 − J y m w 2 m − 1 w 2 m +2 ) − i n X m =1 h m w 2 m − 1 w 2 m (54) L 1 = 1 2 q Γ L 1 ( w 1 − i w 2 ) L 3 = − ( − i) n 2 q Γ R 1 ( w 2 n − 1 − i w 2 n ) W L 2 = 1 2 q Γ L 2 ( w 1 + i w 2 ) L 4 = − ( − i) n 2 q Γ R 2 ( w 2 n − 1 + i w 2 n ) W (55) where W = w 1 w 2 · · · w 2 n is a Casimir op erator whic h comm utes with a ll the elemen ts o f the Clifford algebra generated b y w j and squares to unit y W 2 = 1 . Therefore, W do es not aff e ct the action of ba th op erators (15) where L µ en ter quadratically , so w e find ˆ L 1 + ˆ L 2 = − Γ L + ( ˆ c † 1 ˆ c 1 + ˆ c † 2 ˆ c 2 ) − 2iΓ L − ˆ c † 1 ˆ c † 2 ˆ L 3 + ˆ L 4 = − Γ R + ( ˆ c † 2 n − 1 ˆ c 2 n − 1 + ˆ c † 2 n ˆ c 2 n ) − 2iΓ R − ˆ c † 2 n − 1 ˆ c † 2 n (56) A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 12 leading to the bath shape matrix (50) iden tical to the single fermion case (48). Aga in, carefully following the steps of section 2, w e deriv e the Liouvillean in the form (26) with 4 n × 4 n shap e matrix, whic h w e write in a blo ck tridiag o nal form in terms o f 4 × 4 matrices as A =         B L − h 1 R R 1 0 · · · 0 − R T 1 − h 2 R R 2 . . . 0 0 − R T 2 − h 3 R . . . . . . . . . . . . R n − 1 0 0 · · · − R T n − 1 B R − h n R         (57) and A 0 = Γ L + + Γ R + , where B L := B Γ L + , Γ L − , B R := B Γ R + , Γ R − (in terms of (50)), with Γ L , R ± := Γ L , R 2 ± Γ L , R 1 , and R m :=     0 0 J y m 0 0 0 0 J y m − J x m 0 0 0 0 − J x m 0 0     (58) W e a re not able to p erform a complete diagonalization of the an tisymmetric ma t r ix A (57 ) of t he general XY mo del analytically . F o r example, ev en in the spatially homogeneous case J x , y m ≡ J x , y , h m ≡ h it is not p ossible to pro ceed lik e in the classical harmonic oscillator c hain where the correspo nding matrix is a sum of a T o eplitz and a b ordered mat r ix [30 ]. Namely , in our case A is a sum of a blo ck T o eplitz and blo c k b or der e d matrix and its explicit exact diagonalizatio n remains an op en problem. How ev er, w e stress t ha t ev en relying on n umerical diagonalization of A yielding a set of rapidities β j and pro p erly normalized eigen v ector matrix V , r epresen ts a dramatic progress with resp ect to previously existing nume rical metho ds whic h needed diagonalization of matrices whic h were exponen tially large in n . W e shall later deriv e some exact theoretical and analytical res ults, explaining res ults of exact n umerical computations, in the sp ecial case of a homo gene ous transv erse Ising c hain (subsection 4.1), and the case of a disor der e d XY chain (subsection 4.2) for whic h we shall relate NMM to the problem of Anderson lo calization in o ne dimension, Let us contin ue b y discussing tr a nsp ort observ a bles in the spin chain whose exp ectatio n v alues in NESS are easy to calculate. Note that the bulk Hamiltonian (52) can b e written in terms of the t w o-b o dy e ner gy density op erator H m = − i J x m w 2 m w 2 m +1 + i J y m w 2 m − 1 w 2 m +2 − i h m 2 w 2 m − 1 w 2 m − i h m +1 2 w 2 m +1 w 2 m +2 (59) as H = P m H m . O ne can deriv e the lo cal ener gy curr ent Q m = i[ H m , H m +1 ] fr om the c ontinuity e quation (d / d t ) h H m i = t r H m d ρ/ d t = h i[ H, H m ] i = −h Q m i + h Q m − 1 i (60) where Q m := i[ H m , H m +1 ] Q m = 2i(2 J y m J x m +1 w 2 m − 1 w 2 m +3 + 2 J x m J y m +1 w 2 m w 2 m +4 − A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 13 − J y m h m +1 w 2 m − 1 w 2 m +1 − J x m h m +1 w 2 m w 2 m +2 − − h m +1 J x m +1 w 2 m +1 w 2 m +3 − h m +1 J y m +1 w 2 m +2 w 2 m +4 ) ( 61) The v alidit y of the ab ov e contin uity equation (6 0) dep ends on t w o assumptions only: (i) All Lindblad op erators L µ c om m ute with the densit y H m in the bulk , 2 ≤ m ≤ n − 2 (second equalit y sign), and (ii) [ H m , H m ′ ] = 0 if | m − m ′ | ≥ 2 (third equalit y sign). Using eq. ( 47) o f theorem 3 we can now compute NESS exp ectation v a lues of energy densit y H m and energy curren t Q m , and also of somewhat simpler spin density σ z m = − i w 2 m − 1 w 2 m (62) and spin curr ent S m = σ x m σ y m +1 − σ y m σ x m +1 = − i w 2 m w 2 m +2 − i w 2 m − 1 w 2 m +1 (63) whic h are all quadratic in w j . Note, ho w ev er, that the spin densit y satisfies con tinuit y equation (d / d t ) h σ z m i = −h S m i + h S m − 1 i only in the isotro pic case, when J x m ≡ J y m . 4.1. Homo gene ous tr ansverse I sing ch a in Here w e limit our discuss ion to the spatially homogeneous case J x , y n ≡ J x , y , h n ≡ h . W e shall sho w that in this case the eigen v alue problem A v = β v (64) for (57) can b e most easily and elegan tly treated if f o rm ulated in terms of an abstract inelastic scattering pro blem in one dimension, with asymptotic solutio ns giv en in terms of free normal mo des for the infinite translationally in v arian t chain v = ( . . . , uξ m − 1 , uξ m , uξ m +1 , . . . ) T , where ξ is a complex quasi–momentum (Blo c h) parameter and u is a 4- dimensional amplitude v ector satisfying the condition ( − R T 1 ξ − 1 − h R + R 1 ξ − β I 4 ) u = 0 (65) and the baths pla ying the role of inelastic (absorbing) scatterers at the edges of a finite lattice. The ‘elastic’ (Hamiltonian) v ersion of this tric k has b een used to compute temp oral correlation functions in kic k ed Ising ch ain [31 ]. The singularit y condition for the free mo de equation (65) results, for a general homogeneous X Y model, in eigh t master b ands - different v alues of momenta ξ fo r eac h v alue of the sp ectral parameter (ra pidit y) β . In order to simplify the discussion - whic h will still g et rather in v olv ed - we shall in the fo llowing restrict ourselv es to the tra nsv erse Ising mo del J x = J, J y = 0. In this case w e find just t w o master bands with simple disp ersion relations ξ ± ( β ) = h 2 + J 2 + β 2 ± ω ( β ) 2 hJ ω ( β ) := p ( h 2 + J 2 + β 2 ) 2 − (2 hJ ) 2 (66) but eac h ba nd is doubly-degenerate, since the corresp onding amplitude problem (65) has t w o linearly indep enden t solutions u ± 1 ( β ) =     − h ( h 2 − J 2 + β 2 ± ω ) 0 β ( h 2 + J 2 + β 2 ± ω ) 0     u ± 2 ( β ) =     0 − h ( h 2 − J 2 + β 2 ± ω ) 0 β ( h 2 + J 2 + β 2 ) ± ω )     (67) A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 14 Naiv ely sp eaking, ξ − represen t s left moving and ξ + righ t mo ving free mo des, eac h ha ving t w o p ossible p olarizations. Note that ξ − ξ + = 1. F or a general complex β w e shall c ho o se the branc h of square ro o t ω ( β ) (66) for which | ξ − | ≤ 1. Let us no w write the scattering problem on the le ft bath in terms of an ansatz v =        u L ψ − 1 u − 1 + ψ − 2 u − 2 + ψ + 1 u + 1 + ψ + 2 u + 2 ( ψ − 1 u − 1 + ψ − 2 u − 2 ) ξ − + ( ψ + 1 u + 1 + ψ + 2 u + 2 ) ξ + ( ψ − 1 u − 1 + ψ − 2 u − 2 ) ξ 2 − + ( ψ + 1 u + 1 + ψ + 2 u + 2 ) ξ 2 + . . .        (68) where u L represen t s a 4-dimensional v ector of left-most eigen v ector comp onen ts, ψ − 1 , 2 are the amplitudes of (kno wn) inciden t free mo des, and ψ + 1 , 2 are the amplitudes of the scattered, o utgoing f ree mo des. Plugging t he scattering a nsatz to the eigenproblem (64), the first t w o ro ws of A ( 5 7) result in 6 linearly indep enden t equations for 6 unk nowns ψ + 1 , 2 , u L . Eliminating four v ariables u L w e finally arriv e to the non-unitary 2 × 2 S-matrix  ψ + 1 ψ + 2  = S L  ψ − 1 ψ − 2  (69) with S L 11 = τ − 1 β 2 ( − (Γ L + ) 4 + 4(Γ L + ) 2 ( β 2 − 3 h 2 ) − 16 h ( hJ 2 + iΓ L − ω ) ) S L 12 = τ − 1 β ((Γ L + ) 3 + 8iΓ L − hβ + 4Γ L + ( h 2 − β 2 ))(2i ω ) S L 21 = τ − 1 β ((Γ L + ) 3 − 8iΓ L − hβ + 4Γ L + ( h 2 − β 2 ))( − 2i ω ) S L 22 = τ − 1 β 2 ( − (Γ L + ) 4 + 4(Γ L + ) 2 ( β 2 − 3 h 2 ) − 16 h ( hJ 2 − iΓ L − ω ) ) (70) τ := (Γ L + ) 4 β 2 + 8 β 2 ( h 4 + ( J 2 + β 2 )( J 2 + β 2 − ω ) + h 2 (2 β 2 − ω )) − 2(Γ L + ) 2 ( h 4 + J 4 + 3 β 4 + J 2 (2 β 2 − ω ) − β 2 ω + h 2 ( ω − 2 J 2 − 4 β 2 )) Similarly , one can solv e the scattering problem from t he right bath w ith the scattering ansatz v =        . . . ( ψ + 1 u + 1 + ψ + 2 u + 2 ) ξ − 2 + + ( ψ − 1 u − 1 + ψ − 2 u − 2 ) ξ − 2 − ( ψ + 1 u + 1 + ψ + 2 u + 2 ) ξ − 1 + + ( ψ − 1 u − 1 + ψ − 2 u − 2 ) ξ − 1 − ψ + 1 u + 1 + ψ + 2 u − 2 + ψ − 1 u − 1 + ψ − 2 u − 2 u R        (71) defining the righ t S-matrix  ψ − 1 ψ − 2  = S R  ψ + 1 ψ + 2  (72) Note that since the t w o directions of fr ee mo des (67) do not hav e left-r igh t symmetry a n explicit expression for S R is considerably m ore complicated than (7 0 ) and shall not b e written out here . W e shall now sho w that there exis t t w o qualitativ ely differen t types of NMM - complex rapidities β solving (6 4) for sufficien tly lar ge n . First, we shall discuss the so called evanesc e nt normal master mo des . These a re c haracterized with amplitudes (68) which decay exp onen tially with the distance from A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 15 0.1 0.2 0.3 0.4 Re Β - 2 - 1 0 1 2 Im Β 0.1 0.2 0.3 0.4 Re Β - 2 - 1 0 1 2 Im Β 0.1 0.2 0.3 0.4 Re Β - 2 - 1 0 1 2 Im Β Figure 1. Rapidities β j (black do ts) for a transverse Ising chain with J = 1 . 5 , h = 1 and bath couplings Γ L 1 = 1 , Γ L 2 = 0 . 6, Γ R 1 = 1 , Γ R 2 = 0 . 3, for three different sizes n = 6 (upper ), n = 3 0 (middle), and n = 150 (low er panel). Big blue/ r ed dots indicate po sitions of ev anescent ra pidities (so lutions of eq.(74)) for the left/right bath. – sa y t he left – bath, so t he other – the right b oundary condition b ecomes ph ysically irrelev an t in the limit n → ∞ . Suc h solutions ψ + 1 , 2 = 0 of eq. (69) exist exactly when the determinant of S-matrix v anishes det S L = 0. Using (70) t he determinan t can b e written as det S L = ( β /τ ) 2 p L ( β ) where ¶ p L ( β ) = (Γ L + ) 8 β 2 − 4(Γ L + ) 6 (( h 2 − J 2 ) 2 + (2 J 2 − 4 h 2 ) β 2 + 3 β 4 ) − 16(Γ L + ) 4 (2 h 2 ( h 2 − J 2 ) 2 − (7 h 4 − 6 h 2 J 2 + 2 J 4 ) β 2 + 4( h 2 − J 2 ) β 4 − 3 β 6 ) − 64(Γ L + ) 2 ( h 4 ( h 2 − J 2 ) 2 − 2 h 2 J 4 β 2 − (2 h 4 + 4 h 2 J 2 − J 4 ) β 4 + 2 J 2 β 6 + β 8 ) + 256 h 4 J 4 β 2 (73) ¶ T r ivial zero β = 0 of course doe s no t represent a physical solution since then the whole S-matrix (70) v anishes . A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 16 Th us, for sufficien tly larg e spin c hains we find at most four NMM whose rapidities are giv en as the r o ots of 4-th order p olynomial in β 2 p L ( β ev an ) = 0 (74) that are not sim ultaneously zero es of τ ( β ). Clearly , suc h NMM asy mptotically do not dep end on the c hain size n . In addition, w e find ev anescen t NMM corresp onding to the right bath simply by replacing Γ L + b y Γ R + in (74,73). In fig. 1 w e compare ev anescen t rapidities computed from eq. (74) t o n umerically calculated sp ectrum o f A , at sev eral different sizes n , for a ty pical case o f transv erse Ising c hain, J = 1 . 5 , h = 1 . 0, strongly coupled to tw o baths at considerably differen t temp era t ures, Γ L 1 = 1 , Γ L 2 = 0 . 6, Γ R 1 = 1 , Γ R 2 = 0 . 3 Note that the same pa rameter v alues will b e used for numeric al demonstrations throughout this subsection. Second, w e shall discuss the other extrem e of soft normal master m o des w ith rapidities that are closest to the imaginar y axis, and th us determining the sp ectral gap of the Liouvillean and relaxation time to NESS. Comp osing the scattering from the t w o baths with the free propagation along the chain (back and forth) w e arriv e at the general secular equation for the eigen v alue problem ( 64) in terms of a 2 × 2 determinan t det( ξ 2( n − 3) + S R S L − I 2 ) = 0 (75) In the absence of the baths, Γ L , R ± = 0, t he solutions of t he ab o v e problem exist only fo r real quasi-momenta, namely ξ ± = exp( ± i ϑ ) , ϑ ∈ R . F or suc h extend e d master mo des the lo cal coupling t o the baths can b e considered as a small p erturbation, th us only sligh tly p erturbing the Blo ch-lik e bands β ( e i ϑ ) = ± i ε ( ϑ ) with ‘energy’ ε ( ϑ ) = p h 2 + J 2 − 2 | hJ | cos ϑ (76) The softest NMM, namely the one for which the coupling to the baths is exp ected to b e the weak est, should hav e nearly no des at the ends of the c hain, i.e. ϑ ≈ π /n , or ϑ ≈ π + π /n , and should th us lie near the ba nd edges ± i | h | ± i | J | (see fig. 1 ). In the follo wing w e shall fo cus our calculation on the band edge β ∗ = i( | h | + | J | ) whic h, as can b e chec ke d ap o steriori b y a straig htforw ard but tedious calculation, alwa ys giv es smaller real part of the rapidity than the lo w er edge i( | h | − | J | ), and hence really determines the gap of the Liouvillean. So w e write β = i( | h | + | J | ) + z (77) where z ∈ C is a small parameter, and expand the S-matrices a round the band edge S L , R = − I 2 + 4 g η L , R Z L , R √ − i z + O ( | z | ) (78) where g := q | hJ | 2( | h | + | J | ) , η L , R := ( Γ L , R + ) 4 + 4( Γ L , R + ) 2 (4 h 2 + 2 | hJ | + J 2 ) + 16 h 2 J 2 and Z L 11 = 4 | h | (Γ L + ) 2 + 16 | h | ( | h | + | J | )( | J | − iΓ L − ) Z L 12 = − 2(Γ L + ) 3 − 16Γ L − | h | ( | h | + | J | ) − 8(2 h 2 + 2 | hJ | + J 2 ) Z L 21 = + 2 (Γ L + ) 3 − 16Γ L − | h | ( | h | + | J | ) + 8(2 h 2 + 2 | hJ | + J 2 ) Z L 22 = 4 | h | (Γ L + ) 2 + 16 | h | ( | h | + | J | )( | J | + iΓ L − ) (79) A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 17 and Z R 11 = ( Γ R + ) 4 (2 | h | + | J | ) + 4(Γ R + ) 2 (8 | h | 3 + 9 h 2 | J | + 4 | h | J 2 + | J | 3 ) + 16 h 2 | J | ( | J | (3 | h | + 2 | J | ) − iΓ R − ( | h | + | J | )) Z R 12 = − 2(Γ R + ) 3 − 16Γ R − | h | ( | h | + | J | ) − 8(2 h 2 + 2 | hJ | + J 2 ) Z R 21 = + 2 (Γ R + ) 3 − 16Γ R − | h | ( | h | + | J | ) + 8(2 h 2 + 2 | hJ | + J 2 ) Z R 22 = ( Γ R + ) 4 (2 | h | + | J | ) + 4(Γ R + ) 2 (8 | h | 3 + 9 h 2 | J | + 4 | h | J 2 + | J | 3 ) + 16 h 2 | J | ( | J | (3 | h | + 2 | J | ) + iΓ R − ( | h | + | J | )) (80) Next we expand ξ + (66) in z , yielding ξ + = − 1 − g − 1 √ − i z + O ( | z | ) (81) and so the fr e e pr op agator in (7 5) can b e written as ξ 2( n − 3) + = exp(2 ng − 1 √ − i z ) + O ( | z | ) . (82) In eqs. (78,81,82) the branc h cut along the negativ e real axis has b een chos en for √ − i z . Since the product of S-matrices in (75) is near iden tit y , the free propagator should b e near one a s w ell, hence 2 ng − 1 √ − i z should b e near 2 π i. Let us define z 0 b y setting 2 ng − 1 √ − i z 0 = 2 π i, so z 0 = − i π 2 g 2 n − 2 (83) and write z = z 0 (1 + y ) where | y | ≪ 1 is another small complex parameter. Ho w ev er, since z 0 is purely imaginary , w e need to compute a small but non-v a nishing y whic h will, in the leading order in n , solv e (75) since the real part of the soft mo de’s rapidity is determined as Re β = R e z 0 y = π 2 g 2 n − 2 Im y (84) No w, writing √ − i z = √ − i z 0 √ 1 + y = i π g n − 1 (1 + y / 2 − y 2 / 8) + O ( y 3 ) in (78,82), plugging all that to eq. ( 75) and computing to or der O ( n − 2 ), noting tha t O ( | z | ) = O ( n − 2 ), w e arr iv e to a simple quadratic equation for y , whose solution, plugged to (84), gives the final result, namely the sectral gap of Liouvillean ∆ = 2 Re β ∆ = (2 π hJ ) 2 ( | h | + | J | ) 2 ∆ 1 ∆ 2 n − 3 + O ( n − 4 ) (85) ∆ 1 := 64(Γ L + + Γ R + ) h 2 J 2 (2 h 2 + 2 | hJ | + J 2 ) + 16((Γ L + ) 3 + (Γ R + ) 3 ) h 2 J 2 + 16Γ L + Γ R + (Γ L + + Γ R + )(2 h 2 + 2 | hJ | + J 2 )(4 h 2 + 2 | hJ | + J 2 ) + 4Γ L + Γ R + ((Γ L + ) 3 + (Γ L + ) 2 Γ R + + Γ L + (Γ R + ) 2 + (Γ R + ) 3 )(2 h 2 + 2 | hJ | + J 2 ) + (Γ L + Γ R + ) 3 (Γ L + + Γ R + ) ∆ 2 := ((Γ L + ) 4 + 4( Γ L + ) 2 (4 h 2 + 2 | hJ | + J 2 ) + 16 h 2 J 2 ) × ((Γ R + ) 4 + 4( Γ R + ) 2 (4 h 2 + 2 | hJ | + J 2 ) + 16 h 2 J 2 ) In fig . (2) w e compare this analytical result to exact nume rical calculations o f the eigen v alue of A with minimal real part, confirming b oth, its precise num erical v alue and that the relative scaling of the next order correction is indeed O ( n − 1 ). A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 18 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 - 2 - 1 0 1 Log n Log ¡ n 3 D - C ¥ 0 50 100 150 200 250 300 350 0 5 10 15 n n 3 D Figure 2. Sp ectral ga p ∆ times a third p ow er of the c hain length n for a transverse Ising chain with J = 1 . 5 , h = 1 and bath co uplings Γ L 1 = 1 , Γ L 2 = 0 . 6, Γ R 1 = 1 , Γ R 2 = 0 . 3. Thin horizontal line indica tes the theoretical a symptotic v alue (85). In the inset we show devia tion from asymptotic constant v alue of ∆ n 3 in log-lo g scale and demonstrate that it decays a s ∝ n − 1 (thin line). - 5 - 4 - 3 - 2 - 1 0 - 15 - 10 - 5 0 5 10 15 Re Λ Im Λ Figure 3 . C o mplete sp ectrum of 2 12 complex eigenv alues of Liouville an for a transverse Ising chain with n = 6 spins and J = 1 . 5 , h = 1 and bath couplings Γ L 1 = 1 , Γ L 2 = 0 . 6, Γ R 1 = 1 , Γ R 2 = 0 . 3 (the case of the upp er pane l of fig. 1). Note that, intere stingly , b oth main analytical results of this subsection, namely ev anescen t a nd soft mo de rapidities do not dep end o n Γ L , R − . Ph ysically speaking, they only dep end on the effectiv e strengths of the bath couplings a nd not on t he temp eratures. W e end t his subsection b y presen ting some further n umerical results on heat transp ort in the op en transv erse Ising chain in the Lindblad form. In fig. 3 w e demonstrate expression (42) fo r constructing the full sp ectrum of the Liouvillean in A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 19 0 10 20 30 40 50 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 n < S > , < Q > Figure 4. Energy curr ent (upp er/blue p oints), and average s pin current (lo wer/red po int s), versus chain leng th n for a transverse Is ing chain with J = 1 . 5 , h = 1 a nd bath couplings Γ L 1 = 1 , Γ L 2 = 0 . 6, Γ R 1 = 1 , Γ R 2 = 0 . 3. 0 20 40 60 80 - 0.5 - 0.4 - 0.3 - 0.2 - 0.1 0.0 m , m < > Σ m z 0 2 4 6 8 10 - 25 - 20 - 15 - 10 - 5 Log ∆ H m 70 72 74 76 78 80 - 25 - 20 - 15 - 10 - 5 Log ∆ Σ m z Figure 5. Ener gy density profile (lower, blue p oints), and spin density profile (upper , red p oints), for a transverse Ising chain of n = 80 spins with J = 1 . 5 , h = 1 and bath co upling s Γ L 1 = 1 , Γ L 2 = 0 . 6, Γ R 1 = 1 , Γ R 2 = 0 . 3. The insets display logarithm of the differ e nce to the bulk v alues δ H m := |h H m i − H bulk | (blue points), δ σ z m := |h σ z m i − σ z bulk | (red p oints) in co mparison with ± (4 log ξ − ) m + const with quasi-momentum ξ − = 0 . 58 4692 corre s po nding (66) to the leading ev anes cent r a pidit y β ev an = 0 . 4387 39 (full lines). terms of a set of rapidities, for a short c hain. In fig. 4 w e demonstrate eq. (47) of Theorem 3 b y computing the energy curren t Q m (61), a nd the a v erage spin current S = 1 n − 1 P n − 1 m =1 S m (63) in NESS of a t ypical transv erse Ising c hain. Numerical results giv e a clear indication of b al listic tr ansp ort h Q i = O ( n 0 ) , h S i = O ( n 0 ), ho we ve r its rigorous pro of and analytical calculation of the curren ts w ould re quire full control ov er the complete set of NMM whic h is at presen t not av ailable. In fig. 5 w e plot the energy A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 20 10 20 30 40 50 - 14 - 12 - 10 - 8 - 6 - 4 - 2 n Log Figure 6. Average Lio uvillean sp ectra l g a p h ∆ i versus the chain le ng th n for disordered XY mo dels: (i) J x m = 0 . 5 , J y m = 0 , h m ∈ [1 , 2] (transverse Ising with field disorder, blue p oints), (ii) J x m ∈ [0 . 5 , 2] , J y m = 0 , h m = 1 (transverse Is ing with in teraction disorder, red points), (iii) J x m ∈ [0 . 5 , 1] , J y m ∈ [0 . 5 , 1] , h m = 1 (XY with interaction disorder, golden points), all for ba th couplings Γ L 1 = 1 , Γ L 2 = 0 . 6, Γ R 1 = 1 , Γ R 2 = 0 . 3. F ull lines indicate exp onential fits to right halves of data . Averaging is p er formed ov er 20 00 dis o rder realizations. densit y (59) and spin densit y (62) pro files in NESS. Aga in, w e note flat profiles in the bulk of t he c hain, m, n − m ≫ 1, with exp onential falloff due to adjustmen t to the non-equilibrium bath v alues. Since the densities can b e written, by means of (47), as 4 − point functions in NMM comp onen ts, t he leading falloff exp o nen ts of the profile |h H m i − H bulk | ∼ | ξ − | 4 m is giv en b y the quasi-momen tum ξ − (66) corresp onding to the maximal ev anescen t rapidity β ev an (74). 4.2. Disor d e r e d XY chain In t his subsection we treat the opp osite extreme, a disordered XY chain (52) where three sets of ph ysical parameters are chose n as r andom unc orr elate d v ariables from unifo rm distributions on the interv a ls, J x m ∈ [ J x min , J x max ], J y m ∈ [ J y min , J y max ], h m ∈ [ h min , h max ]. Clearly , the eigen v alue pr o blem (64) for the matrix (57) then b ecomes equiv alen t to the Anderson t ig h t-binding problem in one dimension for a quantum part icle with a 4 − lev el in ternal degree of f reedom. W e do not pursue an y theoretical ana lysis of this problem here, but merely state that nume rical inv estigations indicate existence of exp o nen tial lo calization of al l eigen v ectors (or no rmal master mo des) for disorder of any strength in an y one o f system’s para meters. With the picture of lo calization of NMM in mind, the effect of the couplings to t he heat bat hs at the c hain’s ends on quantum transp o rt can b e predicted b y theoretical argumen ts (see [3 2] for a review of related studies): The sp ectral gap of t he Liouvillean should b e exp onen tially small ∆ ∼ exp ( − n/ℓ ) where ℓ is the lo calization length o f NMM A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 21 10 20 30 40 50 - 6 - 5 - 4 - 3 - 2 - 1 n Log < Q > Figure 7. The scaling of the energy current h Q m i with chain length n of the diso r dered XY model in the same regimes/ parameter s /plot styles as in fig. 6. 0.0 0.2 0.4 0.6 0.8 1.0 - 0.50 - 0.45 - 0.40 - 0.35 - 0.30 m - 1 n - 2 m Figure 8. Scaled ener gy density profile of interaction disorder ed XY c hains (case (iii) of fig. 6) for three chain sizes: n = 2 0 (blue p oints), n = 40 (r ed p oints), n = 60 (golden points). Av era ges over 5 0 000 disorder realizations ha ve b een p erformed. whic h is exp ected to b e prop ortional t o the square of in ve rse disorder strength. This is demonstrated in fig.6. If all NMM are expo nentially lo calized, the curren ts should decrease with the chain size n faster than an y p o we r, p erhaps exponentially , and the system should b eha ve as a n ideal insulator (at all temp eratur es). This is demonstrated b y straigh tforward nume rical calcu latio ns of the heat curren t (61) in fig. 7. In the final figure 8 w e plot the energy densit y profile h H m i ( 5 9) in a t ypical case of disordered XY c hain, v ersus a scaled spatial co o rdinate ( m − 1 ) / ( n − 1) ∈ [0 , 1], for sev eral differen t c hain lengths n , and demonstrate sharping up of energy densit y profiles with increasing n , whic h is again indicating insulating behaviour. A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 22 5. Discussion and conclusions The main result of the pap er is a general metho d of explicit solution of master equations describing dynamics of op en quan tum system, under the condition that the system’s Hamiltonian is quadr atic and all Lindblad op erators are line ar in canonical fermionic op erators (whic h can either represen t real ph ysical fermions o r an y abstract 2-lev el quan tum systems (qubits) thru the Jordan-Wigner tra nsformation). Using a nov el concept of F o ck space o f ph ysical op erators (or densit y op erators o f ph ysical states), and the adjo in t structure of canonical creation and annihilatio n maps o v er this space, the problem can be treated in terms o f a non-Hamilto nian generalization o f the metho d of Lieb, Sc h ultz and Mattis [18] lift ed to an op erator space. W e ha v e explicitly constructed a non-canonical analog of Bog oliub ov transforma t io n of the quantum Liouville map to normal master mo des. Related ideas in the Ha milto nian con text ha v e b een used by the author [31, 33, 3 4] in order to approach t he problem of real time dynamics and ergo dic prop erties of isolate d in teracting ma ny-bo dy quantum syste ms. As an illustratio n of the metho d w e hav e solv ed far fro m equilibrium quan tum heat and spin transp ort problem in Heise nberg XY spin 1/2 c hains whic h are coupled to canonical heat baths only at the t w o ends. Irresp ectiv ely o f the strength of the coupling to the baths and their temp eratures, w e hav e shown a ballistic transp ort in the spatially homogeneous (non- disordered) case, and an ideally insulating b eha viour in the disordere d case asso ciated to lo calization of normal master mo des of the quan tum Liouville op erator. In this contex t the metho d can b e considered as a simple alternativ e to the solution of quan tum L a ngevin equations [24]. Ho w ev er, the metho d should easily b e applicable to v ariety of other ph ysical situations, for example if all fermions are coupled to the baths one c ould mak e a solv able model of gen uine q uantum diffusion, a man y-b o dy generalization of the tight- binding mo del [3 5 ]. W e a lso expect the metho d to b e applicable to the Redfield t yp e of master equations (see e.g. [35]) - which do not conserv e p ositivit y for a short initia l (slippage) time in terv al - prov ided only the system pa r t o f t he Hamiltonian is q uadr atic and system-bath couplings are line ar in fermionic v a riables. F urthermore, extension of the metho d to o p en many-b oson systems should b e straigh tforward, simply b y replacing an ticomm utato r s b y comm utato rs throughout the expo sition of section 2. T reating densit y op erators of NESS a s elemen ts o f a Hilb ert space of op erators one ma y also extend the concept of e n tanglement entr opy , with r esp ect t o a bipa r t ition of a system of many fermions [36], to NESS whic h can in our appro ac h b e view ed as a kind of ground state of the Liouvillean. Saturation of suc h op er ator sp ac e entang lement entr opy [34] (whic h is suggested by n umerical exp erimen ts [3 7]) indicates efficient simulability of NESS b y elab orate n umerical metho ds such as density matrix r enormalization gr oup [38], p erhaps ev en for more general, non-solv able quan tum systems. As last w e mention a more am bitious extension of the presen t w ork: Namely w e prop ose to explore a question, whe ther more inv olved a lgebraic metho ds of solution of in teracting man y-b o dy quan tum systems, lik e e.g. Bethe Ansatz or quantum in ve rse A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 23 scattering [13], could b e generalized to op en quantum systems, e.g. b y means of the prop osed concept of F o c k space of op erators. Could one discuss completely in tegrable op en quan tum systems whic h go b eyond quadratic Liouvilleans? Ac knowle dgemen ts I gratefully ac kno wledge stim ulating discussions with Pierre Gaspard, Keiji Saito and W alter Strunz, thank Carlos Mejia-Monasterio and Thomas H. Seligman for reading the man uscript and many use ful commen ts, and Iztok Pi ˇ zorn and Mark o ˇ Znidari ˇ c for collab oration on related pro jects. The w ork has b een suppo rted by the grants P1- 0044 and J1- 7347 of Slov enian researc h ag ency (ARRS). Explicit analytical calculations rep orted in subsection (4.1) were assisted b y Mathematic a softw a re pac k age. References [1] W. H. Zurek, Rev. Mod. P hys. 75 , 715 (2 003). [2] E. Jo os, H. D. Zeh, C. Kiefer, D. Giulini, J. Kupsch and I.- O. Stamatescu, De c oher enc e and the App e ar anc e of a Classic al World in Qu ant um The ory , (Springer, 2003). [3] J. v o n Neumann, Mathematic al F oundations of Qu ant um Me chanics , T rans. Ro b er t T. Geyer. (Princeton Universit y Pr ess, P rinceton 1955). [4] M. Schlosshauer, Rev. Mo d. Phys. 76 , 1267 (2004). [5] M. A. Nielsen and I. L. Ch uang, Quant u m Computation and Quantum Information (Camb r idg e Univ er sity P ress, Ca mbridge 200 0). [6] G. Benenti, G. Casati and G. Strini, Principles of Quantum Computation and Information. V olume I: Basic Conc epts (W orld Scientific, Singa p o re 200 4); V olume II: Basic T o ols and Sp e cial T opics (W or ld Scientific, Singap ore 2007 ). [7] H. Araki and E. Barouch, J. Sta t. Phys. 31 , 327 (1 983); H. Araki, Publ. RIMS Kyoto Univ. 20 , 277 (198 4 ). [8] D. Ruelle, J. Stat. Phys. 98 , 57 (2000). [9] V. Ja k ˇ siˇ c and C.-A. P illet, J. Stat. P hys. 108 , 787 (2002 ); Commun. Ma th. Phys. 22 6 , 131 (20 02); W. As chbac her, V. Jakˇ siˇ c, Y. Pautrat and C.-A. Pillet, Inr o duction t o non-e quilibrium quantum statistic al me chanics , in O p en Quantu m Systems III. R e c ent D evelopments Lecture Notes in Mathematics, 1882 (2006), 1 -66. [10] P . Gaspar d, Pro g. Theor . P hys. Suppl. 165 , 33 (20 0 6); Ph ys ic a A 369 , 201 (2006). [11] M. H. Lee, Acta Physica Polonica B 38 , 1837 (2007); Phys. Rev. Lett. 87 , 250 601 (2001); Phys. Rev. L ett. 49 , 1072 (19 82). [12] T. Pro sen, J. Phys. A: Math. Theor. 40 , 7881 (2007). [13] V. E. K orepin, N. M. Bo goliub ov, and A. G. Izergin, Quant um Inverse Sc attering and Corr elation functions (Cam bridge Universit y Press , Cambridge 199 7). [14] L. F adee v, P . V an Mo erb eke and F. Lambert (Eds.), Biline ar Int e gr able Syst ems: fr om Classic al to Quantum, Continuous t o Discr ete , (NA TO AR W Pro ceedings), Springer Series: NA TO Science Series II: Ma thematics, Ph ysics and Chemistry , V ol 20 1 (20 0 6). [15] H.-P . Breuer and F. Petruccione, The The ory of Op en Quantum Systems (Oxford University Press, Oxford 2002). [16] F. Haa ke, Q u antum Signatur es of Chaos , 2nd edition (Spring e r, 2001). [17] R. Alicki and K. Lendi, Quant um dynamic al semigr oups and applic ations (Springer, 2007 ). [18] E . H. Lieb, T. D. Sc hultz and D. C. Mattis, Ann. Phys. (New Y ork) 16 , 4 07 (1961). [19] X. Zo tos, F. Naef and P . Prelov ˇ sek, Phys. Rev. B 55 , 110 29 (1997 ). A gene r al metho d to solve ma s ter e quations for quadr atic op en F erm i systems 24 [20] K . Sa ito, S. T akesue and S. Miyashita, Phys. Rev . E 61 , 23 97 (2 000); K. Saito, E urophys. Lett. 61 , 34 (2 003). [21] M. Michel, M. Har tma nn, J. Gemmer and G. Ma hler, Eur. Phys. J. B. 34 , 3 25 (2 003); M. Mic hel, G. Mahler and J. Gemmer, P h ys . Rev. Lett. 95 , 180 602 (200 5); M. Mic hel, J. Gemmer and G. Ma hler, In t. J. Mo d. Phys. B 20 , 4855 (2 0 06); J. Gemmer, R. Steinigeweg a nd M. Michel, Phys. Rev. B 73 , 1 04302 (2006). [22] C. Mejia-Mona sterio, T. P rosen and G. Casa ti, E urophys. Lett. 72 , 520 (20 05); G. Casati a nd C. Mejia-Monas terio, e- print arXiv: 0710. 3500 v1 [cond- mat.st at-mech] . [23] C. Mejia-Mona sterio a nd H. Wich terich, e-pr int ar Xiv:0 709.14 12 v1 [c ond-m at.sta t-mech] . [24] A. Dhar and D. Roy , J . Stat. Phys. 125 , 805 (2006); see als o e-print arXiv:0 711.4 318 v 1 [cond- mat.s tat-mech] . [25] G. Lindblad, Commun. Math. P hys. 48 , 119 (197 6). [26] M. M. W o lf, J. Eisert, T. S. Cubitt a nd J . I. Cir ac, e-print arXiv: 0711. 3172 v 1 [quant-p h] . [27] P . ˇ Semrl, priva te c ommu nic ation . [28] W e note a simila r ity to the for malism of seco nd quantization with no n-orthog o nal o rbitals int ro duced in: M. Moshinsky and T. H. Se lig man, Ann. P h ys . (New Y ork) 66 , 311 (1971). [29] H. Wic hteric h, M. J. Her ich, H. P . B reuer, J. Gemmer and M. Michel, Phys. Rev. E, 76 03 1115 (2007). [30] Z . Rieder, J. L. Leb owitz and E. Lieb, J. Math. Phys. 8 , 10 73 (19 67). [31] T. Pro sen, Pr og. Theo r. Ph ys . Suppl. 139 , 191 (2000). [32] S. Lepr i, R. Livi and A. Politi, Phys. Rep. 377 , 1 (2003). [33] T. Pro sen, Phys. Rev . E 60 , 16 58 (1 999). [34] T. Pro sen and I. Piˇ zorn, Ph ys . Rev. A 76 , 0 32316 (2007). [35] M. Esp osito and P . Gaspard, Ph ys . Rev. B 71 , 214302 (2005); J. Stat. Phys. 121 , 4 6 3 (20 05). [36] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev , P h ys . Rev. Lett. 90 , 22790 2 (2 003); J. I. L a torre, E. Rico , a nd G. Vidal, Qua nt . Inf. Comp. 4 , 48 (2004). [37] T. Pro sen and M. ˇ Znidariˇ c, to b e submitted (2008). [38] S. R. White, Phys. Rev. Lett. 69 , 2863 (1992); U. Schollw¨ oc k a nd S. R. White, in G. G. Ba trouni, and D. P oilbla nc (eds.): E ffective mo dels for low-dimensional strongly correla ted systems, p.155, AIP , Melville, New Y ork (2006); G. Vidal, P h ys . Rev. Lett. 91 , 14790 2 (2003 ); ibid. 93 , 04050 2 (2004).