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

Tomaž Prosen

Department of Physics, FMF, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia

E-mail: tomaz.prosen@fmf.uni-lj.si

Received 7 January 2008
Published 15 April 2008

Contents

• 1. Introduction
• 2. General method of solution for the Lindblad equation
  • 2.1. Fock space of operators
  • 2.2. Reduction to normal master modes
  • 2.3. Non-equilibrium steady states and a complete spectrum of the Liouvillean
  • 2.4. The main general results: uniqueness of NESS, rate of relaxation to NESS and expectation values of physical observables
• 3. Trivial example: a single fermion in a bath
• 4. Nontrivial example: transport in quantum spin chains
  • 4.1. Homogeneous transverse Ising chain
  • 4.2. Disordered XY chain
• 5. Discussion and conclusions
• Acknowledgments
• References

1. Introduction

Understanding the 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 the 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]–[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]–[12], the field of open quantum systems is still lacking nontrivial 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]–[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 (NMM), 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 and currents.

We demonstrate the power of this novel method by applying it to the problem of heat and spin transport far from equilibrium in nearest neighbour 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. [12], [19]–[24] for related recent studies of quantum thermal conductivity in one dimension (1D)), and predict ideally insulating behaviour (at all temperatures) in a disordered case of spatially random interactions/field. Apart from obtaining numerical results which go far beyond what has been accessible so far 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 relaxation to the NESS and show that it scales with the inverse cube of the chain length. (ii) We construct evanescent NMM of the Liouvillean, for long chains, by which we explain quantitatively the exponential fall-off of energy density or temperature profiles near the bath contacts.

The paper is organized as follows. In section 2, we 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.

2. 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  [17, 25] 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₁, x₂,  ...)^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 ,

