lopscience

¡opscience.iop.org

Home Search Collections Journals About Contact us My IOPscience

Exact solution of Markovian master equations for quadratic Fermi systems: thermal baths, open XY spin chains and non-equilibrium phase transition

This content has been downloaded from IOPscience. Please scroll down to see the full text. View the table of contents for this issue, or go to the journal homepage for more

Download details:

IP Address: 128.122.253.212

This content was downloaded on 02/02/2015 at 12:10

Please note that terms and conditions apply.

New Journal of Physics

The open-access journal for physics

Exact solution of Markovian master equations for quadratic Fermi systems: thermal baths, open XY spin chains and non-equilibrium phase transition

Tomaz Prosen and Bojan ZunkoviC

Department of Physics, FMF, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia E-mail: tomaz.prosen@fmf.uni-lj.si

New Journal of Physics 12 (2010) 025016 (23pp)

Received 18 September 2009 Published 26 February 2010 Online at http://www.njp.org/

doi:10.1088/1367-2630/12/2/025016

Abstract. We generalize the method of third quantization to a unified exact treatment of Redfield and Lindblad master equations for open quadratic systems of n fermions in terms of diagonalization of a 4n x 4n matrix. Non-equilibrium thermal driving in terms of the Redfield equation is analyzed in detail. We explain how one can compute all the physically relevant quantities, such as non-equilibrium expectation values of local observables, various entropies or information measures, or time evolution and properties of relaxation. We also discuss how to exactly treat explicitly time-dependent problems. The general formalism is then applied to study a thermally driven open XY spin 1/2 chain. We find that the recently proposed non-equilibrium quantum phase transition in the open XY chain survives the thermal driving within the Redfield model. In particular, the phase of long-range magnetic correlations can be characterized by hypersensitivity of the non-equilibrium steady state to external (bath or bulk) parameters. Studying the heat transport, we find negative differential thermal conductance for sufficiently strong thermal driving as well as non-monotonic dependence of the heat current on the strength of the bath coupling.

New Journal of Physics 12 (2010) 025016 1367-2630/10/025016+23$30.00

© IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

Contents

1. Introduction 2

2. Markovian master equations in non-equilibrium quantum physics 3

3. Diagonalization of quantum Liouvilleans for quadratic Fermi systems 4

3.1. Fock space of operators............................................................5

3.2. Bilinear form of the Liouvillean..................................................5

3.3. Static Liouvillean: normal modes, NESS and decay spectrum....................8

3.4. Static Liouvillean: time-dependent correlation functions..........................9

3.5. Time-dependent Liouvilleans......................................................10

4. XY spin chains 11

4.1. Non-equilibrium phase transition..................................................12

4.2. Heat transport and entropies ......................................................17

5. Discussion 21 Acknowledgments 22 References 22

1. Introduction

One of the main challenges of many-body theory and non-equilibrium statistical mechanics is to understand the properties of relaxation of large interacting quantum systems. There are two common approaches to these types of problems. One important direction is to try to define dynamics in the thermodynamic limit and to investigate its properties with rigorous mathematical methods of operator algebras [1]-[3]. However, in this context explicit results that go beyond existence proofs are quite limited. A second approach is to split a large system into a tensor product of a smaller system of interest, and the rest (environment), and trying to eliminate all the degrees of freedom of the large, macroscopic environment (see e.g. [4, 5]). This approach, although involving a series of approximations, is usually more fruitful for explicit calculations and quantitative analyses. We may be interested in relaxation to either equilibrium or non-equilibrium steady states (NESSs), depending on the equal or non-equal values of thermodynamic potentials assigned to possibly several pieces of environment—which we shall call the baths. Such calculations of the quantitative properties of steady states may be very useful, for example, in the realm of transport theory [6], and may complement the linear response calculations and suggest nonlinear responses or far-from-equilibrium effects.

However, to date we have had very few explicit calculations of the non-equilibrium properties of open many-body quantum systems, and mainly they had to focus on small systems with a single or a pair of degrees of freedom (such as spins or bosons); see for example [7, 8]. The reason is that there has been no theoretical technique to deal with open many-body problems except for the Keldysh formalism of non-equilibrium Green's functions, which however can easily get too involved for explicit calculations. Recently, two new directions have been proposed, both in the direction of numerical simulation and theoretical analysis. Namely, in the context of numerical simulations of open many-body systems, time-dependent density matrix renormalization group techniques [9] have been demonstrated to efficiently simulate relaxation to steady states with the Lindblad master equation [10]. On the other hand, it has been shown [11] that the Lindblad equation for general quadratic fermionic systems, for example, for

XY-like quantum spin chains that are mappable to quadratic fermionic systems, can be solved explicitly with the technique of canonical quantization in the Fock space of operators—third quantization for short.

In this paper, we shall show how the third quantization can be generalized to treat quadratic systems with arbitrary Markovian master equations, which are not necessarily of the Lindblad form. In particular, we shall focus on the Redfield dissipator in terms of which we can simulate simple thermal reservoirs, and thermal driving of the system under non-equilibrium conditions. After giving a short account on mathematical formulation of Markovian master equations and the basic physical assumptions and approximations involved in the derivation (in section 2), we shall in section 3 present a short but self-contained generalization of the theory [11]. In addition, we shall outline the calculation of dynamical correlation functions in Liouvillean dynamics, and formulate an exact treatment of explicitly time-dependent quantum Liouville problems. In section 4, we shall apply our technique to treat an open XY spin chain in the non-equilibrium Redfield model. We shall outline several intriguing exact numerical results on large spin chains. In particular, we show that the recently announced quantum phase transition in the open XY chain in the local Lindblad bath model generalizes also to the non-equilibrium thermal Redfield model with qualitatively identical characteristics. The transition is characterized by the spontaneous emergence of long-range magnetic correlations, and hypersensitivity of the steady state to the external system's parameters, when the transverse magnetic field drops below the critical value |h| < hc = 11 — y21, where y is the anisotropy parameter. Furthermore, we analyze in some detail the heat transport in the XY chain, and find regions of negative differential heat conductance for strong thermal driving, namely non-monotonic dependence of the heat current on the temperature difference between the baths.

2. Markovian master equations in non-equilibrium quantum physics

Decomposing the Hilbert space of the universe into a tensor product H = Hs ® Hb of the central system Hs and the bath (or a set of baths) Hb (environment), one writes the total Hamiltonian as

H = Hs ® 1b + 1s ® Hb + Xm ® YM, (1)

