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

The Lindblad master equation for an arbitrary quadratic system of n fermions is solved explicitly in terms of diagonalization of a 4n x 4n matrix, provided that all Lindblad bath operators are linear in the fermionic variables. The method is applied to the explicit construction of non-equilibrium steady states and the calculation of asymptotic relaxation rates in the far from equilibrium problem of heat and spin transport in a nearest neighbor Heisenberg XY spin 1/2 chain in a transverse magnetic field.


Introduction
Understanding time evolution of an open quantum system of many interacting particles is of primary importance in fundamental problems of quantum physics, such as decoherence [1,2] and closely related quantum measurement problem [3,4], quantum computation [5,6], or the problem of computation of non-equilibrium steady states (NESS) in quantum statistical mechanics [7,8,9]. Even though application of the methods of Hamiltonian dynamical systems and ergodic theory to quantum systems out of equilibrium gives many promising results [10,11,12], the field of open quantum systems is still lacking non-trivial explicitly solvable models, as compared to studies of 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 master equations for open quantum systems are limited to quite restricted models of a single particle, 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 general quadratic system of n interacting fermions which are coupled to a general set of Markovian baths, specified in terms of Lindblad operators which are linear in the fermionic variables -the so called quantum Liouville super-operator (or Liouvillean)can be explicitly diagonalized in terms of 2n normal master modes, i.e. anticommuting super-operators which act on the Fock space of density operators. This can be understood 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 the master 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 real axis. (ii) The covariance matrix of NESS can be computed explicitly in terms of the eigenvectors of 4n × 4n antisymmetric complex matrix. It can be used to completely express 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 of heat and spin transport far from equilibrium in nearest neighbor Heisenberg XY spin 1/2 chains subject to a transverse magnetic field. As a result we reproduce ballistic transport in the integrable spatially homogeneous case (see e.g. [19,20,21,22,23,24,12] for related recent studies of quantum thermal conductivity in one dimension), and predict ideally insulating behaviour (at all temperatures) in a disordered case of spatially random interactions/field. Apart from obtaining numerical results which go by far beyond what was so far accessible by direct numerical solution of the many-particle Lindblad equation, either directly or by means of quantum trajectories [15], we also obtain two notable analytical results in the spatially homogeneous (non-disordered) case: (i) We compute the spectral gap of the Liouvillean i.e. the rate of of relaxation to the NESS and 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 explain quantitatively the exponential falloff of energy density or temperature profiles near the bath contacts.
The paper is organized as follows. In section 2 we shall outline a general method for the diagonalization of the Liouvillean super-operator for finite quadratic open Fermi systems and an explicit construction of NESS. In section 3 we illustrate the method by working out a simple example of a single fermion or a two level quantum system in a bath. In section 4 we demonstrate the usefulness of the new method by applying it to quantum transport in XY spin chains. In section 5 we discuss possible alternative applications and generalizations of the method and reach some conclusions.

