You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: 'The Lindblad master equation for an arbitrary [*quadratic*]{} system of $n$ fermions is solved explicitly in terms of diagonalization of a $4n \times 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.'
address: 'Department of physics, FMF, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia'
author:
- Tomaž Prosen
title: 'Third quantization: a general method to solve master equations for quadratic open Fermi systems'
---

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 [@zurek; @zeh] and closely related quantum measurement problem [@neumann; @schlos], quantum computation [@chuang; @casati], or the problem of computation of [*non-equilibrium steady states*]{} (NESS) in quantum statistical mechanics [@araki; @ruelle; @jaksic]. Even though application of the methods of Hamiltonian dynamical systems and ergodic theory to quantum systems out of equilibrium gives many promising results [@piere2; @lee; @prosenjpa], 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 [@korepin; @fadeev]. 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. [@petruccione; @haake; @alicki]).

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 [@lieb] 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\times 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. [@zotos; @saito; @hartmann; @mejia05; @mejia07; @dhar; @prosenjpa] 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 [@petruccione], 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 \[sect:method\] 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 \[sect:trivia\] 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 \[sect:nontrivia\] we demonstrate the usefulness of the new method by applying it to quantum transport in XY spin chains. In section \[sect:conc\] we discuss possible alternative applications and generalizations of the method and reach some conclusions.