where XM are linear operators over Hs, and YM are linear operators over Hb. Note that XYM can always be chosen to be Hermitian, so this shall be assumed throughout this paper. The Markovian quantum master equation for the time evolution of the central system's density matrix p(t) is derived [4] using three main assumptions: (i) weak coupling (assuming X to be small), (ii) factorizability of the initial density matrix ps(0) ® pb(0) and (iii) Born-Markov approximation, which rests upon the assumption that the bath-correlation functions

r£v(t) := X2tr(Y^(t)Yve—¡Hb)/tre-¡Hb, Y^t) := eitHbYMe—itHb (2)

decay on much shorter timescales than the central system's dynamics XM(t) := eitHsXMe—itHs. We use units in which Planck's constant h = 1, and may use different inverse temperatures ¡5 for different pieces of the environment (for different baths). The resulting master equation is referred to as the Redfield equation

d p(t) = —i[ Hs ,p(t)] + Dp(t), (3)

where the dissipator map has a memoryless kernel with the following general form:

Dp = T dTr^(T)[Xc^-T)p, Xv] + h.c. (4)

If one additionally assumes the so-called rotating wave-approximation, one arrives at the dynamical semi-group that manifestly preserves the positivity of the density matrix at all times1 and can be generally described by the dissipator in the Lindblad form

Dp = £ Yv . m [X¡P. Xv] + h.c. . (5)

where the only condition is that y is a Hermitian yM . v = yv* and a positive definite matrix. The standard Lindblad form is obtained by diagonalizing the matrix y whose eigenvectors yield the usual Lindblad operators. The important property of the bath-correlation functions (2) (which constitute all that we need to know about the baths) is the Kubo-Martin-Schwinger (KMS) condition

rft . v(-t - \P) = r£ ¡(t) . (6)

which is needed to prove that the thermal state pgibbs = e-^Hs/ tre-^Hs is a steady state of the master equation (3), provided that all baths are thermalized to the same inverse temperature2. However, in the case of several thermal baths with possibly different temperatures we may expect that p(t) relaxes to a physically very interesting NESS.

3. Diagonalization of quantum Liouvilleans for quadratic Fermi systems

In this section we give a short account on the general technique of canonical quantization in Liouville space ('third quantization') and complete diagonalization of Markovian master equations (3), with (4) or (5), for quadratic fermionic problems. We treat a finite problem with n fermionic degrees of freedom, described by 2n anti-commuting Hermitian operators w j, j = 1. 2. .... 2n, {wj. wk} = 28jk, in which the Hamiltonian H may take a general quadratic form and the coupling operators may be general linear forms:

Hs = ^ wjHjkwk = w • H w, (7)

Xm = xM.jwj = xi • w. (8)

Thus, 2n x 2n matrix H can be chosen to be antisymmetric HT = —H. Throughout this paper x = (x 1. x2... .)T will designate a vector (column) of appropriate scalar-valued or operator-valued symbols xk. This formalism is immediately applicable for describing either (i) physical fermions cm, m = 1. 2. .... n, where w2m—1 = cm + cm, w2m = i(cm — c^m), or

1 This is not the case for equations (3) and (4) which allow for possible breaking of positivity at an initial short time interval, the so-called sleapage time.

2 With an additional technical condition of neglecting the Cauchy principal value contribution to the time integral (4), see the discussion at the end of section 3.2.

(ii) XY-like systems of spins 1/2 with canonical Pauli operators am, m = 1, 2, ..., n, where the fermionic operators are represented by the famous Jordan-Wigner transformation

W2m-1 = < 0 <' , w2m = < 0 <' ' (9)

3.1. Fock space of operators

The fundamental concept for our analysis is a Fock space structure over the 4n-dimensional Liouville space of operators K, which density matrix p(t) is also a member of. From here on, we shall adopt Dirac bra-ket notation for the operator space K, which is fixed by the following definition of the inner product:

(x |y) = tr xjy, x, y e K. (10)

We note that 22n operator-products | Pa), labeled with a binary multi-index a,

Pai,02,...,«2n := 2-n/2w?1 wa22... wa2nn, aj e {0, 1}, (11)

constitute a complete orthonormal basis of K with respect to an inner product.

In fact it is easy to show that | Pa) is a fermionic Fock basis, and powers 1 in the product (11) can be considered like a sort of fermionic excitation, if we define the following set of linear annihilation maps Cj over3 K,

Cj | P«)= «j |WjPa), (12)

and derive the actions of their Hermitian adjoints—the creation linear maps Ct,

Cj | P„) = (1 -a )|w jP„), (13)

which satisfy canonical anticommutation relations

{Cj, ck } = 0, {Cj, Cj } = j, j, k = 1, 2,..., 2n. (14)

3.2. Bilinear form of the Liouvillean

The aim is now to show that the generator of the master equation (3)

L := —iad H + DD (15)

is in general a quadratic form in these adjoint fermionic maps C j, C j. In order to see that clearly, let us define the left and right multiplication maps over K:

WL|x) := |wjX), wwRR|x) := |xWj). (16)

Inspecting the actions of wLL, wRR on the Fock basis | Pa), one arrives at the following useful identities

wL = Cj + Cj, (17)

wR = P (Cj — Cj) = —(Cj — Cj )PP, (18)

3 We shall use notation where linear maps over the operator space (in physics literature sometimes referred to as 'super-operators') are designated by ~.

P := exp(inNN) and NN := E c)cj (19)

are a parity map and a number map, respectively, which count the parity and number of adjoint fermionic excitations (number of factors in (11)). Note that P anticommutes with all Cj, cj,

hence the second equality of (18), and P2 = i.

The unitary part of the Liouvillean (15) now trivially reads

-i ad Hs = —iwWL • HiWL + iHwR • wR = —4iPt • He . (20)

The dissipator (4) can be represented as a map over K as

2n n TO . .

D = E E Xvk df j (—t) (^»j + rß,ß(T)Lj,k) . (21)

ß,v j,k=1 </0

where f (t) is a (real-valued) propagator of Heisenberg dynamics in the closed system

Xß(t) = xß • exp(—i ad Hst)w =: f ß(t) • w, (22)

which—due to (20)—can be explicitly solved for a quadratic Hamiltonian (7), giving

fß(t) = exp(4iHt) x (23)

Lj,k|x> := |[wjx, wk]>, L^|x) := l[wk, xwj]> (24)

are fundamental basis dissipators which using (17) and (18) evaluate to

Lj ,k = wwL wR — wL wwL = (i + P)(c]ct — C¡ Cj) + (i — P )(CjCk — ckc)), (25)

,k = wwL wR — wR wR = (i + P )(C t Ct — Ct Cj) + (i — P )(C kC j — CkCt). (26)