Here, we are not concerned with physical criteria for the validity of the so-called Markovian approximation under which equation (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 completely general and possibly non-Markovian open system dynamics, is straightforward. See e.g. [26] for a discussion of Markovianity.

2.1. Fock space of operators

We begin by associating a Hilbert space structure x→|x to a linear 2²ⁿ = 4ⁿ dimensional space 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 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  over :

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 map defined by equations (1, 2, 3) can be written as a quadratic form in adjoint Fermi maps  and (or for short, a-fermions)^Note1 .

Firstly, we consider the unitary part of Liouvillean

Since is a Lie algebra, one defines the adjoint representation of a Lie derivative for an arbitrary back on as,  ad A|B:  = |[A,B]. It is now straightforward to 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 , it follows that for an arbitrary quadratic Hamiltonian (2) its Lie derivative has a very similar quadratic form in a-Fermi maps

It is worth stressing here that for an arbitrary (complex) matrix H and (14) conserves the total number of a-fermions , namely .

Secondly, we consider the action of the Lindblad maps

where we write . Again we proceed by computing the actions of on elements of the canonical basis of operator Fock space . 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 in |P_α, namely P_αw_j = (–1)^{|α| + α_j}w_jP_α, and hence

Observing that

one derives from equation (16) the general expression for

Obviously, the maps , and hence also the total Lindblad part of Liouvillean , 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 , with respect to which the operator space can be decomposed into a direct sum , and even/odd operator spaces are orthogonally projected as . Thus, the positive parity subspace is a linear space spanned by |P_α with even |α|. All the maps now act separately on and . For example, the maps defined on even parity subspace are indeed quadratic in a-fermions

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 . The extension to the dynamics of odd observables should be straightforward.

Putting the results of (12), (14), (15) and (21) together we arrive at the final compact quadratic representation of the Liouvillean

where M is a complex Hermitian matrix parametrizing the Lindblad operators

2.2. 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 to 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 :

satisfying the anti-commutation relations

in terms of which the Liouvillean (22) can be rewritten as

where A is an antisymmetric complex 4n×4n matrix with entries

is an identity map over and A₀ 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 –β and A^Tv = –Av = –βv. Hence, eigenvalues always come in pairs β and –β. Let us assume that A can be diagonalized^Note2 , i.e. there exist 4n linearly independent vectors v_r,r = 1,  ..., 4n with the corresponding eigenvalues β₁, –β₁, β₂, –β₂,  ..., β_2n, –β_2n,

ordered such that  Re β₁≥ Re β₂≥≥ 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^Note3 

Let V be 4n×4n matrix whose rth row is given by v_r, V_rs: = v_{r, s}. Then equations (29) and (30) rewrite as

Expressing V^T in terms of equation (32) and plugging the result into equation (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 maps or

We note that due to (25) and (30) NMM maps satisfy almost-canonical anti-commutation relations

namely, could be interpreted as annihilation map and as a creation map of jth NMM, but we should note that is in general not the Hermitian adjoint of .^Note4  In terms of NMM the Liouvillean (34) now achieves a very convenient normal form

where . We shall later show that the constant B₀ is in fact equal to 0.

2.3. 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ⁿ×4ⁿ 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 by , or equivalently . So we have just shown that 0 is always an eigenvalue of , 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 as . From equations (36) it follows that satisfy a projection property , so they are diagonalizable since no nontrivial Jordan block could satisfy the projection property. Furthermore, are mutually commuting , so they can be simultaneously diagonalized and there should exist a vacuum state on which all have value 0. It follows from the stability of completely positive evolution (1) that all eigenvalues λ of 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 also . Thus, we have also shown that the NMM representation (37) is only consistent if B₀ = 0 so we find an interesting sum rule for rapidities

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) ν = (ν₁, ν₂,  ..., ν_2n), ν_j{0, 1},

namely,

where by construction, left and right eigenvectors satisfy the bi-orthonormality relation Θ^L_{ν  '}|Θ^R_ν = δ_{ν  ', ν}.

2.4. The main general results: uniqueness of NESS, rate of relaxation to NESS and expectation values of physical observables

Given a physical observable and an arbitrary initial state with a density operator , the time-dependent expectation value of X can be written in terms of the spectral resolution of the Liouvillean,

namely,

We remind the reader that correctly represents physical Liouvillean only on the subspace of operators with even number of a-fermions. However, since the dynamics is closed on and test physical observable X also belongs to it follows that the component of ρ₀ from 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:

Theorem 1. NESS of Lindblad equation (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 NESS which can be spanned with |Θ^R_{(0,  ..., 0, ν1,  ..., νd)}.

Theorem 2. An arbitrary initial state with a density operator 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_jw_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 equations (42)–(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 through (24) in terms of NMM maps (35) and using the annihilation relations (40).

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 , then one may use a Wick theorem and rewrite it as a sum of products of pair-wise contractions (46).

3. Trivial example: a single fermion in a bath

In order to illustrate the method and demonstrate convenience of the results derived 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₁w₂ + const = 2hc^†c + const  ' and the following bath operators (see e.g. [16, 29])

where the ratio of coupling constants determines the bath temperature T, Γ₂/Γ₁ = 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 Γ_±: = Γ₂±Γ₁. Further, we compute NMM rapidities and the eigenvector matrix

Then, using theorem 3 we compute the expectation value of occupation number which is what we expect in canonical equilibrium.

4. Nontrivial 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}₂/Γ^{L, R}₁ = exp(–2h_1,n/T_{L, R}).

Applying the inverse of Jordan–Wigner transformation (6), σ^x_m = (–i)^m–1∏_j = 1^2m–1w_j and σ^y_m = (–i)^m–1(∏_j = 1^2m–2w_j)w_2m, we rewrite (52) and (53) in terms of Majorana fermions

where W = w₁w₂ ...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² = 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 A₀ = Γ ₊ ^L + Γ ₊ ^R, where , (in terms of (50)), with Γ_±^{L, R}:  = Γ₂^{L, R}±Γ₁^{L, R}, 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} and 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

as . One can derive the local energy current Q_m = i[H_m, H_m + 1] from the continuity equation

where Q_m: = i[H_m, H_m + 1]

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 equation (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)σ_m^z = –S_m + S_m–1 only in the isotropic case, when J^x_m≡J^y_m.

4.1. Homogeneous transverse Ising chain

Here, we limit our discussion to the spatially homogeneous case J_n^{x, y}≡J^{x, y} and h_n≡h. We shall show that in this case the eigenvalue problem

for equation (57) can be most easily and elegantly treated if formulated in terms of an abstract inelastic scattering problem in one dimension, with asymptotic solutions given in terms of free normal modes for the infinite translationally invariant chain v = ( ..., uξ^m–1, uξ^m, uξ^m + 1,  ...)^T, where ξ is a complex quasi–momentum (Bloch) parameter and u is a 4-dimensional amplitude vector satisfying the condition

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 and 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 four-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

with

Similarly, one can solve the scattering problem from the right bath with the scattering ansatz

defining the right S-matrix

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.

Firstly, we shall discuss the so-called evanescent NMM. 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 equation (69) exist exactly when the determinant of S-matrix vanishes det S^L = 0. Using equation (70) the determinant can be written as det S^L = (β/τ)²p^L(β) where^Note5 

Thus, for sufficiently large spin chains we find at most four NMM whose rapidities are given as the roots of fourth-order polynomial in β²

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 equations (73) and (74). In figure 1, we compare evanescent rapidities computed from equation (74) to numerically calculated spectrum of A, at several different sizes n, for a typical case of transverse Ising chain, J = 1.5 and h = 1.0, strongly coupled to two baths at considerably different temperatures, Γ^L₁ = 1, Γ^L₂ = 0.6, Γ^R₁ = 1 and Γ^R₂ = 0.3. Note that the same parameter values will be used for numerical demonstrations throughout this subsection.

Secondly, we shall discuss the other extreme of soft NMM with rapidities that are closest to the imaginary axis, and thus determine 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 . 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ε() 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 figure 1). In the following, we shall focus our calculation on the band edge β* = i(|h| + |J|) which, as can be checked a posteriori 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 is a small parameter, and expand the S-matrices around the band edge

