Third quantization: a general method to solve master equations for quadratic open Fermi systems
aa r X i v : . [ qu a n t - ph ] D ec Third quantization: a general method to solvemaster equations for quadratic open Fermi systems
Tomaˇz Prosen
Department of physics, FMF, University of Ljubljana, Jadranska 19, SI-1000Ljubljana, Slovenia
Abstract.
The Lindblad master equation for an arbitrary quadratic system of n fermions is solved explicitly in terms of diagonalization of a 4 n × n matrix, providedthat all Lindblad bath operators are linear in the fermionic variables. The method isapplied to the explicit construction of non-equilibrium steady states and the calculationof asymptotic relaxation rates in the far from equilibrium problem of heat and spintransport in a nearest neighbor Heisenberg XY spin 1 / general method to solve master equations for quadratic open Fermi systems
1. Introduction
Understanding time evolution of an open quantum system of many interacting particlesis of primary importance in fundamental problems of quantum physics, such asdecoherence [1, 2] and closely related quantum measurement problem [3, 4], quantumcomputation [5, 6], or the problem of computation of non-equilibrium steady states (NESS) in quantum statistical mechanics [7, 8, 9]. Even though application of themethods of Hamiltonian dynamical systems and ergodic theory to quantum systemsout of equilibrium gives many promising results [10, 11, 12], the field of open quantumsystems is still lacking non-trivial explicitly solvable models, as compared to studiesof closed (isolated) quantum systems where we know a large body of the so-called completely integrable systems [13, 14]. Examples of explicitly solvable models of masterequations for open quantum systems are limited to quite restricted models of a singleparticle, single spin or harmonic oscillators (see e.g. [15, 16, 17]).In this paper we show that the generator of the master equation of a generalquadratic system of n interacting fermions which are coupled to a general set ofMarkovian baths, specified in terms of Lindblad operators which are linear in thefermionic variables - the so called quantum Liouville super-operator (or Liouvillean) -can be explicitly diagonalized in terms of 2 n normal master modes , i.e. anticommutingsuper-operators which act on the Fock space of density operators. This can beunderstood as a complex (non-canonical) version of the Bogoliubov transformation [18]lifted on the operator space, and has very powerful consequences: (i) The NESS of themaster equation can be understood as the ‘ground state’ normal mode of the Liouvillean,whereas the long time relaxation rate is given by the eigenmode closest to the realaxis. (ii) The covariance matrix of NESS can be computed explicitly in terms of theeigenvectors of 4 n × n antisymmetric complex matrix. It can be used to completelyexpress physical observables in NESS, such as particle/spin densities, currents, etc.We demonstrate the power of this novel method by applying it to the problem ofheat and spin transport far from equilibrium in nearest neighbor Heisenberg XY spin 1 / ballistic transport in the integrable spatially homogeneous case (see e.g. [19, 20, 21, 22, 23, 24, 12] forrelated recent studies of quantum thermal conductivity in one dimension), and predict ideally insulating behaviour (at all temperatures) in a disordered case of spatially randominteractions/field. Apart from obtaining numerical results which go by far beyondwhat was so far accessible by direct numerical solution of the many-particle Lindbladequation, either directly or by means of quantum trajectories [15], we also obtain twonotable analytical results in the spatially homogeneous (non-disordered) case: (i) Wecompute the spectral gap of the Liouvillean i.e. the rate of of relaxation to the NESSand show that it scales with the inverse cube of the chain length. (ii) We construct evanescent normal master modes of the Liouvillean, for long chains, by which we explainquantitatively the exponential falloff of energy density or temperature profiles near thebath contacts. general method to solve master equations for quadratic open Fermi systems
2. General method of solution for the Lindblad equation
The general master equation governing time evolution of the density matrix ρ ( t ) ofan open quantum system, preserving trace and positivity of ρ , can be written in theLindblad form [25, 17] as (we set ~ = 1)d ρ d t = ˆ L ρ := − i[ H, ρ ] + X µ (cid:0) L µ ρL † µ − { L † µ L µ , ρ } (cid:1) (1)where H is a Hermitian operator (Hamiltonian), [ x, y ] := xy − yx , { x, y } := xy + yx ,and L µ are arbitrary operators representing couplings to different baths (at possiblydifferent values of thermodynamic potentials). We are now going to describe a generalmethod of explicit solution of (1) for a quadratic system of n fermions (or spins 1 / linear bath operators H = n X j,k =1 w j H jk w k = w · H w (2) L µ = n X j =1 l µ,j w j = l µ · w (3)where w j , j = 1 , , . . . , n , are abstract Hermitian
Majorana operators satisfying theanti-commutation relations { w j , w k } = 2 δ j,k j, k = 1 , , . . . , n (4)and generate a Clifford algebra. Thus, 2 n × n matrix H can be chosen to beantisymmetric H T = − H . Throughout this paper x = ( x , x , . . . ) T will designate avector (column) of appropriate scalar valued or operator valued symbols x k .Two notable examples, to which our formalism is immediately applicable, are: (i)canonical fermions c m , m = 1 , , . . . , n , w m − = c m + c † m w m = i( c m − c † m ) (5)or (ii) spins 1 / ~σ m , m = 1 , , . . . , n , w m − = σ x m Y m ′ Hamiltonians H ( t ) and Lindblad operators L µ ( t ), generating more general and possiblynon-Markovian open system dynamics, is straightforward. See e.g. [26] for a discussionof Markovianity . We begin by associating a Hilbert space structure x → | x i to a linear 2 n = 4 n dimensional space K of operators, with a canonical basis | P α i with P α ,α ,...,α n := w α w α · · · w α n n α j ∈ { , } (7) orthonormal with respect to an inner product h x | y i = 2 − n tr x † y (8)The form of the canonical basis of the operator space K suggests that it is just a usualFock space with an unusual physical interpretation. Namely we can define the followingset of 2 n adjoint annihilation linear maps ˆ c j over K ˆ c j | P α i = δ α j , | w j P α i (9)and derive the actions of their Hermitian adjoints - the creation linear maps ˆ c † , h P α ′ | ˆ c † j | P α i = h P α | ˆ c j | P α ′ i ∗ = δ α ′ j , h P α | w j P α ′ i ∗ = δ α j , h P α ′ | w j P α i :ˆ c † j | P α i = δ α j , | w j P α i (10)Straightforward inspection then shows that they satisfy the canonical anticommutationrelations { ˆ c j , ˆ c k } = 0 { ˆ c j , ˆ c † k } = δ j,k j, k = 1 , , . . . , n (11)The key is now to realize that the quantum Liouville map ˆ L defined by eqs. (1,2,3) canbe written as a quadratic form in adjoint Fermi maps ˆ c j , ˆ c † j (or for short, a-fermions ). ‡ First, we consider the unitary part of Liouvilleanˆ L ρ := − i[ H, ρ ] (12)Since K is a Lie algebra, one defines the adjoint representation of a Lie derivative foran arbitrary A ∈ K back on K as, ad A | B i := | [ A, B ] i . It is now straightforward to ‡ Throughout this paper Dirac’s bra-ket notation shall be used only for a Hilbert space K of physicaloperators, including density operators, in a sense of GNS construction although here all spaces willbe finite dimensional. Symbols with a hat shall designate linear maps over the operator space K . Forinstance, we note a key distinction between physical fermions c m (5) and a-fermions ˆ c j (9). general method to solve master equations for quadratic open Fermi systems w j w k | P α i = | w j w k P α i − | P α w j w k i == 2( δ α j , δ α k , + δ α j , δ α k , ) | w j w k P α i == 2(ˆ c j ˆ c † k + ˆ c † j ˆ c k ) | P α i = 2(ˆ c † j ˆ c k − ˆ c † k ˆ c j ) | P α i (13)Extending this relation by linearity to an arbitrary element of K , it follows that for anarbitrary quadratic Hamiltonian (2) its Lie derivative has a very similar quadratic formin a-Fermi mapsˆ L = − i ad H = − n X j,k =1 ˆ c † j H jk ˆ c k = − 4i ˆ c † · H ˆ c (14)It is worth stressing here that for an arbitrary (complex) matrix H , ˆ L (14) conservesthe total number of a-fermions ˆ N := P j ˆ c † j ˆ c j = ˆ c † · ˆ c , namely [ ˆ L , ˆ N ] = 0.Second, we consider the action of the Lindblad mapsˆ L µ ρ := 2 L µ ρL † µ − { L † µ L µ , ρ } = n X j,k =1 l µ,j l ∗ µ,k ˆ L j,k ρ (15)where we write ˆ L j,k ρ := 2 w j ρw k − w k w j ρ − ρw k w j . Again we proceed by computingthe actions of ˆ L j,k on elements of the canonical basis of operator Fock space K . Inorder to do so, it is crucial to observe that the question whether w j commutes oranticommutes with P α depends on the number of a-fermions | α | := P nk =1 α k in | P α i ,namely P α w j = ( − | α | + α j w j P α , and henceˆ L j,k | P α i = (cid:2) − | α | + α k w j w k − w k w j − ( − α j + α k w k w j (cid:3) | P α i (16)Observing that | w j P α i = (ˆ c † j + ˆ c j ) | P α i (17)( − α j | w j P α i = (ˆ c † j − ˆ c j ) | P α i (18)( − | α | | P α i = exp(i π ˆ N ) | P α i (19)one derives from (16) the general expression for ˆ L j,k ˆ L j,k = (cid:16) ˆ + exp(i π ˆ N ) (cid:17) (cid:16) c † j ˆ c † k − ˆ c † j ˆ c k − ˆ c † k ˆ c j (cid:17) + (cid:16) ˆ − exp(i π ˆ N ) (cid:17) (cid:16) c j ˆ c k − ˆ c j ˆ c † k − ˆ c k ˆ c † j (cid:17) (20)Obviously, the maps ˆ L j,k , and hence also the total Lindblad part of Liouvillean P µ ˆ L µ ,do not conserve the number of a-fermions. But they conserve its parity i.e. the productof any two creation/annihilation a-Fermi maps commutes with the parity operationˆ P = exp(i π ˆ N ), with respect to which the operator space can be decomposed into adirect sum K = K + ⊕ K − , and even/odd operator spaces are orthogonally projected as K ± = (ˆ ± exp(i π ˆ N )) K . Thus the positive parity subspace K + is a linear space spanned general method to solve master equations for quadratic open Fermi systems | P α i with even | α | . All the maps ˆ L j,k now act separately on K ± , ˆ L j,k K ± ⊆ K ± . Forexample, the maps defined on even parity subspace are indeed quadratic in a-fermionsˆ L j,k | K + = 4ˆ c † j ˆ c † k − c † j ˆ c k − c † k ˆ c j (21)In this paper we shall focus on physical observables which are products of an even number of Majorana fermions – operators w j – so we shall in the following discussonly Liouville dynamics on the subspace K + . The extension to the dynamics of odd observables should be straightforward.Putting the results (12,14,15,21) together we arrive at the final compact quadraticrepresentation of the Liouvillean ˆ L + := ˆ L| K + ˆ L + = − c † · (2i H + M + M T ) ˆ c + 2 ˆ c † · ( M − M T ) ˆ c † (22)where M is a complex Hermitian matrix parametrizing the Lindblad operators M jk = X µ l µ,j l ∗ µ,k (23) Next we want to show that the representation (22) allows us to reduce it further bya linear transformation of the set of maps { ˆ c j , ˆ c † j ; j = 1 , , . . . , n } to normal mastermodes (NMM) in terms of which the complete spectrum of the Liouvillean, as well asits eigenvectors, can be explicitly constructed; in particular the zero-mode eigenvectorwhich is just the physically relevant NESS.In fact we proceed in analogy to Lieb et al. [18] and define 4 n adjoint HermitianMajorana maps ˆ a r = ˆ a † r , r = 1 , , . . . , n :ˆ a j − = 1 √ c j + ˆ c † j ) ˆ a j = i √ c j − ˆ c † j ) (24)satisfying the anti-commutation relations { ˆ a r , ˆ a s } = δ r,s (25)in terms of which the Liouvillean (22) can be rewritten asˆ L + = ˆ a · A ˆ a − A ˆ (26)where A is an antisymmetric complex 4 n × n matrix with entries A j − , k − = − H jk − M jk + M kj A j − , k = 2i M kj A j, k − = − M jk A j, k = − H jk + M jk − M kj (27)ˆ is an identity map over K and A is a scalar A = 2 n X j =1 M jj = 2 tr M (28) general method to solve master equations for quadratic open Fermi systems bi-linear Liouvillean (26) cannot be brought to a normal form with alinear canonical transformation since the matrix A – which shall in the following bereferred to as a shape matrix of Liouvillean – is not anti-Hermitian like in Hamiltoniansystems. So we should proceed in more general terms.We first recall few facts about complex antisymmetric matrices of even dimension.If v is a right eigenvector A v = βv with complex eigenvalue β , then v is also a left eigenvector with eigenvalue − β , A T v = − A v = − βv . Hence eigenvalues alwayscome in pairs β, − β . Let as assume that A can be diagonalized § , i.e. there exist4 n linearly independent vectors v r , r = 1 , . . . , n with the corresponding eigenvalues β , − β , β , − β , . . . , β n , − β n , A v j − = β j v j − A v j = − β j v j (29)ordered such that Re β ≥ Re β ≥ . . . ≥ Re β n ≥ . The 2 n complex numbers β j shallbe referred to as rapidities . It is easy to check that we can always choose and normalize v r such that k v r · v s = J rs where J := σ x ⊗ I n = · · · · · · · · · · · · ... ... ... ... . . . (30)Let V be 4 n × n matrix whose r th row is given by v r , V rs := v r,s . Then eqs. (29,30)rewrite as AV T = V T D where D := diag { β , − β , β , − β , . . . , β n , − β n } (31) VV T = J (32)Expressing V T in terms of (32) and plugging the result into eq. (31) we arrive at a veryconvenient canonical form of a generic complex antisymmetric matrix AA = V T ΛV where Λ = DJ = β · · ·− β · · · β · · · − β · · · ... ... ... ... . . . (33)Now we apply decomposition (33) to the Liouvillean (26)ˆ L + = ˆ a · V T ΛV ˆ a − A ˆ = ( V ˆ a ) · Λ ( V ˆ a ) − A ˆ (34)Let us define the NMM maps ˆ b := (ˆ b , ˆ b ′ , ˆ b , ˆ b ′ , . . . , ˆ b n , ˆ b ′ n ) := V ˆ a orˆ b j = v j − · ˆ a ˆ b ′ j = v j · ˆ a (35) § It is not known at present whether explict form (27) guarantees diagonalizability of any such A . Notethat one can construct certain types of complex antisymmetric matrices with degenerate eigenvalueswhich cannot be diagonalized [27]. k For a non-degenerate rapidity spectrum { β j } the proof of this statement is a trivial consequence ofantisymmetry A = − A T , whereas in case of degeneracies it can be shown that one can always chooseappropriate linear combinations of eigenvectors. general method to solve master equations for quadratic open Fermi systems almost-canonical anti-commutationrelations { ˆ b j , ˆ b k } = 0 { ˆ b j , ˆ b ′ k } = δ j,k { ˆ b ′ j , ˆ b ′ k } = 0 (36)namely ˆ b j could be interpreted as annihilation map and ˆ b ′ j as a creation map of j thNMM, but we should note that ˆ b ′ j is in general not the Hermitian adjoint of ˆ b j [28]. Interms of NMM the Liouvillean (34) now achieves a very convenient normal formˆ L + = − n X j =1 β j ˆ b ′ j ˆ b j − B ˆ (37)where B = A − P nj =1 β j . We shall later show that the constant B is in fact equal to0. The Liouvillean can always be represented in terms of a large but finite 4 n × n matrix.We shall now outline the procedure of complete construction of its spectrum in termsof NMM which are easy to calculate in terms of diagonalization of 4 n × n matrix A asdescribed in the previous subsection.We proceed by constructing the Liouvillean ‘vacuum’. From the representation(22) it follows immediately that h | = h P (0 , ..., | is left-annihilated by ˆ L + , h | ˆ L + = 0,or equivalently ˆ L † + | i = 0. So we have just shown that 0 is always an eigenvalue of ˆ L + ,hence there should also exist the corresponding right eigenvector | NESS i , normalized as h | NESS i = tr ρ NESS = 1, which represents physical NESS, i.e. stationary solutions ofthe Lindblad equation (1)ˆ L + | NESS i = 0 (38)Let us define NMM number maps as ˆ N j := ˆ b ′ j ˆ b j . From eqs. (36) it follows that ˆ N j satisfya projection property ˆ N j = ˆ N j , so they are diagonalizable since no nontrivial Jordanblock could satisfy the projection property. Furthermore, ˆ N j are mutually commuting[ ˆ N j , ˆ N k ] = 0, so they can be simultaneously diagonalized and there should exist avacuum state on which all ˆ N j have value 0. It follows from the stability of completelypositive evolution (1) that all eigenvalues λ of ˆ L + should obey Re λ ≤ 0, and since byassumption Re β j ≥ h | and | NESS i should be the left and right vacua which aresimultaneously annihilated by NMM maps h | ˆ b ′ j = 0 ˆ b j | NESS i = 0 (39)and hence also ˆ N j | NESS i = 0. Thus we have also shown that the NMM representation(37) is only consistent if B = 0 so we find an interesting sum rule for rapidities n X j =1 β j = 2 tr M (40) general method to solve master equations for quadratic open Fermi systems n binary integers (NMM occupationnumbers) ν = ( ν , ν , . . . , ν n ), ν j ∈ { , } , h Θ L ν | ˆ L + = λ ν h Θ L ν | ˆ L + | Θ R ν i = λ ν | Θ R ν i (41)namely λ ν = − n X j =1 β j ν j (42) h Θ L ν | = h | ˆ b ν n n · · · ˆ b ν ˆ b ν | Θ R ν i = ˆ b ′ ν ˆ b ′ ν · · · ˆ b ′ ν n n | NESS i (43)where by construction, left and right eigenvectors satisfy the bi-orthonormality relation h Θ L ν ′ | Θ R ν i = δ ν ′ ,ν . Given a physical observable X ∈ K + and an arbitrary initial state with a density operator ρ ∈ K , the time dependent expectation value of X can be written in terms of thespectral resolution of the Liouvillean,exp( t ˆ L + ) = X ν exp( tλ ν ) | Θ R ν ih Θ L ν | (44)namely h X ( t ) i = tr Xρ ( t ) = tr h X exp( t ˆ L + ) ρ i = X ν exp( tλ ν ) h Θ L ν | ρ X | Θ R ν i (45)We remind the reader that ˆ L + correctly represents physical Liouvillean only on thesubspace K + of operators with even number of a-fermions. However, since the dynamicsis closed on K + and test physical observable X also belongs to K + it follows that thecomponent of ρ from K − does not contribute to the expectation value h X ( t ) i .Given the exact and explicit constructions developed in this section we can nowmake the following rigorous and useful conclusions, assuming throughout that Liouvil-lean shape matrix A is diagonalizable: Theorem 1: NESS of Lindblad equation (1) is unique if and only if the rapidity spec-trum { β j } does not contain 0, in our ordering convention, if β n = 0. In the oppositecase, if we have d ≥ d dimensional convex set ofnon-equilibrium steady states which can be spanned with | Θ R(0 ,..., ,ν ,...,ν d ) i . Theorem 2: An arbitrary initial state with a density operator ρ ∈ K converges withtime to NESS if and only if all rapidities have strictly positive real parts, Re β j > spectral gap ∆ of theLiouvillean which equals ∆ = 2 Re β n . general method to solve master equations for quadratic open Fermi systems Theorem 3: Assume that the rapidity spectrum does not contain 0, i.e. β n = 0.Then the expectation value of any quadratic observable w j w k in a (unique) NESS canbe explicitly computed as h w j w k i NESS = δ j,k + h | ˆ c j ˆ c k | NESS i = (46)= δ j,k + 12 n X m =1 (cid:16) v m, j − v m − , k − − v m, j v m − , k − i v m, j v m − , k − − i v m, j − v m − , k (cid:17) (47)The statements of theorems 1 and 2 simply follow from exact and explicit spectraldecomposition (42,43,44).The proof of theorem 3 is also straightforward: The first expression (46) followsfrom the definition of the annihilation maps (9) and the explicit representation of thedensity operator of NESS, ρ NESS , in the canonical basis P α . The second, very usefulequality (47) is then obtained by expressing ˆ c j thru (24) in terms of NMM maps (35)and using the annihilation relations (39).The quadratic correlator of theorem 3 covers many physically interestingobservables such as densities or currents. However if one needs expectation valuesof more general observables, e.g. an expectation value of a high order monomial h P α i NESS = h | ˆ c α ˆ c α · · · ˆ c α n n | NESS i , then one may use a Wick theorem and rewriteit as a sum of products of pair-wise contractions (46). 3. Trivial example: A single fermion in a bath In order to illustrate the method and demonstrate convenience of the results derived inthe previous section we first work out a simple example of a single fermion n = 1 (orequivalently, an arbitrary qubit, a two-level quantum system), in a thermal bath. Wetake the most general single fermion Hamiltonian H = − i hw w + const = 2 hc † c + const ′ and the following bath operators (see e.g. [16, 29]) L = 12 p Γ ( w − i w ) = p Γ c L = 12 p Γ ( w + i w ) = p Γ c † (48)where the ratio of coupling constants determine the bath temperature T , Γ / Γ =exp( − h/T ). Leaving out the details of a straightforward calculation, simply followingthe steps of the previous section, we arrive at the following shape matrix of theLiouvillean (26) A = − h R + B Γ + , Γ − A = Γ + (49)where R := − − B Γ + , Γ − := i2 Γ + − i2 Γ − Γ − − i2 Γ + Γ − i2 Γ − i2 Γ − − Γ − i2 Γ + − Γ − − i2 Γ − − i2 Γ + (50) general method to solve master equations for quadratic open Fermi systems ± := Γ ± Γ . Further, we compute NMM rapidities β , = Γ + ± i h and theeigenvector matrix V = Γ − Γ + − Γ − Γ + + i − i Γ − Γ + + i Γ − Γ + + 1 − − i4 − i4 14Γ − Γ + + 1 i Γ − Γ + − i i Γ − Γ + + i − Γ − Γ + + 1 14 i4 − i4 14 (51)Then, using theorem 3 we compute the expectation value of occupation number h c † c i = − i2 h w w i = Γ / (Γ + Γ ) which is what we expect in canonical equilibrium. 4. Non-trivial example: transport in quantum spin chains Here we work out a physically more interesting example, namely we construct NESSfor the magnetic and heat transport of a Heisenberg XY spin 1 / H = n − X m =1 (cid:0) J x m σ x m σ x m +1 + J y m σ y m σ y m +1 (cid:1) + n X m =1 h m σ z m (52)which is coupled to two thermal/magnetic baths at the ends of the chain, generated bytwo pairs of canonical Lindblad operators [29] (similar to (48)) L = 12 q Γ L1 σ − L = 12 q Γ R1 σ − n L = 12 q Γ L2 σ +1 L = 12 q Γ R2 σ + n (53)where σ ± m = σ x m ± i σ y m and Γ L , R1 , are positive coupling constants related to bathtemperatures/magnetizations, for example if spins were non-interacting the bathtemperatures T L , R would be given with Γ L , R2 / Γ L , R1 = exp( − h , n /T L , R ).Applying the inverse of Jordan-Wigner transformation (6), σ x m = ( − i) m − Q m − j =1 w j , σ y m = ( − i) m − ( Q m − j =1 w j ) w m , we rewrite (52,53) in terms of Majorana fermions H = − i n − X m =1 ( J x m w m w m +1 − J y m w m − w m +2 ) − i n X m =1 h m w m − w m (54) L = 12 q Γ L1 ( w − i w ) L = − ( − i) n q Γ R1 ( w n − − i w n ) WL = 12 q Γ L2 ( w + i w ) L = − ( − i) n q Γ R2 ( w n − + i w n ) W (55)where W = w w · · · w n is a Casimir operator which commutes with all the elements ofthe Clifford algebra generated by w j and squares to unity W = 1. Therefore, W doesnot affect the action of bath operators (15) where L µ enter quadratically, so we findˆ L + ˆ L = − Γ L+ (ˆ c † ˆ c + ˆ c † ˆ c ) − L − ˆ c † ˆ c † ˆ L + ˆ L = − Γ R+ (ˆ c † n − ˆ c n − + ˆ c † n ˆ c n ) − R − ˆ c † n − ˆ c † n (56) general method to solve master equations for quadratic open Fermi systems n × n shape matrix, which we write in a block tridiagonal form in terms of 4 × A = B L − h R R · · · − R T − h R R . . . − R T − h R ...... . . . . . . R n − · · · − R Tn − B R − h n R (57)and A = Γ L+ + Γ R+ , where B L := B Γ L+ , Γ L − , B R := B Γ R+ , Γ R − (in terms of (50)), withΓ L , R ± := Γ L , R2 ± Γ L , R1 , and R m := J y m 00 0 0 J y m − J x m − J x m (58)We are not able to perform a complete diagonalization of the antisymmetric matrix A (57) of the general XY model analytically. For example, even in the spatiallyhomogeneous case J x , y m ≡ J x , y , h m ≡ h it is not possible to proceed like in theclassical harmonic oscillator chain where the corresponding matrix is a sum of aToeplitz and a bordered matrix [30]. Namely, in our case A is a sum of a blockToeplitz and block bordered matrix and its explicit exact diagonalization remains anopen problem. However, we stress that even relying on numerical diagonalization of A yielding a set of rapidities β j and properly normalized eigenvector matrix V , representsa dramatic progress with respect to previously existing numerical methods which neededdiagonalization of matrices which were exponentially large in n . We shall later derivesome exact theoretical and analytical results, explaining results of exact numericalcomputations, in the special case of a homogeneous transverse Ising chain (subsection4.1), and the case of a disordered XY chain (subsection 4.2) for which we shall relateNMM to the problem of Anderson localization in one dimension,Let us continue by discussing transport observables in the spin chain whoseexpectation values in NESS are easy to calculate. Note that the bulk Hamiltonian(52) can be written in terms of the two-body energy density operator H m = − i J x m w m w m +1 + i J y m w m − w m +2 − i h m w m − w m − i h m +1 w m +1 w m +2 (59)as H = P m H m . One can derive the local energy current Q m = i[ H m , H m +1 ] from the continuity equation (d / d t ) h H m i = tr H m d ρ/ d t = h i[ H, H m ] i = −h Q m i + h Q m − i (60)where Q m := i[ H m , H m +1 ] Q m = 2i(2 J y m J x m +1 w m − w m +3 + 2 J x m J y m +1 w m w m +4 − general method to solve master equations for quadratic open Fermi systems − J y m h m +1 w m − w m +1 − J x m h m +1 w m w m +2 −− h m +1 J x m +1 w m +1 w m +3 − h m +1 J y m +1 w m +2 w m +4 ) (61)The validity of the above continuity equation (60) depends on two assumptions only:(i) All Lindblad operators L µ commute with the density H m in the bulk , 2 ≤ m ≤ n − H m , H m ′ ] = 0 if | m − m ′ | ≥ H m and energy current Q m , and also of somewhat simpler spin density σ z m = − i w m − w m (62)and spin current S m = σ x m σ y m +1 − σ y m σ x m +1 = − i w m w m +2 − i w m − w m +1 (63)which are all quadratic in w j . Note, however, that the spin density satisfies continuityequation (d / d t ) h σ z m i = −h S m i + h S m − i only in the isotropic case, when J x m ≡ J y m . Here we limit our discussion to the spatially homogeneous case J x , y n ≡ J x , y , h n ≡ h . Weshall show that in this case the eigenvalue problem A v = βv (64)for (57) can be most easily and elegantly treated if formulated in terms of anabstract inelastic scattering problem in one dimension, with asymptotic solutionsgiven in terms of free normal modes for the infinite translationally invariant chain v = ( . . . , uξ m − , uξ m , uξ m +1 , . . . ) T , where ξ is a complex quasi–momentum (Bloch)parameter and u is a 4-dimensional amplitude vector satisfying the condition( − R T ξ − − h R + R ξ − β I ) u = 0 (65)and the baths playing the role of inelastic (absorbing) scatterers at the edges of a finitelattice. The ‘elastic’ (Hamiltonian) version of this trick has been used to computetemporal correlation functions in kicked Ising chain [31].The singularity condition for the free mode equation (65) results, for a generalhomogeneous XY model, in eight master bands - different values of momenta ξ for eachvalue of the spectral parameter (rapidity) β . In order to simplify the discussion - whichwill still get rather involved - we shall in the following restrict ourselves to the transverseIsing model J x = J, J y = 0. In this case we find just two master bands with simpledispersion relations ξ ± ( β ) = h + J + β ± ω ( β )2 hJ ω ( β ) := p ( h + J + β ) − (2 hJ ) (66)but each band is doubly-degenerate, since the corresponding amplitude problem (65)has two linearly independent solutions u ± ( β ) = − h ( h − J + β ± ω )0 β ( h + J + β ± ω )0 u ± ( β ) = − h ( h − J + β ± ω )0 β ( h + J + β ) ± ω ) (67) general method to solve master equations for quadratic open Fermi systems ξ − represents left moving and ξ + right moving free modes, each havingtwo possible polarizations. Note that ξ − ξ + = 1. For a general complex β we shall choosethe branch of square root ω ( β ) (66) for which | ξ − | ≤ 1. Let us now write the scatteringproblem on the left bath in terms of an ansatz v = u L ψ − u − + ψ − u − + ψ +1 u +1 + ψ +2 u +2 ( ψ − u − + ψ − u − ) ξ − + ( ψ +1 u +1 + ψ +2 u +2 ) ξ + ( ψ − u − + ψ − u − ) ξ − + ( ψ +1 u +1 + ψ +2 u +2 ) ξ ... (68)where u L represents a 4-dimensional vector of left-most eigenvector components, ψ − , are the amplitudes of (known) incident free modes, and ψ +1 , are the amplitudes of thescattered, outgoing free modes. Plugging the scattering ansatz to the eigenproblem (64),the first two rows of A (57) result in 6 linearly independent equations for 6 unknowns ψ +1 , , u L . Eliminating four variables u L we finally arrive to the non-unitary 2 × (cid:18) ψ +1 ψ +2 (cid:19) = S L (cid:18) ψ − ψ − (cid:19) (69)with S L11 = τ − β ( − (Γ L+ ) + 4(Γ L+ ) ( β − h ) − h ( hJ + iΓ L − ω )) S L12 = τ − β ((Γ L+ ) + 8iΓ L − hβ + 4Γ L+ ( h − β ))(2i ω ) S L21 = τ − β ((Γ L+ ) − L − hβ + 4Γ L+ ( h − β ))( − ω ) S L22 = τ − β ( − (Γ L+ ) + 4(Γ L+ ) ( β − h ) − h ( hJ − iΓ L − ω )) (70) τ := (Γ L+ ) β + 8 β ( h + ( J + β )( J + β − ω ) + h (2 β − ω )) − L+ ) ( h + J + 3 β + J (2 β − ω ) − β ω + h ( ω − J − β ))Similarly, one can solve the scattering problem from the right bath with the scatteringansatz v = ...( ψ +1 u +1 + ψ +2 u +2 ) ξ − + ( ψ − u − + ψ − u − ) ξ − − ( ψ +1 u +1 + ψ +2 u +2 ) ξ − + ( ψ − u − + ψ − u − ) ξ − − ψ +1 u +1 + ψ +2 u − + ψ − u − + ψ − u − u R (71)defining the right S-matrix (cid:18) ψ − ψ − (cid:19) = S R (cid:18) ψ +1 ψ +2 (cid:19) (72)Note that since the two directions of free modes (67) do not have left-right symmetry anexplicit expression for S R is considerably more complicated than (70) and shall not bewritten out here. We shall now show that there exist two qualitatively different typesof NMM - complex rapidities β solving (64) for sufficiently large n .First, we shall discuss the so called evanescent normal master modes . These arecharacterized with amplitudes (68) which decay exponentially with the distance from general method to solve master equations for quadratic open Fermi systems Β- - Β Β- - Β Β- - Β Figure 1. Rapidities β j (black dots) for a transverse Ising chain with J = 1 . , h = 1and bath couplings Γ L1 = 1 , Γ L2 = 0 . 6, Γ R1 = 1 , Γ R2 = 0 . 3, for three different sizes n = 6(upper), n = 30 (middle), and n = 150 (lower panel). Big blue/red dots indicatepositions of evanescent rapidities (solutions of eq.(74)) for the left/right bath. – say the left – bath, so the other – the right boundary condition becomes physicallyirrelevant in the limit n → ∞ . Such solutions ψ +1 , = 0 of eq. (69) exist exactly whenthe determinant of S-matrix vanishes det S L = 0. Using (70) the determinant can bewritten as det S L = ( β/τ ) p L ( β ) where ¶ p L ( β ) = (Γ L+ ) β − L+ ) (( h − J ) + (2 J − h ) β + 3 β ) − L+ ) (2 h ( h − J ) − (7 h − h J +2 J ) β + 4( h − J ) β − β ) − L+ ) ( h ( h − J ) − h J β − (2 h +4 h J − J ) β + 2 J β + β )+ 256 h J β (73) ¶ Trivial zero β = 0 of course does not represent a physical solution since then the whole S-matrix(70) vanishes. general method to solve master equations for quadratic open Fermi systems β p L ( β evan ) = 0 (74)that are not simultaneously zeroes of τ ( β ). Clearly, such NMM asymptotically do notdepend on the chain size n . In addition, we find evanescent NMM correspondingto the right bath simply by replacing Γ L+ by Γ R+ in (74,73). In fig. 1 we compareevanescent rapidities computed from eq. (74) to numerically calculated spectrum of A ,at several different sizes n , for a typical case of transverse Ising chain, J = 1 . , h = 1 . L1 = 1 , Γ L2 = 0 . R1 = 1 , Γ R2 = 0 . soft normal master modes withrapidities that are closest to the imaginary axis, and thus determining the spectralgap of the Liouvillean and relaxation time to NESS. Composing the scattering from thetwo baths with the free propagation along the chain (back and forth) we arrive at thegeneral secular equation for the eigenvalue problem (64) in terms of a 2 × ξ n − S R S L − I ) = 0 (75)In the absence of the baths, Γ L , R ± = 0, the solutions of the above problem exist only forreal quasi-momenta, namely ξ ± = exp( ± i ϑ ) , ϑ ∈ R . For such extended master modesthe local coupling to the baths can be considered as a small perturbation, thus onlyslightly perturbing the Bloch-like bands β ( e i ϑ ) = ± i ε ( ϑ ) with ‘energy’ ε ( ϑ ) = p h + J − | hJ | cos ϑ (76)The softest NMM, namely the one for which the coupling to the baths is expected tobe the weakest, should have nearly nodes at the ends of the chain, i.e. ϑ ≈ π/n , or ϑ ≈ π + π/n , and should thus lie near the band edges ± i | h | ± i | J | (see fig. 1). In thefollowing we shall focus our calculation on the band edge β ∗ = i( | h | + | J | ) which, as canbe checked aposteriori by a straightforward but tedious calculation, always gives smallerreal part of the rapidity than the lower edge i( | h | − | J | ), and hence really determinesthe gap of the Liouvillean. So we write β = i( | h | + | J | ) + z (77)where z ∈ C is a small parameter, and expand the S-matrices around the band edge S L , R = − I + 4 gη L , R Z L , R √− i z + O ( | z | ) (78)where g := q | hJ | | h | + | J | ) , η L , R := (Γ L , R+ ) + 4(Γ L , R+ ) (4 h + 2 | hJ | + J ) + 16 h J and Z L11 = 4 | h | (Γ L+ ) + 16 | h | ( | h | + | J | )( | J | − iΓ L − ) Z L12 = − L+ ) − L − | h | ( | h | + | J | ) − h + 2 | hJ | + J ) Z L21 = + 2(Γ L+ ) − L − | h | ( | h | + | J | ) + 8(2 h + 2 | hJ | + J ) Z L22 = 4 | h | (Γ L+ ) + 16 | h | ( | h | + | J | )( | J | + iΓ L − ) (79) general method to solve master equations for quadratic open Fermi systems Z R11 = (Γ R+ ) (2 | h | + | J | ) + 4(Γ R+ ) (8 | h | + 9 h | J | + 4 | h | J + | J | )+ 16 h | J | ( | J | (3 | h | + 2 | J | ) − iΓ R − ( | h | + | J | )) Z R12 = − R+ ) − R − | h | ( | h | + | J | ) − h + 2 | hJ | + J ) Z R21 = + 2(Γ R+ ) − R − | h | ( | h | + | J | ) + 8(2 h + 2 | hJ | + J ) Z R22 = (Γ R+ ) (2 | h | + | J | ) + 4(Γ R+ ) (8 | h | + 9 h | J | + 4 | h | J + | J | )+ 16 h | J | ( | J | (3 | h | + 2 | J | ) + iΓ R − ( | h | + | J | )) (80)Next we expand ξ + (66) in z , yielding ξ + = − − g − √− i z + O ( | z | ) (81)and so the free propagator in (75) can be written as ξ n − = exp(2 ng − √− i z ) + O ( | z | ) . (82)In eqs. (78,81,82) the branch cut along the negative real axis has been chosen for √− i z .Since the product of S-matrices in (75) is near identity, the free propagator should benear one as well, hence 2 ng − √− i z should be near 2 π i. Let us define z by setting2 ng − √− i z = 2 π i, so z = − i π g n − (83)and write z = z (1 + y ) where | y | ≪ z is purely imaginary, we need to compute a small but non-vanishing y whichwill, in the leading order in n , solve (75) since the real part of the soft mode’s rapidityis determined asRe β = Re z y = π g n − Im y (84)Now, writing √− i z = √− i z √ y = i πgn − (1 + y/ − y / 8) + O ( y ) in(78,82), plugging all that to eq. (75) and computing to order O ( n − ), noting that O ( | z | ) = O ( n − ), we arrive to a simple quadratic equation for y , whose solution, pluggedto (84), gives the final result, namely the sectral gap of Liouvillean ∆ = 2 Re β ∆ = (2 πhJ ) ( | h | + | J | ) ∆ ∆ n − + O ( n − ) (85)∆ := 64(Γ L+ + Γ R+ ) h J (2 h + 2 | hJ | + J )+ 16((Γ L+ ) + (Γ R+ ) ) h J + 16Γ L+ Γ R+ (Γ L+ + Γ R+ )(2 h + 2 | hJ | + J )(4 h + 2 | hJ | + J )+ 4Γ L+ Γ R+ ((Γ L+ ) + (Γ L+ ) Γ R+ + Γ L+ (Γ R+ ) + (Γ R+ ) )(2 h + 2 | hJ | + J )+ (Γ L+ Γ R+ ) (Γ L+ + Γ R+ )∆ := ((Γ L+ ) + 4(Γ L+ ) (4 h + 2 | hJ | + J ) + 16 h J ) × ((Γ R+ ) + 4(Γ R+ ) (4 h + 2 | hJ | + J ) + 16 h J )In fig. (2) we compare this analytical result to exact numerical calculations of theeigenvalue of A with minimal real part, confirming both, its precise numerical valueand that the relative scaling of the next order correction is indeed O ( n − ). general method to solve master equations for quadratic open Fermi systems - - 101 Log n L og ¡ n D - C ¥ n n D Figure 2. Spectral gap ∆ times a third power of the chain length n for a transverseIsing chain with J = 1 . , h = 1 and bath couplings Γ L1 = 1 , Γ L2 = 0 . 6, Γ R1 = 1 , Γ R2 = 0 . n in log-log scale and demonstratethat it decays as ∝ n − (thin line). - - - - - - - - Λ I m Λ Figure 3. Complete spectrum of 2 complex eigenvalues of Liouvillean for atransverse Ising chain with n = 6 spins and J = 1 . , h = 1 and bath couplingsΓ L1 = 1 , Γ L2 = 0 . 6, Γ R1 = 1 , Γ R2 = 0 . Note that, interestingly, both main analytical results of this subsection, namelyevanescent and soft mode rapidities do not depend on Γ L , R − . Physically speaking, theyonly depend on the effective strengths of the bath couplings and not on the temperatures.We end this subsection by presenting some further numerical results on heattransport in the open transverse Ising chain in the Lindblad form. In fig. 3 wedemonstrate expression (42) for constructing the full spectrum of the Liouvillean in general method to solve master equations for quadratic open Fermi systems n < S > , < Q > Figure 4. Energy current (upper/blue points), and average spin current (lower/redpoints), versus chain length n for a transverse Ising chain with J = 1 . , h = 1 and bathcouplings Γ L1 = 1 , Γ L2 = 0 . 6, Γ R1 = 1 , Γ R2 = 0 . - - - - - m < H > , m < > Σ m z - - - - - L og ∆ H m 70 72 74 76 78 80 - - - - - L og ∆ Σ m z Figure 5. Energy density profile (lower, blue points), and spin density profile(upper, red points), for a transverse Ising chain of n = 80 spins with J = 1 . , h = 1and bath couplings Γ L1 = 1 , Γ L2 = 0 . 6, Γ R1 = 1 , Γ R2 = 0 . 3. The insets displaylogarithm of the difference to the bulk values δH m := |h H m i − H bulk | (blue points), δσ z m := |h σ z m i − σ zbulk | (red points) in comparison with ± (4 log ξ − ) m + const withquasi-momentum ξ − = 0 . β evan = 0 . terms of a set of rapidities, for a short chain. In fig. 4 we demonstrate eq. (47) ofTheorem 3 by computing the energy current Q m (61), and the average spin current S = n − P n − m =1 S m (63) in NESS of a typical transverse Ising chain. Numerical resultsgive a clear indication of ballistic transport h Q i = O ( n ) , h S i = O ( n ), however itsrigorous proof and analytical calculation of the currents would require full control overthe complete set of NMM which is at present not available. In fig. 5 we plot the energy general method to solve master equations for quadratic open Fermi systems 10 20 30 40 50 - - - - - - - n L og < D > Figure 6. Average Liouvillean spectral gap h ∆ i versus the chain length n fordisordered XY models: (i) J x m = 0 . , J y m = 0 , h m ∈ [1 , 2] (transverse Ising withfield disorder, blue points), (ii) J x m ∈ [0 . , , J y m = 0 , h m = 1 (transverse Isingwith interaction disorder, red points), (iii) J x m ∈ [0 . , , J y m ∈ [0 . , , h m = 1 (XYwith interaction disorder, golden points), all for bath couplings Γ L1 = 1 , Γ L2 = 0 . R1 = 1 , Γ R2 = 0 . 3. Full lines indicate exponential fits to right halves of data. Averagingis performed over 2000 disorder realizations. density (59) and spin density (62) profiles in NESS. Again, we note flat profiles in thebulk of the chain, m, n − m ≫ 1, with exponential falloff due to adjustment to thenon-equilibrium bath values. Since the densities can be written, by means of (47),as 4 − point functions in NMM components, the leading falloff exponents of the profile |h H m i − H bulk | ∼ | ξ − | m is given by the quasi-momentum ξ − (66) corresponding to themaximal evanescent rapidity β evan (74). In this subsection we treat the opposite extreme, a disordered XY chain (52) where threesets of physical parameters are chosen as random uncorrelated variables from uniform distributions on the intervals, J x m ∈ [ J xmin , J xmax ], J y m ∈ [ J ymin , J ymax ], h m ∈ [ h min , h max ].Clearly, the eigenvalue problem (64) for the matrix (57) then becomes equivalent to theAnderson tight-binding problem in one dimension for a quantum particle with a 4 − levelinternal degree of freedom. We do not pursue any theoretical analysis of this problemhere, but merely state that numerical investigations indicate existence of exponentiallocalization of all eigenvectors (or normal master modes) for disorder of any strength inanyone of system’s parameters.With the picture of localization of NMM in mind, the effect of the couplings to theheat baths at the chain’s ends on quantum transport can be predicted by theoreticalarguments (see [32] for a review of related studies): The spectral gap of the Liouvilleanshould be exponentially small ∆ ∼ exp( − n/ℓ ) where ℓ is the localization length of NMM general method to solve master equations for quadratic open Fermi systems 10 20 30 40 50 - - - - - - n L og < Q > Figure 7. The scaling of the energy current h Q m i with chain length n of the disorderedXY model in the same regimes/parameters/plot styles as in fig. 6. - - - - - m - n - < H > m Figure 8. Scaled energy density profile of interaction disordered XY chains (case(iii) of fig. 6) for three chain sizes: n = 20 (blue points), n = 40 (red points), n = 60(golden points). Averages over 50000 disorder realizations have been performed. which is expected to be proportional to the square of inverse disorder strength. Thisis demonstrated in fig.6. If all NMM are exponentially localized, the currents shoulddecrease with the chain size n faster than any power, perhaps exponentially, and thesystem should behave as an ideal insulator (at all temperatures). This is demonstratedby straightforward numerical calculations of the heat current (61) in fig. 7. In the finalfigure 8 we plot the energy density profile h H m i (59) in a typical case of disordered XYchain, versus a scaled spatial coordinate ( m − / ( n − ∈ [0 , n , and demonstrate sharping up of energy density profiles with increasing n , which is again indicating insulating behaviour. general method to solve master equations for quadratic open Fermi systems 5. Discussion and conclusions The main result of the paper is a general method of explicit solution of master equationsdescribing dynamics of open quantum system, under the condition that the system’sHamiltonian is quadratic and all Lindblad operators are linear in canonical fermionicoperators (which can either represent real physical fermions or any abstract 2-levelquantum systems (qubits) thru the Jordan-Wigner transformation). Using a novelconcept of Fock space of physical operators (or density operators of physical states), andthe adjoint structure of canonical creation and annihilation maps over this space, theproblem can be treated in terms of a non-Hamiltonian generalization of the method ofLieb, Schultz and Mattis [18] lifted to an operator space. We have explicitly constructeda non-canonical analog of Bogoliubov transformation of the quantum Liouville map tonormal master modes. Related ideas in the Hamiltonian context have been used by theauthor [31, 33, 34] in order to approach the problem of real time dynamics and ergodicproperties of isolated interacting many-body quantum systems.As an illustration of the method we have solved far from equilibrium quantumheat and spin transport problem in Heisenberg XY spin 1/2 chains which are coupledto canonical heat baths only at the two ends. Irrespectively of the strength of thecoupling to the baths and their temperatures, we have shown a ballistic transport inthe spatially homogeneous (non-disordered) case, and an ideally insulating behaviour inthe disordered case associated to localization of normal master modes of the quantumLiouville operator. In this context the method can be considered as a simple alternativeto the solution of quantum Langevin equations [24].However, the method should easily be applicable to variety of other physicalsituations, for example if all fermions are coupled to the baths one could make asolvable model of genuine quantum diffusion, a many-body generalization of the tight-binding model [35]. We also expect the method to be applicable to the Redfield typeof master equations (see e.g. [35]) - which do not conserve positivity for a short initial(slippage) time interval - provided only the system part of the Hamiltonian is quadratic and system-bath couplings are linear in fermionic variables. Furthermore, extension ofthe method to open many-boson systems should be straightforward, simply by replacinganticommutators by commutators throughout the exposition of section 2.Treating density operators of NESS as elements of a Hilbert space of operators onemay also extend the concept of entanglement entropy , with respect to a bipartition of asystem of many fermions [36], to NESS which can in our approach be viewed as a kind ofground state of the Liouvillean. Saturation of such operator space entanglement entropy [34] (which is suggested by numerical experiments [37]) indicates efficient simulability of NESS by elaborate numerical methods such as density matrix renormalization group [38], perhaps even for more general, non-solvable quantum systems.As last we mention a more ambitious extension of the present work: Namely wepropose to explore a question, whether more involved algebraic methods of solution ofinteracting many-body quantum systems, like e.g. Bethe Ansatz or quantum inverse general method to solve master equations for quadratic open Fermi systems Acknowledgements I gratefully acknowledge stimulating discussions with Pierre Gaspard, Keiji Saito andWalter Strunz, thank Carlos Mejia-Monasterio and Thomas H. Seligman for readingthe manuscript and many useful comments, and Iztok Piˇzorn and Marko ˇZnidariˇc forcollaboration on related projects. The work has been supported by the grants P1-0044 and J1-7347 of Slovenian research agency (ARRS). Explicit analytical calculationsreported in subsection (4.1) were assisted by Mathematica software package. References [1] W. H. Zurek, Rev. Mod. Phys. , 715 (2003).[2] E. Joos, H. D. Zeh, C. Kiefer, D. Giulini, J. Kupsch and I.-O. Stamatescu, Decoherence and theAppearance of a Classical World in Quantum Theory , (Springer, 2003).[3] J. von Neumann, Mathematical Foundations of Quantum Mechanics , Trans. Robert T. Geyer.(Princeton University Press, Princeton 1955).[4] M. Schlosshauer, Rev. Mod. Phys. , 1267 (2004).[5] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (CambridgeUniversity Press, Cambridge 2000).[6] G. Benenti, G. Casati and G. Strini, Principles of Quantum Computation and Information. VolumeI: Basic Concepts (World Scientific, Singapore 2004); Volume II: Basic Tools and Special Topics (World Scientific, Singapore 2007).[7] H. Araki and E. Barouch, J. Stat. Phys. , 327 (1983);H. Araki, Publ. RIMS Kyoto Univ. , 277 (1984).[8] D. Ruelle, J. Stat. Phys. , 57 (2000).[9] V. Jakˇsiˇc and C.-A. Pillet, J. Stat. Phys. , 787 (2002); Commun. Math. Phys. , 131 (2002);W. Aschbacher, V. Jakˇsiˇc, Y. Pautrat and C.-A. Pillet, Inroduction to non-equilibrium quantumstatistical mechanics , in Open Quantum Systems III. Recent Developments Lecture Notes inMathematics, (2006), 1-66.[10] P. Gaspard, Prog. Theor. Phys. Suppl. , 33 (2006); Physica A , 201 (2006).[11] M. H. Lee, Acta Physica Polonica B , 1837 (2007); Phys. Rev. Lett. , 250601 (2001); Phys.Rev. Lett. , 1072 (1982).[12] T. Prosen, J. Phys. A: Math. Theor. , 7881 (2007).[13] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, Quantum Inverse Scattering and Correlationfunctions (Cambridge University Press, Cambridge 1997).[14] L. Fadeev, P. Van Moerbeke and F. Lambert (Eds.), Bilinear Integrable Systems: from Classical toQuantum, Continuous to Discrete , (NATO ARW Proceedings), Springer Series: NATO ScienceSeries II: Mathematics, Physics and Chemistry, Vol 201 (2006).[15] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press,Oxford 2002).[16] F. Haake, Quantum Signatures of Chaos , 2nd edition (Springer, 2001).[17] R. Alicki and K. Lendi, Quantum dynamical semigroups and applications (Springer, 2007).[18] E. H. Lieb, T. D. Schultz and D. C. Mattis, Ann. Phys. (New York) , 407 (1961).[19] X. Zotos, F. Naef and P. Prelovˇsek, Phys. Rev. B , 11029 (1997). general method to solve master equations for quadratic open Fermi systems [20] K. Saito, S. Takesue and S. Miyashita, Phys. Rev. E , 2397 (2000);K. Saito, Europhys. Lett. , 34 (2003).[21] M. Michel, M. Hartmann, J. Gemmer and G. Mahler, Eur. Phys. J. B. , 325 (2003);M. Michel, G. Mahler and J. Gemmer, Phys. Rev. Lett. , 180602 (2005);M. Michel, J. Gemmer and G. Mahler, Int. J. Mod. Phys. B , 4855 (2006);J. Gemmer, R. Steinigeweg and M. Michel, Phys. Rev. B , 104302 (2006).[22] C. Mejia-Monasterio, T. Prosen and G. Casati, Europhys. Lett. , 520 (2005);G. Casati and C. Mejia-Monasterio, e-print arXiv:0710.3500v1 [cond-mat.stat-mech] .[23] C. Mejia-Monasterio and H. Wichterich, e-print arXiv:0709.1412v1 [cond-mat.stat-mech] .[24] A. Dhar and D. Roy, J. Stat. Phys. , 805 (2006); see also e-print arXiv:0711.4318v1[cond-mat.stat-mech] .[25] G. Lindblad, Commun. Math. Phys. , 119 (1976).[26] M. M. Wolf, J. Eisert, T. S. Cubitt and J. I. Cirac, e-print arXiv:0711.3172v1 [quant-ph] .[27] P. ˇSemrl, private communication .[28] We note a similarity to the formalism of second quantization with non-orthogonal orbitalsintroduced in: M. Moshinsky and T. H. Seligman, Ann. Phys. (New York) , 311 (1971).[29] H. Wichterich, M. J. Herich, H. P. Breuer, J. Gemmer and M. Michel, Phys. Rev. E, , 1073 (1967).[31] T. Prosen, Prog. Theor. Phys. Suppl. , 191 (2000).[32] S. Lepri, R. Livi and A. Politi, Phys. Rep. , 1 (2003).[33] T. Prosen, Phys. Rev. E , 1658 (1999).[34] T. Prosen and I. Piˇzorn, Phys. Rev. A , 032316 (2007).[35] M. Esposito and P. Gaspard, Phys. Rev. B , 214302 (2005); J. Stat. Phys. , 463 (2005).[36] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. , 227902 (2003);J. I. Latorre, E. Rico, and G. Vidal, Quant. Inf. Comp. , 48 (2004).[37] T. Prosen and M. ˇZnidariˇc, to be submitted (2008).[38] S. R. White, Phys. Rev. Lett. , 2863 (1992);U. Schollw¨ock and S. R. White, in G. G. Batrouni, and D. Poilblanc (eds.): Effective modelsfor low-dimensional strongly correlated systems, p.155, AIP, Melville, New York (2006);G. Vidal, Phys. Rev. Lett. , 147902 (2003); ibid.93