It will prove useful if we express the internal dynamics (23) explicitly in terms of eigenvalues and eigenvectors of the Hamiltonian matrix H. Since a 2n x 2n matrix is antisymmetric and Hermitian, its real eigenvalues come in pairs 6m, —em, j = 1,..., n, with the corresponding eigenvectors um, u_m, namely Hum = emum and Hum = —'€mU^ since H* = —H. The eigenvectors may and should always be chosen orthonormal (even in the case of degeneracies), meaning

U - Um = 0, u • um = S¡,m. (27)

Then the spectral decomposition of the Heisenberg dynamics reads

fß(t) = E (e_4l€mt(£ß • UmU +e4lfmt(xß • Um)Um). (28)

Note that P± = (i ± P)/2 are orthogonal projectors that commute with all the terms (20), (25) and (26) that constitute the Liouvillean (15), [P±, L] = 0, and hence the dynamics (3) does not mix the operator subspaces K± = P±K composed of an even/odd number of fermionic operators. Since we are mainly interested in the expectation values of even observables, such as

currents and densities, we shall in the present paper focus on the dynamics in subspace K+ only, and consider the corresponding Liouvillean L|K+:

L+ = PP+ LPP+. (29)

The extension to the odd parity subspace is straightforward. Collecting the results (20), (21), (25) and (26), it is now obvious that L + is a bilinear form in ct and cj. For convenience, we define 4n Hermitian Majorana maps ar, r = 1, ... 4n,

a2j-i = (cj+c5), a2j = (cj- c5), (30)

and express the Liouvillean as

L+ = a- A a- A0 i, (31)

where the 4n x 4n complex antisymmetric matrix A, later referred to as a structure matrix, and a scalar A0, can be expressed as

A2j-1,2k-1 = —2iHjk - Mjk + Mk, j,

A2 j-1,2k = i Mk, j +i Mlk,

A2 j ,2k-1 = -i Mjk - i Ml j, (32)

A2j,2k = -2iHj k - Ml,k + Mkl j, A0 = tr M + tr Ml where M is a 2n x 2n bath matrix that can be compactly written as

M := £ Xv ® V (33)

dTOr)/M(-r). (34)

Defining the bath-spectral functions := ^ fZcdt flv(t)e—iwt for which the KMS

condition reads

f£,v(—= e^f (35)

and extending the range of integration in (34) to [ — to, to], or better still, neglecting the Cauchy principal value parts in the integrals—which exactly amounts to neglecting the Lamb-shift Hamiltonian term [4] in the master equation—we obtain a very simple expression (involving only finite sums) for the bath vectors:

^=* E Ef ) «*, ■ um +e4em ■ um )«m). (36)

, m= 1

At this point a remark on neglecting the Lamb-shift term is in order. As the Redfield model already involves a series of physical assumptions and approximations, it is somewhat difficult to argue under what conditions these terms can be dropped on the same level of approximations. However, one can straightforwardly show, using the KMS condition (6) and hermiticity (r, v(r))* = r^,(r), that only if the Cauchy principal value terms are dropped (i.e. if the range of integration in (4) is extended to [ — to, to]) the Redfield dissipator annihilates

the Gibbs state DD|e—3Hs > = 0, and hence the Gibbs state is the steady state of the equilibrium thermal Redfield model.

Note again that the inverse temperature in (36) could, in principle, be a function of the bath index /3 = /3v in case one would be interested in the non-equilibrium situation with couplings to several different temperatures. But we should stress that different temperatures only make sense among uncorrelated baths labelled by n, v, for which Tj// v = 0 for any /3.

We note also that the present formalism uniformly covers both the Redfield and the Lindblad master equations, as the Lindblad dissipator (5) is obtained from (4) by simply taking the limit T/v(t) = yu,v5(t + 0), and then the bath matrix reduces to a Hermitian form M = v u Yv,mXv ® Xu = Mf, which is equivalent to the one used in [11].

3.3. Static Liouvillean: normal modes, NESS and decay spectrum

Having the compact form of the Liouvillean (31)—and assuming for the time being that the structure matrix A is static, i.e. there is no explicit time dependence in the matrix H or coupling vectors Xu —we follow [11] and explicitly construct its normal form, the NESS, which is

exactly the right-vacuum state of (31) L + |NESS> = 0, the spectral gap, and the full spectrum of Liouvillean decay modes, all in terms of spectral decomposition of 4n x 4n matrix A. We state the main results here in a compact form.

Assuming the structure matrix is diagonalizable, its eigenvalues can be paired as /, —/, j = 1,..., 2n, assuming Re/ > 0, and its eigenvectors v2j-1 (corresponding to /) and v2j (corresponding to — /) can always be normalized—irrespective of possible degeneracies among / j , which shall be called rapidities

VVT = J, J := ax ® 12„ =

/0 1 0 0

-such that

0 ... \

Vr,s := vr,s. Thus the structure matrix

where V is a 4n x 4n matrix whose rth row is given by v allows the following decomposition,

A = VTdiag{/1, —/1,..., /2n, —/2n}JV, which after plugging into the Liouvillean (31) immediately brings it to a normal form

2X ß Â' bj

bj := V2j-i ■a,

b'j := v2j ■a,

are the normal-master-mode (NMM) maps, satisfying almost canonical anti-commutation relations

{bj, bk} = 0, {bj, b[} = Sj.k, {b;, b'k} = 0.

The map bj could be interpreted as an annihilation map and bj as a creation map of the j th NMM, but we should note that bj is in general not the Hermitian adjoint of bj. The right vacuum

is now essentially defined by bj |NESS) = 0, whereas the left vacuum is trivial (1|L + = 0 and satisfies (1 |bj = 0 . Note that the left vacuum simply corresponds to the identity operator 1 e K, i.e. an empty state 2n/2P0,0,...,0 with respect to a-fermions Cj.

Assuming that the whole rapidity spectrum is strictly away from the imaginary line Re ¡5j > 0, we state the following exact results:

1. |NESS) is unique.

2. Almost any initial density matrix relaxes to NESS with an exponential rate A = 2 min Re 5j (the spectral gap of the Liouvillean). The complete spectrum of 4n eigenvalues of L+ is obtained by all possible binary linear combinations = —2v • 5, Vj e {0, 1}.

3. The expectation value of any quadratic observable wj wk in a (unique) NESS can be explicitly computed as

(wj Wk)NESS := trwj WkPness = 2(1 |a2j—1 a2k—1 |NESS) (42)

= 2 ^^ v2m,2j — 1 v2m — 1,2k —1 (43)

=-- du G2j-i,2k —l (u),

n J-CO

G(o) := (A — i^1)—1 (45)

is a matrix of the non-equilibrium Green's function. The first equality is proven in [11]4 whereas the last equality requires a simple contour integration on the spectral decomposition of the resolvent (45).

4. The Wick theorem may be used for the calculation of expectation values of arbitrary higher order (even!) observables by the sums of all possible pairwise contractions of the form (42).

Note that as soon as some of the rapidities condense to the imaginary axis, or vanish, NESS typically becomes non-unique (see [13] for a detailed discussion of Liouvillean degeneracies).

3.4. Static Liouvillean: time-dependent correlation functions The complete Liouvillean propagator can be written explicitly as

exp(tL+) = ^ exp(—2tv • £)(£;)Vl ... (b2n)V2n|NESS)<1|(¿2n)V2n - - (Pi)Vl. (46)

V e|0,1|2n

It may be of some physical interest to evaluate dynamical response after perturbing the NESS by multiplying it with some local observable. In order to avoid discussion of negative parity dynamics , we take a pair of simplest even-order, quadratic observables, and define

4 Small simplification has been made with respect to the statement of theorem 3 of [11] which has been pointed out by Pizorn [12].

the corresponding non-equilibrium time-dependent correlation function—or non-equilibrium response function—as

,k ),(l ,m

)(t) :_ j(t)wk(t)wi(0)wm(0))NESS