General method of solution for the Lindblad equation {#sect:method}
====================================================

The general master equation governing time evolution of the density matrix $\rho(t)$ of an open quantum system, preserving trace and positivity of $\rho$, can be written in the Lindblad form [@lindblad; @alicki] as (we set $\hbar=1$) $$\frac{{ {\rm d} }\rho }{{ {\rm d} }t} = {{\hat {\cal L}}}\rho :=
-{ {\rm i} }[H,\rho] + \sum_{\mu} \left(2 L_\mu \rho L_\mu^\dagger - \{L_\mu^\dagger L_\mu,\rho\} \right)
\label{eq:lind}$$ where $H$ is a Hermitian operator (Hamiltonian), $[x,y]:=xy-yx$, $\{x,y\}:=xy+yx$, and $L_\mu$ 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 (\[eq:lind\]) for a [*quadratic*]{} system of $n$ fermions (or spins $1/2$) with [*linear*]{} bath operators $$\begin{aligned}
H &=& \sum_{j,k=1}^{2n} w_j H_{jk} w_k = {{\underline{w}}} \cdot {{\mathbf{H}}}\, {{\underline{w}}} \label{eq:hamil}  \\
L_\mu &=& \sum_{j=1}^{2n} l_{\mu,j} w_j = {{\underline{l}}}_\mu \cdot {{\underline{w}}}\  \label{eq:lindb}\end{aligned}$$ where $w_j$, $j=1,2,\ldots,2n$, are abstract [*Hermitian*]{} Majorana operators satisfying the anti-commutation relations $$\{w_j,w_k\} = 2\delta_{j,k} \qquad j,k =1,2,\ldots, 2n$$ and generate a Clifford algebra. Thus, $2n \times 2n$ matrix ${{\mathbf{H}}}$ can be chosen to be antisymmetric ${{\mathbf{H}}}^T = -{{\mathbf{H}}}$. Throughout this paper ${{\underline{x}}}=(x_1,x_2,\ldots)^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,\ldots, n$, $$w_{2m-1}= c_m + c_m^\dagger \qquad
w_{2m} = { {\rm i} }(c_m- c^\dagger_m)
\label{eq:physfermi}$$ or (ii) spins $1/2$ with canonical Pauli operators $\vec{\sigma}_m$, $m=1,2,\ldots, n$, $$w_{2m-1} = \sigma^{{\rm x}}_m \prod_{m'<m} \sigma^{{\rm z}}_{m'} \qquad
w_{2m} = \sigma^{{\rm y}}_m \prod_{m'<m} \sigma^{{\rm z}}_{m'}
\label{eq:jordan}$$

Here we are not concerned with physical criteria for the validity of the so-called Markovian approximation under which eq. (\[eq:lind\]) is derived, so we shall make no assumptions on the smallness of the bath coupling constants $l_{\mu,j}$. We merely consider the Lindblad equation (\[eq:lind\]) 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_\mu(t)$, generating more general and possibly non-Markovian open system dynamics, is straightforward. See e.g. [@wolf] for a discussion of [*Markovianity*]{}.

Fock space of operators
-----------------------

We begin by associating a Hilbert space structure $x \to {{\vert x \rangle}}$ to a linear $2^{2n}=4^n$ dimensional space ${\cal K}$ of operators, with a canonical basis ${{\vert P_{{{\underline{\alpha}}}} \rangle}}$ with $$P_{\alpha_1,\alpha_2,\ldots,\alpha_{2n}} := w_1^{\alpha_1}w_2^{\alpha_2}\cdots w_{2n}^{\alpha_{2n}} \qquad
\alpha_j\in\{ 0,1\}$$ [*orthonormal*]{} with respect to an inner product $${\langle x \vert y \rangle} = 2^{-n} \tr x^\dagger y$$ The form of the canonical basis of the operator space ${\cal 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*]{} $\hat{c}_j$ over ${\cal K}$ $${ {\hat c} }_j {{\vert P_{{{\underline{\alpha}}}} \rangle}} = \delta_{\alpha_j,1} {{\vert w_j P_{{{\underline{\alpha}}}} \rangle}}
\label{eq:defanih}$$ and derive the actions of their Hermitian adjoints - the [*creation linear maps*]{} $\hat{c}^\dagger$, ${{\langle P_{{{\underline{\alpha}}}'} \vert}}{ {\hat c} }^\dagger_j{{\vert P_{{{\underline{\alpha}}}} \rangle}} =
{{\langle P_{{{\underline{\alpha}}}} \vert}}{ {\hat c} }_j{{\vert P_{{{\underline{\alpha}}}'} \rangle}}^* =
\delta_{\alpha'_j,1} {\langle P_{{{\underline{\alpha}}}} \vert w_j P_{{{\underline{\alpha}}}'} \rangle}^*=
\delta_{\alpha_j,0} {\langle P_{{{\underline{\alpha}}}'} \vert w_j P_{{{\underline{\alpha}}}} \rangle} $: $${ {\hat c} }^\dagger_j {{\vert P_{{{\underline{\alpha}}}} \rangle}} = \delta_{\alpha_j,0} {{\vert w_j P_{{{\underline{\alpha}}}} \rangle}}$$ Straightforward inspection then shows that they satisfy the canonical anticommutation relations $$\{{ {\hat c} }_j,{ {\hat c} }_k\} = 0 \qquad \{{ {\hat c} }_j,{ {\hat c} }_k^\dagger\} = \delta_{j,k} \qquad j,k=1,2,\ldots, 2n$$ The key is now to realize that the quantum Liouville map ${{\hat {\cal L}}}$ defined by eqs. (\[eq:lind\],\[eq:hamil\],\[eq:lindb\]) can be written as a quadratic form in [*adjoint Fermi maps*]{} ${ {\hat c} }_j,{ {\hat c} }^\dagger_j$ (or for short, [*a-fermions*]{}). [^1]\
\
First, we consider the unitary part of Liouvillean $${{\hat {\cal L}}}_0 \rho := -{ {\rm i} }[H,\rho]
\label{eq:exprunitary}$$ Since ${\cal K}$ is a Lie algebra, one defines the [*adjoint representation*]{} of a Lie derivative for an arbitrary $A\in {\cal K}$ back on ${\cal K}$ as, $\ad A {{\vert B \rangle}} := {{\vert [A,B] \rangle}}$. It is now straightforward to compute the action of a Lie derivative of a product of two Majorana operators on an arbitrary basis element $$\begin{aligned}
\ad w_j w_k {{\vert P_{{\underline{\alpha}}} \rangle}} &=& {{\vert w_j w_k P_{{\underline{\alpha}}} \rangle}} - {{\vert P_{{\underline{\alpha}}} w_j w_k \rangle}} = \nonumber \\
&=& 2 (\delta_{\alpha_j,1}\delta_{\alpha_k,0} + \delta_{\alpha_j,0}\delta_{\alpha_k,1}){{\vert w_j w_k P_{{\underline{\alpha}}} \rangle}} = \nonumber \\
&=& 2 ({ {\hat c} }_j { {\hat c} }^\dagger_k + { {\hat c} }^\dagger_j { {\hat c} }_k) {{\vert P_{{\underline{\alpha}}} \rangle}} =
2 ({ {\hat c} }^\dagger_j { {\hat c} }_k - { {\hat c} }^\dagger_k{ {\hat c} }_j) {{\vert P_{{\underline{\alpha}}} \rangle}}\end{aligned}$$ Extending this relation by linearity to an arbitrary element of ${\cal K}$, it follows that for an arbitrary quadratic Hamiltonian (\[eq:hamil\]) its Lie derivative has a very similar quadratic form in a-Fermi maps $${{\hat {\cal L}}}_0 = -{ {\rm i} }\ad H = -4{ {\rm i} }\sum_{j,k=1}^{2n} { {\hat c} }^\dagger_j H_{jk} { {\hat c} }_k =
-4{ {\rm i} }\, {{\underline{{ {\hat c} }}}}^\dagger\cdot {{\mathbf{H}}}\,{{\underline{{ {\hat c} }}}}
\label{eq:unitary}$$ It is worth stressing here that for an arbitrary (complex) matrix ${{\mathbf{H}}}$, ${{\hat {\cal L}}}_0$ (\[eq:unitary\]) conserves the total number of a-fermions ${{\hat {\cal N}}}:= \sum_{j} { {\hat c} }^\dagger_j { {\hat c} }_j = {{\underline{{ {\hat c} }}}}^\dagger\cdot{{\underline{{ {\hat c} }}}}$, namely $[{{\hat {\cal L}}}_0,{{\hat {\cal N}}}] = 0$.\
\
Second, we consider the action of the Lindblad maps $${{\hat {\cal L}}}_\mu \rho := 2 L_\mu \rho L_\mu^\dagger - \{L_\mu^\dagger L_\mu,\rho\} =
\sum_{j,k=1}^{2n} l_{\mu,j} l_{\mu,k}^* {{\hat {\cal L}}}_{j,k}\rho
\label{eq:exprlindb}$$ where we write ${{\hat {\cal L}}}_{j,k}\rho := 2 w_j \rho w_k - w_k w_j \rho - \rho w_k w_j$. Again we proceed by computing the actions of ${{\hat {\cal L}}}_{j,k}$ on elements of the canonical basis of operator Fock space ${\cal K}$. In order to do so, it is crucial to observe that the question whether $w_j$ commutes or anticommutes with $P_{{{\underline{\alpha}}}}$ depends on the number of a-fermions $|{{\underline{\alpha}}}|:=\sum_{k=1}^{2n} \alpha_k$ in ${{\vert P_{{{\underline{\alpha}}}} \rangle}}$, namely $P_{{{\underline{\alpha}}}}w_j = (-1)^{|{{\underline{\alpha}}}|+\alpha_j} w_j P_{{{\underline{\alpha}}}}$, and hence $${{\hat {\cal L}}}_{j,k} {{\vert P_{{{\underline{\alpha}}}} \rangle}}
= \left[2(-1)^{|{{\underline{\alpha}}}|+\alpha_k} w_j w_k - w_k w_j - (-1)^{\alpha_j+\alpha_k} w_k w_j\right]
{{\vert P_{{{\underline{\alpha}}}} \rangle}}
\label{eq:Ljk1}$$ Observing that $$\begin{aligned}
\phantom{(-1)^{\alpha_j}}{{\vert w_j P_{{{\underline{\alpha}}}} \rangle}} &=& ({ {\hat c} }^\dagger_j + { {\hat c} }_j){{\vert P_{{{\underline{\alpha}}}} \rangle}} \\
(-1)^{\alpha_j} {{\vert w_j P_{{{\underline{\alpha}}}} \rangle}} &=& ({ {\hat c} }^\dagger_j - { {\hat c} }_j){{\vert P_{{{\underline{\alpha}}}} \rangle}} \\
(-1)^{|{{\underline{\alpha}}}|} {{\vert P_{{{\underline{\alpha}}}} \rangle}} &=& \exp({ {\rm i} }\pi{{\hat {\cal N}}}) {{\vert P_{{{\underline{\alpha}}}} \rangle}}\end{aligned}$$ one derives from (\[eq:Ljk1\]) the general expression for ${{\hat {\cal L}}}_{j,k}$ $$\begin{aligned}
{{\hat {\cal L}}}_{j,k} &=& \left(\hat{\one}+\exp({ {\rm i} }\pi{{\hat {\cal N}}})\right)\left(2 { {\hat c} }^\dagger_j { {\hat c} }^\dagger_k - { {\hat c} }^\dagger_j { {\hat c} }_k -
{ {\hat c} }^\dagger_k { {\hat c} }_j\right) \nonumber \\
&+& \left(\hat{\one}-\exp({ {\rm i} }\pi{{\hat {\cal N}}})\right)\left(2 { {\hat c} }_j { {\hat c} }_k - { {\hat c} }_j { {\hat c} }^\dagger_k -
{ {\hat c} }_k { {\hat c} }^\dagger_j\right)\end{aligned}$$ Obviously, the maps ${{\hat {\cal L}}}_{j,k}$, and hence also the total Lindblad part of Liouvillean $\sum_\mu{{\hat {\cal L}}}_\mu$, 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 ${{\hat {\cal P}}}= \exp({ {\rm i} }\pi{{\hat {\cal N}}})$, with respect to which the operator space can be decomposed into a direct sum ${\cal K} = {\cal K}^+ \oplus {\cal K}^-$, and even/odd operator spaces are orthogonally projected as ${\cal K}^{\pm} = (\hat{\one} \pm \exp({ {\rm i} }\pi{{\hat {\cal N}}})){\cal K}$. Thus the positive parity subspace ${\cal K}^+$ is a linear space spanned by ${{\vert P_{{{\underline{\alpha}}}} \rangle}}$ with [*even*]{} $|{{\underline{\alpha}}}|$. All the maps ${{\hat {\cal L}}}_{j,k}$ now act separately on ${\cal K}^\pm$, ${{\hat {\cal L}}}_{j,k} {\cal K}^\pm
\subseteq {\cal K}^\pm$. For example, the maps defined on even parity subspace are indeed quadratic in a-fermions $${{\hat {\cal L}}}_{j,k}\vert_{{\cal K}^+} = 4 { {\hat c} }^\dagger_j { {\hat c} }^\dagger_k - 2 { {\hat c} }^\dagger_j { {\hat c} }_k -
2 { {\hat c} }^\dagger_k { {\hat c} }_j
\label{eq:lindbpart}$$ 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 ${\cal K}^+$. The extension to the dynamics of [*odd*]{} observables should be straightforward.

Putting the results (\[eq:exprunitary\],\[eq:unitary\],\[eq:exprlindb\],\[eq:lindbpart\]) together we arrive at the final compact quadratic representation of the Liouvillean ${{\hat {\cal L}}}_+ := {{\hat {\cal L}}}\vert_{{\cal K}^+}$ $${{\hat {\cal L}}}_+ = -2\, {{\underline{{ {\hat c} }}}}^\dagger\cdot(2{ {\rm i} }{{\mathbf{H}}} + {{\mathbf{M}}} + {{\mathbf{M}}}^T)\,{{\underline{{ {\hat c} }}}}
+ 2\,{{\underline{{ {\hat c} }}}}^\dagger\cdot ({{\mathbf{M}}}-{{\mathbf{M}}}^T)\,{{\underline{{ {\hat c} }}}}^\dagger
\label{eq:liouv1}$$ where ${{\mathbf{M}}}$ is a complex Hermitian matrix parametrizing the Lindblad operators $$M_{jk} = \sum_{\mu} l_{\mu,j} l_{\mu,k}^*$$

Reduction to normal master modes
--------------------------------

Next we want to show that the representation (\[eq:liouv1\]) allows us to reduce it further by a linear transformation of the set of maps $\{{ {\hat c} }_j,{ {\hat c} }_j^\dagger; j=1,2,\ldots,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.*]{} [@lieb] and define $4n$ adjoint Hermitian Majorana maps ${ {\hat a} }_r = { {\hat a} }_r^\dagger$, $r=1,2,\ldots,4n$: $${ {\hat a} }_{2j-1} = \frac{1}{\sqrt{2}}({ {\hat c} }_j + { {\hat c} }^\dagger_j) \qquad
{ {\hat a} }_{2j} = \frac{{ {\rm i} }}{\sqrt{2}}({ {\hat c} }_j - { {\hat c} }^\dagger_j)
\label{eq:aexpr}$$ satisfying the anti-commutation relations $$\{{ {\hat a} }_r,{ {\hat a} }_s\} = \delta_{r,s}
\label{eq:clifa}$$ in terms of which the Liouvillean (\[eq:liouv1\]) can be rewritten as $${{\hat {\cal L}}}_+ = {{\underline{{ {\hat a} }}}}\cdot{{\mathbf{A}}}\,{{\underline{{ {\hat a} }}}} - A_0 \hat{\one}
\label{eq:liouv2}$$ where ${{\mathbf{A}}}$ is an antisymmetric complex $4n\times 4n$ matrix with entries $$\begin{aligned}
A_{2j-1,2k-1} &=&-2{ {\rm i} }H_{jk}-M_{jk}+M_{kj} \nonumber \\
A_{2j-1,2k}   &=&  \phantom{wj}2{ {\rm i} }M_{kj} \nonumber \\
A_{2j,2k-1}   &=& -2{ {\rm i} }M_{jk} \nonumber \\
A_{2j,2k}     &=&-2{ {\rm i} }H_{jk}+M_{jk}-M_{kj}
\label{eq:explA}\end{aligned}$$ $\hat{\one}$ is an identity map over ${\cal K}$ and $A_0$ is a scalar $$A_0 = 2\sum_{j=1}^{2n} M_{jj} = 2\tr{{\mathbf{M}}}$$ Obviously, the [*bi-linear*]{} Liouvillean (\[eq:liouv2\]) cannot be brought to a normal form with a linear [*canonical*]{} transformation since the matrix ${{\mathbf{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 ${{\underline{v}}}$ is a [*right*]{} eigenvector ${{\mathbf{A}}}{{\underline{v}}} = \beta {{\underline{v}}}$ with complex eigenvalue $\beta$, then ${{\underline{v}}}$ is also a [*left*]{} eigenvector with eigenvalue $-\beta$, ${{\mathbf{A}}}^T {{\underline{v}}} = -{{\mathbf{A}}} {{\underline{v}}} = -\beta {{\underline{v}}}$. Hence eigenvalues always come in pairs $\beta,-\beta$. Let as assume that ${{\mathbf{A}}}$ can be [*diagonalized*]{}[^2], i.e. there exist $4n$ linearly independent vectors ${{\underline{v}}}_r,r=1,\ldots,4n$ with the corresponding eigenvalues $\beta_1,-\beta_1,\beta_2,-\beta_2,\ldots,\beta_{2n},-\beta_{2n}$, $${{\mathbf{A}}}{{\underline{v}}}_{2j-1} = \beta_j {{\underline{v}}}_{2j-1} \qquad
{{\mathbf{A}}}{{\underline{v}}}_{2j} = -\beta_j {{\underline{v}}}_{2j}
\label{eq:eigv}$$ ordered such that $\re\beta_1\ge \re\beta_2\ge \ldots \ge \re\beta_{2n} \ge 0.$ The $2n$ complex numbers $\beta_j$ shall be referred to as [*rapidities*]{}. It is easy to check that we can always choose and normalize ${{\underline{v}}}_r$ such that[^3] $${{\underline{v}}}_r \cdot {{\underline{v}}}_s = J_{rs} \quad \textrm{where} \quad
{{\mathbf{J}}} := \sigma^{{{\rm x}}}\otimes {{\mathbf{I}}}_{2n} =
\pmatrix{0 & 1 & 0 & 0 & \cdots \cr
                 1 & 0 & 0 & 0 & \cdots \cr
                 0 & 0 & 0 & 1 & \cdots \cr
                 0 & 0 & 1 & 0 & \cdots \cr
\vdots &\vdots & \vdots & \vdots& \ddots \cr}
\label{eq:norm}$$ Let ${{\mathbf{V}}}$ be $4n\times 4n$ matrix whose $r$th row is given by ${{\underline{v}}}_r$, $V_{rs} := v_{r,s}$. Then eqs. (\[eq:eigv\],\[eq:norm\]) rewrite as $$\begin{aligned}
{{\mathbf{A}}}{{\mathbf{V}}}^T &=& {{\mathbf{V}}}^T{{\mathbf{D}}}\quad \textrm{where}\quad {{\mathbf{D}}} := {\rm diag}\{\beta_1,-\beta_1,
\beta_2,-\beta_2,\ldots,\beta_{2n},-\beta_{2n}\} \label{eq:AVVD} \\
{{\mathbf{V}}} {{\mathbf{V}}}^T &=& {{\mathbf{J}}}
\label{eq:VVJ}\end{aligned}$$ Expressing ${{\mathbf{V}}}^T$ in terms of (\[eq:VVJ\]) and plugging the result into eq. (\[eq:AVVD\]) we arrive at a very convenient canonical form of a generic complex antisymmetric matrix ${{\mathbf{A}}}$ $${{\mathbf{A}}} = {{\mathbf{V}}}^T {{\mathbf{\Lambda}}} {{\mathbf{V}}}\quad\textrm{where}\quad
{{\mathbf{\Lambda}}} = {{\mathbf{D}}}{{\mathbf{J}}} =
\pmatrix{0 & \beta_1 & 0 & 0 & \cdots \cr
                 -\beta_1 & 0 & 0 & 0 & \cdots \cr
                 0 & 0 & 0 & \beta_2 & \cdots \cr
                 0 & 0 & -\beta_2 & 0 & \cdots \cr
\vdots &\vdots & \vdots & \vdots& \ddots \cr}
\label{eq:canform}$$

Now we apply decomposition (\[eq:canform\]) to the Liouvillean (\[eq:liouv2\]) $${{\hat {\cal L}}}_+ = {{\underline{{ {\hat a} }}}}\cdot {{\mathbf{V}}}^T {{\mathbf{\Lambda}}}{{\mathbf{V}}}{{\underline{{ {\hat a} }}}}  - A_0\hat{\one} =
({{\mathbf{V}}}{{\underline{{ {\hat a} }}}})\cdot{{\mathbf{\Lambda}}}({{\mathbf{V}}}{{\underline{{ {\hat a} }}}}) - A_0\hat{\one}
\label{eq:liouv4}$$ Let us define the NMM maps ${{\underline{{ {\hat b} }}}} := ({ {\hat b} }_1,{ {\hat b} }'_1,{ {\hat b} }_2,{ {\hat b} }'_2,\ldots,{ {\hat b} }_{2n},{ {\hat b} }'_{2n}) :=
{{\mathbf{V}}}{{\underline{{ {\hat a} }}}}$ or $${ {\hat b} }_j = {{\underline{v}}}_{2j-1}\cdot{{\underline{{ {\hat a} }}}} \qquad
{ {\hat b} }'_j = {{\underline{v}}}_{2j}\cdot{{\underline{{ {\hat a} }}}}
\label{eq:bexpr}$$ We note that due to (\[eq:clifa\],\[eq:norm\]) NMM maps satisfy [*almost-canonical*]{} anti-commutation relations $$\{{ {\hat b} }_j,{ {\hat b} }_k\} = 0 \qquad \{{ {\hat b} }_j,{ {\hat b} }'_k\} = \delta_{j,k} \qquad \{{ {\hat b} }'_j,{ {\hat b} }'_k\}=0
\label{acar}$$ namely ${ {\hat b} }_j$ could be interpreted as annihilation map and ${ {\hat b} }'_j$ as a creation map of $j$th NMM, but we should note that ${ {\hat b} }'_j$ is in general [*not*]{} the Hermitian adjoint of ${ {\hat b} }_j$ [@thomas]. In terms of NMM the Liouvillean (\[eq:liouv4\]) now achieves a very convenient normal form $${{\hat {\cal L}}}_+ = -2\sum_{j=1}^{2n} \beta_j { {\hat b} }'_j { {\hat b} }_j - B_0 \hat{\one}
\label{eq:liouv3}$$ where $B_0 = A_0 - \sum_{j=1}^{2n}\beta_j$. 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 \times 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 \times 4n$ matrix ${{\mathbf{A}}}$ as described in the previous subsection.

We proceed by constructing the Liouvillean ‘vacuum’. From the representation (\[eq:liouv1\]) it follows immediately that ${{\langle 1 \vert}} = {{\langle P_{(0,0\ldots,0)} \vert}}$ is left-annihilated by ${{\hat {\cal L}}}_+$, ${{\langle 1 \vert}}{{\hat {\cal L}}}_+ = 0$, or equivalently ${{\hat {\cal L}}}_+^\dagger {{\vert 1 \rangle}} = 0$. So we have just shown that $0$ is always an eigenvalue of ${{\hat {\cal L}}}_+$, hence there should also exist the corresponding right eigenvector $\ness$, normalized as ${\langle 1 \vert {\rm NESS} \rangle} = \tr \rho_{\rm NESS} = 1$, which represents physical NESS, i.e. stationary solutions of the Lindblad equation (\[eq:lind\]) $${{\hat {\cal L}}}_+\ness = 0$$ Let us define NMM number maps as ${{\hat {\cal N}}}_j := { {\hat b} }'_j { {\hat b} }_j$. From eqs. (\[acar\]) it follows that ${{\hat {\cal N}}}_j$ satisfy a projection property ${{\hat {\cal N}}}_j^2 = {{\hat {\cal N}}}_j$, so they are diagonalizable since no nontrivial Jordan block could satisfy the projection property. Furthermore, ${{\hat {\cal N}}}_j$ are mutually commuting $[{{\hat {\cal N}}}_j,{{\hat {\cal N}}}_k]=0$, so they can be simultaneously diagonalized and there should exist a vacuum state on which all ${{\hat {\cal N}}}_j$ have value 0. It follows from the stability of completely positive evolution (\[eq:lind\]) that all eigenvalues $\lambda$ of ${{\hat {\cal L}}}_+$ should obey $\re\lambda \le 0$, and since by assumption $\re\beta_j \ge 0$, ${{\langle 1 \vert}}$ and $\ness$ should be the left and right vacua which are simultaneously annihilated by NMM maps $${{\langle 1 \vert}} { {\hat b} }'_j = 0 \qquad { {\hat b} }_j \ness = 0
\label{eq:anihb}$$ and hence also ${{\hat {\cal N}}}_j\ness = 0$. Thus we have also shown that the NMM representation (\[eq:liouv3\]) is only consistent if $B_0=0$ so we find an interesting sum rule for rapidities $$\sum_{j=1}^{2n} \beta_j = 2\tr{{\mathbf{M}}}$$

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) ${{\underline{\nu}}}=(\nu_1,\nu_2,\ldots,\nu_{2n})$, $\nu_j\in\{0,1\}$, $${{\langle \Theta^{\rm L}_{{\underline{\nu}}} \vert}}{{\hat {\cal L}}}_+ = \lambda_{{{\underline{\nu}}}} {{\langle \Theta^{\rm L}_{{\underline{\nu}}} \vert}}
\qquad
{{\hat {\cal L}}}_+ {{\vert \Theta^{\rm R}_{{\underline{\nu}}} \rangle}} = \lambda_{{{\underline{\nu}}}} {{\vert \Theta^{\rm R}_{{\underline{\nu}}} \rangle}}$$ namely $$\begin{aligned}
\lambda_{{\underline{\nu}}} &=& -2\sum_{j=1}^{2n} \beta_j \nu_j \label{eq:eval}\\
{{\langle \Theta^{\rm L}_{{\underline{\nu}}} \vert}} &=& {{\langle 1 \vert}}{{ {\hat b} }_{2n}}^{\nu_{2n}}\cdots {{ {\hat b} }_2}^{\nu_2}\,{{ {\hat b} }_1}^{\nu_1}
\qquad
{{\vert \Theta^{\rm R}_{{\underline{\nu}}} \rangle}} = {{ {\hat b} }_1}^{'\nu_1}\,{{ {\hat b} }_2}^{'\nu_2} \cdots {{ {\hat b} }_{2n}}^{' \nu_{2n}} \ness \label{eq:evec}\end{aligned}$$ where by construction, left and right eigenvectors satisfy the bi-orthonormality relation ${\langle \Theta^{\rm L}_{{\underline{\nu'}}} \vert \Theta^{\rm R}_{{\underline{\nu}}} \rangle} = \delta_{{{\underline{\nu}}}',{{\underline{\nu}}}}$.

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

Given a [*physical observable*]{} $X\in {\cal K}^+$ and an [*arbitrary*]{} initial state with a density operator $\rho_0 \in {\cal K}$, the time dependent expectation value of $X$ can be written in terms of the spectral resolution of the Liouvillean, $$\exp(t {{\hat {\cal L}}}_+) = \sum_{{{\underline{\nu}}}} \exp(t \lambda_{{{\underline{\nu}}}})
{{\vert \Theta^{\rm R}_{{{\underline{\nu}}}} \rangle}}{{\langle \Theta^{\rm L}_{{{\underline{\nu}}}} \vert}}
\label{eq:propag}$$ namely $${{\langle X(t)\rangle}} = \tr X\rho(t) = \tr\!\!\left[X\exp(t{{\hat {\cal L}}}_+)\rho_0\right] =
\sum_{{{\underline{\nu}}}} \exp(t\lambda_{{{\underline{\nu}}}}){{\langle \Theta^{\rm L}_{{{\underline{\nu}}}} \vert}}\rho_0 X{{\vert \Theta^{\rm R}_{{{\underline{\nu}}}} \rangle}}
\label{eq:timedeptexp}$$ We remind the reader that ${{\hat {\cal L}}}_+$ correctly represents physical Liouvillean only on the subspace ${\cal K}^+$ of operators with [*even*]{} number of a-fermions. However, since the dynamics is closed on ${\cal K}^+$ and test physical observable $X$ also belongs to ${\cal K}^+$ it follows that the component of $\rho_0$ from ${\cal K}^-$ does not contribute to the expectation value ${{\langle X(t)\rangle}}$.

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 ${{\mathbf{A}}}$ is diagonalizable:\
\
[**Theorem 1:**]{} NESS of Lindblad equation (\[eq:lind\]) is [*unique*]{} if and only if the rapidity spectrum $\{\beta_j\}$ does not contain $0$, in our ordering convention, if $\beta_{2n} \ne 0$. In the opposite case, if we have $d \ge 1$ vanishing rapidities, then we have a $2^d$ dimensional convex set of non-equilibrium steady states which can be spanned with ${{\vert \Theta^{\rm R}_{(0,\ldots,0,\nu_{1},\ldots,\nu_{d})} \rangle}}$.\
\
[**Theorem 2:**]{} An arbitrary initial state with a density operator $\rho_0 \in {\cal K}$ converges with time to NESS if and only if all rapidities have [*strictly positive*]{} real parts, $\re \beta_j > 0$. Then, the rate of exponential relaxation to NESS is given by the [*spectral gap*]{} $\Delta$ of the Liouvillean which equals $\Delta=2\re\beta_{2n}$.\
\
[**Theorem 3:**]{} Assume that the rapidity spectrum does not contain 0, i.e. $\beta_{2n}\ne 0$. Then the expectation value of any quadratic observable $w_j w_k$ in a (unique) NESS can be explicitly computed as $$\begin{aligned}
{{\langle w_j w_k\rangle}}_{\rm NESS} &=& \delta_{j,k} + {{\langle 1 \vert}}{ {\hat c} }_j { {\hat c} }_k \ness = \label{eq:pairwise} \\
&=& \delta_{j,k} + \frac{1}{2}\sum_{m=1}^{2n} \biggl(
v_{2m,2j-1} v_{2m-1,2k-1} - v_{2m,2j} v_{2m-1,2k}  \nonumber \\
&& \qquad\qquad- { {\rm i} }v_{2m,2j} v_{2m-1,2k-1} - { {\rm i} }v_{2m,2j-1} v_{2m-1,2k}\biggr)
\label{eq:pairwise2}\end{aligned}$$

The statements of theorems 1 and 2 simply follow from exact and explicit spectral decomposition (\[eq:eval\],\[eq:evec\],\[eq:propag\]).

The proof of theorem 3 is also straightforward: The first expression (\[eq:pairwise\]) follows from the definition of the annihilation maps (\[eq:defanih\]) and the explicit representation of the density operator of NESS, $\rho_{\rm NESS}$, in the canonical basis $P_{{{\underline{\alpha}}}}$. The second, very useful equality (\[eq:pairwise2\]) is then obtained by expressing ${ {\hat c} }_j$ thru (\[eq:aexpr\]) in terms of NMM maps (\[eq:bexpr\]) and using the annihilation relations (\[eq:anihb\]).

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 ${{\langle P_{{{\underline{\alpha}}}}\rangle}}_{\rm NESS} =
{{\langle 1 \vert}}{ {\hat c} }^{\alpha_1}_{1}{ {\hat c} }^{\alpha_2}_2\cdots{ {\hat c} }^{\alpha_{2n}}_{2n}\ness
$, then one may use a [*Wick theorem*]{} and rewrite it as a sum of products of pair-wise contractions (\[eq:pairwise\]).

Trivial example: A single fermion in a bath {#sect:trivia}
===========================================

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 = -{ {\rm i} }h w_1 w_2 + {\rm const} =
2h c^\dagger c + {\rm const}'$ and the following bath operators (see e.g. [@haake; @wich]) $$L_1 = \frac{1}{2}\sqrt{\Gamma_1} (w_1 - { {\rm i} }w_2)  =
\sqrt{\Gamma_1} c \qquad
L_2 = \frac{1}{2}\sqrt{\Gamma_2} (w_1 + { {\rm i} }w_2) =
\sqrt{\Gamma_2} c^\dagger
\label{eq:canlind}$$ where the ratio of coupling constants determine the bath temperature $T$, $\Gamma_2/\Gamma_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 (\[eq:liouv2\]) $${{\mathbf{A}}} = -h{{\mathbf{R}}} + {{\mathbf{B}}}_{\Gamma_+,\Gamma_-} \qquad A_0 = \Gamma_+$$ where $${{\mathbf{R}}} := \pmatrix{ 0 & 0 & 1 & 0 \cr
                  0 & 0 & 0 & 1 \cr
                 \!\!-1 & 0 & 0 & 0 \cr
                  0 &\!\!\!-1 & 0 & 0}
\quad
{{\mathbf{B}}}_{\Gamma_+,\Gamma_-} :=
\pmatrix{
0 & \frac{{ {\rm i} }}{2} \Gamma_+ &\!\!\!-\frac{{ {\rm i} }}{2}\Gamma_- & \frac{1}{2}\Gamma_- \cr
\!\!-\frac{{ {\rm i} }}{2}\Gamma_+ & 0 & \frac{1}{2}\Gamma_- & \frac{{ {\rm i} }}{2}\Gamma_- \cr
\frac{{ {\rm i} }}{2}\Gamma_- &\!\!\!-\frac{1}{2}\Gamma_- & 0 & \frac{{ {\rm i} }}{2}\Gamma_+ \cr
\!\!-\frac{1}{2}\Gamma_- &\!\!\!-\frac{{ {\rm i} }}{2}\Gamma_- &\!\!\!-\frac{{ {\rm i} }}{2}\Gamma_+ & 0}
\label{eq:canbath}$$ and $\Gamma_\pm := \Gamma_2 \pm \Gamma_1$. Further, we compute NMM rapidities $\beta_{1,2} = \frac{1}{2}\Gamma_+ \pm { {\rm i} }h$ and the eigenvector matrix $${{\mathbf{V}}} = \pmatrix{
\frac{\Gamma_-}{\Gamma_+} - 1 & { {\rm i} }\frac{\Gamma_-}{\Gamma_+} + { {\rm i} }&
 -{ {\rm i} }\frac{\Gamma_-}{\Gamma_+} + { {\rm i} }&  \frac{\Gamma_-}{\Gamma_+} + 1 \cr
 -\frac{1}{4} & -\frac{{ {\rm i} }}{4} & -\frac{{ {\rm i} }}{4} & \frac{1}{4} \cr
 \frac{\Gamma_-}{\Gamma_+} + 1 &  { {\rm i} }\frac{\Gamma_-}{\Gamma_+} - { {\rm i} }&
  { {\rm i} }\frac{\Gamma_-}{\Gamma_+} + { {\rm i} }&  -\frac{\Gamma_-}{\Gamma_+} + 1 \cr
  \frac{1}{4} & \frac{{ {\rm i} }}{4} & -\frac{{ {\rm i} }}{4} & \frac{1}{4}}\qquad$$ Then, using theorem 3 we compute the expectation value of occupation number ${{\langle c^\dagger c\rangle}} = \frac{1}{2} - \frac{{ {\rm i} }}{2} {{\langle w_1 w_2\rangle}} = \Gamma_2/(\Gamma_1 + \Gamma_2)$ which is what we expect in canonical equilibrium.

Non-trivial example: transport in quantum spin chains {#sect:nontrivia}
=====================================================

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 $$H =
\sum_{m=1}^{n-1} \left( J^{{\rm x}}_m \sigma^{{\rm x}}_m \sigma^{{\rm x}}_{m+1} + J^{{\rm y}}_m \sigma^{{\rm y}}_m \sigma^{{\rm y}}_{m+1}\right)
+ \sum_{m=1}^n h_m \sigma^{{\rm z}}_m
\label{eq:hamsc}$$ which is coupled to [*two*]{} thermal/magnetic baths [*at the ends*]{} of the chain, generated by two pairs of canonical Lindblad operators [@wich] (similar to (\[eq:canlind\])) $$\begin{aligned}
L_1 &=& \frac{1}{2}\sqrt{\Gamma_1^{\rm L}} \sigma^{-}_1 \qquad
L_3 = \frac{1}{2}\sqrt{\Gamma_1^{\rm R}} \sigma^{-}_n \nonumber \\
 L_2 &=& \frac{1}{2}\sqrt{\Gamma_2^{\rm L}} \sigma^{+}_1 \qquad
 L_4 = \frac{1}{2}\sqrt{\Gamma_2^{\rm R}} \sigma^{+}_n
\label{eq:bathsc}\end{aligned}$$ where $\sigma^{\pm}_m = \sigma^{{\rm x}}_m \pm { {\rm i} }\sigma^{{\rm y}}_m$ and $\Gamma^{{\rm L},{\rm R}}_{1,2}$ are positive coupling constants related to bath temperatures/magnetizations, for example if spins were non-interacting the bath temperatures $T_{{\rm L},{\rm R}}$ would be given with $\Gamma^{{\rm L},{\rm R}}_2/\Gamma^{{\rm L},{\rm R}}_1 =
\exp(-2h_{\rm 1,n}/T_{{\rm L},{\rm R}})$.

Applying the inverse of Jordan-Wigner transformation (\[eq:jordan\]), $\sigma^{{\rm x}}_m = (-{ {\rm i} })^{m-1}\prod_{j=1}^{2m-1} w_j$, $\sigma^{{\rm y}}_m = (-{ {\rm i} })^{m-1}(\prod_{j=1}^{2m-2} w_j)w_{2m}$, we rewrite (\[eq:hamsc\],\[eq:bathsc\]) in terms of Majorana fermions $$\begin{aligned}
H &=& -{ {\rm i} }\sum_{m=1}^{n-1}\left(J^{{\rm x}}_m w_{2m} w_{2m+1} - J^{{\rm y}}_m w_{2m-1}w_{2m+2}\right)
     -{ {\rm i} }\sum_{m=1}^n h_m w_{2m-1}w_{2m} \\
L_1 &=& \frac{1}{2}\sqrt{\Gamma_1^{\rm L}} (w_1 - { {\rm i} }w_2)
\qquad
L_3 = -\frac{(-{ {\rm i} })^n}{2}\sqrt{\Gamma_1^{\rm R}} (w_{2n-1} - { {\rm i} }w_{2n})W \nonumber \\
 L_2 &=& \frac{1}{2}\sqrt{\Gamma_2^{\rm L}}  (w_1 + { {\rm i} }w_2)
\qquad
L_4 = -\frac{(-{ {\rm i} })^n}{2}\sqrt{\Gamma_2^{\rm R}} (w_{2n-1} + { {\rm i} }w_{2n})W\end{aligned}$$ where $W=w_{1}w_{2}\cdots 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 (\[eq:exprlindb\]) where $L_\mu$ enter quadratically, so we find $$\begin{aligned}
{{\hat {\cal L}}}_1 + {{\hat {\cal L}}}_2 &=& -\Gamma_+^{\rm L} ({ {\hat c} }^\dagger_1{ {\hat c} }_1 + { {\hat c} }^\dagger_2{ {\hat c} }_2) -
2{ {\rm i} }\Gamma_-^{\rm L} { {\hat c} }^\dagger_1{ {\hat c} }^\dagger_2 \nonumber \\
{{\hat {\cal L}}}_3 + {{\hat {\cal L}}}_4 &=& -\Gamma_+^{\rm R} ({ {\hat c} }^\dagger_{2n-1}{ {\hat c} }_{2n-1} + { {\hat c} }^\dagger_{2n}{ {\hat c} }_{2n}) -
2{ {\rm i} }\Gamma_-^{\rm R} { {\hat c} }^\dagger_{2n-1}{ {\hat c} }^\dagger_{2n}\end{aligned}$$ leading to the bath shape matrix (\[eq:canbath\]) identical to the single fermion case (\[eq:canlind\]). Again, carefully following the steps of section \[sect:method\], we derive the Liouvillean in the form (\[eq:liouv2\]) with $4n\times 4n$ shape matrix, which we write in a [*block tridiagonal*]{} form in terms of $4\times 4$ matrices as $${{\mathbf{A}}}=
\pmatrix{
{{\mathbf{B}}}_{\rm L} - h_1 {{\mathbf{R}}} & {{\mathbf{R}}}_1 & {{\mathbf{0}}} & \cdots & {{\mathbf{0}}} \cr
-{{\mathbf{R}}}^T_1 & - h_2 {{\mathbf{R}}} & {{\mathbf{R}}}_2 &  \ddots & {{\mathbf{0}}} \cr
{{\mathbf{0}}} & -{{\mathbf{R}}}^T_2 & -h_3 {{\mathbf{R}}} &  & \vdots \cr
\vdots & \ddots & & \ddots & {{\mathbf{R}}}_{n-1} \cr
{{\mathbf{0}}} & {{\mathbf{0}}} & \cdots & -{{\mathbf{R}}}^T_{n-1} & {{\mathbf{B}}}_{\rm R} - h_n {{\mathbf{R}}} }
\label{eq:bigA}$$ and $A_0= \Gamma_+^{\rm L} + \Gamma_+^{\rm R}$, where ${{\mathbf{B}}}_{\rm L}:={{\mathbf{B}}}_{\Gamma_+^{\rm L},\Gamma_-^{\rm L}}$, ${{\mathbf{B}}}_{\rm R}:={{\mathbf{B}}}_{\Gamma_+^{\rm R},\Gamma_-^{\rm R}}$ (in terms of (\[eq:canbath\])), with $\Gamma_{\pm}^{{\rm L},{\rm R}} := \Gamma_2^{{\rm L},{\rm R}}\pm\Gamma_1^{{\rm L},{\rm R}}$, and $${{\mathbf{R}}}_m := \pmatrix{ 0 & 0 & J^{{\rm y}}_m & 0 \cr
                    0 & 0 & 0 & J^{{\rm y}}_m \cr
                    -J_m^{{\rm x}}& 0 & 0 & 0 \cr
                    0 & -J_m^{{\rm x}}& 0 & 0 }$$ We are not able to perform a complete diagonalization of the antisymmetric matrix ${{\mathbf{A}}}$ (\[eq:bigA\]) of the general XY model analytically. For example, even in the spatially homogeneous case $J^{{{\rm x}},{{\rm y}}}_m \equiv J^{{{\rm x}},{{\rm y}}}, h_m\equiv 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 [@rieder]. Namely, in our case ${{\mathbf{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 ${{\mathbf{A}}}$ yielding a set of rapidities $\beta_j$ and properly normalized eigenvector matrix ${{\mathbf{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 \[sect:ti\]), and the case of a [*disordered*]{} XY chain (subsection \[sect:disord\]) 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 (\[eq:hamsc\]) can be written in terms of the two-body [*energy density*]{} operator $$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
H_m =  -{ {\rm i} }J^{{\rm x}}_m w_{2m} w_{2m+1} + { {\rm i} }J^{{\rm y}}_m w_{2m-1}w_{2m+2}
     -\frac{{ {\rm i} }h_m}{2}w_{2m-1}w_{2m}-\frac{{ {\rm i} }h_{m+1}}{2}w_{2m+1}w_{2m+2}
     \label{eq:hdens}$$ as $H=\sum_m H_m$. One can derive the local [*energy current*]{} $Q_m={ {\rm i} }[H_m,H_{m+1}]$ from the [*continuity equation*]{} $$({ {\rm d} }/{ {\rm d} }t){{\langle H_m\rangle}} = \tr H_m { {\rm d} }\rho/{ {\rm d} }t = {{\langle { {\rm i} }[H,H_m]\rangle}}
= -{{\langle Q_{m}\rangle}} + {{\langle Q_{m-1}\rangle}}
\label{eq:cont}$$ where $Q_m := { {\rm i} }[H_{m},H_{m+1}] $ $$\begin{aligned}
Q_m =
&&2{ {\rm i} }(2J^{{\rm y}}_m J^{{\rm x}}_{m+1} w_{2m-1}w_{2m+3} + 2J^{{\rm x}}_mJ^{{\rm y}}_{m+1} w_{2m}w_{2m+4} -
\nonumber\\
  &&-J^{{\rm y}}_m h_{m+1} w_{2m-1}w_{2m+1} - J^{{\rm x}}_m h_{m+1} w_{2m}w_{2m+2}- \nonumber\\
  &&-h_{m+1}J^{{\rm x}}_{m+1}w_{2m+1}w_{2m+3}-h_{m+1}J^{{\rm y}}_{m+1}w_{2m+2}w_{2m+4}) \label{eq:hcurr}\end{aligned}$$ The validity of the above continuity equation (\[eq:cont\]) depends on two assumptions only: (i) All Lindblad operators $L_\mu$ [*commute*]{} with the density $H_m$ in the [*bulk*]{}, $2 \le m \le n-2$ (second equality sign), and (ii) $[H_m,H_{m'}] = 0$ if $|m-m'| \ge 2$ (third equality sign).

Using eq. (\[eq:pairwise2\]) 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*]{} $$\sigma^{{\rm z}}_m = -{ {\rm i} }w_{2m-1}w_{2m}
\label{eq:sdens}$$ and [*spin current*]{} $$S_m = \sigma^{{\rm x}}_{m}\sigma^{{\rm y}}_{m+1}-\sigma^{{\rm y}}_{m}\sigma^{{\rm x}}_{m+1} = -{ {\rm i} }w_{2m}w_{2m+2} -{ {\rm i} }w_{2m-1}w_{2m+1}
\label{eq:scurr}$$ which are all quadratic in $w_j$. Note, however, that the spin density satisfies continuity equation $({ {\rm d} }/{ {\rm d} }t){{\langle \sigma^{{\rm z}}_m\rangle}} = -{{\langle S_{m}\rangle}} + {{\langle S_{m-1}\rangle}}$ only in the isotropic case, when $J^{{\rm x}}_m \equiv J^{{\rm y}}_m$.

Homogeneous transverse Ising chain {#sect:ti}
----------------------------------

Here we limit our discussion to the spatially homogeneous case $J_n^{{{\rm x}},{{\rm y}}}\equiv J^{{{\rm x}},{{\rm y}}},h_n\equiv h$. We shall show that in this case the eigenvalue problem $${{\mathbf{A}}}{{\underline{v}}} =
\beta{{\underline{v}}}
\label{eq:evA}$$ for (\[eq:bigA\]) 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 ${{\underline{v}}}=(\ldots, {{\underline{u}}} \xi^{m-1},{{\underline{u}}} \xi^{m}, {{\underline{u}}} \xi^{m+1},\ldots)^T$, where $\xi$ is a complex [*quasi–momentum*]{} (Bloch) parameter and ${{\underline{u}}}$ is a 4-dimensional amplitude vector satisfying the condition $$(-{{\mathbf{R}}}^T_1 \xi^{-1} - h {{\mathbf{R}}} + {{\mathbf{R}}}_1 \xi - \beta {{\mathbf{I}}}_4) {{\underline{u}}} = 0
\label{eq:freemode}$$ 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 [@ProsenPTPS].

The singularity condition for the free mode equation (\[eq:freemode\]) results, for a general homogeneous $XY$ model, in eight [*master bands*]{} - different values of momenta $\xi$ for each value of the spectral parameter (rapidity) $\beta$. 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^{{\rm x}}= J, J^{{\rm y}}=0$. In this case we find just two master bands with simple dispersion relations $$\xi_{\pm}(\beta) = \frac{h^2 + J^2 + \beta^2 \pm \omega(\beta)}{2 h J}
\quad
\omega(\beta):=\sqrt{(h^2+J^2+\beta^2)^2 - (2 h J)^2}
\label{eq:dispersion}$$ but each band is doubly-degenerate, since the corresponding amplitude problem (\[eq:freemode\]) has two linearly independent solutions $$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
{{\underline{u}}}^\pm_{1}(\beta) = \pmatrix{
-h (h^2 - J^2+\beta^2 \pm \omega) \cr
0 \cr
\beta ( h^2 + J^2 + \beta^2 \pm \omega) \cr
 0
 } \quad
{{\underline{u}}}^\pm_{2}(\beta) = \pmatrix{
0 \cr
-h (h^2 - J^2+\beta^2 \pm \omega)\cr
0 \cr
\beta ( h^2 + J^2 + \beta^2) \pm \omega)
\label{eq:freemodes}
 }$$ Naively speaking, $\xi_-$ represents left moving and $\xi_+$ right moving free modes, each having two possible polarizations. Note that $\xi_- \xi_+ = 1$. For a general complex $\beta$ we shall choose the branch of square root $\omega(\beta)$ (\[eq:dispersion\]) for which $|\xi_-| \le 1$. Let us now write the scattering problem on the [*left*]{} bath in terms of an ansatz $${{\underline{v}}} = \pmatrix{ {{\underline{u}}}_{\rm L} \cr
 \psi_1^- {{\underline{u}}}^-_1 + \psi_2^- {{\underline{u}}}^-_2  + \psi_1^+ {{\underline{u}}}^+_1 + \psi_2^+ {{\underline{u}}}^+_2 \cr
 (\psi_1^- {{\underline{u}}}^-_1 + \psi_2^- {{\underline{u}}}^-_2)\xi_-  + (\psi_1^+ {{\underline{u}}}^+_1 + \psi_2^+ {{\underline{u}}}^+_2)\xi_+ \cr
 (\psi_1^- {{\underline{u}}}^-_1 + \psi_2^- {{\underline{u}}}^-_2)\xi_-^2  + (\psi_1^+ {{\underline{u}}}^+_1 + \psi_2^+ {{\underline{u}}}^+_2)\xi_+^2 \cr
 \vdots}
 \label{eq:scatL}$$ where ${{\underline{u}}}_{\rm L}$ represents a $4$-dimensional vector of left-most eigenvector components, $\psi^-_{1,2}$ are the amplitudes of (known) incident free modes, and $\psi^+_{1,2}$ are the amplitudes of the scattered, outgoing free modes. Plugging the scattering ansatz to the eigenproblem (\[eq:evA\]), the first two rows of ${{\mathbf{A}}}$ (\[eq:bigA\]) result in 6 linearly independent equations for 6 unknowns $\psi^+_{1,2},{{\underline{u}}}_{\rm L}$. Eliminating four variables ${{\underline{u}}}_{\rm L}$ we finally arrive to the non-unitary $2 \times 2$ S-matrix $$\pmatrix{
 \psi^+_1 \cr
 \psi^+_2} = {{\mathbf{S}}}^{\rm L} \pmatrix{
 \psi^-_1 \cr
 \psi^-_2} \label{eq:SL} \\$$ with $$\begin{aligned}
 S^{\rm L}_{11} &=& \tau^{-1} \beta^2 (-(\GLp)^4+4 (\GLp)^2 (\beta^2-3h^2) -16 h (h J^2 + { {\rm i} }\GLm \omega)) \nonumber \\
 S^{\rm L}_{12} &=& \tau^{-1} \beta ( (\GLp)^3 + 8{ {\rm i} }\GLm h \beta + 4 \GLp (h^2 - \beta^2)) (2{ {\rm i} }\omega) \nonumber \\
 S^{\rm L}_{21} &=& \tau^{-1} \beta ( (\GLp)^3 - 8{ {\rm i} }\GLm h \beta + 4 \GLp (h^2 - \beta^2)) (-2{ {\rm i} }\omega) \nonumber \\
 S^{\rm L}_{22} &=& \tau^{-1} \beta^2 (-(\GLp)^4+4 (\GLp)^2 (\beta^2-3h^2) -16 h (h J^2 - { {\rm i} }\GLm \omega)) \label{eq:Sexplicit} \\
\tau &:=& (\GLp)^4 \beta^2 + 8 \beta^2 (h^4 + (J^2 + \beta^2)(J^2 + \beta^2 - \omega) + h^2 (2 \beta^2 - \omega)) \nonumber \\
&-& 2(\GLp)^2 (h^4 + J^4 + 3\beta^4 + J^2 (2 \beta^2-\omega) - \beta^2\omega + h^2 (\omega - 2 J^2 - 4\beta^2)) \nonumber
 \end{aligned}$$ Similarly, one can solve the scattering problem from the [*right*]{} bath with the scattering ansatz $${{\underline{v}}} = \pmatrix{
 \vdots \cr
 (\psi_1^+ {{\underline{u}}}^+_1 + \psi_2^+ {{\underline{u}}}^+_2)\xi_+^{-2}  + (\psi_1^- {{\underline{u}}}^-_1 + \psi_2^- {{\underline{u}}}^-_2)\xi_-^{-2} \cr
 (\psi_1^+ {{\underline{u}}}^+_1 + \psi_2^+ {{\underline{u}}}^+_2)\xi_+^{-1}  + (\psi_1^- {{\underline{u}}}^-_1 + \psi_2^- {{\underline{u}}}^-_2)\xi_-^{-1} \cr
 \psi_1^+ {{\underline{u}}}^+_1 + \psi_2^+ {{\underline{u}}}^-_2  + \psi_1^- {{\underline{u}}}^-_1 + \psi_2^- {{\underline{u}}}^-_2 \cr
 {{\underline{u}}}_{\rm R}}$$ defining the right S-matrix $$\pmatrix{
 \psi^-_1 \cr
 \psi^-_2} = {{\mathbf{S}}}^{\rm R} \pmatrix{
 \psi^+_1 \cr
 \psi^+_2}$$ Note that since the two directions of free modes (\[eq:freemodes\]) do not have left-right symmetry an explicit expression for $S^{\rm R}$ is considerably more complicated than (\[eq:Sexplicit\]) and shall not be written out here. We shall now show that there exist two qualitatively different types of NMM - complex rapidities $\beta$ solving (\[eq:evA\]) for sufficiently [*large*]{} $n$.

![ Rapidities $\beta_j$ (black dots) for a transverse Ising chain with $J=1.5,h=1$ and bath couplings $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$, for three different sizes $n=6$ (upper), $n=30$ (middle), and $n=150$ (lower panel). Big blue/red dots indicate positions of evanescent rapidities (solutions of eq.(\[eq:4order\])) for the left/right bath. \[fig:rapid\]](rapid.eps)

First, we shall discuss the so called [*evanescent normal master modes*]{}. These are characterized with amplitudes (\[eq:scatL\]) 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\to \infty$. Such solutions $\psi^{+}_{1,2} = 0$ of eq. (\[eq:SL\]) exist exactly when the determinant of S-matrix vanishes $\det S^{\rm L} = 0$. Using (\[eq:Sexplicit\]) the determinant can be written as $\det S^{\rm L} = (\beta/\tau)^2 p^{\rm L}(\beta)$ where[^4] $$\begin{aligned}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!p^{\rm L}(\beta) &=& (\GLp)^8 \beta^2 - 4(\GLp)^6 ((h^2\!-\!J^2)^2  + (2J^2\!-\!4h^2)\beta^2 + 3\beta^4) \nonumber \\
&-& 16 (\GLp)^4(2 h^2 (h^2\!-\!J^2)^2 - (7h^4\!-\!6 h^2J^2\!+\!2 J^4)\beta^2 + 4(h^2\!-\!J^2)\beta^4 - 3 \beta^6)   \nonumber \\
&-& 64 (\GLp)^2(h^4 (h^2\!-\!J^2)^2-2h^2 J^4\beta^2 -(2h^4\!+\!4 h^2 J^2\!-\!J^4)\beta^4 + 2 J^2 \beta^6 + \beta^8) \nonumber \\
&+& 256 h^4 J^4 \beta^2 \label{eq:p}\end{aligned}$$ 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 $\beta^2$ $$p^{\rm L}(\beta_{\rm evan}) = 0 \label{eq:evan}
\label{eq:4order}$$ that are [*not*]{} simultaneously zeroes of $\tau(\beta)$. 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 $\GLp$ by $\GRp$ in (\[eq:evan\],\[eq:p\]). In fig. \[fig:rapid\] we compare evanescent rapidities computed from eq. (\[eq:evan\]) to numerically calculated spectrum of ${{\mathbf{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, $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm 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 (\[eq:evA\]) in terms of a $2\times 2$ determinant $$\det ( \xi_+^{2(n-3)}  {{\mathbf{S}}}^{\rm R} {{\mathbf{S}}}^{\rm L}- {{\mathbf{I}}}_2) = 0
\label{eq:sec}$$ In the absence of the baths, $\Gamma^{\rm L,R}_{\pm} = 0$, the solutions of the above problem exist only for real quasi-momenta, namely $\xi_{\pm} = \exp(\pm{ {\rm i} }\vartheta), \vartheta \in {\mathbb{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 $\beta(e^{{ {\rm i} }\vartheta}) = \pm { {\rm i} }\varepsilon(\vartheta)$ with ‘energy’ $$\varepsilon(\vartheta) = \sqrt{h^2 + J^2 - 2 |h J| \cos\vartheta}$$ 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. $\vartheta \approx \pi/n$, or $\vartheta \approx \pi + \pi/n$, and should thus lie near the band edges $\pm { {\rm i} }|h| \pm { {\rm i} }|J|$ (see fig. \[fig:rapid\]). In the following we shall focus our calculation on the band edge $\beta^* = { {\rm 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 ${ {\rm i} }(|h|-|J|)$, and hence really determines the gap of the Liouvillean. So we write $$\beta = { {\rm i} }(|h| + |J|) + z$$ where $z\in{\mathbb{C}}$ is a small parameter, and expand the S-matrices around the band edge $${{\mathbf{S}}}^{\rm L,R} = -{{\mathbf{I}}}_2 + \frac{4g}{\eta^{\rm L,R}}
{{\mathbf{Z}}}^{\rm L,R} \sqrt{-{ {\rm i} }z} + {\cal O}(|z|)
\label{eq:SZ}$$ where $g:=\sqrt{\frac{|hJ|}{2(|h|+|J|)}}$, $\eta^{\rm L,R}:=(\Gamma^{\rm L,R}_+)^4 + 4 (\Gamma^{\rm L,R}_+)^2(4 h^2 + 2 |h J| + J^2) + 16 h^2 J^2$ and $$\begin{aligned}
Z^{\rm L}_{11} &=& 4 |h|(\GLp)^2 + 16|h|(|h|+|J|)(|J|-{ {\rm i} }\GLm) \nonumber \\
Z^{\rm L}_{12} &=& -2 (\GLp)^3 - 16 \GLm |h|(|h|+|J|)-8(2h^2+2|hJ|+J^2) \nonumber\\
Z^{\rm L}_{21} &=& +2 (\GLp)^3 - 16 \GLm |h|(|h|+|J|)+8(2h^2+2|hJ|+J^2) \nonumber\\
Z^{\rm L}_{22} &=& 4 |h|(\GLp)^2 + 16|h|(|h|+|J|)(|J|+{ {\rm i} }\GLm) \label{eq:ZL}\end{aligned}$$ and $$\begin{aligned}
Z^{\rm R}_{11} &=& (\GRp)^4(2|h|+|J|)+4(\GRp)^2(8|h|^3+9h^2|J|+4|h|J^2+|J|^3) \nonumber \\
                          &+& 16 h^2|J|(|J|(3|h|+2|J|)-{ {\rm i} }\GRm(|h|+|J|)) \nonumber \\
Z^{\rm R}_{12} &=& -2 (\GRp)^3 - 16 \GRm |h|(|h|+|J|)-8(2h^2+2|hJ|+J^2) \nonumber\\
Z^{\rm R}_{21} &=& +2 (\GRp)^3 - 16 \GRm |h|(|h|+|J|)+8(2h^2+2|hJ|+J^2) \nonumber\\
Z^{\rm R}_{22} &=& (\GRp)^4(2|h|+|J|)+4(\GRp)^2(8|h|^3+9h^2|J|+4|h|J^2+|J|^3) \nonumber \\
                          &+& 16 h^2|J|(|J|(3|h|+2|J|)+{ {\rm i} }\GRm(|h|+|J|)) 
\label{eq:ZR}\end{aligned}$$ Next we expand $\xi_+$ (\[eq:dispersion\]) in $z$, yielding $$\xi_+ = -1 - g^{-1} \sqrt{-{ {\rm i} }z} + {\cal O}(|z|)
\label{eq:xi1}$$ and so the [*free propagator*]{} in (\[eq:sec\]) can be written as $$\xi_+^{2(n-3)} = \exp(2n g^{-1}\sqrt{-{ {\rm i} }z}) + {\cal O}(|z|).
\label{eq:xiexp}$$ In eqs. (\[eq:SZ\],\[eq:xi1\],\[eq:xiexp\]) the branch cut along the negative real axis has been chosen for $\sqrt{-{ {\rm i} }z}$. Since the product of S-matrices in (\[eq:sec\]) is near identity, the free propagator should be near one as well, hence $2n g^{-1}\sqrt{-{ {\rm i} }z}$ should be near $2\pi{ {\rm i} }$. Let us define $z_0$ by setting $2n g^{-1}\sqrt{-{ {\rm i} }z_0} = 2\pi{ {\rm i} }$, so $$z_0 = -{ {\rm i} }\pi^2 g^2 n^{-2}$$ and write $z = z_0 (1 + y) $ where $|y| \ll 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 (\[eq:sec\]) since the real part of the soft mode’s rapidity is determined as $$\re \beta = \re z_0 y = \pi^2 g^2 n^{-2} \im y
\label{eq:y}$$

Now, writing $\sqrt{-{ {\rm i} }z} = \sqrt{-{ {\rm i} }z_0}\sqrt{1+y} = { {\rm i} }\pi g n^{-1} (1 + y/2 - y^2/8) + {\cal O}(y^3)$ in (\[eq:SZ\],\[eq:xiexp\]), plugging all that to eq. (\[eq:sec\]) and computing to order ${\cal O}(n^{-2})$, noting that ${\cal O}(|z|)={\cal O}(n^{-2})$, we arrive to a simple quadratic equation for $y$, whose solution, plugged to (\[eq:y\]), gives the final result, namely the sectral gap of Liouvillean $\Delta= 2\re \beta $ $$\begin{aligned}
\Delta &=& \frac{(2\pi h J)^2}{(|h|+|J|)^2} \frac{\Delta_1}{\Delta_2} n^{-3} + {\cal O}(n^{-4}) \label{eq:delta} \\
\Delta_1 &:=& 64 (\GLp + \GRp) h^2 J^2 (2 h^2 + 2|h J| + J^2)  \nonumber \\
                &+& 16 ((\GLp)^3 + (\GRp)^3) h^2 J^2  \nonumber \\
                &+& 16 \GLp \GRp (\GLp + \GRp) (2 h^2 + 2|hJ| + J^2)(4 h^2 + 2|hJ| + J^2) \nonumber \\
                &+& 4 \GLp \GRp ((\GLp)^3 + (\GLp)^2 \GRp + \GLp (\GRp)^2 + (\GRp)^3) (2 h^2+ 2|h J| + J^2) \nonumber \\
                 &+& (\GLp \GRp)^3 (\GLp+\GRp) \nonumber \\
\Delta_2 &:=&  ((\GLp)^4 + 4(\GLp)^2 (4 h^2 + 2|hJ| + J^2)+ 16 h^2 J^2 ) \nonumber \\
               &\times& ((\GRp)^4 + 4(\GRp)^2 (4 h^2 + 2|hJ| + J^2)+ 16 h^2 J^2) \nonumber
 \end{aligned}$$ In fig. (\[fig:delta\]) we compare this analytical result to exact numerical calculations of the eigenvalue of ${{\mathbf{A}}}$ with minimal real part, confirming both, its precise numerical value and that the relative scaling of the next order correction is indeed ${\cal O}(n^{-1})$.

Note that, interestingly, both main analytical results of this subsection, namely evanescent and soft mode rapidities [*do not depend*]{} on $\Gamma^{\rm L,R}_-$. Physically speaking, they only depend on the effective strengths of the bath couplings and not on the temperatures.

![ Spectral gap $\Delta$ times a third power of the chain length $n$ for a transverse Ising chain with $J=1.5,h=1$ and bath couplings $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$. Thin horizontal line indicates the theoretical asymptotic value (\[eq:delta\]). In the inset we show deviation from asymptotic constant value of $\Delta n^3$ in log-log scale and demonstrate that it decays as $\propto n^{-1}$ (thin line). \[fig:delta\]](deltafine.eps)

![ Spectral gap $\Delta$ times a third power of the chain length $n$ for a transverse Ising chain with $J=1.5,h=1$ and bath couplings $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$. Thin horizontal line indicates the theoretical asymptotic value (\[eq:delta\]). In the inset we show deviation from asymptotic constant value of $\Delta n^3$ in log-log scale and demonstrate that it decays as $\propto n^{-1}$ (thin line). \[fig:delta\]](delta.eps)

![ Complete spectrum of $2^{12}$ complex eigenvalues of Liouvillean for a transverse Ising chain with $n=6$ spins and $J=1.5,h=1$ and bath couplings $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$ (the case of the upper panel of fig. \[fig:rapid\]). \[fig:spc\]](spc.eps)

![ Energy current (upper/blue points), and average spin current (lower/red points), versus chain length $n$ for a transverse Ising chain with $J=1.5,h=1$ and bath couplings $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$. \[fig:current\]](current.eps)

![ 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 $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$. The insets display logarithm of the difference to the bulk values $\delta H_m :=  |{{\langle H_m\rangle}} - H_{\rm bulk}|$ (blue points), $\delta \sigma^{{\rm z}}_m :=  |{{\langle \sigma^{{\rm z}}_m\rangle}} - \sigma^{{\rm z}}_{\rm bulk}|$ (red points) in comparison with $\pm (4 \log \xi_- ) m + {\rm const}$ with quasi-momentum $\xi_- = 0.584692$ corresponding (\[eq:dispersion\]) to the leading evanescent rapidity $\beta_{\rm evan} = 0.438739$ (full lines). \[fig:profile\]](profile1.eps)

![ 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 $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$. The insets display logarithm of the difference to the bulk values $\delta H_m :=  |{{\langle H_m\rangle}} - H_{\rm bulk}|$ (blue points), $\delta \sigma^{{\rm z}}_m :=  |{{\langle \sigma^{{\rm z}}_m\rangle}} - \sigma^{{\rm z}}_{\rm bulk}|$ (red points) in comparison with $\pm (4 \log \xi_- ) m + {\rm const}$ with quasi-momentum $\xi_- = 0.584692$ corresponding (\[eq:dispersion\]) to the leading evanescent rapidity $\beta_{\rm evan} = 0.438739$ (full lines). \[fig:profile\]](profile2.eps "fig:")![ 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 $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$. The insets display logarithm of the difference to the bulk values $\delta H_m :=  |{{\langle H_m\rangle}} - H_{\rm bulk}|$ (blue points), $\delta \sigma^{{\rm z}}_m :=  |{{\langle \sigma^{{\rm z}}_m\rangle}} - \sigma^{{\rm z}}_{\rm bulk}|$ (red points) in comparison with $\pm (4 \log \xi_- ) m + {\rm const}$ with quasi-momentum $\xi_- = 0.584692$ corresponding (\[eq:dispersion\]) to the leading evanescent rapidity $\beta_{\rm evan} = 0.438739$ (full lines). \[fig:profile\]](profile3.eps "fig:")

We end this subsection by presenting some further numerical results on heat transport in the open transverse Ising chain in the Lindblad form. In fig. \[fig:spc\] we demonstrate expression (\[eq:eval\]) for constructing the full spectrum of the Liouvillean in terms of a set of rapidities, for a short chain. In fig. \[fig:current\] we demonstrate eq. (\[eq:pairwise2\]) of Theorem 3 by computing the energy current $Q_m$ (\[eq:hcurr\]), and the average spin current $S=\frac{1}{n-1}\sum_{m=1}^{n-1} S_m$ (\[eq:scurr\]) in NESS of a typical transverse Ising chain. Numerical results give a clear indication of [*ballistic transport*]{} ${{\langle Q\rangle}} = {\cal O}(n^0), {{\langle S\rangle}}={\cal O}(n^0)$, 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 fig. \[fig:profile\] we plot the energy density (\[eq:hdens\]) and spin density (\[eq:sdens\]) profiles in NESS. Again, we note flat profiles in the bulk of the chain, $m,n-m \gg 1$, with exponential falloff due to adjustment to the non-equilibrium bath values. Since the densities can be written, by means of (\[eq:pairwise2\]), as $4-$point functions in NMM components, the leading falloff exponents of the profile $|{{\langle H_m\rangle}} - H_{\rm bulk}| \sim
|\xi_-|^{4m} $ is given by the quasi-momentum $\xi_-$ (\[eq:dispersion\]) corresponding to the maximal evanescent rapidity $\beta_{\rm evan}$ (\[eq:evan\]).

Disordered XY chain {#sect:disord}
-------------------

In this subsection we treat the opposite extreme, a disordered XY chain (\[eq:hamsc\]) where three sets of physical parameters are chosen as [*random uncorrelated*]{} variables from [*uniform*]{} distributions on the intervals, $J^{{\rm x}}_m \in [J^{{\rm x}}_{\rm min},J^{{\rm x}}_{\rm max}]$, $J^{{\rm y}}_m \in [J^{{\rm y}}_{\rm min},J^{{\rm y}}_{\rm max}]$, $h_m \in [h_{\rm min},h_{\rm max}]$. Clearly, the eigenvalue problem (\[eq:evA\]) for the matrix (\[eq:bigA\]) 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.

![ Average Liouvillean spectral gap ${{\langle \Delta\rangle}}$ versus the chain length $n$ for disordered XY models: (i) $J^{{\rm x}}_m = 0.5, J^{{\rm y}}_m = 0, h_m \in [1,2]$ (transverse Ising with field disorder, blue points), (ii) $J^{{\rm x}}_m \in [0.5,2], J^{{\rm y}}_m = 0, h_m = 1$ (transverse Ising with interaction disorder, red points), (iii) $J^{{\rm x}}_m \in [0.5,1],J^{{\rm y}}_m \in [0.5,1], h_m = 1$ (XY with interaction disorder, golden points), all for bath couplings $\Gamma^{\rm L}_1=1,\Gamma^{\rm L}_2=0.6$, $\Gamma^{\rm R}_1=1,\Gamma^{\rm R}_2=0.3$. Full lines indicate exponential fits to right halves of data. Averaging is performed over $2000$ disorder realizations. \[fig:disdelta\]](disdelta.eps)

![ The scaling of the energy current ${{\langle Q_m\rangle}}$ with chain length $n$ of the disordered XY model in the same regimes/parameters/plot styles as in fig. \[fig:disdelta\]. \[fig:discurrent\]](discurrent.eps)

![ Scaled energy density profile of interaction disordered XY chains (case (iii) of fig. \[fig:disdelta\]) for three chain sizes: $n=20$ (blue points), $n=40$ (red points), $n=60$ (golden points). Averages over 50000 disorder realizations have been performed. \[fig:disprofile\]](disprofile.eps)

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 [@livi] for a review of related studies): The spectral gap of the Liouvillean should be exponentially small $\Delta \sim \exp(-n/\ell)$ where $\ell$ is the localization length of NMM which is expected to be proportional to the square of inverse disorder strength. This is demonstrated in fig.\[fig:disdelta\]. 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 (\[eq:hcurr\]) in fig. \[fig:discurrent\]. In the final figure \[fig:disprofile\] we plot the energy density profile ${{\langle H_m\rangle}}$ (\[eq:hdens\]) in a typical case of disordered XY chain, versus a scaled spatial coordinate $(m-1)/(n-1) \in [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 {#sect:conc}
==========================

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 [@lieb] 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 [@ProsenPTPS; @prosen; @pizorn] 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 [@dhar].

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 tight-binding model [@piere]. We also expect the method to be applicable to the Redfield type of master equations (see e.g. [@piere]) - 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 \[sect:method\].

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 [@latorre], 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*]{} [@pizorn] (which is suggested by numerical experiments [@PZ08]) indicates [*efficient simulability*]{} of NESS by elaborate numerical methods such as [*density matrix renormalization group*]{} [@dmrg], 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 [@korepin], 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?

Acknowledgements {#acknowledgements .unnumbered}
================

I gratefully acknowledge stimulating discussions with Pierre Gaspard, Keiji Saito and Walter Strunz, 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 Slovenian research agency (ARRS). Explicit analytical calculations reported in subsection (\[sect:ti\]) were assisted by [*Mathematica*]{} software package.

References {#references .unnumbered}
==========

[10]{}

W. H. Zurek, Rev. Mod. Phys. [**75**]{}, 715 (2003).

E. Joos, H. D. Zeh, C. Kiefer, D. Giulini, J. Kupsch and I.-O. Stamatescu, [*Decoherence and the Appearance of a Classical World in Quantum Theory*]{}, (Springer, 2003).

J. von Neumann, [*Mathematical Foundations of Quantum Mechanics*]{}, Trans. Robert T. Geyer. (Princeton University Press, Princeton 1955).

M. Schlosshauer, Rev. Mod. Phys. [**76**]{}, 1267 (2004).

M. A. Nielsen and I. L. Chuang, [*Quantum Computation and Quantum Information*]{} (Cambridge University Press, Cambridge 2000).

G. Benenti, G. Casati and G. Strini, [*Principles of Quantum Computation and Information. Volume I: Basic Concepts*]{} (World Scientific, Singapore 2004); [*Volume II: Basic Tools and Special Topics*]{} (World Scientific, Singapore 2007).

H. Araki and E. Barouch, J. Stat. Phys. [**31**]{}, 327 (1983);\
H. Araki, Publ. RIMS Kyoto Univ. [**20**]{}, 277 (1984).

D. Ruelle, J. Stat. Phys. [**98**]{}, 57 (2000).

V. Jakšič and C.-A. Pillet, J. Stat. Phys. [**108**]{}, 787 (2002); Commun. Math. Phys. [**226**]{}, 131 (2002);\
W. Aschbacher, V. Jakšič, Y. Pautrat and C.-A. Pillet, [*Inroduction to non-equilibrium quantum statistical mechanics*]{}, in [*Open Quantum Systems III. Recent Developments*]{} Lecture Notes in Mathematics, [**1882**]{} (2006), 1-66.

P. Gaspard, Prog. Theor. Phys. Suppl. [**165**]{}, 33 (2006); Physica A [**369**]{}, 201 (2006).

M. H. Lee, Acta Physica Polonica B [**38**]{}, 1837 (2007); Phys. Rev. Lett. [**87**]{}, 250601 (2001); Phys. Rev. Lett. [**49**]{}, 1072 (1982).

T. Prosen, J. Phys. A: Math. Theor. [**40**]{}, 7881 (2007).

V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, [*Quantum Inverse Scattering and Correlation functions*]{} (Cambridge University Press, Cambridge 1997).

L. Fadeev, P. Van Moerbeke and F. Lambert (Eds.), [ *Bilinear Integrable Systems: from Classical to Quantum, Continuous to Discrete*]{}, (NATO ARW Proceedings), Springer Series: NATO Science Series II: Mathematics, Physics and Chemistry, Vol 201 (2006).

H.-P. Breuer and F. Petruccione, [*The Theory of Open Quantum Systems*]{} (Oxford University Press, Oxford 2002).

F. Haake, [*Quantum Signatures of Chaos*]{}, 2nd edition (Springer, 2001).

R. Alicki and K. Lendi, [*Quantum dynamical semigroups and applications*]{} (Springer, 2007).

E. H. Lieb, T. D. Schultz and D. C. Mattis, Ann. Phys. (New York) [**16**]{}, 407 (1961).

X. Zotos, F. Naef and P. Prelovšek, Phys. Rev. B [**55**]{}, 11029 (1997).

K. Saito, S. Takesue and S. Miyashita, Phys. Rev. E [**61**]{}, 2397 (2000);\
K. Saito, Europhys. Lett. [**61**]{}, 34 (2003).

M. Michel, M. Hartmann, J. Gemmer and G. Mahler, Eur. Phys. J. B. [**34**]{}, 325 (2003);\
M. Michel, G. Mahler and J. Gemmer, Phys. Rev. Lett. [**95**]{}, 180602 (2005);\
M. Michel, J. Gemmer and G. Mahler, Int. J. Mod. Phys. B [**20**]{}, 4855 (2006);\
J. Gemmer, R. Steinigeweg and M. Michel, Phys. Rev. B [**73**]{}, 104302 (2006).

C. Mejia-Monasterio, T. Prosen and G. Casati, Europhys. Lett. [**72**]{}, 520 (2005);\
G. Casati and C. Mejia-Monasterio, e-print [arXiv:0710.3500v1 \[cond-mat.stat-mech\]]{}.

C. Mejia-Monasterio and H. Wichterich, e-print [arXiv:0709.1412v1 \[cond-mat.stat-mech\]]{}.

A. Dhar and D. Roy, J. Stat. Phys. [**125**]{}, 805 (2006); see also e-print [arXiv:0711.4318v1 \[cond-mat.stat-mech\]]{}.

G. Lindblad, Commun. Math. Phys. [**48**]{}, 119 (1976).

M. M. Wolf, J. Eisert, T. S. Cubitt and J. I. Cirac, e-print [arXiv:0711.3172v1 \[quant-ph\]]{}.

P. [Šemrl]{}, [*private communication*]{}.

We note a similarity to the formalism of second quantization with non-orthogonal orbitals introduced in: M. Moshinsky and T. H. Seligman, Ann. Phys. (New York) [**66**]{}, 311 (1971).

H. Wichterich, M. J. Herich, H. P. Breuer, J. Gemmer and M. Michel, Phys. Rev. E, [**76**]{} 031115 (2007).

Z. Rieder, J. L. Lebowitz and E. Lieb, J. Math. Phys. [**8**]{}, 1073 (1967).

T. Prosen, Prog. Theor. Phys. Suppl. [**139**]{}, 191 (2000).

S. Lepri, R. Livi and A. Politi, Phys. Rep. [**377**]{}, 1 (2003).

T. Prosen, Phys. Rev. E **60**, 1658 (1999).

T. Prosen and I. Pižorn, Phys. Rev. A [**76**]{}, 032316 (2007).

M. Esposito and P. Gaspard, Phys. Rev. B [**71**]{}, 214302 (2005); J. Stat. Phys. [**121**]{}, 463 (2005).

G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. [**90**]{}, 227902 (2003);\
J. I. Latorre, E. Rico, and G. Vidal, Quant. Inf. Comp.[**4**]{}, 48 (2004).

T. Prosen and M. Žnidarič, to be submitted (2008).

S. R. White, Phys. Rev. Lett. **69**, 2863 (1992);\
U. Schollw[" ock]{} and S. R. White, in G. G. Batrouni, and D. Poilblanc (eds.): Effective models for low-dimensional strongly correlated systems, p.155, AIP, Melville, New York (2006);\
G. Vidal, Phys. Rev. Lett. [**91**]{}, 147902 (2003); [*ibid.*]{} [**93**]{}, 040502 (2004).

[^1]: Throughout this paper Dirac’s bra-ket notation shall be used only for a Hilbert space ${\cal 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 ${\cal K}$. For instance, we note a key distinction between [*physical fermions*]{} $c_m$ (\[eq:physfermi\]) and a-fermions $\hat{c}_j$ (\[eq:defanih\]).

[^2]: It is not known at present whether explict form (\[eq:explA\]) guarantees diagonalizability of any such ${{\mathbf{A}}}$. Note that one can construct certain types of complex antisymmetric matrices with degenerate eigenvalues which cannot be diagonalized [@semrl].

[^3]: For a non-degenerate rapidity spectrum $\{\beta_j\}$ the proof of this statement is a trivial consequence of antisymmetry ${{\mathbf{A}}}=-{{\mathbf{A}}}^T$, whereas in case of degeneracies it can be shown that one can always choose appropriate linear combinations of eigenvectors.

[^4]: Trivial zero $\beta=0$ of course does not represent a physical solution since then the whole S-matrix (\[eq:Sexplicit\]) vanishes.