General method of solution for the Lindblad equation
The general master equation governing time evolution of the density matrix ρ(t) of an open quantum system, preserving trace and positivity of ρ, can be written in the Lindblad form [25,17] as (we set = 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 possibly different values of thermodynamic potentials). We are now going to describe a general method of explicit solution of (1) for a quadratic system of n fermions (or spins 1/2) with linear bath operators where w j , j = 1, 2, . . . , 2n, are abstract Hermitian Majorana operators satisfying the anti-commutation relations and generate a Clifford algebra. Thus, 2n × 2n matrix H can be chosen to be antisymmetric H T = −H. Throughout this paper x = (x 1 , x 2 , . . .) T will designate a vector (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, 2, . . . , n, or (ii) spins 1/2 with canonical Pauli operators σ m , m = 1, 2, . . . , n, Here we are not concerned with physical criteria for the validity of the so-called Markovian approximation under which eq. (1) is derived, so we shall make no assumptions on the smallness of the bath coupling constants l µ,j . We merely consider the Lindblad equation (1) as a possible parametrization of an important subset of Markovian completely positive quantum channels and demonstrate its complete solvability for quadratic systems. Note that generalization of our formalism to explicitly time dependent Hamiltonians H(t) and Lindblad operators L µ (t), generating more general and possibly non-Markovian open system dynamics, is straightforward. See e.g. [26] for a discussion of Markovianity.

Fock space of operators
We begin by associating a Hilbert space structure x → |x to a linear 2 2n = 4 n dimensional space K of operators, with a canonical basis |P α with orthonormal with respect to an inner product The form of the canonical basis of the operator space K suggests that it is just a usual Fock space with an unusual physical interpretation. Namely we can define the following set of 2n adjoint annihilation linear mapsĉ j over K c j |P α = δ α j ,1 |w j P α (9) and derive the actions of their Hermitian adjoints -the creation linear mapsĉ † , Straightforward inspection then shows that they satisfy the canonical anticommutation relations The key is now to realize that the quantum Liouville mapL defined by eqs. (1,2,3) can be written as a quadratic form in adjoint Fermi mapsĉ j ,ĉ † j (or for short, a-fermions). ‡ First, we consider the unitary part of Liouvillean Since K is a Lie algebra, one defines the adjoint representation of a Lie derivative for an arbitrary A ∈ K back on K as, ad A|B := |[A, B] . It is now straightforward to ‡ Throughout this paper Dirac's bra-ket notation shall be used only for a Hilbert space K of physical operators, including density operators, in a sense of GNS construction although here all spaces will be finite dimensional. Symbols with a hat shall designate linear maps over the operator space K. For instance, we note a key distinction between physical fermions c m (5) and a-fermionsĉ j (9).
compute the action of a Lie derivative of a product of two Majorana operators on an arbitrary basis element Extending this relation by linearity to an arbitrary element of K, it follows that for an arbitrary quadratic Hamiltonian (2) its Lie derivative has a very similar quadratic form in a-Fermi mapŝ It is worth stressing here that for an arbitrary (complex) matrix H,L 0 (14) conserves the total number of a-fermionsN := jĉ † jĉ j =ĉ † ·ĉ, namely [L 0 ,N ] = 0.
Second, we consider the action of the Lindblad mapŝ where we writeL j,k ρ := 2w j ρw k − w k w j ρ − ρw k w j . Again we proceed by computing the actions ofL j,k on elements of the canonical basis of operator Fock space K. In order to do so, it is crucial to observe that the question whether w j commutes or anticommutes with P α depends on the number of a-fermions |α| := 2n k=1 α k in |P α , namely P α w j = (−1) |α|+α j w j P α , and hencê Observing that one derives from (16) the general expression forL j,k Obviously, the mapsL j,k , and hence also the total Lindblad part of Liouvillean µL µ , do not conserve the number of a-fermions. But they conserve its parity i.e. the product of 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 a direct sum K = K + ⊕ K − , and even/odd operator spaces are orthogonally projected as K ± = (1 ± exp(iπN ))K. Thus the positive parity subspace K + is a linear space spanned by |P α with even |α|. All the mapsL j,k now act separately on K ± ,L j,k K ± ⊆ K ± . For example, the maps defined on even parity subspace are indeed quadratic in a-fermionŝ 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 discuss only 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 quadratic representation of the LiouvilleanL + : where M is a complex Hermitian matrix parametrizing the Lindblad operators

Reduction to normal master modes
Next we want to show that the representation (22) allows us to reduce it further by a linear transformation of the set of maps {ĉ j ,ĉ † j ; j = 1, 2, . . . , 2n} to normal master modes (NMM) in terms of which the complete spectrum of the Liouvillean, as well as its eigenvectors, can be explicitly constructed; in particular the zero-mode eigenvector which is just the physically relevant NESS.
In fact we proceed in analogy to Lieb et al. [18] and define 4n adjoint Hermitian Majorana mapsâ r =â † r , r = 1, 2, . . . , 4n: satisfying the anti-commutation relations in terms of which the Liouvillean (22) can be rewritten aŝ where A is an antisymmetric complex 4n × 4n matrix with entries ½ is an identity map over K and A 0 is a scalar Obviously, the bi-linear Liouvillean (26) cannot be brought to a normal form with a linear canonical transformation since the matrix A -which shall in the following be referred to as a shape matrix of Liouvillean -is not anti-Hermitian like in Hamiltonian systems. 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 Av = βv with complex eigenvalue β, then v is also a left eigenvector with eigenvalue −β, A T v = −Av = −βv. Hence eigenvalues always come in pairs β, −β. Let as assume that A can be diagonalized §, i.e. there exist 4n linearly independent vectors v r , r = 1, . . . , 4n with the corresponding eigenvalues ordered such that Re β 1 ≥ Re β 2 ≥ . . . ≥ Re β 2n ≥ 0. The 2n complex numbers β j shall be referred to as rapidities. It is easy to check that we can always choose and normalize v r such that Let V be 4n × 4n matrix whose rth row is given by v r , V rs := v r,s . Then eqs. (29,30) rewrite as Expressing V T in terms of (32) and plugging the result into eq. (31) we arrive at a very convenient canonical form of a generic complex antisymmetric matrix A Now we apply decomposition (33) to the Liouvillean (26) Let us define the NMM mapsb : is not known at present whether explict form (27) guarantees diagonalizability of any such A. Note that one can construct certain types of complex antisymmetric matrices with degenerate eigenvalues which cannot be diagonalized [27].
For a non-degenerate rapidity spectrum {β j } the proof of this statement is a trivial consequence of antisymmetry A = −A T , whereas in case of degeneracies it can be shown that one can always choose appropriate linear combinations of eigenvectors.
We note that due to (25,30) NMM maps satisfy almost-canonical anti-commutation relations namelyb j could be interpreted as annihilation map andb ′ j as a creation map of jth NMM, but we should note thatb ′ j is in general not the Hermitian adjoint ofb j [28]. In terms of NMM the Liouvillean (34) now achieves a very convenient normal form where We shall later show that the constant B 0 is in fact equal to 0.

Non-equilibrium steady states and a complete spectrum of the Liouvillean
The Liouvillean can always be represented in terms of a large but finite 4 n × 4 n matrix.
We shall now outline the procedure of complete construction of its spectrum in terms of NMM which are easy to calculate in terms of diagonalization of 4n × 4n matrix A as described in the previous subsection.
We proceed by constructing the Liouvillean 'vacuum'. From the representation (22) it follows immediately that 1| = P (0,0...,0) | is left-annihilated byL + , 1|L + = 0, or equivalentlyL † + |1 = 0. So we have just shown that 0 is always an eigenvalue ofL + , hence there should also exist the corresponding right eigenvector |NESS , normalized as 1|NESS = tr ρ NESS = 1, which represents physical NESS, i.e. stationary solutions of the Lindblad equation (1) Let us define NMM number maps asN j :=b ′ jb j . From eqs. (36) it follows thatN j satisfy a projection propertyN 2 j =N j , so they are diagonalizable since no nontrivial Jordan block 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 a vacuum state on which allN j have value 0. It follows from the stability of completely positive evolution (1) that all eigenvalues λ ofL + should obey Re λ ≤ 0, and since by assumption Re β j ≥ 0, 1| and |NESS should be the left and right vacua which are simultaneously annihilated by NMM maps and hence alsoN j |NESS = 0. Thus we have also shown that the NMM representation (37) is only consistent if B 0 = 0 so we find an interesting sum rule for rapidities The complete excitation spectrum and the corresponding left/right eigenvectors of the Liouvillean are given in terms of a sequence of 2n binary integers (NMM occupation numbers) ν = (ν 1 , ν 2 , . . . , ν 2n ), ν j ∈ {0, 1}, where by construction, left and right eigenvectors satisfy the bi-orthonormality relation

The main general results: uniqueness of NESS, rate of relaxation to NESS, and expectation values of physical observables
Given a physical observable X ∈ K + and an arbitrary initial state with a density operator ρ 0 ∈ K, the time dependent expectation value of X can be written in terms of the spectral resolution of the Liouvillean, We remind the reader thatL + correctly represents physical Liouvillean only on the subspace K + of operators with even number of a-fermions. However, since the dynamics is closed on K + and test physical observable X also belongs to K + it follows that the component of ρ 0 from K − does not contribute to the expectation value X(t) . Given the exact and explicit constructions developed in this section we can now make the following rigorous and useful conclusions, assuming throughout that Liouvillean shape matrix A is diagonalizable: (1) is unique if and only if the rapidity spectrum {β j } does not contain 0, in our ordering convention, if β 2n = 0. In the opposite case, if we have d ≥ 1 vanishing rapidities, then we have a 2 d dimensional convex set of non-equilibrium steady states which can be spanned with |Θ R (0,...,0,ν 1 ,...,ν d ) .
Theorem 2: An arbitrary initial state with a density operator ρ 0 ∈ K converges with time to NESS if and only if all rapidities have strictly positive real parts, Re β j > 0. Then, the rate of exponential relaxation to NESS is given by the spectral gap ∆ of the Liouvillean which equals ∆ = 2 Re β 2n .
Theorem 3: Assume that the rapidity spectrum does not contain 0, i.e. β 2n = 0. Then the expectation value of any quadratic observable w j w k in a (unique) NESS can be explicitly computed as The statements of theorems 1 and 2 simply follow from exact and explicit spectral decomposition (42,43,44).
The proof of theorem 3 is also straightforward: The first expression (46) follows from the definition of the annihilation maps (9) and the explicit representation of the density operator of NESS, ρ NESS , in the canonical basis P α . The second, very useful equality (47) is then obtained by expressingĉ j thru (24) in terms of NMM maps (35) and using the annihilation relations (39).
The quadratic correlator of theorem 3 covers many physically interesting observables such as densities or currents. However if one needs expectation values of more general observables, e.g. an expectation value of a high order monomial 2n |NESS , then one may use a Wick theorem and rewrite it as a sum of products of pair-wise contractions (46).

Trivial example: A single fermion in a bath
In order to illustrate the method and demonstrate convenience of the results derived in the previous section we first work out a simple example of a single fermion n = 1 (or equivalently, an arbitrary qubit, a two-level quantum system), in a thermal bath. We take the most general single fermion Hamiltonian H = −ihw 1 w 2 +const = 2hc † c+const ′ and the following bath operators (see e.g. [16,29]) where the ratio of coupling constants determine the bath temperature T , Γ 2 /Γ 1 = exp(−2h/T ). Leaving out the details of a straightforward calculation, simply following the steps of the previous section, we arrive at the following shape matrix of the Liouvillean (26) where and Γ ± := Γ 2 ± Γ 1 . Further, we compute NMM rapidities β 1,2 = 1 2 Γ + ± ih and the eigenvector matrix Then, using theorem 3 we compute the expectation value of occupation number which is what we expect in canonical equilibrium.

Non-trivial example: transport in quantum spin chains
Here we work out a physically more interesting example, namely we construct NESS for the magnetic and heat transport of a Heisenberg XY spin 1/2 chain, with arbitrary -either homogeneous or positionally dependent (e.g. disordered) -nearest neighbour interaction which is coupled to two thermal/magnetic baths at the ends of the chain, generated by two pairs of canonical Lindblad operators [29] (similar to (48)) where σ ± m = σ x m ± iσ y m and Γ L,R 1,2 are positive coupling constants related to bath temperatures/magnetizations, for example if spins were non-interacting the bath temperatures T L,R would be given with Γ L,R 2 /Γ L,R 1 = exp(−2h 1,n /T L,R ). Applying the inverse of Jordan-Wigner transformation (6), where W = w 1 w 2 · · · w 2n is a Casimir operator which commutes with all the elements of the Clifford algebra generated by w j and squares to unity W 2 = 1. Therefore, W does not affect the action of bath operators (15) where L µ enter quadratically, so we find leading to the bath shape matrix (50) identical to the single fermion case (48). Again, carefully following the steps of section 2, we derive the Liouvillean in the form (26) with 4n × 4n shape matrix, which we write in a block tridiagonal form in terms of 4 × 4 matrices as and − (in terms of (50)), with Γ L,R ± := Γ L,R 2 ± Γ L,R 1 , and 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 spatially homogeneous case J x,y m ≡ J x,y , h m ≡ h it is not possible to proceed like in the classical harmonic oscillator chain where the corresponding matrix is a sum of a Toeplitz and a bordered matrix [30]. Namely, in our case A is a sum of a block Toeplitz and block bordered matrix and its explicit exact diagonalization remains an open problem. However, we stress that even relying on numerical diagonalization of A yielding a set of rapidities β j and properly normalized eigenvector matrix V, represents a dramatic progress with respect to previously existing numerical methods which needed diagonalization of matrices which were exponentially large in n. We shall later derive some exact theoretical and analytical results, explaining results of exact numerical computations, in the special case of a homogeneous transverse Ising chain (subsection 4.1), and the case of a disordered XY chain (subsection 4.2) for which we shall relate NMM to the problem of Anderson localization in one dimension, Let us continue by discussing transport observables in the spin chain whose expectation 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 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 − 2 (second equality sign), and (ii) [H m , H m ′ ] = 0 if |m − m ′ | ≥ 2 (third equality sign). Using eq. (47) of theorem 3 we can now compute NESS expectation values of energy density H m and energy current Q m , and also of somewhat simpler spin density and spin current which are all quadratic in w j . Note, however, that the spin density satisfies continuity equation (d/dt) σ z m = − S m + S m−1 only in the isotropic case, when J x m ≡ J y m .

Homogeneous transverse Ising chain
Here we limit our discussion to the spatially homogeneous case J x,y n ≡ J x,y , h n ≡ h. We shall show that in this case the eigenvalue problem and the baths playing the role of inelastic (absorbing) scatterers at the edges of a finite lattice. The 'elastic' (Hamiltonian) version of this trick has been used to compute temporal correlation functions in kicked Ising chain [31]. The singularity condition for the free mode equation (65) results, for a general homogeneous XY model, in eight master bands -different values of momenta ξ for each value of the spectral parameter (rapidity) β. In order to simplify the discussion -which will still get rather involved -we shall in the following restrict ourselves to the transverse Ising model J x = J, J y = 0. In this case we find just two master bands with simple dispersion relations but each band is doubly-degenerate, since the corresponding amplitude problem (65) has two linearly independent solutions Naively speaking, ξ − represents left moving and ξ + right moving free modes, each having two possible polarizations. Note that ξ − ξ + = 1. For a general complex β we shall choose the branch of square root ω(β) (66) for which |ξ − | ≤ 1. Let us now write the scattering problem on the left bath in terms of an ansatz where u L represents a 4-dimensional vector of left-most eigenvector components, ψ − 1,2 are the amplitudes of (known) incident free modes, and ψ + 1,2 are the amplitudes of the scattered, 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,2 , u L . Eliminating four variables u L we finally arrive to the non-unitary 2×2 S-matrix ) Similarly, one can solve the scattering problem from the right bath with the scattering ansatz Note that since the two directions of free modes (67) do not have left-right symmetry an explicit expression for S R is considerably more complicated than (70) and shall not be written out here. We shall now show that there exist two qualitatively different types of NMM -complex rapidities β solving (64) for sufficiently large n. First, we shall discuss the so called evanescent normal master modes. These are characterized with amplitudes (68) which decay exponentially with the distance from -say the left -bath, so the other -the right boundary condition becomes physically irrelevant in the limit n → ∞. Such solutions ψ + 1,2 = 0 of eq. (69) exist exactly when the determinant of S-matrix vanishes det S L = 0. Using (70) the determinant can be written as det Thus, for sufficiently large spin chains we find at most four NMM whose rapidities are given as the roots of 4-th order polynomial in β 2 p L (β evan ) = 0 (74) that are not simultaneously zeroes of τ (β). Clearly, such NMM asymptotically do not depend on the chain size n. In addition, we find evanescent NMM corresponding to the right bath simply by replacing Γ L + by Γ R + in (74,73). In fig. 1 we compare evanescent 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.5, h = 1.0, strongly coupled to two baths at considerably different temperatures, Γ L 1 = 1, Γ L 2 = 0.6, Γ R 1 = 1, Γ R 2 = 0.3 Note that the same parameter values will be used for numerical demonstrations throughout this subsection.
Second, we shall discuss the other extreme of soft normal master modes with rapidities that are closest to the imaginary axis, and thus determining the spectral gap of the Liouvillean and relaxation time to NESS. Composing the scattering from the two baths with the free propagation along the chain (back and forth) we arrive at the general secular equation for the eigenvalue problem (64) in terms of a 2 × 2 determinant In the absence of the baths, Γ L,R ± = 0, the solutions of the above problem exist only for real quasi-momenta, namely ξ ± = exp(±iϑ), ϑ ∈ R. For such extended master modes the local coupling to the baths can be considered as a small perturbation, thus only slightly perturbing the Bloch-like bands β(e iϑ ) = ±iε(ϑ) with 'energy' The softest NMM, namely the one for which the coupling to the baths is expected to be 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 the following we shall focus our calculation on the band edge β * = i(|h| + |J|) which, as can be checked aposteriori by a straightforward but tedious calculation, always gives smaller real part of the rapidity than the lower edge i(|h| − |J|), and hence really determines the gap of the Liouvillean. So we write where z ∈ C is a small parameter, and expand the S-matrices around the band edge where g := |hJ| 2(|h|+|J|) , η L,R := (Γ L,R + ) 4 + 4(Γ L,R + ) 2 (4h 2 + 2|hJ| + J 2 ) + 16h 2 J 2 and and Next we expand ξ + (66) in z, yielding and so the free propagator in (75) can be written as In eqs. (78,81,82) the branch cut along the negative real axis has been chosen for √ −iz. Since the product of S-matrices in (75) is near identity, the free propagator should be near one as well, hence 2ng −1 √ −iz should be near 2πi. Let us define z 0 by setting 2ng −1 √ −iz 0 = 2πi, so and write z = z 0 (1 + y) where |y| ≪ 1 is another small complex parameter. However, since z 0 is purely imaginary, we need to compute a small but non-vanishing y which will, in the leading order in n, solve (75) since the real part of the soft mode's rapidity is determined as Re β = Re z 0 y = π 2 g 2 n −2 Im y Now, writing √ −iz = √ −iz 0 √ 1 + y = iπgn −1 (1 + y/2 − y 2 /8) + O(y 3 ) in (78,82), plugging all that to eq. (75) and computing to order O(n −2 ), noting that O(|z|) = O(n −2 ), we arrive to a simple quadratic equation for y, whose solution, plugged to (84), gives the final result, namely the sectral gap of Liouvillean ∆ = 2 Re β fig. (2) we compare this analytical result to exact numerical calculations of the eigenvalue of A with minimal real part, confirming both, its precise numerical value and that the relative scaling of the next order correction is indeed O(n −1 ).  terms of a set of rapidities, for a short chain. In fig. 4 we demonstrate eq. density (59) and spin density (62) profiles in NESS. Again, we note flat profiles in the bulk of the chain, m, n − m ≫ 1, with exponential falloff due to adjustment to the non-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 m − H bulk | ∼ |ξ − | 4m is given by the quasi-momentum ξ − (66) corresponding to the maximal evanescent rapidity β evan (74).