= 4(1|fl2j-1 a2k-i exp(tC+)a1i- a2m-1 |NESS). (47)

Expressing the multiplication maps a2j-1 _ J]r=1 (V2r,2j-1 + V2r-12j-1 b'r) and plugging in the propagator (46), while noting that only the terms with 0 or 2 Liouvillean excitations contribute, we obtain a simple expression

(2n \ / 2n

^^ v2r,2j-1 v2r-1,2k-1 II ^^ v2r;,2l-1 v2r'-1,2m-1

4 e 2t(ßr +ßr'4V2r>,2j-1 V2r,2k-1 -V2r,2j-1 V2r',2k-l)

l^r <r;

(v2r'-1,2!-1 v2r-1,2m-1 - v2r-1,2/-1 v2r'-1,2m-l) . (48)

3.5. Time-dependent Liouvilleans

In this section, we indicate how to efficiently treat explicitly time-dependent master equations, written in third quantized form as

-|p(t)) = L+ (t )|p(t)), L + (t) = a> A(t )a_ — A0 (t)1, (49)

where explicit time-dependence of the structure matrix A(t) may physically arise due to driving by means of an external time-dependent force (time-dependent matrix H(t)) or time-dependent coupling operators (time-dependent vectors x^(t)). In this situation NESS cannot exist, but we shall show that one may still efficiently evaluate the propagator

|p(t)> = W|p(0)>, U := Texp / dr£+(r) , (50)

where T indicates a time-ordered product.

The procedure is the following. Note that the space of all anti-symmetric complex structure matrices forms a Lie algebra so(4n, C). The following straightforward identity

[2a • Aa, 2a • Ba] = 2a • [A, B]a, (51)

holding for any pair of complex 4n x 4n matrices A, B, indicates that Liouvilleans (31) and (49) generate a 4n-dimensional representation of so(4n, C). Thus, the time-ordered product (50) can be evaluated within a Lie group SO (4n, C) of 4n x 4n matrices,

U = Texp jf drA(r)^ , (52)

and5 then the full Liouvillean propagator is written as

U = exp(a_- Ca- Co 1), C = 1ln U, Co = j dr Ao (r). (53)

5 Even if this has to be done numerically, using Trotter-Suzuki decomposition schemes, the computational complexity is only polynomial in n.

11 IOP Institute of Physics <J>deutsche physikalische Gesellschaft

The logarithm of U can now be considered as a 'static' Liouvillean, so we can diagonalize it by the methods of section 3.3, leading to spectral decomposition of the form (46).

4. XY spin chains

The theory of the previous sections shall now be applied to investigate a homogeneous,

finite XY chain of n spins, described by Pauli matrices aj, y ,z, j = 1,... n, with the Hamiltonian

n — I / * * \ n

H = E -X J + ^--J+ E h-J, (54)

which is described by two real parameters, anisotropy y and transverse magnetic field h. Without loss of generality we may assume that y, h e [0, to). We decide to couple the XY chain thermally only at its ends, so we consider the most general four coupling operators that allow for an explicit solution,

X1 = k1 (aj cos + a1y sin ), X3 = k3 (aN cos 03 + aN sin 03),

y y (55)

