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

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 4n x 4n matrix. Non-equilibrium thermal driving in terms of the Redfield equation is analyzed in detail. We explain how to compute all 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 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 thermal conductance for sufficiently strong thermal driving, as well as non-monotonic dependence of the heat current on the strength of the bath coupling.


Introduction
One of the main challenges of the 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 this type 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,2,3]. However, in this context explicit results which 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 either in relaxation to equilibrium or nonequilibrium steady states, 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] as may complement the linear response calculations and suggest non-linear response or far-from-equilibrium effects.
However, to date we have had a very few explicit calculations of 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 freedoms (such as spins, or bosons), see for example [7,8]. The reason is that there has been no theoretical techniques 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 manybody 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 which 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 is 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 nonequilibrium Redfield model. We shall outline several intriguing exact numerical results on large spin chains. In particular, we show that recently announced quantum phase transition in the open XY chain in the local Lindblad bath model, generalizes also to non-equilibrium thermal Redfield model with qualitatively identical characteristics. The transition is characterized by spontaneous emergence of long range magnetic correlations, and hypersensitivity of the steady state to external system's parameters, when the transverse magnetic field drops bellow the critical value |h| < h c = |1 − γ 2 | where γ is the anisotropy parameter. Furthermore, we analyze in some detail the heat transport in 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.

Markovian master equations in non-equilibrium quantum physics
Decomposing the Hilbert space of the universe into a tensor product H = H s ⊗ H b of the central system H s and the bath (or a set of baths) H b (environment), one writes the total Hamiltonian as where X µ , are linear operators over H s , and Y µ linear operators over H b . Note that X µ , Y µ 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 systems's density matrix ρ(t) is derived [4] using three main assumptions: (i) weak coupling (assuming λ to be small), (ii) factorizability of the initial density matrix ρ s (0) ⊗ ρ b (0), and (iii) Born-Markov approximation which rests upon the assumption that the bath-correlation functions decay on much shorter time scale than the central systems dynamicsX µ (t) := e itHs X µ e −itHs . We use units in which Planck's constant = 1, and may use different inverse temperatures β for different pieces of the environment (for different baths). The resulting master equation is referred to as the Redfield equation d dt where the dissipator-map has a memoryless kernel with the following general form If one additionally assumes the so-called rotating wave-approximation, one arrives at the dynamical semi-group which manifestly preserves the positivity of density matrix at all times ‡ and can be generally described by the dissipator in the Lindblad form where the only condition is that γ is a Hermitian γ µ,ν = γ * ν,µ and positive definite matrix. The standard Lindblad form is obtained by diagonalizing the matrix γ whose eigenvectors yield the usual Lindblad operators. The important property of the bathcorrelation functions (2) (which constitute all that we need to know about the baths) is the Kubo-Martin-Schwinger(KMS) condition which is needed to prove that the thermal state ρ gibbs = e −βHs / tr e −βHs is a steady state of the master equation (3), provided that all baths are thermalized to the same inverse temperature §. However, in case of several thermal baths with possibly different temperatures we may expect that ρ(t) relaxes to a physically very interesting nonequilibrium-steady-state (NESS).

Diagonalization of quantum Liouvilleans for quadratic fermi systems
In this section we give a short account on the general technique of canonical quantization in the Liouvile 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-comuting Hermitian operators w j , j = 1, 2, . . . , 2n, {w j , w k } = 2δ j,k , in which the Hamiltonian H may take a general quadratic form and the coupling operators may be general linear forms: 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 . This formalism is immedately applicable either for describing, (i) physical fermions c m , m = 1, 2, . . . , n, where