where and

and

Next, we expand ξ ₊  (66) in z, yielding

and so the free propagator in (75) can be written as

In equations (78), (81) and (82) the branch cut along the negative real axis has been chosen for . Since the product of S-matrices in (75) is near identity, the free propagator should be near one as well, hence should be near 2πi. Let us define z₀ by setting , so

and write z = z₀(1 + y) where |y|  1 is another small complex parameter. However, since z₀ 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

Now, writing in (78) and (82), plugging all that to equation (75) and computing to order , noting that , we arrive at a simple quadratic equation for y, whose solution, plugged to (84), gives the final result, namely the spectral gap of Liouvillean Δ = 2 Re β

In figure 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 .

Figure 2. Spectral gap Δ times a third power of the chain length n for a transverse Ising chain with J = 1.5 and h = 1, and bath couplings ΓL1 = 1, ΓL2 = 0.6, ΓR1 = 1 and ΓR2 = 0.3. Thin horizontal line indicates the theoretical asymptotic value (85). In the inset, we show deviation from asymptotic constant value of Δn3 in log–log scale and demonstrate that it decays as n- 1 (thin line).

Note that, interestingly, both main analytical results of this subsection, namely evanescent and soft mode rapidities do not depend on Γ^{L, R}_–. Physically speaking, they only depend on the effective strengths of the bath couplings and not on the temperatures.

We end this subsection by presenting some further numerical results on heat transport in the open transverse Ising chain in the Lindblad form. In figure 3, we demonstrate expression (42) for constructing the full spectrum of the Liouvillean in terms of a set of rapidities, for a short chain. In figure 4, we demonstrate equation (47) of theorem 3 by computing the energy current Q_m (61), and the average spin current  (63) in NESS of a typical transverse Ising chain. Numerical results give a clear indication of ballistic transport  and , however, its rigorous proof and analytical calculation of the currents would require full control over the complete set of NMM which is at present not available. In figure 5, we plot the energy 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).

Figure 5. Energy density profile (lower, blue points), and spin density profile (upper, red points), for a transverse Ising chain of n = 80 spins with J = 1.5, h = 1 and bath couplings ΓL1 = 1, ΓL2 = 0.6, ΓR1 = 1 and ΓR2 = 0.3. The insets display logarithm of the difference to the bulk values δHm:  = |Hm-Hbulk| (blue points), δσzm:  = |σzm-σzbulk| (red points) in comparison with ±(4log ξ-)m + const with quasi-momentum ξ- = 0.584692 corresponding (66) to the leading evanescent rapidity βevan = 0.438739 (full lines).

4.2. 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, J^x_m[J^x_min, J^x_max], J^y_m[J^y_min, J^y_max] and h_m[h_min, h_max]. 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 NMM) 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 figure 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 figure 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 again indicates insulating behaviour.