X2 = k2 (aX cos 02 + a{ sin 02), X4 = k4(a^ cos 04 + aN sin 04),

and fully decorrelated baths T^ v = We take standard baths of harmonic oscillators at

two ends with possibly different inverse temperatures, and Ohmic spectral functions

TA*» = X -, 01,2 = ^L, ft,4 = ^R. (56)

exp(^ßj — I'

Note that the frequency cutoff in the spectral function is irrelevant as we neglect the Lamb-shift term in the master equation.

The entire problem can be fermionized by means of Jordan-Wigner transformation (9), namely the Hamiltonian and the coupling operators transform to

H = —iW2j W2j+I--^^ W2j — I W2j+2) — i^ h W2j — I W2j ,

I — Y I+Y

~Y~ W2 J w2 j+I — ~Y~

J=I j=I

XI = KI (wI cos 0I + w2 sin 0I), X3 = Wk3(w2n cos 03 — w2n—I sin 03),

X2 = K2(wI cos 02 + w2 sin 02), X4 = Wk4(w2n cos 04 — w2n_I sin 04),

where W = (—i)n I wI w2 ... w2n is an operator that commutes with all the elements of K+ (or anti-commutes with all the elements of K—) and satisfies WWf = Wf W = I, hence it has no effect on the dissipator (4) in L+. We note, however, that the commutation of W through p in (4) for the dynamics in K— produces a minus sign in all the bath terms, i.e. it changes the sign of DD, with respect to a pure fermionic problem.

The 4n x 4n structure matrix now has a specific block-tridiagonal + block-bordered form,

A = A7 + B, (58)

/a b 0 cab Oca

0 0 ...

\0 0 ...

in—I

In—I 0

rn I rn

where a, b, c are 4 x 4 matrices:

a = — ih 12

b = 212 ® (i-y — Y-x), c = —bT.

The sequences of 4 x 4 matrices lj, lj, rj, rj, which form the block-bordered part B, can be straightforwardly computed (see (32)) from the form of the coupling vectors x 1 2 = (k12 cos 01,2, k12 sin 01,2, 0, ... , 0)T, x3 4 = (0, ..., 0, —k3,4 sin 03,4, k3,4 cos 03,4)T, and their bath-transformations (36) with (56). Although we are unable to give closed-form general expressions, we can make an asymptotic estimate—for large n—on the decay of these matrices with their distance from the diagonal

Ill J ||~||1J ll~l|rn+I.

n+I—J ||a exp(—KJ).

The coefficient K > 0 in general depends only on y, h and ft (for lj) or ft (for r j). Note that for the special case of local Lindblad coupling (5) with the same local coupling operators (57), the only non-vanishing blocks that remain are the diagonal ones l1 and rn, given explicitly in [11].

Below we shall present some intriguing numerical results of the non-equilibrium thermal Redfield equation (3) and (4) for the open XY chain given by (54)-(56), in comparison with the local non-equilibrium Lindblad model (5) where a suitable set of coupling operators of the form (55) and 4 x 4 coupling matrix yM,v can be chosen to parametrize the Lindblad operators

L12 = y rL2af, L3 4 = y rR2anT, parametrized in exactly the same way as in [11, 14]. For all the numerical results reported for the thermal Redfield model, we consider the bath parameter values k1 = k3 = 1, k2 = k4 = 0, = 03 = n/6, and ft = 0.3, ft = 5.2 unless ^s are varying, and X = 0.1 unless X is varying, whereas for the Lindblad model we always take the bath parameters TL = 0.5, TL = 0.3, TR = 0.5, TR = 0.1.

4.1. Non-equilibrium phase transition

In [14] an intriguing suggestion of a quantum phase transition far from equilibrium in the steady state of an open boundary driven XY spin chain has been put forward. Numerical and heuristic theoretical evidence has been given for the spontaneous emergence of long-range magnetic order in NESS as soon as the magnetic field drops below the critical value |h| < hc,

h c = |1 — y 21. (62)

However, that study was done with local Lindblad reservoirs, so the question remained whether the effect persists in the presence of local thermal reservoirs satisfying KMS conditions for non-vanishing temperatures. It is an easy task now to follow the recipes of section 3.3 and numerically evaluate the spin-spin correlator (note the use of the Wick theorem as the spin-spin correlator is of fourth order in w j):

Cl m = tr (azamPness) — tr (azPness) tr (amPness)

= <W2/- 1 w2m — 1) NESS < W2l w2m ) NESS < W2l—1 w2m ) NESS < W2l w2m — 1) NESS.

Figure 1. Correlation matrices Cl,m, l horizontal axis (left to right), m vertical axis (bottom to top), of the non-equilibrium thermal Redfield model of an open XY chain for y = 0.5 and different field strength h indicated on the figures (note that hc = 0.75) and two different system sizes n (indicated). Bath parameters are specified in the text.

First, we use efficient prescription (43) to compute correlation matrices at non-equilibrium conditions pL = 0.3 = = 5.2 and plot them for two different system sizes and five different values of h around hc in figure 1. The results look qualitatively identical to those for Lindblad driving, even for other quantities that were investigated numerically in detail in [14].

For example, in figure 2, we plot the phase diagram of the residual correlator Cres = Ellmm'>n/2 IC,mLm'>n/2 1, which also reveals possible criticality in the region of a large anisotropy y > 1 previously not discussed. We note that the size dependence of the residual magnetic correlator Cres shows very characteristic behavior: namely

Cres a exp(-nn) with n > 0 for |h| > hc or h = 0, (64)

Cres a 1/n for 0 < |h| < hc. (65)

Thus we shall refer to the regime with 0 < |h | < hc as long-range magnetic correlation (LRMC) phase6, the regime with |h | > hc, or h = 0, as non-LRMC phase, and the regime with |h | = hc as critical. Scaling (64) and (65) is illustrated in figure 3. Exponential decay of the Cres (n) in non-LRMC phase (64) is consistent with the exponential decay of a two-point correlator with the distance between sites C(r) = j-i=r Ci,j/J2j-i=r 1 ~ exp(-§r), as can be qualitatively noted already in figure 1. However, we demonstrate in figure 4 that the exponent § could, in principle, be very different between the Redfield and local Lindblad models. Furthermore, as for the Lindblad model the exponents § and n (of (64)) appear to be equal, for the Redfield

6 Note, interestingly, that unlike for the local Lindblad driving [14] the XX line y = 0, 0 < |h | < 1, also exhibits long-range magnetic correlations for the thermal Redfield driving.

Figure 2. Phase diagram for the non-equilibrium thermal Redfield model of an open XY chain. We plot the residual correlator Cres against the bulk parameters y, h. The dashed curve indicates the critical line hc(y) (62). The system size is fixed to n = 100 and bath parameters are specified in the text.

-10"13 o

500 1000

Figure 3. Residual correlator Cres as a function of system size n for the LRMC phase (y = 0.5, h = 0.2, left plot) and the non-LRMC phase (y = 0.5, h = 0.9, right plot), where we compare the non-equilibrium thermal Redfield model (red squares) and the non-equilibrium Lindblad model (blue circles) with bath parameters as specified in the text. The thin lines indicate the suggested behavior 1/n (on the left) and exp(-nn) (on the right) (with the numerical best fit n = 1.192 for the Redfield model and n = 0.937 for the Lindblad model).

model they do not seem to be simply related. Analytical estimation of these exponents presents a challenge for future theoretical work.

However, we note that with the thermal driving with Redfield dissipators, the long-range magnetic order disappears when the temperatures of the baths become equal, = and there we recover, consistently, all the properties of the thermal state [15] that are most easily numerically reproduced by the method of [16], i.e. fast decay of correlations for any h and absence of long-range order. For example, it is interesting to note how the residual correlator Cres

Figure 4. Comparison of the decay of the two-point spin-spin correlator C(r) = —i=r Ci j/Y]j-i =r 1 ~ exp(-§r) between the non-equilibrium thermal Red-field model (red squares) and the non-equilibrium Lindblad model (blue circles) for the same values of bulk parameters in the non-LRMC phase (h = 1.05, y = 0.2, n = 200) and bath parameters specified in the text. The thin lines indicate suggested exponential decays a exp(-§r) with the exponents § = 1.635 (fitting the Redfield model) and § = 0.937 (fitting the Lindblad model).