Disordered XY chain
In this subsection we treat the opposite extreme, a disordered XY chain (52) where three sets of physical parameters are chosen as random uncorrelated variables from uniform distributions on the intervals, . Clearly, the eigenvalue problem (64) for the matrix (57) then becomes equivalent to the Anderson tight-binding problem in one dimension for a quantum particle with a 4−level internal degree of freedom. We do not pursue any theoretical analysis of this problem here, but merely state that numerical investigations indicate existence of exponential localization of all eigenvectors (or normal master modes) for disorder of any strength in anyone of system's parameters.
With the picture of localization of NMM in mind, the effect of the couplings to the heat baths at the chain's ends on quantum transport can be predicted by theoretical arguments (see [32] for a review of related studies): The spectral gap of the Liouvillean should be exponentially small ∆ ∼ exp(−n/ℓ) where ℓ is the localization length of NMM which is expected to be proportional to the square of inverse disorder strength. This is demonstrated in fig.6. If all NMM are exponentially localized, the currents should decrease with the chain size n faster than any power, perhaps exponentially, and the system should behave as an ideal insulator (at all temperatures). This is demonstrated by straightforward numerical calculations of the heat current (61) in fig. 7. In the final figure 8 we plot the energy density profile H m (59) in a typical case of disordered XY chain, versus a scaled spatial coordinate (m − 1)/(n − 1) ∈ [0, 1], for several different chain lengths n, and demonstrate sharping up of energy density profiles with increasing n, which is again indicating insulating behaviour.