Figure 6. Average Liouvillean spectral gap Δ versus the chain length n for disordered XY models: (i) Jxm = 0.5, Jym = 0 and hm[1, 2] (transverse Ising with field disorder, blue points), (ii) Jxm[0.5, 2], Jym = 0 and hm = 1 (transverse Ising with interaction disorder, red points) and (iii) Jxm[0.5, 1], Jym[0.5, 1] and hm = 1 (XY with interaction disorder, golden points), all for bath couplings ΓL1 = 1, ΓL2 = 0.6, ΓR1 = 1 and ΓR2 = 0.3. Full lines indicate exponential fits to right halves of data. Averaging is performed over 2000 disorder realizations.

Figure 7. The scaling of the energy current Qm with chain length n of the disordered XY model in the same regimes/parameters/plot styles as in figure 6.

5. Discussion and conclusions

The main result of the paper is a general method of explicit solution of master equations describing the dynamics of an 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) through 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 et al  [18] lifted to an operator space. We have explicitly constructed a non-canonical analog of Bogoliubov transformation of the quantum Liouville map to NMM. 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 problems in Heisenberg XY spin 1/2 chains which are coupled to canonical heat baths only at the two ends. Irrespective 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 with localization of NMM 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 a 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 tight-binding 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.

Finally, we mention a more ambitious extension of the present work: namely, we propose to explore the question of whether more involved algebraic methods of solution of interacting many-body quantum systems, 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?

Acknowledgments

I gratefully acknowledge stimulating discussions with Pierre Gaspard, Keiji Saito and Walter Strunz, I thank Carlos Mejia-Monasterio and Thomas H Seligman for reading the manuscript and many useful comments, and Iztok Pižorn and Marko Žnidarič for collaboration on related projects. The work has been supported by the grants P1-0044 and J1-7347 of the Slovenian research agency (ARRS). Explicit analytical calculations reported in subsection (4.1) were assisted by the Mathematica software package.

References

[1]

Zurek W H 2003 Rev. Mod. Phys. 75 715 
CrossRef

[2]

Joos E, Zeh H D, Kiefer C, Giulini D, Kupsch J and Stamatescu I-O 2003 Decoherence and the Appearance of a Classical World in Quantum Theory  (Berlin: Springer) 

[3]

von Neumann J 1955 Mathematical Foundations of Quantum Mechanics transl. R T Geyer  (Princeton: Princeton University Press) 

[4]

Schlosshauer M 2004 Rev. Mod. Phys. 76 1267 
CrossRef

[5]

Nielsen M A and Chuang I L 2000 Quantum Computation and Quantum Information  (Cambridge: Cambridge University Press) 

[6]

Benenti G, Casati G and Strini G 2004 Principles of Quantum Computation and Information Basic Concepts vol 1  (Singapore: World Scientific) 
CrossRef
Benenti G, Casati G and Strini G 2007 Principles of Quantum Computation and Information Basic Tools and Special Topics vol 2  (Singapore: World Scientific) 

[7]

Araki H and Barouch E 1983 J. Stat. Phys. 31 327 
CrossRef
Araki H 1984 Publ. Res. Inst. Math. Sci. Kyoto Univ. 20 277 
CrossRef

[8]

Ruelle D 2000 J. Stat. Phys. 98 57 
CrossRef

[9]

Jakšič V and Pillet C-A 2002 J. Stat. Phys. 108 787 
CrossRef
Jakšič V and Pillet C-A 2002 Commun. Math. Phys. 226 131 
CrossRef
Aschbacher W, Jakšič 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 

[10]

Gaspard P 2006 Prog. Theor. Phys. Suppl. 165 33 
CrossRef
Gaspard P 2006 Physica A 369 201 
CrossRef

[11]

Lee M H 2007 Acta Phys. Pol. B 38 1837 

Lee M H 2001 Phys. Rev. Lett. 87 250601 
CrossRefPubMed
Lee M H 1982 Phys. Rev. Lett. 49 1072 
CrossRef

[12]

Prosen T 2007 J. Phys. A: Math. Theor. 40 7881 
IOPscience

[13]

Korepin V E, Bogoliubov N M and Izergin A G 1997 Quantum Inverse Scattering and Correlation Functions  (Cambridge: Cambridge University Press) 

[14]

 Fadeev L, Van Moerbeke P and Lambert F (ed) 2006 Bilinear Integrable Systems: from Classical to Quantum, Continuous to Discrete (NATO ARW Proc.) (Springer Series: NATO Science Series II: Mathematics, Physics and Chemistry) vol 201 (Berlin: Springer) 