iogio(A®

Figure 5. Residual correlation Cres versus (inverse) temperature drop Ay0 between the left and the right bath, ft = 2 - A^/2, ft = 2 + A^/2, for the non-equilibrium thermal Redfield model of an open XY chain in LRMC phase (Y = 0.5, h = 0.3), system size n = 100, and the bath parameters specified in the text. The thin line indicates suggested | A012 behavior.

(for large n in the LRMC phase) decreases as a function of the difference of inverse temperatures A;0 = ft — ft, namely the numerics of figure 5 suggests clearly that Cres a (A^)2.

The heuristic explanation of this non-equilibrium phase transition is rather straightforward [14], however, its exact proof and also the quantitative dependence of the decay exponent n(Y, h) are still lacking. We note that the transition point h = hc is characterized by a simple

10 20 50 100 200 500 1000

Figure 6. Liouvillean spectral gap A for the non-equilibrium thermal Redfield model of an open XY chain. We plot three different cases with: y = 0.5, h = 0.8 > hc (non-LRMC phase, light blue circles), y = 0.5, h = 0.75 = hc (critical regime, dark blue squares), y = 0.5, h = 0.3 < hc (LRMC phase, black diamonds), whereas bath parameters are specified in the text. Suggested power law decays n-3 and n-5 are indicated with thin lines.

property of the XY spin chain quasi-particle dispersion relation

) = y^(cos q — h)2 + y2 sin2 q, (66)

where ej = ^(2nj/n) would be exactly the (positive) eigenvalues of matrix H if periodic boundary conditions were imposed on the closed system. Namely, in non-LRMC phase |h | > hc there exists only a single pair of trivial stationary points q* = 0, n, whereas in LRMC phase |h | < hc there exists another pair of non-trivial stationary points ±q* = 0, n, d^/dq |q =q* = 0, which introduces a new non-trivial length scale 1 /q* that determines typical sizes of correlated regions in the matrix C/,m (see figure 1). Therefore, this simple non-equilibrium quasi-particle picture predicts mean-field critical exponent 1 /q* ~ |hc — h|—1/2 as h f hc (confirmed in [14]).

The non-equilibrium quantum phase transition can also be characterized by the scaling of the Liouvillean spectral gap A(n), namely in the critical regime one expects a qualitative increase in the relaxation time 1/A to NESS. Numerical results (see figure 6) suggest that the spectral gap of the Liouvillean remains like in the local Lindblad case [11],

A a n—3 for h = hc, A a n—5 for h = hc, (67)

although we are at the moment unable to prove this conjecture. Also note the slight fluctuations of A(n) in the LRMC phase as opposed to a smooth power law in the non-LRMC phase.

Long-range correlations for |h | < hc naturally imply sensitivity of NESS to tiny variations in the system's parameters. For example, one may also expect that local observables in NESS will then be sensitive functions of the bath-driving or even bulk parameters, such as the magnetic field h. In figure 7, we plot local magnetization in the center of the chain sz = (a^)NESS versus field strength h. Indeed, we notice that for |h| > hc, sz(h) is a smooth function whereas for |h| < hc, sz(h) becomes a rapidly oscillating or, better still, fluctuating function. Even though the amplitude of these oscillations decreases with n, the scale of h on which sz (h) fluctuates decreases with n even faster, so we predict that in the thermodynamic limit n ^ to, in the

Figure 7. Hypersensitivity of NESS to magnetic field strength h. We plot local magnetization sz(h) = (0^/2)NESS for the non-equilibrium thermal Redfield model of an open XY chain with y = 0.5 and bath parameters as written in the text. Big blue (small red) circles represent data for n = 50 (n = 100), whereas the vertical line denotes the critical value h = hc.

LRMC phase the local susceptibility dsz/dh would be ill defined. In summary, the LRMC phase can be characterized by hypersensitivity of NESS to external parameters.

4.2. Heat transport and entropies

An important non-equilibrium physical effect that one can investigate more deeply in an open XY chain is heat transport, which has recently been intensively studied in quantum spin chains, see e.g. [17]-[21] or [22] for a recent review on the topic.

Writing the Hamiltonian (54) in the bulk as a sum H = m Hm with a two-body energy density operator

tj _ -1+ Y , -1 - Y hm .hm+1 _

Hm = —i—2— W2m w2m+1 + i-2-W2m-1 W2m +2 — W2m-1 W2m — i 2 W2m+1 W2m+2, (68)

one can derive the local energy current Qm = i[ Hm, Hm+1]

= i(1 — Y )(W2m — 1 W2m+3 + W2mW2m +4) — 2ih(1 — Y)(W2m — 1 W2m+1 + W2m +2W2m +4)

— 2ih (1 + Y)(W2m W2m +2 + W2m+1 W2m+3 ), (69)

which, by construction, satisfies the continuity equation

(d/dt)( Hm ) = (i[ H, Hm ]) + tr Hm D p(t) = —( Qm ) + ( Qm — 1). (70)

The two terms between the two equality signs above correspond to the unitary and dissipative term in the master equation (3). The unitary term has already been transformed to a simple expectation value using cyclicity of the trace tr x[y, z] = tr y [z, x], while the dissipative term can be further shown to vanish in the bulk 2 ^ m ^ n — 2 by exercising the cyclicity of the trace again and transforming the integrand of (4) to terms of the form trXM(—t)p[Xv, Hm] = 0. The rhs expression of equation (70) then follows from the nearest-neighbor locality of the

3 -1 0 1 2

fog10 (Pr)