Fock space of operators
The fundamental concept for our analysis is a Fock space structure over the 4 n dimensional Liouville space of operators K, which density matrix ρ(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 x † y, x, y ∈ K.
We note that 2 2n operator-products |P α , labelled with a binary multi-index α P α 1 ,α 2 ,...,α 2n : constitute a complete orthonormal basis of K with respect to an inner product.
In fact it is easy to show that |P α is a fermionic Fock basis, and powers 1 in the product (11) can be considered like a sort of Fermionic excitations, if we define the following set of linear annihilation mapsĉ j over K and derive the actions of their Hermitian adjoints -the creation linear mapsĉ † , which satisfy canonical anticommutation relations

Bilinear form of the Liouvillean
The aim is now to show that the generator of the master equation (3) is in general a quadratic form in these adjoint fermionic mapsĉ j ,ĉ † j . In order to see that clearly, let us define the left and right multiplication maps over K Inspecting the actions ofŵ L j ,ŵ R j on the Fock basis |P α one arrives at the following useful identitieŝ We shall use notation where linear maps over the operator space (in physics literature sometimes referred to as "super-operators") are designated byˆ. are a parity map, and a number map, respectively, which count the parity and number of the adjoint fermionic excitations (number of factors in (11)). Note thatP, anticommutes with allĉ j ,ĉ † j , hence the second equality of (18), andP 2 =1. The unitary part of the Liouvillean (15) now trivially reads The dissipator (4) can be represented as a map over K aŝ where f µ (t) is a (real-valued) propagator of Heisenberg dynamics in the closed system which -due to (20) -can be explicitly solved for a quadratic Hamiltonian (7), giving are fundamental basis dissipators which using (17,18) evaluate tô It will prove useful if we express the internal dynamics (23) explicitly in terms of eigenvalues and eigenvectors of the Hamiltonian matrix H. Since 2n × 2n matrix is anti-symmetric and Hermitian, its real eigenvalues come in pairs m , − m , j = 1, . . . , n, with the corresponding eigenvectors u m , u * m , namely Hu m = m u m and Hu * m = − m u * m since H * = −H. The eigenvectors may and should always be chosen orthonormal (even in the case of degeneracies), meaning Then the spectral decomposition of the Heisenberg dynamics reads Note thatP ± = (1 ±P)/2 are orthogonal projectors which commute with all the terms (20,25,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 even/odd number of fermionic operators. Since we are mainly interested in expectation values of even observables, such as currents and densities, we shall in the present paper focus on the dynamics in the subspace K + only, and consider the corresponding The extension to the odd parity subspace is straightforward. Collecting the results (20,21,25,26) it is now obvious thatL + is a bilinear form inĉ † j andĉ j . For convenience, we define 4n Hermitian Majorana mapsâ r , r = 1, . . . 4n and express the Liouvillean aŝ where the 4n × 4n complex antisymmetrix matrix A, later referred to as a structure matrix, and a scalar A 0 , can be expressed as where M is a 2n × 2n bath-matrix which can be compactly written as Defining the bath-spectral functionsΓ β µ,ν (ω) : and extending the range of integration in (34) to [−∞, ∞], or better to say, 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 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 (Γ β µ,ν (τ )) * = Γ β ν,µ (τ ) that only if the Cauchy principal value terms are dropped (i.e. if the range of integration in (4) is extended to [−∞, ∞]) the Redfield dissipator annihilates the Gibbs stateD|e −βHs = 0, and hence Gibbs state is the steady state of equilibrium thermal Redfield model.
Note again that the inverse temperature in (36) could in principle be a function of the bath-index β = β ν in case one would be interested in non-equilibrium situation with couplings to several different temperatures. But we should stress that different temperatures only make sense among uncorrelated baths for which Γ β µ,ν ≡ 0 for any β. 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) = γ µ,ν δ(t + 0), and then the bath-matrix reduces to a Hermitian form M = ν,µ γ ν,µ x ν ⊗ x µ = M † which is equivalent to the one used in [11].

Static Liouvillean: normal modes, non-equilibrium steady state 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 x µ -we follow Ref. [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 × 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 , −β j , j = 1, . . . , 2n, assuming Re β j ≥ 0, and its eigenvectors v 2j−1 (corresponding to β j ), and v 2j (corresponding to −β j ) can always be normalized -irrespective of possible degeneracies of among β j , which shall be called rapidities -such that where V is 4n × 4n matrix whose rth row is given by v r , V r,s := v r,s . Thus the structure matrix allows the following decomposition which after plugging into the Liouvillean (31) immediately brings it to a normal form are the normal-master-mode (NMM) maps, satisfying almost canonical anticommutation relations The mapb j could be interpreted as an 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 . The right-vacuum is now essentially defined byb j |NESS = 0, whereas the left-vacuum is trivial 1|L + = 0 and satisfies 1|b j = 0.
Assuming that the whole rapidity spectrum is strictly away from the real line Re β j > 0, we state the following exact results: (i) |NESS is unique.
(ii) Almost any initial density matrix relaxes to NESS with an exponential rate ∆ = 2 min Re β j (the spectral gap of the Liouvillean). The complete spectrum of 4 n eigenvalues ofL + is obtained by all possible binary linear combinations λ ν = −2ν·β, ν j ∈ {0, 1}.
(iii) The expectation value of any quadratic observable w j w k in a (unique) NESS can be explicitly computed as where is a matrix of the non-equilibrium Green's function. The first equality is proven in [11] ¶ whereas the last equality requires a simple contour integration on the spectral decomposition of the resolvent (45).
(iv) The Wick theorem may be used for calculation of expectation values of arbitrary higher order (even!) observables by 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 Ref. [13] for a detailed discussion of Liouvillean degeneracies).

Static Liouvillean: time-dependent correlation functions
The complete Liouvillean propagator can be written explicitly as 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 dynamicsL − we take a pair of simplest even-order, quadratic observables, and define the corresponding non-equilibrum time-dependent correlation function -or non-equlibrium response function -as (47) ¶ Small simplification has been made with respect to the statement of Theorem 3 of Ref. [11] which has been pointed out by I. Pižorn [12].

Time-dependent Liouvilleans
In this subsecton we indicate how to efficiently treat explicitly time-dependent master equations, written in third quantized form as d dt where explicit time-dependece 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 whereT indicates a time-ordered product. The procedure is the following. Note that the space of all anti-symmetric complex structure matrices form a Lie algebra so(4n, C). The following straightforward identity holding for any pair of complex 4n×4n matrices A, B, indicates that Liouvilleans (31,49) generate 4 n dimensional representation of so(4n, C). Thus, the time-ordered product (50) can be evaluated within a Lie group SO(4n, C) of 4n × 4n matrices, and + then full Liouvillean propagator is written aŝ The logarithm ofÛ can be now considered as a 'static' Liouvillean, so we can diagonalize it by the methods of subsection (3.3), leading to spectral decomposition of the form (46). + Even if this has to be done numerically, using Trotter-Suzuki decomposition schemes, the computational compexity is only polynomial in n.

XY spin chains
The theory of the previous two sections shall now be applied to investigate a homogeneous, finite XY chain of n spins, described by Pauli matrices σ x,y,z j , j = 1, . . . n with the Hamiltonian which is described by two real parameters, anisotropy γ and transverse magnetic field h.
Without loss of generality we may assume that γ, h ∈ [0, ∞). We decide to couple XY chain thermally only at its ends, so we consider the most general four coupling operators which allow for an explicit solution and fully decorrelated baths Γ β µ,ν = δ µ,ν Γ β µ . We take standard baths of harmonic oscillators at two ends with possibly different inverse temperatures, and Ohmic spectral functionsΓ Note that frequency cutoff in the spectral function is irrelevant as we neglect the Lamb shift term in the master equation.
The enitre problem can be fermionized by means of Jordan-Wigner transformation (9), namely the Hamiltonian and the coupling operators transform to where W = (−i) n−1 w 1 w 2 · · · w 2n is an operator which commutes with all the elements of K + (or anti-commutes with all the elements of K − ) and satisfies W W † = W † W = 1, hence it has no effect on the dissipator (4) inL + . We note however, that the commutation of W thru ρ in (4) for the dynamics in K − produces a minus sign in all the bath terms, i.e. it changes the sign ofP −DP− , with respect to a pure fermionic problem. The 4n × 4n structure matrix has now a specific block-tridiagonal + block-bordered form, where a, b, c are 4 × 4 matrices The sequences of 4 × 4 matrices l j , l j , r j , r j which form the block-bordered part B can be straightforwardly computed [seeing (32)] from the form of the coupling vectors x 1,2 = (κ 1,2 cos θ 1,2 , κ 1,2 sin θ 1,2 , 0, . . . 0) T , x 3,4 = (0, . . . , 0, −κ 3,4 sin θ 3,4 , κ 3,4 cos θ 3,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 The coefficient K > 0 in general depends only on γ, h, and β L (for l j ) or β R (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 which remain are the diagonal ones l 1 and r n , given explicitly in Ref. [11]. Below we shall present some intriguing numerical results of the non-equilibrium thermal Redfield equation (3,4) for the open XY chain given by (54,55,56), in comparison with the local non-equilibrium Lindblad model (5) where a suitable set of coupling operators of the form (55) and 4 × 4 coupling matrix γ µ,ν can be chosen to parametrize the Lindblad operators L 1,2 = Γ L 1,2 σ ∓ 1 , L 3,4 = Γ R 1,2 σ ∓ n , parametrized exactly in the same way as in Refs. [11,14]. For all the numerical results reported for the thermal Redfield model we consider the bath parameter values κ 1 = κ 3 = 1, κ 2 = κ 4 = 0, θ 1 = θ 3 = π/6, and β L = 0.3, β R = 5.2 unless β's are varying, and λ = 0.1 unless λ is varying, whereas for the Lindblad model we always take the bath parameters Γ L 1 = 0.5, Γ L 2 = 0.3, Γ R 1 = 0.5, Γ R 2 = 0.1.

Non-equilibrium phase transition
In Ref. [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| < h c , However, that study was done with local Lindblad reservoirs, so the questions 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 subsection 3.3 and numerically evaluate the spin-spin correlator (note the use of Wick theorem as the spin-spin correlator is of fourth order in w j ): First, we use efficient prescription (43) to compute correlation matrices at nonequilibrium conditions β L = 0.3 = β R = 5.2 and plot them for two different system sizes and five different values of h around h c in figure 1. Results look qualitatively identical to those for the 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 C res = |l−m|>n/2 l,m |C l,m |/ |l−m|>n/2 l,m 1, which also reveals possible criticality in the region of a large anisotropy γ > 1 previously not discussed. We note that size dependence of the residual magnetic correlator C res shows a very characteristic behaviour: namely C res ∝ exp(−ηn) with η > 0 for |h| > h c or h = 0 (64) Thus we shall refer to the regime with 0 < |h| < h c as long range magnetic correlation (LRMC) phase * , the regime with |h| > h c , or h = 0, as non-LRMC phase, and the regime with |h| = h c as critical. Scaling (64,65) is illustrated in figure 3. Exponential decay of the C res (n) in non-LRMC phase (64) is consistent with the exponential decay of 2-point correlator with the distance between sites C(r) = j−i=r C i,j / j−i=r 1 ∼ exp(−ξr), as can be qualitatively noted already in the figure 1. However, we demonstrate in figure 4 that the exponents ξ could in principle be very different between the Redfield and local Lindblad models. Futhermore, as for the Linbdlad model the exponents ξ and η [of (64)] appear to be equal, for the Redfield model they don't seem to be simply related. Analytical estimation of these exponents present a challenge for future theoretical work. However, we note that with the thermal driving with Redfield dissipators, the longrange-magnetic order disappears when the temperatures of the baths become equal, β L = β R , and there we recover, consistently, all the properties of the thermal state [15] which are most easily numerically reproduced by the method of Ref. [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 C res (for large n in the LRMC phase) decreases as a function of the difference of inverse temperatures ∆β = β R − β L , namely numerics of figure 5 suggests clearly that C res ∝ (∆β) 2 .
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 η(γ, h) are still lacking. We note that the transition point h = h c is charac- * Note, interestingly, that unlike for the local Lindblad driving [14] the XX line γ = 0, 0 < |h| < 1, also exhibits long range magnetic correlations for the thermal Redfield driving.  terized by a simple property of the XY spin chain quasiparticle dispersion relation where j = ω(2πj/n) would be exactly the (positive) eigenvalues of matrix H if periodic boundary conditions would be imposed on the closed system. Namely, in non-LRMC phase |h| > h c there exist only a single pair of trivial stationary points q * = 0, π, whereas in LRMC phase |h| < h c there exist another pair of nontrivial stationary points ±q * = 0, π, dω/dq| q=q * = 0, which introduces a new non-trivial length scale 1/q * which determines typical sizes of correlated regions in the matrix C l,m (see figure 1). Therefore this simple non-equilibrium quasi-particle picture predicts mean-field critical exponent 1/q * ∼ |h c − h| −1/2 as h ↑ h c (confirmed in Ref. [14]). The non-equilibrium quantum phase transition can also be characterized by the scaling of the Liouvillean spectral gap ∆(n), namely in the critical regime one expects a qualitative increase in the relaxation time 1/∆ to NESS. Numerical results (see figrure 6) suggest the spectral gap of the Liouvillean remains like in the local Lindblad case [11] although we are at the moment unable to prove this conjecture. Also note slight fluctuations of ∆(n) in the LRMC phase as opposed to a smooth power law in the non-LRMC phase. Long range correlations for |h| < h c naturally imply sensitivity of NESS to tiny variations in system's parameters. For example, one may expect also that local observables in NESS will be then 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 s z = σ z n/2 NESS versus the field strength h. Indeed, we notice that for |h| > h c , s z (h) is a smooth function wheres for |h| < h c , s z (h) becomes rapidly oscillating or better to say, fluctuating, function. Even though the amplitude of these oscillations decreases with n, the scale of h on which s z (h) fluctuates decreases with n even much faster, so we predict that in the thermodynamic limit n → ∞, in LRMC phase the local susceptibility ds z /dh would be ill defined. In summary, LRMC phase can be characterized by hypersensitivity of NESS to external parameters.

Heat transport and entropies
An important non-equilibrium physical effect which one can investigate more deeply in an open XY chain is the heat transport, which has been recently intensively studied in quantum spin chains, see e.g. [17,18,19,20,21] or [22] for a recent review on the topic.
Writing the Hamiltonian (54) in the bulk as a sum H = m H m with a two body energy density operator one can derive the local energy current which, by construction, satisfies the continuity equation 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 been already 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 excercising the cyclicity of the trace again and transforming the integrand of (4) to terms of the form trX µ (−τ )ρ[X ν , H m ] ≡ 0. The RHS expression of eq. (70) then follows from the nearest-neighbour locality of the Hamiltonian. Therefore, in NESS the expectation value of the current Q m NESS should be independent of the position m. By looking at the dependence of the steady-state current on the system size we clearly find ballistic transport, namely Q m NESS = O(n 0 ), irrespectively of the temperature differences between the baths and bulk parameters of the model (i.e. whether being in the LRMC phase, non-LRMC phase, or critical). However, we find very interesting dependence of the heat current on the temperature driving, i.e. on the two temperatures of the thermal baths. In figure 8 we plot Q m NESS versus β L and β R and find a maximum of the current for intermediate driving, namely when one of the temperatures is less than one 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-equilibroum 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 w j w k NESS and the Wick theorem, one can adopt the recipe which has been proposed in Ref. [24] for computing block entropies (or entanglement entropies) applied to the entire lattice. In fact, taking an arbitrary block region A ⊆ {1, . . . , n}, one can compute Von Neumann entropy S A (ρ) = −tr A ρ A log 2 ρ A (in base 2), where ρ A = trĀρ is a reduced density matrix andĀ denotes the complement of A, as and ±iν j are the eigenvalues of the 2#(A) × 2#(A) part of the correlation matrix B j,k defined by w j w k NESS =: δ j,k + iB j,k , restricted to Majorana operators w j , w k corresponding to spins from the block A. The same general procedure has been applied to thermal (Gibbs) states in Ref. [16]. When taking the maximal block A = {1, . . . , n} we obtain exactly the standard Gibbs entropy of NESS. In figure 9 we plot the Gibbs entropy S {1,...,n} 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 ∆β may be understood as the entropy production rate.
Calculation of Gibbs entropy of NESS provides also 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 β's), or for very strong bath coupling λ, the positivity of NESS might be slightly violated (red region in figure 9), namely some of the correlation matrix eigenvalues ν j become slightly larger than 1 (but in our numerical experience never by more than 10 −7 or so). Figure 8. NESS expectation value of the heat current Q m NESS versus two inverse temperatures β L and β R for the non-equilibrium thermal Redfield model of an open XY chain with γ = 0.5, h = 0.9, sistem size n = 53, and bath parameters given in the text. Note that the 'shoulders' of maxima, around β L ≈ 0.05, β R > 1, and with L and R exchanged, could be interpreted as negative differential heat conductance.
We can use the concept of block entropy of NESS to further characterize the nonequilibrium 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} . Figure 9. Gibbs entropy of NESS versus two inverse temperatures β L and β R for the same parameters as in the 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.
Interestingly, we find (see figure 10) that QMI saturates I(n) = O(n 0 ) in the non-LRMC phase (for |h| > h c ), whereas in LRMC phase (for 0 < |h| < h c ) QMI becomes extensive I(n) = O(n) indicating a drastic enhancement of correlations in NESS. This is again very similar to the behaviour of operator space entanglement entropy (OSEE) (analized for the Lindblad model in [14]), so one may extend the relationship between QMI and OSEE which has been conjectured for thermal states in Ref. [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) = | H m NESS −H|/|H| whereH = (n−3) −1 n−2 m=2 H m NESS is the averaged energy density. Quite strikingly, we observe a big variation of f (m) from site to site for LRMC phase and very smooth (non-fluctuating) behaviour for the non-LRMC phase which is characterized with a bulk-constant f (m) which is exponentially small in n. This behaviour can again be considered as a manifestation of hypersensitivity of NESS and LRMC. At last we check the dependence of the heat current Q m NESS on the system-bath coupling strength λ. It was recently reported by Karevski and Platini [25] that the spin current J m in the local Lindblad model of an open isotropic XX chain γ = 0 has a non-monotonic dependence on λ which can be universally described by a formula J m NESS = a λ 2 /(b + λ 4 ) where a , b are 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 a very similar behaviour for the heat current where a, b are again some constants which may depend on all system's parameters except λ. This is particulary interesting as in the anisotropic XY model the spin current is not even well defined as there is no corresponding conservation law. This behaviour is demonstrated in figure 12 where on 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.

Discussion
The purpose of the present paper was three-fold. Firstly, we have outlined a general method for 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 obsevables 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 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 non-equilibrium steady states, for example such states which carry non-vanishing heat current. We wish to stress this physically obvious but mathematically delicate point with a particular care, as we have found a qualitatively different result for Lindblad-Davies dissipators which generate proper dynamical semigroups and satisfy detailed balance condition with respect to Gibbs states [26,5]. 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 have found that the resulting steady state (fixed point of the Liouvillean dynamics) is simply some convex combination of two Gibbs states corresponding to the bath temperatues, and as such has always zero heat current and cannot represent 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 which can describe non-equilibrium thermal driving of a (non-self-thermalizing, e.g. integrable) open quantum system.