Discussion and conclusions
The main result of the paper is a general method of explicit solution of master equations describing dynamics of open quantum system, under the condition that the system's Hamiltonian is quadratic and all Lindblad operators are linear in canonical fermionic operators (which can either represent real physical fermions or any abstract 2-level quantum systems (qubits) thru the Jordan-Wigner transformation). Using a novel concept of Fock space of physical operators (or density operators of physical states), and the adjoint structure of canonical creation and annihilation maps over this space, the problem can be treated in terms of a non-Hamiltonian generalization of the method of Lieb, Schultz and Mattis [18] lifted to an operator space. We have explicitly constructed a non-canonical analog of Bogoliubov transformation of the quantum Liouville map to normal master modes. Related ideas in the Hamiltonian context have been used by the author [31,33,34] in order to approach the problem of real time dynamics and ergodic properties of isolated interacting many-body quantum systems.
As an illustration of the method we have solved far from equilibrium quantum heat and spin transport problem in Heisenberg XY spin 1/2 chains which are coupled to canonical heat baths only at the two ends. Irrespectively of the strength of the coupling to the baths and their temperatures, we have shown a ballistic transport in the spatially homogeneous (non-disordered) case, and an ideally insulating behaviour in the disordered case associated to localization of normal master modes of the quantum Liouville operator. In this context the method can be considered as a simple alternative to the solution of quantum Langevin equations [24].
However, the method should easily be applicable to variety of other physical situations, for example if all fermions are coupled to the baths one could make a solvable model of genuine quantum diffusion, a many-body generalization of the tightbinding model [35]. We also expect the method to be applicable to the Redfield type of 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 of the method to open many-boson systems should be straightforward, simply by replacing anticommutators by commutators throughout the exposition of section 2.
Treating density operators of NESS as elements of a Hilbert space of operators one may also extend the concept of entanglement entropy, with respect to a bipartition of a system of many fermions [36], to NESS which can in our approach be viewed as a kind of ground 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 we propose to explore a question, whether more involved algebraic methods of solution of interacting many-body quantum systems, like e.g. Bethe Ansatz or quantum inverse scattering [13], could be generalized to open quantum systems, e.g. by means of the proposed concept of Fock space of operators. Could one discuss completely integrable open quantum systems which go beyond quadratic Liouvilleans?