Figure 8. NESS expectation value of heat current (Qm>ness versus two inverse temperatures and f3R for the non-equilibrium thermal Redfield model of an open XY chain with y = 0.5, h = 0.9, system size n = 53, and bath parameters given in the text. Note that the 'shoulders' of maxima, around ¡3L ^ 0.05, > 1, and with L and R exchanged, could be interpreted as negative differential heat conductance.

Hamiltonian. Therefore, in NESS the expectation value of the current (Qm>NESS should be independent of the position m. By looking at the dependence of the steady-state current on system size we clearly find ballistic transport, namely (Qm>NESS = O(n0), irrespective of the temperature differences between the baths and bulk parameters of the model (i.e. whether in the LRMC phase, non-LRMC phase, or critical). However, we find a very interesting dependence of heat current on temperature driving, i.e. on the two temperatures of the thermal baths. In figure 8, we plot (Qm >NESS versus ¡3L and and find a maximum of the current for intermediate driving, namely when one of the temperatures is less than 1/^L < 1 and the other temperature is about 1/^R ^ 20. This is a clear signature of negative differential heat conductance, which could perhaps be related to similar far-from-equilibrium effects recently observed in spin and charge transport [23].

This behavior can be nicely characterized by computing the Gibbs entropy of NESS. Since NESS is completely characterized by quadratic correlations (wj wk>NESS and the Wick theorem, one can adopt the recipe that has been proposed in [24] for computing block entropies (or entanglement entropies) applied to the entire lattice. In fact, taking an arbitrary block region A c {1,..., n}, one can compute Von Neumann entropy SA(p) = -trApA log2pA (in base 2), where pA = trap is a reduced density matrix and A denotes the complement of A, as

SA = ^ H2((1 + vj)/2) with H2(x) :=-x log2 x - (1 - x) log2 x, (71)

and ±ivj are the eigenvalues of the 2#(A) x 2#(A) part of the correlation matrix Bjk defined by (wjwk>NESS =: 8jk + iBjk, restricted to Majorana operators wj, wk corresponding to spins from block A. The same general procedure has been applied to thermal (Gibbs) states in [16]. When taking the maximal block A = {1,..., n}, we obtain exactly the standard Gibbs entropy

-3-2-10 1 2

log10 iAr)

Figure 9. Gibbs entropy of NESS versus two inverse temperatures ¡L and ft for the same parameters as in figure 8. Note that in the red region (of both large inverse temperatures), NESS is no longer a density matrix (at least one of the eigenvalues becomes slightly negative); hence the Gibbs entropy is strictly no longer defined there.

of NESS. In figure 9, we plot the Gibbs entropy 5{i,...,„} as a function of two bath temperatures and show that, quite remarkably, the regions of large (maximal) heat current correspond to regions of large (locally maximal) Gibbs entropy. This is not unexpected as the product of the heat current and the inverse temperature difference A5 may be understood as the entropy production rate.

Calculation of Gibbs entropy of NESS also provides a nice way of controlling the positivity of NESS as a density matrix, since this is by no means guaranteed by the Redfield master equation. Indeed we find that for very small temperatures (large ¡5s), or for very strong bath coupling X, the positivity of NESS might be slightly violated (red region in figure 9), namely some of the correlation matrix eigenvalues v j become slightly larger than 1 (but in our numerical experience never by more than 10-7 or so).

We can use the concept of block entropy of NESS to further characterize the non-equilibrium phase transition. For example, we may compute the total (quantum plus classical) correlations between two halves of the spin chain in NESS as given by quantum mutual information (QMI) I (n) = S{1,...,n/2} + S{n/2+1,...,n} - S{1,...,n}.

Interestingly, we find (see figure 10) that QMI saturates I (n) = O(n0) in the non-LRMC phase (for |h| > hc), whereas in the LRMC phase (for 0 < |h| < hc) QMI becomes extensive I(n) = O(n), indicating a drastic enhancement of correlations in NESS. This is again very similar to the behavior of operator space entanglement entropy (OSEE) (analyzed for the Lindblad model in [14]), so one may extend the relationship between QMI and OSEE which has been conjectured for thermal states in [16] to NESS.

In the context of energy transport, it is interesting to look at the energy density profiles in NESS. In figure 11, we plot the relative spatial fluctuation of the energy density f (m) = kHm)ness - H\/\H|, where H = (n - 3)-1 ¿=2 (Hm)ness is the averaged energy density. Quite strikingly, we observe a big variation of f (m) from site to site for the LRMC phase and

S 3 2 1 0

400 600 n

Figure 10. Another manifestation of the non-equilibrium phase transition: QMI of NESS for the non-equilibrium thermal Redfield model of an open XY spin chain. The bulk parameters are y = 0.5 and h = 0.9 > hc = 0.75 (lightest blue, saturated curve), h = 0.7, h = 0.5 and h = 0.3 (from lighter to darker blue curves). Thin red lines indicate the linear growth of QMI for h < hc.

0.001 -

Figure 11. Another manifestation of non-equilibrium phase transition: positional fluctuations in energy density in NESS of the non-equilibrium thermal Redfield model of an openXY chain. We plot the relative fluctuation f (m) = |(Hm)NESS — H/|/|H/|, where H is the bulk average of energy density (Hm)NESS. Three curves correspond to y = 0.5 and h = 0.7 < hc (black curve), h = 0.75 = hc (dark blue curve) and h = 0.8 > hc (light blue curve), while the system size is n = 253.

very smooth (non fluctuating) behavior for the non-LRMC phase, which is characterized with a bulk-constant f (m) that is exponentially small in n. This behavior can again be considered as a manifestation of hypersensitivity of NESS and LRMC.

Finally, we check the dependence of the heat current (Qm )NESS on the system-bath coupling strength A. It was recently reported by Karevski and Platini [25] that the spin current Jm in the local Lindblad model of an open isotropic XX chain y = 0 has non-monotonic dependence on A which can be universally described by a formula (Jm)NESS = a'A2/(b' + A4), where a', b' are

0.5 0.4 0.3

0.2 0.1 0.0

0.0 0.5 1.0 1.5 2.0

Figure 12. NESS expectation value of the heat current (Qm>NEss versus the coupling strength X for the non-equilibrium thermal Redfield model of an open XY chain with h = 0.5 < hc (black curve), h = 0.75 = hc (dark blue curve) and h = 1.0 > hc (light blue curve), system size n = 200 and other bath parameters as given in the text. Note that the full line gives numerically excellent fit to the Karevski and Platini formula [25] (Qm > = aX2/(b + X4) for best fitted parameters a = 0.040, b = 0.0070, a = 0.066, b = 0.0076, and a = 0.088, b = 0.0071 for the three cases of h = 0.5, 0.75, 1.0, respectively.

some constants. For the anisotropic XY model and general non-equilibrium thermal Redfield driving, we are unable to derive an exact analytic result; however, our numerical simulations suggest very similar behavior for the heat current

(Qm>ness ~ aX2/(b + X4), (72)

where a and b are again some constants that may depend on all system parameters except X. This is particularly interesting because in the anisotropic XY model the spin current is not even well defined as there is no corresponding conservation law. This behavior is demonstrated in figure 12, where one may also notice small but detectable deviations between numerics and the best fit to (72). We note that the error of the fit does not decrease but is roughly constant when we increase the system size n.

5. Discussion

The purpose of the present paper was three-fold. Firstly, we have outlined a general method for the exact treatment of quadratic many-body Markovian master equations. Our formalism, which rests upon treating density operators as elements of a suitable operator Fock space (or Liouville-Fock space), is quite flexible and allows for explicit solution of static and time-dependent quantum many-body Liouvillean problems, for example, computation of arbitrary physical observables in the non-equilibrium steady state, decay rates of approach to the steady state, or even time evolution of the density matrix of externally forced systems described by explicitly time-dependent Liouvilleans, all with polynomial computation complexity in number of particles (fermionic degrees of freedom).

Secondly, we have analyzed in detail the Redfield model of thermal baths within our framework. In spite of the fact that the Redfield model does not define a proper dynamical semigroup,

namely it is not guaranteed to preserve positivity of the density operator, we have confirmed that steady states typically correspond to proper (positive) density operators. Tiny deviations from positivity have only been observed in some test cases for very small temperatures or very large couplings to the baths (which anyway violate the weak coupling assumption). Furthermore, we have shown that coupling the central system with several thermal baths of the Redfield type at different temperatures produces physically interesting NESSs, for example, such states that carry non-vanishing heat current. We wish to stress this physically obvious but mathematically delicate point with particular care, as we have found qualitatively different results for Lindblad-Davies dissipators which generate proper dynamical semigroups and satisfy the detailed balance condition with respect to Gibbs states [5, 26]. Namely, when we constructed a Lindblad-Davies dissipator with respect to two baths with two different temperatures coupled to two ends of the system (spin chain), we found that the resulting steady state (fixed point of Liouvillean dynamics) is simply some convex combination of two Gibbs states corresponding to the bath temperatures, and as such has always zero heat current and cannot represent the physical steady state. This implies that the secular approximation (sometimes called the rotating wave approximation), which is the one-step from the Redfield to the Lindblad-Davies bath model, prohibits the emergence of the physical out-of-equilibrium steady states with currents; therefore, the seemingly harmless rapidly oscillating terms in the Redfield dissipator may be absolutely essential for non-equilibrium physics. Thus we conjecture that the thermal Redfield model is somehow a minimal mathematical model that can describe non-equilibrium thermal driving of a non-self-thermalizing (e.g. integrable) open quantum system.

Thirdly, we have applied our theory to analyze non-equilibrium quantum phase transition and heat transport in an open XY spin 1/2 chain. We have carefully compared numerical results for the non-equilibrium thermal Redfield model and the local Lindblad model, which has been discussed before [11, 14]. We have found that the phase diagram of the non-equilibrium XY model is insensitive to the theory with which we describe the baths, and the differences were only quantitative. In particular, we wish to stress that thermally driven heat current in the XY chains exhibits non-monotonic dependence on the temperature difference, which may be interpreted as negative differential heat conductance. We believe that our numerical results on a non-equilibrium open XY chain provide a strong motivation for further analytical work. In particular, we believe that the block-tridiagonal plus block-bordered structure of the Liouvillean structure matrix (58) and (59) could be explored in combination with the non-equilibrium Green function formula for the observables (44) and (45) to yield explicit asymptotic results for large n.

Note added: formally quite a similar approach to non-equilibrium quasi-particles has recently been developed independently by Kosov [27].

Acknowledgments

We acknowledge financial support by Programme P1-0044 and Grant J1-2208 of the Slovenian Research Agency (ARRS).

References

[1] Araki H and Barouch E 1983 J. Stat. Phys. 31 327

Araki H 1984 Publ. RIMS Kyoto Univ. 20 277

[2] Ruelle D 2000 J. Stat. Phys. 98 57

[3] Jaksic V and Pillet C-A 2002 J. Stat. Phys. 108 787

Jaksic V and Pillet C-A 2002 Commun. Math. Phys. 226 131

Aschbacher W, Jaksic V, Pautrat Y and Pillet C-A 2006 Introduction to non-equilibrium quantum statistical mechanics Open Quantum Systems III. Recent Developments (Lecture Notes in Mathematics vol 1882) pp 1-66

[4] Breuer H-P and Petruccione F 2002 The Theory of Open Quantum Systems (Oxford: Oxford University Press)

[5] Alicki R and Lendi K 2007 Quantum Dynamical Semigroups and Applications (Berlin: Springer)

[6] Wichterich H, Herich M J, Breuer H P, Gemmer J and Michel M 2007 Phys. Rev. E 76 031115

[7] Sinaysky I, Petruccione F and Burgarth D 2008 Phys. Rev. A 78 062301

[8] Ban M 2009 J. Mod. Opt. 56 577

[9] White S R 1992 Phys. Rev. Lett. 69 2863 Schollwock U 2005 Rev. Mod. Phys. 77 259 Vidal G 2003 Phys. Rev. Lett. 91 147902 Vidal G 2004 Phys. Rev. Lett. 93 040502

[10] Prosen T and Znidaric M 2009 J. Stat. Mech. P02035

[11] Prosen T 2008 New J. Phys. 10 040326

[12] PiZorn 12009 Entanglement and quantum critical phenomena: operator spaces and random matrix theory PhD

Thesis University of Ljubljana

[13] Prosen T and PiZorn I 2009 to be published

[14] Prosen T and PiZorn I 2008 Phys. Rev. Lett. 101 105701

[15] Barouch E and McCoy B M 1971 Phys. Rev. A 3 786 Barouch E and McCoy B M 1971 Phys. Rev. A 3 2140

[16] Znidaric M, Prosen T and PiZorn I 2008 Phys. Rev. A 78 022103

[17] Saito K, Takesue S and Miyashita S 2000 Phys. Rev. E 61 2397 Saito K 2003 Europhys. Lett. 61 34

[18] Heidrich-Meisner F, Honecker A, Cabra D C and Brenig W 2002 Phys. Rev. B 66 140406 Heidrich-Meisner F, Honecker A and Brenig W 2005 Phys. Rev. B 71 184415

[19] Mejia-Monasterio C, Prosen T and Casati G 2005 Europhys. Lett. 72 520

[20] Michel M, Hartmann M, Gemmer J and Mahler G 2003 Eur. Phys. J. B 34 325 Michel M, Mahler G and Gemmer J 2005 Phys. Rev. Lett. 95 180602 Gemmer J, Steinigeweg R and Michel M 2006 Phys. Rev. B 73 104302 Steinigeweg R, Ogiewa M and Gemmer J 2009 Europhys. Lett. 87 10002

[21] Arrachea L, Lozano G S and Aligia A A 2009 Phys. Rev. B 80 014425

[22] Dhar A 2008 Adv. Phys. 57 457

[23] Benenti G, Casati G, Prosen T and Rossini D 2009 Europhys. Lett. 85 37001 Benenti G, Casati G, Prosen T, Rossini D and Znidaric M 2009 Phys. Rev. B 80 035110

[24] Vidal G, Latorre J I, Rico E and Kitaev A 2003 Phys. Rev. Lett. 90 227902 Latorre J I, Rico E and Vidal G 2004 Quantum Inform. Comput. 4 48

[25] Karevski D and Platini T 2009 Phys. Rev. Lett. 102 207207

[26] Davies E B 1974 Comm. Math. Phys. 39 91 Alicki R 1976 Rep. Math. Phys. 10 249

Kossakovski A, Frigerio A, Gorini V and Verri M 1977 Comm. Math. Phys. 57 97 Spohn H 1977 Lett. Math. Phys. 2 33 Frigerio A 1977 Lett. Math. Phys. 2 79

[27] Kosov D S 2009 J. Chem. Phys. 131 171102