CrossRef

[15]

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

[16]

Haake F 2001 Quantum Signatures of Chaos 2nd edn  (Berlin: Springer) 

[17]

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

[18]

Lieb E H, Schultz T D and Mattis D C 1961 Ann. Phys. (NY) 16 407 
CrossRef

[19]

Zotos X, Naef F and Prelovšek P 1997 Phys. Rev. B 55 11029 
CrossRef

[20]

Saito K, Takesue S and Miyashita S 2000 Phys. Rev. E 61 2397 
CrossRef
Saito K 2003 Europhys. Lett. 61 34 
IOPscience

[21]

Michel M, Hartmann M, Gemmer J and Mahler G 2003 Eur. Phys. J. B 34 325 
CrossRef
Michel M, Mahler G and Gemmer J 2005 Phys. Rev. Lett. 95 180602 
CrossRefPubMed
Michel M, Gemmer J and Mahler G 2006 Int. J. Mod. Phys. B 20 4855 
CrossRef
Gemmer J, Steinigeweg R and Michel M 2006 Phys. Rev. B 73 104302 
CrossRef

[22]

Mejia-Monasterio C, Prosen T and Casati G 2005 Europhys. Lett. 72 520 
IOPscience
Casati G and Mejia-Monasterio C 2007 Preprint 0710.3500v1
Preprint

[23]

Mejia-Monasterio C and Wichterich H 2007 Preprint 0709.1412v1
Preprint

[24]

Dhar A and Roy D 2006 J. Stat. Phys. 125 805 (Preprint 0711.4318v1)
CrossRefPreprint

[25]

Lindblad G 1976 Commun. Math. Phys. 48 119 
CrossRef

[26]

Wolf M M, Eisert J, Cubitt T S and Cirac J I 2007 Preprint 0711.3172v1
Preprint

[27]

Šemrl P 2008 private communication

[28]

Moshinsky M and Seligman T H 1971 Ann. Phys., NY 66 311 
CrossRef

[29]

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

[30]

Rieder Z, Lebowitz J L and Lieb E 1967 J. Math. Phys. 8 1073 
CrossRef

[31]

Prosen T 2000 Prog. Theor. Phys. Suppl. 139 191 
CrossRef

[32]

Lepri S, Livi R and Politi A 2003 Phys. Rep. 377 1 
CrossRef

[33]

Prosen T 1999 Phys. Rev. E 60 1658 
CrossRef

[34]

Prosen T and Pižorn I 2007 Phys. Rev. A 76 032316 
CrossRef

[35]

Esposito M and Gaspard P 2005 Phys. Rev. B 71 214302 
CrossRef
Esposito M and Gaspard P 2005 J. Stat. Phys. 121 463 
CrossRef

[36]

Vidal G, Latorre J I, Rico E and Kitaev A 2003 Phys. Rev. Lett. 90 227902 
CrossRefPubMed
Latorre J I, Rico E and Vidal G 2004 Quantum Inf. Comput. 4 48 

[37]

Prosen T and Žnidarič M 2008 to be submitted

[38]

White S R 1992 Phys. Rev. Lett. 69 2863 
CrossRefPubMed
Schollwöck U and White S R 2006 Effective Models for Low-Dimensional Strongly Correlated Systems ed G G Batrouni and D Poilblanc (Melville, NY: AIP) pp 155 

Vidal G 2003 Phys. Rev. Lett. 91 147902 
CrossRefPubMed
Vidal G 2004 Phys. Rev. Lett. 93 040502 
CrossRefPubMed

Notes

Note1  Throughout this paper Dirac's bra-ket notation shall be used only for a Hilbert space of physical operators, including density operators, in a sense of Gelfand–Naimark–Segal construction, although here all spaces will be finite dimensional. Symbols with a hat shall designate linear maps over the operator space . For instance, we note a key distinction between physical fermions c_m (5) and a-fermions (9).

Note2  It is not known at present whether the explicit 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].

Note3  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.

Note4  We note a similarity to the formalism of second quantization with non-orthogonal orbitals introduced in [28].

Note5  Trivial zero β = 0 of course does not represent a physical solution since then the whole S-matrix (70) vanishes.