# Exact steady state of the open $XX$-spin chain:
entanglement and transport properties

F. Benatti  Department of Physics, University of Trieste, I-34151 Trieste,
Italy  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Trieste,
I-34151 Trieste, Italy  R. Floreanini  Istituto Nazionale di Fisica Nucleare
(INFN), Sezione di Trieste, I-34151 Trieste, Italy  f@infn L. Memarzadeh
<EMAIL_ADDRESS>
The authors contributed equally to this study and are listed in alphabetical
order. Department of Physics, Sharif University of Technology, Tehran, Iran

###### Abstract

We study the reduced dynamics of open quantum spin chains of arbitrary length
$N$ with nearest neighbour $XX$ interactions, immersed within an external
constant magnetic field along the $z$ direction, whose end spins are weakly
coupled to heat baths at different temperatures, via energy preserving
couplings. We find the analytic expression of the unique stationary state of
the master equation obtained in the so-called global approach based on the
spectralization of the full chain Hamiltonian. Hinging upon the explicit
stationary state, we reveal the presence of sink and source terms in the spin-
flow continuity equation and compare their behaviour with that of the
stationary heat flow. Moreover, we also obtain analytic expressions for the
steady state two-spin reduced density matrices and for their concurrence. We
then set up an algorithm suited to compute the stationary bipartite
entanglement along the chain and to study its dependence on the Hamiltonian
parameters and on the bath temperatures.

##  I Introduction

Transport phenomena in open interacting quantum spin chains have recently
received an increasing attention as instances of many-body systems driven by
intrinsic inter-spin interactions and coupled to external heat baths at the
two ends of the chain. Specific experimental realizations have been reported
in scenarios involving ultracold-atoms, light-harvesting complexes and quantum
thermodynamics at large Datta -Prosen5 .

In presence of external baths, the reduced dynamics of any open quantum system
is obtained by tracing over the baths’ degrees of freedom. When the strength
of the system-baths interaction is small, applying the weak-coupling limit
techniques yields a dissipative irreversible time evolution that is generated
by a master equation in Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form
Alicki-Lendi -Merkli .

The derivation of the GKSL master equation requires the diagonalization of the
full spin-chain Hamiltonian. Due to the many degrees of freedom of the quantum
spin chain and their mutual interactions, dissipative effects then arise
involving all spins in the chain together with environment-induced excitation
transfer between different sites (e.g. see Davies4 -Rivas2 ). These gives rise
to new, global effects in transport phenomena that can not be captured using
other, simplified approaches to the chain open dynamics.111 As finding
eigenvalues and eigenvectors of the system Hamiltonian might in general be
laboriously difficult, an alternative approach has been often advocated,
consisting in neglecting the inter-spin interaction in the derivation of the
master equa-tion (e.g. see Michel -Hovhannisyan ). Although the two
approaches, named global and local, are regularly adopted in applications and
compared Rivas1 ; Guimaraes ; Werlang ; Santos ; Migliore ; Zoubi , Levy
-Giovannetti20 , it turns out that the local approach might not be able to
capture all the correct system transport properties BFM .

In the following, we focus on the study of the stationary transport and
bipartite entanglement properties of open $XX$-chains with energy conserving
couplings to external baths. We provide an explicit analytic form for the
chain stationary state by means of which we obtain analytic expressions for
the spin flow, revealing the presence of sink and source terms, and for the
heat flow. Remarkably, we have also been able to explicitly compute the
reduced two-spin density matrices resulting from the stationary state and
study the corresponding bipartite entanglement along the chain. For the latter
task, we develop a suitable algorithmic representation of the stationary state
in the spin representation.
The structure of the paper is as follows: in section II we set the framework
for the derivation of the open chain dynamics and diagonalize the chain
Hamiltonian by turning the spin representation into a Fermionic one. In
Section III we derive the Lindblad operators yielding the dissipative
contribution to the GKSL master equation, prove that the latter has a unique
stationary state and explicitly derive its expression in the Fermionic
representation. Based on it, we then discuss the ensuing transport properties
in terms of spin and heat flows. In section IV we rewrite the stationary state
in the spin representation and show that it provides reduced two-spin density
matrices in the so-called $X$ form which allows for a simple analytical
expression of the concurrence. Then, we set up a representation of the
stationary state that is best suited for the study of bipartite entanglement
and its dependence on the various parameters of the chain and on the
temperatures of the baths coupled to it. We conclude by summarizing and
discussing the results, while the more technical issues are presented in
various Appendices.

##  II Open $XX$ spin chain of length $N$

As mentioned in the Introduction, in the following we address an open quantum
chain consisting of $N$ spins at sites $1,2,\ldots,N$, immersed in a constant
magnetic field along the $z$ direction, with $XX$ nearest neighbour
interactions among themselves. The ensuing closed chain dynamics is thus
generated by the following Hamiltonian:

$H=g\sum_{\ell=1}^{N-1}\left(\sigma_{x}^{(\ell)}\sigma_{x}^{(\ell+1)}+\sigma_{y}^{(\ell)}\sigma_{y}^{(\ell+1)}\right)+\Delta\sum_{\ell=1}^{N}\sigma_{z}^{(\ell)}\
,$ (1)

with free boundary ocnditions, where $\Delta>0$ is the intensity of the
constant transverse magnetic field, $\sigma_{x,y,z}^{(\ell)}$ are the Pauli
matrices at site $\ell$, and $g>0$ is the strength of the nearest neighbour
interaction. Throughout the paper we work in natural units where both Planck
and Boltzmann constants are set to 1, $\hbar=\kappa_{B}=1$.

The spin chain is then turned into an open many-body quantum system by
coupling the two end spins, at site 1 on the left end, $L$, and at site $N$ on
the right end, $R$, to two independent free Bosonic thermal baths with
Hamiltonians

$H_{\alpha}=\int_{0}^{+\infty}{\rm
d}\nu\,\nu\,\mathfrak{b}^{\dagger}_{\alpha}(\nu)\,\mathfrak{b}_{\alpha}(\nu)\
,\qquad\alpha=L,R\ ,$ (2)

where $\mathfrak{b}_{\alpha}(\nu)$, $\mathfrak{b}^{\dagger}_{\alpha}(\nu)$ are
Bosonic operators satisfying the canonical commutation relations

$[\mathfrak{b}_{\alpha}(\nu),\,\mathfrak{b}_{\beta}^{\dagger}(\nu^{\prime})]=\delta_{\alpha\beta}\,\delta(\nu-\nu^{\prime})\
.$

The coupling of the baths to the left and right spins are described by the
interaction Hamiltonian

$H^{\prime}=\lambda\,\sum_{\alpha=L,R}\Big{(}\sigma_{+}^{(\alpha)}\mathfrak{B}_{\alpha}+\sigma_{-}^{(\alpha)}\mathfrak{B}_{\alpha}^{\dagger}\Big{)}\
,$ (3)

where $\lambda<<1$ is a dimensionless coupling constant,

$\sigma_{\pm}^{(\ell)}\equiv\frac{1}{2}\big{(}\sigma^{(\ell)}_{x}\pm
i\sigma^{(\ell)}_{y}\big{)}\ ,$ (4)

are spin ladder operators at site $\ell$, whence
$\sigma_{\pm}^{(L)}=\sigma_{\pm}^{(1)}$ and
$\sigma_{\pm}^{(R)}=\sigma_{\pm}^{(N)}$, while

$\mathfrak{B}_{\alpha}=\int_{0}^{\infty}{\rm
d}\nu\,h_{\alpha}(\nu)\,\mathfrak{b}_{\alpha}(\nu)\
,\quad[h_{\alpha}(\nu)]^{*}=h_{\alpha}(\nu)\ ,$ (5)

are bath operators, with $*$ meaning complex conjugation and $h_{L,R}(\nu)$
are suitable smearing functions. Referring to BFM  for more details, we begin
by shortly reviewing the rigorous weak-coupling limit derivation of the open
chain master equation of GKSL type in the so-called global approach. As we
shall show, the resulting dissipative dynamics of the $N$ spins in the
presence of the two baths involves the full inter-spin interactions.

Assuming the free Boson baths to be in their equilibrium Gibbs states at
temperatures $T_{L}\equiv 1/\beta_{L}$ and $T_{R}\equiv 1/\beta_{R}$, the
state of the environment is then given by

$\rho_{env}=\frac{{\rm e}^{-\beta_{L}\,H_{L}}}{{\rm Tr}\Big{(}{\rm
e}^{-\beta_{L}\,H_{L}}\Big{)}}\,\otimes\,\frac{{\rm
e}^{-\beta_{R}\,H_{R}}}{{\rm Tr}\Big{(}{\rm e}^{-\beta_{R}\,H_{R}}\Big{)}}\ .$
(6)

It is invariant under the bath dynamics generated by
$H_{env}=\sum_{\alpha=L,R}H_{\alpha}$ and exhibits thermal expectations of the
form

$\displaystyle\hskip-5.69046pt{\rm
Tr}_{B}\Big{(}\rho_{env}b^{\dagger}_{\alpha}(\nu)b_{\alpha^{\prime}}(\nu^{\prime})\Big{)}=\delta_{\alpha\alpha^{\prime}}\delta(\nu-\nu^{\prime})\,n_{\alpha}(\nu)$
(7) $\displaystyle\hskip-5.69046pt{\rm
Tr}_{B}\Big{(}\rho_{env}b_{\alpha}(\nu)b^{\dagger}_{\alpha^{\prime}}(\nu^{\prime})\Big{)}=\delta_{\alpha\alpha^{\prime}}\delta(\nu-\nu^{\prime})\,(1+n_{\alpha}(\nu)),$
(8)

with thermal mean occupation numbers

$n_{\alpha}(\nu)=\frac{1}{{\rm e}^{\beta_{\alpha}\nu}-1}\ ,\qquad\nu\geq 0\ .$
(9)

Finally, choosing the initial state of the compound system chain plus baths of
the form $\rho_{\rm tot}(0)=\rho(0)\otimes\rho_{env}$, with $\rho(0)$ an
initial state of the $N$ spins of the chain, in presence of a fast decay of
the thermal correlation functions, one applies the weak-coupling limit
techniques and obtains a fully physically consistent dissipative chain
dynamics Alicki-Lendi -Merkli . In practice, the initial state of the compound
system spin-chain plus baths evolves into

$\rho_{\rm tot}(t)={\rm e}^{-itH_{\rm tot}}\,\rho_{\rm tot}(0)\,{\rm
e}^{itH_{\rm tot}}\ ,$ (10)

where $H_{\rm tot}=H+H_{env}+H^{\prime}$ is the total system Hamiltonian. The
state of the open chain at time $t$, $\rho(t)$, is then retrieved by tracing
over the baths degrees of freedom, $\rho(t)={\rm Tr}_{env}\Big{(}\rho_{\rm
tot}(t)\Big{)}$. Then, one rescales the physical time variable $t$ to
$\tau=t\,\lambda^{2}$ and takes the limit $\lambda\to 0$ in $\partial_{t}{\rm
Tr}_{env}\Big{(}\rho_{\rm tot}(t)\Big{)}$. In doing so, too fast oscillations
with respect to the chain transition frequencies $\omega=E_{i}-E_{j}$ are
suppressed, where $E_{i}$ and $|E_{i}\rangle$ solve the spin Hamiltonian
eigenvalue equation $H|E_{i}\rangle=E_{i}\,|E_{i}\rangle$. This procedure
corresponds to the so-called rotating wave approximation leading to a master
equation of the GKSL form

${\partial\rho(t)\over\partial
t}=-i\big{[}H+\lambda^{2}\,H_{LS},\rho(t)]+{\mathbb{D}}[\rho(t)]=\mathbb{L}[\rho(t)]\
.$ (11)

On the right hand side of the above time-evolution equation, one distinguishes
a Hamiltonian term $\lambda^{2}\,H_{LS}$ which provides a Lamb-shift
correction to the spin-chain Hamiltonian $H$ and a purely dissipative term
${\mathbb{D}}[\rho(t)]$. As we shall see, in the specific physical context
here considered, ${\mathbb{D}}[\rho(t)]$ consists of contributions
${\mathbb{D}}^{(\alpha)}_{\omega}[\rho(t)]$ resulting from positive transition
frequencies $\omega\geq 0$, only:

${\mathbb{D}}[\rho(t)]=\lambda^{2}\sum_{\alpha=L,R}\ \sum_{\omega\geq
0}{\mathbb{D}}^{(\alpha)}_{\omega}[\rho(t)]\ .$ (12)

Their explicit form reads

$\displaystyle\hskip-28.45274pt{\mathbb{D}}^{(\alpha)}_{\omega}[\rho(t)]=C^{(\alpha)}_{\omega}\bigg{[}A_{\alpha}(\omega)\rho(t)A_{\alpha}^{\dagger}(\omega)$
$\displaystyle\hskip
56.9055pt-\,\frac{1}{2}\bigg{\\{}A_{\alpha}^{\dagger}(\omega)A_{\alpha}(\omega),\rho(t)\bigg{\\}}\bigg{]}$
(13) $\displaystyle\hskip
28.45274pt+\widetilde{C}^{(\alpha)}_{\omega}\bigg{[}A_{\alpha}^{\dagger}(\omega)\rho(t)A_{\alpha}(\omega)$
$\displaystyle\hskip
56.9055pt-\,\frac{1}{2}\bigg{\\{}A_{\alpha}(\omega)A_{\alpha}^{\dagger}(\omega),\rho(t)\bigg{\\}}\bigg{]}\
,$ (14)

whose coefficients

$\displaystyle C^{(\alpha)}_{\omega}$ $\displaystyle=$ $\displaystyle
2\pi\,[h_{\alpha}(\omega)]^{2}\,\big{(}n_{\alpha}(\omega)+1\big{)}\ ,$ (15)
$\displaystyle\widetilde{C}^{(\alpha)}_{\omega}$ $\displaystyle=$
$\displaystyle 2\pi\,[h_{\alpha}(\omega)]^{2}\,n_{\alpha}(\omega)\ ,$ (16)

with $n_{\alpha}(\omega)$ as in (9), come from the real parts of the half-
Fourier transforms of the bath correlation functions.

Instead, the Lamb-shift correction amounts to the Hamiltonian

$\displaystyle H_{LS}=\sum_{\alpha=L,R}\sum_{\omega}\
\bigg{[}S^{(\alpha)}_{\omega}A_{\alpha}^{\dagger}(\omega)A_{\alpha}(\omega)$
$\displaystyle\hskip
56.9055pt+\,\widetilde{S}^{(\alpha)}_{\omega}A_{\alpha}(\omega)A_{\alpha}^{\dagger}(\omega)\bigg{]}\
,$ (17)

where, unlike the dissipative term, the sum runs now over all positive and
negative transition frequencies and whose coefficients read

$\displaystyle S^{(\alpha)}_{\omega}$ $\displaystyle=$ $\displaystyle
P\,\int_{0}^{+\infty}{\rm
d}\nu\,[h_{\alpha}(\nu)]^{2}\frac{1+n_{\alpha}(\nu)}{\omega-\nu}\ ,$ (18)
$\displaystyle\widetilde{S}^{(\alpha)}_{\omega}$ $\displaystyle=$
$\displaystyle P\,\int_{0}^{+\infty}{\rm
d}\nu\,[h_{\alpha}(\nu)]^{2}\frac{n_{\alpha}(\nu)}{\nu-\omega}\ ,$ (19)

with $P$ denoting the principal value.

In all the previous expressions there appear Lindblad operators of the form

$A^{\dagger}_{\alpha}(\omega)=\sum_{E_{i}-E_{j}=\omega}|E_{i}\rangle\langle
E_{i}|\,\sigma_{+}^{(\alpha)}\,|E_{j}\rangle\langle E_{j}|\ ,$ (20)

together with their Hermitean conjugates

$A_{\alpha}(\omega)=\sum_{E_{i}-E_{j}=\omega}|E_{j}\rangle\langle
E_{j}|\,\sigma_{-}^{(\alpha)}\,|E_{i}\rangle\langle E_{i}|\ .$ (21)

In order to obtain explicit expressions for the elements of the master
equation (11), one needs to work with eigenvalues and eigenvectors of the full
spin Hamiltonian $H$: this point of view is known as global approach to open
quantum spin chains. This way of proceeding is in contrast with the so-called
local approach where the weak-coupling limit is implemented by switching off
the spin interactions, thus obtaining strictly local dissipative terms that
involve only the left and the right spins. The spin interactions are then
reinserted at the end of the weak-coupling procedure.

######  Remark 1.

The fact that the dissipative contribution to the generator,
${\mathbb{D}}[\rho(t)]$, involves only transition frequencies $\omega\geq 0$
is due the thermal bath energies being positive and to the form of the
interaction in (3). Indeed, in the interaction representation, terms as
$A^{\dagger}_{\alpha}(\omega)\mathfrak{b}_{\alpha}(\nu)$ contribute with time
oscillations $\exp(\pm i(\nu-\omega)t)$. On the time scale
$\tau=t/\lambda^{2}$ and in the weak-coupling limit when $\lambda\to 0$, fast
oscillations select contributions with $\omega=\nu\geq 0$. Negative
transitions frequencies, $\omega\leq 0$, would also be selected if in (3)
there were interaction terms of the form
$\sigma^{(\alpha)}_{+}\,\mathfrak{b}^{\dagger}_{\alpha}(\nu)$ which, together
with their Hermitian conjugates, would correspond to the presence of terms of
the form $A^{\dagger}_{\alpha}(\omega)\mathfrak{b}^{\dagger}_{\alpha}(\nu)$,
and Hermitian conjugates, contributing with time oscillations $\exp(\pm
i(\nu+\omega)t)$ .

###  II.1 Spin-chain Hamiltonian: eigenvalues and eigenvectors

In order to address how the presence of the baths modifies the chain dynamics
in the weak-coupling limit and within the global approach, we first need
diagonalize the chain Hamiltonian in (1). By means of the $\ell$-th spin
ladder operators in (4) one rewrites

$H=\Delta\,\sum_{\ell=1}^{N}\sigma^{(\ell)}_{z}\,+\,2g\,\sum_{\ell=1}^{N-1}\Big{(}\sigma^{(\ell)}_{+}\sigma^{(\ell+1)}_{-}\,+\,\sigma^{(\ell)}_{-}\sigma^{(\ell+1)}_{+}\Big{)}\
.$ (22)

By means of the Jordan-Wigner transformation Coleman , one introduces
Fermionic annihilation and creation operators

$a_{j}:=\prod_{k=1}^{j-1}(-\sigma_{z}^{(k)})\,\sigma_{-}^{(j)}\ ,\quad
a^{\dagger}_{j}=\prod_{k=1}^{j-1}(-\sigma_{z}^{(k)})\,\sigma_{+}^{(j)}\ ,$
(23)

with the convention that $\prod_{k=1}^{j-1}(-\sigma^{(k)}_{z})=1$ for $j=1$,
satisfying the anti-commutation relations

$\left\\{a_{j}\,,\,a^{\dagger}_{k}\right\\}=\delta_{jk}\ .$ (24)

Let $|\uparrow\rangle$ and $|\downarrow\rangle$ be the eigenvectors of
$\sigma_{z}$, $\sigma_{z}|\uparrow\rangle=|\uparrow\rangle$,
$\sigma_{z}|\downarrow\rangle=-|\downarrow\rangle$. Since
$\sigma_{-}|\downarrow\rangle=0$, the vacuum vector such that
$a_{j}|vac\rangle=0$, for all $j=1,2,\ldots,N$, amounts to

$|vac\rangle=|\downarrow\rangle^{\otimes N}\ .$ (25)

Using that $\displaystyle
a_{j}^{\dagger}\,a_{j}=\sigma^{(j)}_{+}\sigma^{(j)}_{-}=\frac{1^{(j)}+\sigma^{(j)}_{z}}{2}$,
one inverts the transformation (23):

$\sigma^{(j)}_{-}=(\sigma^{(j)}_{+})^{\dagger}=\prod_{k=1}^{j-1}\Big{(}1-2a^{\dagger}_{k}\,a_{k}\Big{)}\,a_{j}\
,\quad\sigma^{(j)}_{z}=2\,a_{j}^{\dagger}\,a_{j}\,-\,1\ ,$ (26)

finally turning the spin Hamiltonian into a Fermionic one,
$H=-\,N\,\Delta\,+\,2g\,\widetilde{H}$, where

$\widetilde{H}=\gamma\,\sum_{j=1}^{N}a^{\dagger}_{j}a_{j}\,+\,\sum_{j=1}^{N-1}\Big{(}a^{\dagger}_{j}\,a_{j+1}\,+\,a^{\dagger}_{j+1}a_{j}\Big{)}\
,\ \gamma:=\frac{\Delta}{g}\ .$ (27)

As shown in Appendix A, $H$ can then be diagonalized,

$H=-N\,\Delta+\sum_{\ell=1}^{N}\left(2\,\Delta+4g\cos\left(\frac{\ell\pi}{N+1}\right)\right)b^{\dagger}_{\ell}\,b_{\ell}\
,$ (28)

where the operators

$b_{\ell}:=\sum_{j=1}^{N}u_{\ell j}\,a_{j}\ ,\quad
b^{\dagger}_{\ell}:=\sum_{j=1}^{N}u_{\ell j}\,a^{\dagger}_{j}\ ,$ (29)

are also Fermionic, $\left\\{b_{j}\,,\,b^{\dagger}_{k}\right\\}=\delta_{jk}$
with the same vacuum as the operators $a_{j}$: $b_{\ell}|vac\rangle=0$ for all
$\ell=1,2,\ldots,N$, while the coefficients

$u_{\ell k}=\sqrt{\frac{2}{N+1}}\sin\left(\frac{\ell k\pi}{N+1}\right)\ .$
(30)

form an orthogonal and symmetric matrix $U=[u_{k\ell}]$.

In the following, we shall denote by $\mathbf{n}$ the $N$-tuple
$n_{1},n_{2},\ldots,n_{N}$, where $n_{j}=0,1$ is the occupation number of the
$j$-th mode relative to the operators $b_{j}$ and $b_{j}^{\dagger}$. The
eigenvectors of the Hamiltonian (1) have thus the form

$|\mathbf{n}\rangle=(b^{\dagger}_{1})^{n_{1}}(b_{2}^{\dagger})^{n_{2}}\cdots(b_{N}^{\dagger})^{n_{N}}\,|vac\rangle\
.$ (31)

Indeed, according to (31),

$\displaystyle b_{\ell}|\mathbf{n}\rangle$ $\displaystyle=$
$\displaystyle\delta_{n_{\ell},1}\,(-1)^{\sum_{j=1}^{\ell-1}n_{j}}\sqrt{n_{\ell}}\,|\mathbf{n}_{\ell}^{-}\rangle\
,$ (32) $\displaystyle b^{\dagger}_{\ell}|\mathbf{n}\rangle$ $\displaystyle=$
$\displaystyle\delta_{n_{\ell},0}\,(-1)^{\sum_{j=1}^{\ell-1}n_{j}}\sqrt{1-n_{\ell}}\,|\mathbf{n}^{+}_{\ell}\rangle\
,$ (33) $\displaystyle b^{\dagger}_{\ell}\,b_{\ell}|\mathbf{n}\rangle$
$\displaystyle=$ $\displaystyle n_{\ell}\,|\mathbf{n}\rangle\ ,$ (34)

where, $\mathbf{n}_{\ell}^{\pm}$ denote the $N$-tuples
$n_{1},\ldots,n_{\ell}\pm 1,\ldots,n_{N}$. Then, one verifies that
$H\,|\mathbf{n}\rangle=E_{\mathbf{n}}\,|\mathbf{n}\rangle$, where

$E_{\mathbf{n}}=\Delta\Big{(}2\sum_{\ell=1}^{N}n_{\ell}\,-\,N\Big{)}+4g\,\sum_{\ell=1}^{N}n_{\ell}\,\cos\left(\frac{\ell\pi}{N+1}\right)\
.$ (35)

######  Remark 2.

Since the matrix $U=[u_{k\ell}]$ with entries $u_{k\ell}$ as in (30) is real
and symmetric, from (29) one obtains

$a_{j}=\sum_{\ell=1}^{N}u_{j\ell}\,b_{\ell}\ ,\quad
a^{\dagger}_{j}=\sum_{\ell=1}^{N}u_{j\ell}\,b^{\dagger}_{\ell}\ .$ (36)

Then, by using (26), one expresses the Fermionic operators of type $b$ in
terms of spin operators:

$\displaystyle b_{\ell}$ $\displaystyle=$ $\displaystyle\sum_{j=1}^{N}u_{\ell
j}\,\prod_{k=1}^{j-1}(-\sigma^{(k)}_{z})\,\sigma^{(j)}_{-}\ ,$ (37)
$\displaystyle b^{\dagger}_{\ell}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N}u_{\ell
j}\,\prod_{k=1}^{j-1}(-\sigma^{(k)}_{z})\,\sigma^{(j)}_{+}\ .$ (38)

A comparison with known results is provided in Appendix B.

##  III Coupling to external baths

With the notation of the previous section, the Lindblad operators (21) now
read

$A^{\dagger}_{\alpha}(\omega)=\sum_{E_{\mathbf{m}}-E_{\mathbf{n}}=\omega}|\mathbf{m}\rangle\langle\mathbf{m}|\,\sigma_{+}^{(\alpha)}\,|\mathbf{n}\rangle\langle\mathbf{n}|\
.$ (39)

Their explicit form can be derived by expressing the spin operators
$\sigma^{(\alpha)}_{+}$ first in terms of the Fermionic operators
$a_{j},a^{\dagger}_{j}$,

$\sigma^{(L)}_{+}=a^{\dagger}_{1}\
,\quad\sigma^{(R)}_{+}=\prod_{j=1}^{N-1}(1-2\,a^{\dagger}_{j}\,a_{j})\,a_{N}^{\dagger}\
,$ (40)

and then in terms of the operators $b_{\ell},b^{\dagger}_{\ell}$. Using (36),
one immediately derives

$\sigma^{(L)}_{+}=\sum_{\ell=1}^{N}u_{1\ell}\,b^{\dagger}_{\ell}\ ,$ (41)

while the presence of $\prod_{j=1}^{N-1}(1-2\,a^{\dagger}_{j}\,a_{j})$ in the
expression for $\sigma^{(R)}_{+}$ requires some preliminary manipulation.
Firstly, using that
$1-2\,a^{\dagger}_{j}\,a_{j}=\exp(i\pi\,a^{\dagger}_{j}a_{j})$ and that the
relations (36) yield

$\sum_{j=1}^{N-1}\,a^{\dagger}_{j}\,a_{j}=\sum_{j=1}^{N}a^{\dagger}_{j}\,a_{j}\,-\,a^{\dagger}_{N}\,a_{N}=\sum_{\ell=1}^{N}b^{\dagger}_{\ell}\,b_{\ell}\,-\,a^{\dagger}_{N}\,a_{N}\
,$ (42)

one gets

$\displaystyle{\rm
e}^{i\,\pi\,\sum_{j=1}^{N-1}\,a^{\dagger}_{j}\,a_{j}}\,a^{\dagger}_{N}$
$\displaystyle=$ $\displaystyle{\rm
e}^{i\,\pi\,\sum_{\ell=1}^{N}\,b^{\dagger}_{\ell}\,b_{\ell}}\,{\rm
e}^{-i\,\pi\,a^{\dagger}_{N}\,a_{N}}\,a^{\dagger}_{N}$ (43) $\displaystyle=$
$\displaystyle{\rm
e}^{i\,\pi\,\sum_{\ell=1}^{N}\,b^{\dagger}_{\ell}\,b_{\ell}}\,(1\,-\,2a^{\dagger}_{N}\,a_{N})\,a^{\dagger}_{N}$
$\displaystyle=$ $\displaystyle-\,{\rm
e}^{i\,\pi\,\sum_{\ell=1}^{N}\,b^{\dagger}_{\ell}\,b_{\ell}}\,a^{\dagger}_{N}\
,$

whence, finally,

$\sigma^{(R)}_{+}=\,-\,\left({\rm
e}^{i\,\pi\,\sum_{\ell=1}^{N}\,b^{\dagger}_{\ell}\,b_{\ell}}\right)\,\sum_{\ell=1}^{N}u_{N\ell}\,b^{\dagger}_{\ell}\
.$ (44)

By means of (41) and (33), one then computes the transition amplitudes

$\displaystyle\hskip-19.91684pt\langle\mathbf{m}|\sigma^{(L)}_{+}|\mathbf{n}\rangle=\sum_{\ell=1}^{N}(-1)^{\sum_{j=1}^{\ell-1}n_{j}}\,\sqrt{1-n_{\ell}}\,u_{1\ell}\,\delta_{\mathbf{m}\mathbf{n}^{+}_{\ell}}\
.$ (45)

Let $\mathbf{n}_{0_{\ell}}$ respectively $\mathbf{n}_{1_{\ell}}$ denote the
$N$-tuples with fixed digits $n_{\ell}=0$, respectively $n_{\ell}=1$, at site
$\ell$. Then, the only contributing transition amplitudes are

$\displaystyle\langle\mathbf{n}_{1_{\ell}}|\sigma^{(L)}_{+}|\mathbf{n}_{0_{\ell}}\rangle=(-1)^{\sum_{j=1}^{\ell-1}n_{j}}\,\sqrt{1-n_{\ell}}\,u_{1\ell}\
,$ (46)

with $\ell=1,2,\ldots,N$. Also, from (35) and (21), the transition frequencies
associated with such amplitudes are

$\omega_{\ell}=E_{\mathbf{n}_{1_{\ell}}}\,-\,E_{\mathbf{n}_{0_{\ell}}}=2\,\Delta\,+\,4\,g\,\cos\left(\frac{\ell\pi}{N+1}\right)\
,$ (47)

while the corresponding Lindblad operators in (39) read

$A_{L}^{\dagger}(\omega_{\ell})=u_{1\ell}\,\sum_{\widehat{\mathbf{n}}_{\ell}}(-1)^{\sum_{j=1}^{\ell-1}n_{j}}\,|\mathbf{n}_{1_{\ell}}\rangle\langle\mathbf{n}_{0_{\ell}}|\
,$ (48)

where the symbol $\sum_{\widehat{\mathbf{n}}_{\ell}}$ means that the summation
is performed over all binary $2^{N-1}$-tuples of indices $n_{j}=0,1$ with
$j\neq\ell$. It thus follows thats

$\displaystyle A^{\dagger}_{L}(\omega_{\ell})\,A_{L}(\omega_{\ell})$
$\displaystyle=$ $\displaystyle
u^{2}_{1\ell}\,\sum_{\widehat{\mathbf{n}}_{\ell}}\
|\mathbf{n}_{1_{\ell}}\rangle\langle\mathbf{n}_{1_{\ell}}|\ ,$ (49)
$\displaystyle A_{L}(\omega_{\ell})\,A^{\dagger}_{L}(\omega_{\ell})$
$\displaystyle=$ $\displaystyle
u^{2}_{1\ell}\,\sum_{\widehat{\mathbf{n}}_{\ell}}\
|\mathbf{n}_{0_{\ell}}\rangle\langle\mathbf{n}_{0_{\ell}}|\ .$ (50)

In a similar way, from (44), one obtains that the only contributing transition
amplitudes associated to the right bath are

$\langle\mathbf{n}_{1_{\ell}}|\sigma^{(R)}_{+}|\mathbf{n}_{0_{\ell}}\rangle=(-1)^{\sum_{j=1}^{\ell+1}n_{j}}\,\sqrt{1-n_{\ell}}\,u_{N\ell}\
,$ (51)

with Lindblad operators

$A_{R}^{\dagger}(\omega_{\ell})=\,u_{N\ell}\,\sum_{\widehat{\mathbf{n}}_{\ell}}(-1)^{\sum_{j=\ell+1}^{N}n_{j}}\,|\mathbf{n}_{1_{\ell}}\rangle\langle\mathbf{n}_{0_{\ell}}|\
,$ (52)

whence

$\displaystyle A^{\dagger}_{R}(\omega_{\ell})\,A_{R}(\omega_{\ell})$
$\displaystyle=$ $\displaystyle
u^{2}_{N\ell}\,\sum_{\widehat{\mathbf{n}}_{\ell}}\
|\mathbf{n}_{1_{\ell}}\rangle\langle\mathbf{n}_{1_{\ell}}|\ ,$ (53)
$\displaystyle A_{R}(\omega_{\ell})\,A^{\dagger}_{R}(\omega_{\ell})$
$\displaystyle=$ $\displaystyle
u^{2}_{N\ell}\,\sum_{\widehat{\mathbf{n}}_{\ell}}\
|\mathbf{n}_{0_{\ell}}\rangle\langle\mathbf{n}_{0_{\ell}}|\ .$ (54)

Notice that, in the spin representation, all Lindblad operators
$A^{\dagger}_{\alpha}(\omega_{\ell})$ involve, through the relations (26) and
(36), products of all on-site spin operators. This structure is typical of the
global approach to open spin chains and strikingly differs from the local one
which yields Lindblad operators involving only spin operators pertaining to
the first and last spin of the chain.

The operators $A_{\alpha}(\omega_{\ell})$ and
$A^{\dagger}_{\alpha}(\omega_{\ell})$, $\alpha=L,R$ , have to be inserted into
the expressions (13) and (14) when $\omega=\omega_{\ell}\geq 0$ together with
$C^{(\alpha)}_{\omega_{\ell}}$ and $\widetilde{C}^{(\alpha)}_{\omega_{\ell}}$
as in (15) and (16). Instead, the Lamb-shift Hamiltonian (17) contributing to
(11) requires the operators
$A_{\alpha}^{\dagger}(\omega_{\ell})A_{\alpha}(\omega_{\ell})$ and
$A_{\alpha}(\omega_{\ell})A^{\dagger}_{\alpha}(\omega_{\ell})$, $\alpha=L,R$,
with both positive and negative $\omega_{\ell}$. The Hamiltonian
$H+\lambda^{2}\,H_{LS}$ in (11) is thus diagonal in the energy eigenbasis
$\\{|\mathbf{n}\rangle\\}$.

The $N=2$ and $N=3$ cases are explicitly worked out in Appendix C.

######  Remark 3.

Some observations are in order at this point: the first one is that, as a
consequence of the fact that the transition frequencies contributing to the
dissipative generator in (12) are positive, not all those corresponding to the
non-vanishing Lindblad operators $A^{\dagger}_{\alpha}(\omega_{\ell})$ in (48)
and (52) need be such. This means that, the Lindblad operators
$A^{\dagger}_{\alpha}(\omega_{\ell})$ with $\omega_{\ell}<0$ can only
contribute to the Lamb-shift Hamiltonian and not to the dissipative part of
the generator. The sign of $\omega_{\ell}$ depends on the strength of the
inter-spin coupling constant $g$; indeed,

$\cos\left(\frac{\pi\ell}{N+1}\right)<0\quad\hbox{for}\quad
N\geq\ell>\frac{N+1}{2}\ .$

Correspondingly $\omega_{\ell}<0$ for $\displaystyle
g>\frac{\Delta}{2\,\left|\cos\left(\frac{\ell\pi}{N+1}\right)\right|}$.
The second observation is that, should any of the transition frequencies
$\omega_{\ell}$ in the list (47) be negative, the opposite one,
$-\omega_{\ell}=E_{\mathbf{n}_{0_{\ell}}}-E_{\mathbf{n}_{1_{\ell}}}$, not
being in the list, would not give rise to a dissipative contribution of the
form $\mathbb{D}_{-\omega_{\ell}}[\rho(t)]$, as
$A^{\dagger}_{\alpha}(-\omega_{\ell})\equiv 0$. In the following we shall
assume

$g\leq\frac{\Delta}{2\,\cos\left(\frac{\pi}{N+1}\right)}\ ,$ (55)

so that $\omega_{\ell}\geq 0$, $\ell=1,2,\ldots,N$, and leave the study of the
presence of negative transition frequencies for future investigations.

##  IV Stationary state

The master equation (11) possesses a unique stationary state left invariant by
the generated reduced dynamics, namely such that
${\mathbb{L}}[\rho_{\infty}]=0$. This follows from the fact that, as shown in
Appendix D, the only operator commuting with all Lindblad operators
$A_{L,R}(\omega_{\ell})$ and $A^{\dagger}_{L,R}(\omega_{\ell})$ and with the
Hamiltonian must be multiples of the identity Spohn2 – Fagnola2 .

Because of the diagonal form of $H+\lambda^{2}H_{LS}$ in the energy
eigenbasis, $\Big{[}H+\lambda^{2}H_{LS}\,,\,P_{\mathbf{k}}\Big{]}=0$, for all
energy eigenprojections
$P_{\mathbf{k}}:=|\mathbf{k}\rangle\langle\mathbf{k}|$. On the other hand, by
inserting $\rho(t)=P_{\mathbf{k}}$ into (13) and (14), one obtains

$\displaystyle\mathbb{D}^{(L)}_{\omega_{\ell}}[P_{\mathbf{k}}]=u^{2}_{1\ell}\,\Big{(}C^{(L)}_{\omega_{\ell}}\delta_{k_{\ell}1}-\widetilde{C}^{(L)}_{\omega_{\ell}}\delta_{k_{\ell}0}\Big{)}\times$
$\displaystyle\hskip
28.45274pt\times\,\Big{(}P_{\mathbf{k}_{0_{\ell}}}-P_{\mathbf{k}_{1_{\ell}}}\Big{)}$
(56)
$\displaystyle\mathbb{D}^{(R)}_{\omega_{\ell}}[P_{\mathbf{k}}]=u^{2}_{N\ell}\,\Big{(}C^{(R)}_{\omega_{\ell}}\delta_{k_{\ell}1}-\widetilde{C}^{(R)}_{\omega_{\ell}}\delta_{k_{\ell}0}\Big{)}\,\times$
$\displaystyle\hskip
28.45274pt\times\,\Big{(}P_{\mathbf{k}_{0_{\ell}}}-P_{\mathbf{k}_{1_{\ell}}}\Big{)}\
.$ (57)

From (12), using the two previous expressions one finds

$\mathbb{D}\left[P_{\mathbf{k}}\right]=\lambda^{2}\,\sum_{\ell=1}^{N}\Big{(}\delta_{k_{\ell}1}\,d_{\ell}\,-\,\delta_{k_{\ell}0}\,\widetilde{d}_{\ell}\Big{)}\,\Big{(}P_{\mathbf{k}_{0_{\ell}}}\,-\,P_{\mathbf{k}_{1_{\ell}}}\Big{)}\
,$ (58)

where, using (15) and (16),

$\displaystyle\hskip-14.22636ptd_{\ell}$ $\displaystyle=$ $\displaystyle
C^{(L)}_{\omega_{\ell}}\,u^{2}_{1\ell}\,+\,C^{(R)}_{\omega_{\ell}}\,u^{2}_{N\ell}$
(59) $\displaystyle=$
$\displaystyle\frac{4\pi\,}{N+1}\sin^{2}\left(\frac{\pi\ell}{N+1}\right)\,\Big{[}[h_{L}(\omega_{\ell})]^{2}\Big{(}1+n_{L}(\omega_{\ell})\Big{)}$
$\displaystyle+$
$\displaystyle[h_{R}(\omega_{\ell})]^{2}\Big{(}1+n_{R}(\omega_{\ell})\Big{)}\Big{]}$
$\displaystyle\hskip-28.45274pt\widetilde{d}_{\ell}$ $\displaystyle=$
$\displaystyle\widetilde{C}^{(L)}_{\omega_{\ell}}\,u^{2}_{1\ell}\,+\,\widetilde{C}^{(R)}_{\omega_{\ell}}\,u^{2}_{N\ell}$
(60) $\displaystyle=$
$\displaystyle\frac{4\pi\,}{N+1}\sin^{2}\left(\frac{\pi\ell}{N+1}\right)\,\Big{[}[h_{L}(\omega_{\ell})]^{2}\,n_{L}(\omega_{\ell})\,+$
$\displaystyle+$
$\displaystyle[h_{R}(\omega_{\ell})]^{2}\,n_{R}(\omega_{\ell})\Big{]}\ .$

Consider the diagonal expression
$X_{diag}=\sum_{\mathbf{n}}x_{\mathbf{n}}\,P_{\mathbf{n}}$; then, the
dissipator maps it into

$\mathbb{D}[X_{diag}]=\sum_{\mathbf{n}}\widetilde{x}_{\mathbf{n}}\,P_{\mathbf{n}}\
,$ (61)

with

$\widetilde{x}_{\mathbf{n}}={\rm
Tr}\Big{(}P_{\mathbf{n}}\,\mathbb{D}[X_{diag}]\Big{)}=\sum_{\mathbf{m}}\,x_{\mathbf{m}}\,{\rm
Tr}\Big{(}P_{\mathbf{n}}\mathbb{D}\Big{[}P_{\mathbf{m}}\Big{]}\Big{)}\ .$ (62)

From (58) it follows that

$\displaystyle\hskip-5.69046pt\mathbb{D}[X_{diag}]=\lambda^{2}\,\sum_{\mathbf{n}}\,x_{\mathbf{n}}\sum_{\ell=1}^{N}\Big{(}\delta_{n_{\ell}1}\,d_{\ell}\,-\,\delta_{n_{\ell}0}\,\widetilde{d}_{\ell}\Big{)}\,\times$
$\displaystyle\hskip
28.45274pt\times\Big{(}P_{\mathbf{n}_{0_{\ell}}}\,-\,P_{\mathbf{n}_{1_{\ell}}}\Big{)}\,=$
(63)
$\displaystyle\hskip-5.69046pt=\sum_{\ell=1}^{N}\,\sum_{\hat{\mathbf{n}}_{\ell}}\Big{(}d_{\ell}\,x_{\mathbf{n}_{1_{\ell}}}\,-\,\widetilde{d}_{\ell}\,x_{\mathbf{n}_{0_{\ell}}}\Big{)}\,\Big{(}P_{\mathbf{n}_{0_{\ell}}}\,-\,P_{\mathbf{n}_{1_{\ell}}}\Big{)}\
,$ (64)

whence

$\widetilde{x}_{\mathbf{n}}=\lambda^{2}\,\sum_{\ell=1}^{N}\Big{(}\delta_{n_{\ell}0}\,-\,\delta_{n_{\ell}1}\Big{)}\,\Big{(}d_{\ell}\,x_{\mathbf{n}_{1_{\ell}}}\,-\,\widetilde{d}_{\ell}\,x_{\mathbf{n}_{0_{\ell}}}\Big{)}\
.$ (65)

Therefore, $\mathbb{D}[X_{diag}]=0$ is obtained by the factorized expressions

$x_{\mathbf{n}}=\prod_{\ell=1}^{N}\,x^{(\ell)}_{n_{\ell}}\ ,\qquad
x^{(\ell)}_{0}=d_{\ell}\ ,\ x^{(\ell)}_{1}=\widetilde{d}_{\ell}\ .$ (66)

All $x^{(\ell)}\geq 0$ and, after normalization, the uniqueness of the
stationary state together with the expressions (59) and (60) yield

$\displaystyle\rho_{\infty}=\sum_{\mathbf{n}}\Lambda_{\mathbf{n}}\,P_{\mathbf{n}}\
,\quad\Lambda_{\mathbf{n}}=\prod_{\ell=1}^{N}\lambda^{(\ell)}_{n_{\ell}}\ ,$
(67)
$\displaystyle\lambda^{(\ell)}_{n_{\ell}}:=\frac{x^{(\ell)}_{n_{\ell}}}{x^{(\ell)}_{0}+x^{(\ell)}_{1}}=\frac{R^{(\ell)}_{n_{\ell}}}{R_{\ell}}\
,$ (68) $\displaystyle
R^{(\ell)}_{n_{\ell}}:=[h_{L}(\omega_{\ell})]^{2}\Big{(}1-n_{\ell}+n_{L}(\omega_{\ell})\Big{)}\
+$ $\displaystyle\hskip
56.9055pt+\,[h_{R}(\omega_{\ell})]^{2}\Big{(}1-n_{\ell}+n_{R}(\omega_{\ell})\Big{)}\
,$ (69) $\displaystyle
R_{\ell}:=[h_{L}(\omega_{\ell})]^{2}\Big{(}1+2n_{L}(\omega_{\ell})\Big{)}\ +$
$\displaystyle\hskip
56.9055pt+\,[h_{R}(\omega_{\ell})]^{2}\Big{(}1+2n_{R}(\omega_{\ell})\Big{)}\
.$ (70)

With the simplifying assumption $h_{L,R}(\omega_{\ell})=h$, for each
$\ell=1,2,\ldots,N$, one retrieves

$\lambda^{(\ell)}_{n_{\ell}}=\frac{1}{2}\left[1+\frac{(-1)^{n_{\ell}}}{1+n_{L}(\omega_{\ell})+n_{R}(\omega_{\ell})}\right]\
.$ (71)

If we further restrict to identical baths, by imposing equal temperatures and
thus $\beta_{L}=\beta_{R}=\beta$, one computes

$\lambda^{(\ell)}_{n_{\ell}}=\frac{{\rm
e}^{\beta(1-n_{\ell})\omega_{\ell}}}{{\rm e}^{\beta\omega_{\ell}}+1}\ ,$ (72)

so that

$\rho_{\infty}=\sum_{\mathbf{n}}\prod_{\ell=1}^{N}\frac{{\rm
e}^{\beta(1-n_{\ell})\omega_{\ell}}}{{\rm
e}^{\beta\omega_{\ell}}+1}\,|\mathbf{n}\rangle\langle\mathbf{n}|\ .$ (73)

On the other hand, using (35),

$\frac{{\rm e}^{-\beta\,E_{\mathbf{n}}}}{\sum_{\mathbf{n}}{\rm
e}^{-\beta\,E_{\mathbf{n}}}}=\prod_{\ell=1}^{N}\frac{{\rm
e}^{2\beta(1-n_{\ell})\big{(}\Delta+2g\cos(\frac{\pi\ell}{N+1})\big{)}}}{{\rm
e}^{2\beta(1-n_{\ell})\big{(}\Delta+2g\cos(\frac{\pi\ell}{N+1})\big{)}}+1}\ .$
(74)

Then, (47) implies that the open chain stationary state $\rho_{\infty}$ is the
Gibbs state at inverse temperature $\beta$ with Hamiltonian $H$ as given in
(1):

$\rho_{\infty}=\frac{{\rm e}^{-\beta\,H}}{{\rm Tr}\big{(}{\rm
e}^{-\beta\,H}\big{)}}\ .$ (75)

##  V Transport properties

Having determined the explicit, analytic form of the stationary state, we can
now study its transport properties by analyzing the spin and heat flows along
he chain, driven by the two external baths.

###  V.1 Stationary spin flow: sinks and sources

The spin flow at site $k=1,2,\ldots,N$ along the spin chain corresponds to the
rate of change in time of the average of $\sigma^{(k)}_{z}$ given by the
quantity

$\displaystyle\frac{d}{dt}\text{Tr}\big{[}\sigma_{z}^{(k)}\rho(t)\big{]}=\text{Tr}\big{[}\sigma_{z}^{(k)}\mathbb{L}[\rho(t])\big{]}$
$\displaystyle\hskip
85.35826pt=\text{Tr}\big{[}\widetilde{\mathbb{L}}[\sigma^{(k)}_{z}]\,\rho(t)\big{]}\
.$ (76)

In the first equality $\mathbb{L}$ is the generator at the right hand side of
(11), while the second equality follows from the cyclicity of the trace,
$\text{Tr}(XY)=\text{Tr}(YX)$ and defines the so-called dual generator
$\widetilde{\mathbb{L}}$

$\widetilde{\mathbb{L}}[\sigma^{(k)}_{z}]=i\big{[}H+\lambda^{2}H_{LS},\sigma^{(k)}_{z}\big{]}+\widetilde{\mathbb{D}}[\sigma^{(k)}_{z}]\
,$ (77)

with
$\widetilde{\mathbb{D}}[\sigma^{(k)}_{z}]=\lambda^{2}\sum_{\alpha=L,R}\sum_{\omega_{\ell}\geq
0}^{N}\ \widetilde{\mathbb{D}}^{(\alpha)}_{\omega_{\ell}}[\sigma^{(k)}_{z}]$,
where

$\displaystyle\widetilde{\mathbb{D}}^{(\alpha)}_{\omega_{\ell}}[\sigma^{(k)}_{z}]$
$\displaystyle=C^{(\alpha)}_{\omega_{\ell}}\bigg{[}A_{\alpha}^{\dagger}(\omega_{\ell})\sigma^{(k)}_{z}A_{\alpha}(\omega_{\ell})$
$\displaystyle\hskip
28.45274pt-\frac{1}{2}\bigg{\\{}A^{\dagger}_{\alpha}(\omega_{\ell})A_{\alpha}(\omega_{\ell}),\sigma^{(k)}_{z}\bigg{\\}}\bigg{]}$
(78)
$\displaystyle+\widetilde{C}^{(\alpha)}_{\omega_{\ell}}\bigg{[}A_{\alpha}(\omega_{\ell})\sigma^{(k)}_{z}A_{\alpha}^{\dagger}(\omega_{\ell})$
$\displaystyle\hskip
28.45274pt-\frac{1}{2}\bigg{\\{}A_{\alpha}(\omega_{\ell})A_{\alpha}^{\dagger}(\omega_{\ell}),\sigma^{(k)}_{z}\bigg{\\}}\bigg{]}.$
(79)

The Hamiltonian contribution to the rate of change in time of the average of
$\sigma_{z}$ can be expressed in terms of the dimensionless spin currents:

$\displaystyle
J^{(k,k+1)}=4i\Big{(}\sigma_{-}^{(k)}\sigma_{+}^{(k+1)}-\sigma_{+}^{(k)}\sigma_{-}^{(k+1)}\Big{)}\
,$ (80)

as

$i\Big{[}H+\lambda^{2}H_{LS},\sigma_{z}^{(k)}\Big{]}=(g+\kappa)\,\Big{(}J^{(k-1,k)}-J^{(k,k+1)}\Big{)}\
,$ (81)

where the Lamb-shift contribution is characterized by a constant

$\kappa=\frac{i\lambda^{2}}{8\sqrt{2}}\sum_{\alpha=L,R}\sum_{\ell=1}^{N}\,\Big{(}S_{\omega_{\ell}}^{(\alpha)}-\widetilde{S}_{\omega_{\ell}}^{(\alpha)}\Big{)}\
.$ (82)

The operator differences in (81) thus contribute to the continuity equation
(76) as current divergence terms with the right dimension of energy. Since we
are interested in the stationary transport properties, we set
$\rho(t)=\rho_{\infty}$ in the right hand side of (76) and find $\langle
J^{(k,k+1)}\rangle_{\infty}:=\text{Tr}\big{(}\rho_{\infty}\,J^{(k,k+1)}\big{)}=0$.
Indeed, passing from spin to Fermionic operators, by (26) and (36), one finds

$\displaystyle J^{(k,k+1)}$ $\displaystyle=$
$\displaystyle-4i\big{(}a_{k}a^{\dagger}_{k+1}+a_{k}^{\dagger}a_{k+1}\big{)}$
(83) $\displaystyle=$
$\displaystyle-4i\sum_{j\,\ell=1}^{N}u_{kj}u_{k+1\ell}(b_{j}b^{\dagger}_{\ell}+b_{j}^{\dagger}b_{\ell})\
.$

Hence, all their averages with respect to the energy eigenstates vanish,

$\hskip-5.69046pt\langle\mathbf{n}|J^{(k,k+1)}|\mathbf{n}\rangle=-4i\sum_{j=1}^{N}u_{kj}u_{k+1j}=\langle
u_{k}|u_{k+1}\rangle=0\ .$ (84)

Indeed, (32)-(33) yield
$\langle\mathbf{n}|b_{j}b^{\dagger}_{\ell}|\mathbf{n}\rangle=\delta_{j\ell}\,(1-n_{j})$
and
$\langle\mathbf{n}|b^{\dagger}_{j}b_{\ell}|\mathbf{n}\rangle=\delta_{j\ell}\,n_{j}$,
while the columns $|u_{k}\rangle$ of the orthogonal and symmetric matrix $U$
(see Remark 2) are orthogonal. Thus the stationary left, $\langle
J^{(k-1,k)}\rangle_{\infty}$, and right, $\langle
J^{(k,k+1)}\rangle_{\infty}$, spin currents through site $k$ both vanish.

Clearly, being $\rho_{\infty}$ time-independent, the right hand side of (76)
then yields
$\text{Tr}\big{(}\rho_{\infty}\,\widetilde{\mathbb{D}}[\sigma^{(k)}]\big{)}=0$.
However, the left and right purely dissipative contributions,

$\mathfrak{Q}^{(k)}_{\alpha}:=\lambda^{2}\sum_{\ell=1}^{N}\text{Tr}\big{(}\rho_{\infty}\,\widetilde{\mathbb{D}}^{(\alpha)}_{\omega_{\ell}}[\sigma^{(k)}_{z}]\big{)}\
,$ (85)

do not separately vanish; indeed, as shown in Appendix E,

$\displaystyle\mathfrak{Q}^{(k)}_{L}=2\pi\lambda^{2}\sum_{\ell=1}^{N}\,u^{2}_{k\ell}\,u^{2}_{1\ell}\,[h_{L}(\omega_{\ell})]^{2}\,[h_{R}(\omega_{\ell})]^{2}\,\times$
$\displaystyle\hskip
56.9055pt\times\,\frac{n_{L}(\omega_{\ell})-n_{R}(\omega_{\ell})}{R_{\ell}}\
,$ (86)
$\displaystyle\mathfrak{Q}^{(k)}_{R}=2\pi\lambda^{2}\sum_{\ell=1}^{N}\,u^{2}_{k\ell}\,u^{2}_{N\ell}\,[h_{L}(\omega_{\ell})]^{2}\,[h_{R}(\omega_{\ell})]^{2}\,\times$
$\displaystyle\hskip
56.9055pt\times\,\frac{n_{R}(\omega_{\ell})-n_{L}(\omega_{\ell})}{R_{\ell}}\
,$ (87)

with $R_{\ell}$ as in (70). Furthermore, since from (30) one finds that
$u_{1\ell}=(-)^{\ell}u_{N\ell}$ for all $\ell=1,2,\ldots,N$, it follows that
$\mathfrak{Q}^{(k)}_{L}=-\mathfrak{Q}^{(k)}_{R}$ as it should physically be.

Also, assuming $h_{L}(\omega)=h_{R}(\omega)=h$, one gets

$\displaystyle\hskip-14.22636pt\mathfrak{Q}^{(k)}_{L}=\pi\lambda^{2}\sum_{\ell=1}^{N}\,u^{2}_{k\ell}\,u^{2}_{1\ell}\,\frac{n_{L}(\omega_{\ell})-n_{R}(\omega_{\ell})}{1+n_{L}(\omega_{\ell})+n_{R}(\omega_{\ell})}\
,$ (88)
$\displaystyle\hskip-14.22636pt\mathfrak{Q}^{(k)}_{R}=\pi\lambda^{2}\sum_{\ell=1}^{N}\,u^{2}_{k\ell}\,u^{2}_{N\ell}\,\frac{n_{R}(\omega_{\ell})-n_{L}(\omega_{\ell})}{1+n_{L}(\omega_{\ell})+n_{R}(\omega_{\ell})}\
.$ (89)

One thus sees that, while the continuity equation (76) in the stationary case
does not contain any current divergence at site $k$, it does however contain
terms of a different origin that are due to the presence of the two baths.
These terms vanish only if the temperatures are the same so that
$n_{L}(\omega_{\ell})=n_{R}(\omega_{\ell})$ and are thus interpretable as spin
flow source or sink contributions, depending on whether they are positive or
negative.

Figure 1: Behaviour of the source term $\mathfrak{Q}^{(4)}_{R}$ as functions
of $T_{R}$ in a $N=8$ spin chain with left bath temperature $T_{L}=10$,
$\lambda=1$, $\Delta=15,30,50$ and $g$ close to the saturation values in (55).

Note that, due to the scaling as $1/N^{2}$ of the products
$u^{2}_{k\ell}\,u^{2}_{N\ell}$ (see (30)) and the presence of $N$ of them in
(88) and (89), the source and sink terms scale as $1/N$ with increasing number
of spins. In Figure 1 we consider a chain with $N=8$ spins, set $T_{L}=~{}0$
so that $n_{L}(\omega_{\ell})=0$ and show the dependence of the source term

$\mathfrak{Q}^{(4)}_{R}=\pi\lambda^{2}\sum_{\ell=1}^{8}\,u^{2}_{4\ell}\,u^{2}_{8\ell}\,{\rm
e}^{-\beta_{R}\omega_{\ell}}\ ,$ (90)

in the middle of the chain as a function of the right temperature $T_{R}$ and
various values of the transverse magnetic field $\Delta$ and the inter-spin
coupling strength $g$. The values of $g$ associated with $\Delta$ are chosen
close to the bound (55), for reasons that will become clear later when we
discuss the bipartite stationary entanglement.

######  Remark 4.

The presence of sink and source contributions at sites $k\neq 1,N$ is strictly
related to the global structure of the Lindblad operators in (48) and (52)
that involve all spins of the chain. Should the Lindblad operators depend only
on the leftmost and rightmost spin operators as in the local approach to open
spin chains (see Remark 3), sink and source terms would disappear as is the
case for the two spins in Levy . Notice that in the global approach developed
before sinks and sources are present even in the limit where the inter-spin
coupling $g\to 0$; indeed, $g$ appears in the thermal factors
$n_{L,R}(\omega_{\ell})$ through the transition frequencies $\omega_{\ell}$
(see (9)). These terms remain different and non zero whenever
$\beta_{L}\neq\beta_{R}$, even for $g=0$.

###  V.2 Stationary heat flow

Beside the spin flow, the presence of the two baths at the far ends of the
chain also establishes heat flows in and out of the chain. According to
standard quantum thermodynamics arguments Alicki1 ; Spohn1 , the heat flow
through an open quantum system due to its weak coupling to a thermal bath, is
measured by

$\mathfrak{H}(t):=\text{Tr}\left(\frac{{\rm d}\rho(t)}{{\rm
d}t}\,H\right)=\text{Tr}\left(\mathbb{L}[\rho(t)]\,H\right)\ ,$ (91)

where $\rho\mapsto\rho(t)$ is the dissipative evolution due to the bath,
generated by $\mathbb{L}$, while $H$ is the open system time-independent
Hamiltonian. Because of the structure of the GKSL equation as in (11), only
the dissipative term of the generator contributes to the heat flow; therefore,
in the spin chain stationary state, the heat flow due to the left,
respectively right bath is given by

$\mathfrak{H}^{st}_{\alpha}=\sum_{\ell=1}^{N}\text{Tr}\left(\mathbb{D}_{\omega_{\ell}}^{(\alpha)}[\rho_{\infty}]\,H\right)\
,\quad\alpha=L,R\ .$ (92)

Certainly, $\mathbb{D}[\rho_{\infty}]=0$ implies
$\mathfrak{H}^{st}_{L}+\mathfrak{H}^{st}_{R}=0$; however, as for the spin
flow, the single bath contributions to the heat flow need not separately
vanish and their sign, if positive, corresponds to heat flowing into the chain
from the bath, or to heat flowing out of the chain and into the bath.

Using (47), (56), (67)-(70), (15) and (16) one computes

$\displaystyle\mathfrak{H}^{st}_{L}=\sum_{\ell=1}^{N}\text{Tr}\left(\mathbb{D}^{(L)}_{\omega_{\ell}}[\rho_{\infty}]\,H\right)=\lambda^{2}\sum_{\ell=1}^{N}\sum_{\mathbf{n},\mathbf{k}}\,E_{\mathbf{n}}\,\Lambda_{\mathbf{k}}\,u^{2}_{1\ell}\,\times$
$\displaystyle\hskip
28.45274pt\times\,\Big{(}C^{(L)}_{\omega_{\ell}}\delta_{1k_{\ell}}-\widetilde{C}^{(L)}_{\omega_{\ell}}\delta_{0k_{\ell}}\Big{)}\,\langle\mathbf{n}|\Big{(}P_{\mathbf{k}_{0_{\ell}}}-P_{\mathbf{k}_{1_{\ell}}}\Big{)}|\mathbf{n}\rangle$
$\displaystyle\hskip
14.22636pt=\lambda^{2}\sum_{\ell=1}^{N}\omega_{\ell}\,u^{2}_{1\ell}\,\Big{(}C^{(L)}_{\omega_{\ell}}\,\lambda^{(\ell)}_{1}\,-\,\widetilde{C}^{(L)}_{\omega_{\ell}}\,\lambda^{(\ell)}_{0}\Big{)}$
$\displaystyle\hskip
14.22636pt=\lambda^{2}\sum_{\ell=1}^{N}\omega_{\ell}\,u^{2}_{1\ell}\,[h_{L}(\omega_{\ell})]^{2}[h_{R}(\omega_{\ell})]^{2}\,\times$
$\displaystyle\hskip
28.45274pt\times\,\frac{n_{L}(\omega_{\ell})-n_{R}(\omega_{\ell})}{R_{\ell}}\
.$ (93)

Notice that the heat flow is positive, namely it flows from the left bath into
the chain if $n_{L}(\omega)>n_{R}(\omega)$, that is (see (9)) if the left bath
is at higher temperature than the right one. Furthermore, the simplifying
assumption $h_{L}(\omega)=h_{R}(\omega)=h$ yields

$\mathfrak{H}^{st}_{R}=\pi\,\lambda^{2}\,\sum_{\ell=1}^{N}\omega_{\ell}\,u^{2}_{1\ell}\,\frac{n_{R}(\omega_{\ell})-n_{L}(\omega_{\ell})}{1+n_{L}(\omega_{\ell})+n_{R}(\omega_{\ell})}\
.$ (94)

Furthermore, the transition frequencies $\omega_{\ell}$ in (47) are of order
$1$ with respect to increasing $N$, whence each of the $N$ contributions
$\omega_{\ell}\,u^{2}_{1\ell}$ to the heat flow scales as $1/N$ due to (30).
Thus, unlike the sink and source terms in (88) and (89) that scale as $1/N$,
the heat flow does not vanish with increasing $N$. Setting $N=8$ and $T_{L}=0$
as for the source terms in (90), and choosing the same set of parameters
$\Delta$ and $g$ as in Figure 1, the behaviour of the heat flow
$\mathfrak{H}^{st}_{R}$ as a function of $T_{R}$ is reported in Figure 2.

Figure 2: Behaviour of the heat flow $\mathfrak{H}^{st}_{R}$ as a function of
the right bath temperature $T_{R}$ in a $N=8$ spin chain with left bath
temperature $T_{L}=0$, $\lambda=1$, $\Delta=15,30,50$ and $g$ close to the
corresponding saturation values.

##  VI Entanglement properties

Besides transport phenomena, open spin chains represent attractive models of
many-body systems due to their entanglement properties. Indeed, although the
external, transverse magnetic field tends to align all spins in a separable
state, the inter-spin interaction instead is able to generate quantum
correlations among all spins. The presence of the external baths at the chain
end points constitutes interesting additional driving factors influencing the
behaviour of the spin entanglement.

In what follows we shall focus upon the entanglement between any two spins,
$r$ and $s>r$, in the stationary state via the concurrence of the reduced
bipartite density matrix $\rho_{r,s}$ obtained from (67) by tracing over the
spins at sites different from $r$ and $s$. In order to achieve this goal, one
needs to re-express the stationary state in (67) in terms of spin operators,
rather than Fermionic ones.

###  VI.1 Stationary state: spin representation

In this respect, instead of the standard lexicograhic ordering, it proves
convenient to regroup the binary strings $\mathbf{n}$ in terms of the number
of ones they contain. We then introduce the enumeration of the $2^{N}$ binary
$N$-tuples $\mathbf{n}$ known as combinatorial numbering of degree $p$
combinadic , that we shall refer to as combinadic ordering for sake of
shortness. For fixed $p=0,1,2,\ldots,N$, one bijiectively associates to each
of the $N\choose p$ $N$-tuples with $n_{i_{\ell}}=1$ at sites
$i_{1}<\cdots<i_{p}$ the integers

$1\leq\mathcal{N}_{p}=\,1\,+\,\sum_{\ell=1}^{p}{i_{\ell}-1\choose\ell}\leq{N\choose
p}\ ,$ (95)

where the binomial coefficients ${i_{\ell}-1\choose\ell}$ are set to vanish if
$i_{\ell}-1<\ell$. According to such a numbering, we identify $\mathbf{n}$
with a unique $\mathcal{N}_{p}$ for some $p=0,1,\ldots,N$; then, the
stationary state $\rho_{\infty}$ may be written as

$\rho_{\infty}=\sum_{p=0}^{N}\sum_{\mathcal{N}_{p}=1}^{{N\choose
p}}\mathcal{L}^{(p)}_{\mathcal{N}_{p}}\,|\mathcal{N}_{p}\rangle\langle\mathcal{N}_{p}|\
,$ (96)

where $\displaystyle\mathcal{L}^{(p)}_{\mathcal{N}_{p}}$ denotes the
eigenvalue $\Lambda_{\mathbf{n}}$ in (67) corresponding to the binary
$N$-tuple $\mathbf{n}$ with $p$ 1’s, indexed by the combinadic integer
$\mathcal{N}_{p}$. Applications of the above formalism to the $N=2$ and $N=3$
cases can be found in Appendix F.

Notice that, for any fixed $p=0,1,\cdots,N$, the integers $\mathcal{N}_{p}$ in
(95) correspond to the Fermionic excitations of the modes $i_{1}<\cdots<i_{p}$
of type $b$. Indeed, $\mathcal{N}_{p}$ identifies the binary $N$-tuple
$\mathbf{n}$, where $n_{i_{1}}=\cdots=n_{i_{p}}=1$ while the remaining $n_{j}$
vanish. We can thus consistently label:

$|\mathcal{N}_{p}\rangle=b_{i_{1}}^{\dagger}\cdots
b^{\dagger}_{i_{p}}|vac\rangle\ .$ (97)

Then, using (29) one writes

$|\mathcal{N}_{p}\rangle=\sum_{j_{1},j_{2},\ldots,j_{p}}u_{i_{1}j_{1}}u_{i_{2}j_{2}}\cdots
u_{i_{p}j_{p}}\,a^{\dagger}_{j_{1}}a^{\dagger}_{j_{2}}\cdots
a^{\dagger}_{j_{p}}|vac\rangle\ .$ (98)

Notice that, unlike the indices $i_{1},\ldots,i_{p}$, the indices
$j_{1},j_{2},\ldots,j_{p}$ are in general not ordered; their reordering such
that $j_{1}<j_{2}<\cdots<j_{p}$ yields

$|\mathcal{N}_{p}\rangle=\sum_{j_{1}<j_{2}<\ldots<j_{p}}D^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}\,a^{\dagger}_{j_{1}}a^{\dagger}_{j_{2}}\cdots
a^{\dagger}_{j_{p}}|vac\rangle\ ,$ (99)

where

$D^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}=\hbox{det}\Big{(}U^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}\Big{)}\
,$ (100)

is the determinant of the $p\times p$ sub-matrix
$U^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}$ of the orthogonal and symmetric
matrix $U$ in Remark 2 with $p$ rows indexed by $i_{1}<\cdots<i_{p}$ and $p$
columns by $j_{1}<\cdots<j_{p}$. Its entries are thus given by

$\left(U^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}\right)_{ab}=u_{i_{a}j_{b}}=\sqrt{\frac{2}{N+1}}\,\sin\Big(\frac{\pi\,i_{a}j_{b}}{N+1}\Big{missing})\
.$ (101)

######  Remark 5.

For $p=0$, all $n_{j}=0$ whence there are no $i_{\ell}$ and $j_{\ell}$ to
choose and one sets $U^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}=1$. Instead,
if only $n_{i_{\ell}}=1$ then the matrices
$U^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}$ reduce to the scalars
$U^{i_{\ell}}_{j_{\ell}}=u_{i_{\ell}j_{\ell}}$. Finally, if all $n_{j}=1$,
then there is only one contributing matrix, $U^{12\cdots N}_{12\cdots N}=U$,
and ${\rm Det}(U)=-1$. Unlike the matrix $U$ the sub-matrices
$U^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}$ are not symmetric; however,

$D^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}=D^{j_{1}<\cdots<j_{p}}_{i_{1}<\cdots<i_{p}}\
.$ (102)

It is convenient to introduce the ${N\choose p}\times{N\choose p}$ matrices
$\mathcal{D}^{(p)}$, where $\mathcal{D}^{(0)}=1$ and $\mathcal{D}^{(N)}=-1$
are scalars, otherwise $\mathcal{D}^{(p)}$ has entries
$\mathcal{D}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N}^{\prime\prime}_{p}}$
corresponding to the determinants
$D^{i^{\prime}_{1}<\cdots<i^{\prime}_{p}}_{i^{\prime\prime}_{1}<\cdots<i^{\prime\prime}_{p}}$,
where $\mathcal{N}^{\prime}_{p}$ identifies an $N$-tuple with $1$’s at sites
$i^{\prime}_{1}<\cdots<i^{\prime}_{p}$, $\mathcal{N}^{\prime\prime}_{p}$
identifies an $N$-tuple with $1$’s at sites
$i^{\prime\prime}_{1}<\cdots<i^{\prime\prime}_{p}$. Because of (102), the
matrices $\mathcal{D}^{(p)}$ are symmetric for all $p=0,1,\ldots,N$.

In the spin representation $|vac\rangle=|\downarrow\rangle^{\otimes N}$;
therefore, as shown in Appendix G, setting $|\uparrow\rangle=|1\rangle_{S}$
and $|\downarrow\rangle=|0\rangle_{S}$ so that
$\sigma_{+}|0\rangle_{S}=|1\rangle_{S}$, one can express the Fermionic states
$|\mathcal{N}_{p}\rangle$ with $p$ excitations at sites $i_{1}<\cdots<i_{p}$
as linear combinations of the spin states
$|\mathcal{N}^{\prime}_{p}\rangle_{S}$ with $p$ spins flipped up at the sites
$i^{\prime}_{1}<\cdots<i^{\prime}_{p}$ identified by the combinadic index
$\mathcal{N}^{\prime}_{p}$. It follows that, with respect to the standard spin
basis, the stationary state $\rho_{\infty}$ in (67) can be recast as

$\displaystyle\rho_{\infty}$ $\displaystyle=$
$\displaystyle\sum_{p=0}^{N}\,\sum_{\mathcal{N}^{\prime}_{p},\mathcal{N}^{\prime\prime}_{p}}\,\mathcal{S}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N^{\prime}}^{\prime}_{p}}\,|\mathcal{N}^{\prime}_{p}\rangle_{S}{}_{S}\langle\mathcal{N}^{\prime\prime}_{p}|\
,$ (103)
$\displaystyle\mathcal{S}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N}^{\prime\prime}_{p}}$
$\displaystyle:=$
$\displaystyle\sum_{\mathcal{N}_{p}}\,\mathcal{L}_{\mathcal{N}_{p}}\,\mathcal{D}^{(p)}_{\mathcal{N}_{p}\mathcal{N}^{\prime}_{p}}\,\mathcal{D}^{(p)}_{\mathcal{N}_{p}\mathcal{N}^{\prime\prime}_{p}}\
,$ (104)

where $\displaystyle\mathcal{L}_{\mathcal{N}_{p}}$ are the eigenvalues of
$\rho_{\infty}$ as in (67).

From (102) and (104) it follows that, in the standard spin basis, the
stationary state is represented by the block-diagonal matrix

$\mathcal{S}=\bigoplus_{p=0}^{N}\mathcal{S}^{(p)}\ ,\
\mathcal{S}^{(p)}={\mathcal{D}}^{(p)}\,{\mathcal{L}}^{(p)}\,{\mathcal{D}}^{(p)}\
,\ \mathcal{D}=\bigoplus_{p=0}^{N}\mathcal{D}^{(p)}\ ,$ (105)

where $\mathcal{L}^{(p)}$ is the diagonal matrix whose entries are the
eigenvalues in (67) labelled by the combinadic integers $\mathcal{N}_{p}$
while the matrices $\mathcal{D}^{(p)}$ are as in the previous remark.

Finally, again in Appendix G it is shown that the stationary state
$\rho_{\infty}$ has the following structure in terms of spin operators

$\rho_{\infty}=\sum_{p=0}^{N}\,\sum_{\mathcal{N}^{\prime}_{p},\mathcal{N}^{\prime\prime}_{p}}\,\mathcal{S}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N}^{\prime\prime}_{p}}\,\prod_{\ell=1}^{N}\left(X_{n^{\prime}_{\ell}}^{(\ell)}\left(X_{n^{\prime\prime}_{\ell}}^{(\ell)}\right)^{\dagger}\right)\
,$ (106)

where $\displaystyle X_{0}^{(\ell)}=\frac{1-\sigma^{(\ell)}_{z}}{2}$ and
$X^{(\ell)}_{1}=\sigma_{+}^{(\ell)}$, while $n^{\prime}_{\ell}$ and
$n^{\prime\prime}_{\ell}$ are the digits of the binary $N$-tuples with
combinadic indices $\mathcal{N}^{\prime}_{p}$ and
$\mathcal{N}^{\prime\prime}_{p}$.

The above expression of the stationary state $\rho_{\infty}$ can thus be
algorithmically computed for any $N$; the cases $N=2,3$ provide concrete and
informative analytical instances of the above structure and are dealt with in
Appendix H.

###  VI.2 Two-spin entanglement

The spin-operator expression of the stationary state $\rho_{\infty}$ is useful
to investigate the entanglement content of any pair of spins along the chain
and its dependence on their positions $N\geq s>r\geq 1$. We shall quantify the
two-spin entanglement by means of the concurrence Wootters  of the two-spin
reduced density matrix which is obtained by tracing over the spins at sites
different from $r$ and $s$, operation that will be denoted as ${\rm
Tr}_{(r,s)}$. Considering the expression (106) one readily computes

$\displaystyle\hskip-14.22636pt{\rm
Tr}_{(r,s)}\left(\prod_{\ell=1}^{N}\left(X_{n^{\prime}_{\ell}}^{(\ell)}\left(X_{n^{\prime\prime}_{\ell}}^{(\ell)}\right)^{\dagger}\right)\right)=$
$\displaystyle\hskip 14.22636pt=\left(\prod_{\ell\neq r,\ell\neq s}{\rm
Tr}\left(X_{n^{\prime}_{\ell}}^{(\ell)}\,\left(X_{n^{\prime\prime}_{\ell}}^{(\ell)}\right)^{\dagger}\right)\right)\times$
$\displaystyle\hskip-14.22636pt\times\,\left(X_{n^{\prime}_{r}}^{(r)}\,\left(X^{(r)}_{n^{\prime\prime}_{r}}\right)^{\dagger}\right)\,\left(X_{n^{\prime}_{s}}^{(s)}\left(X^{(s)}_{n^{\prime\prime}_{s}}\right)^{\dagger}\right)$
$\displaystyle\hskip-14.22636pt=\left(\prod_{\ell\neq r,\ell\neq
s}\delta_{n^{\prime}_{\ell}n^{\prime\prime}_{\ell}}\right)\
\left(X_{n^{\prime}_{r}}\,X^{\dagger}_{n^{\prime\prime}_{r}}\right)\,\otimes\left(X_{n^{\prime}_{s}}\,X^{\dagger}_{n^{\prime\prime}_{s}}\right)\
,$ (107)

where, in the final two-spin expression, the reference to the spin sites has
safely been neglected. We thus see that the partial trace reduces the double
sum over all possible binary strings $\mathbf{n}^{\prime}$ and
$\mathbf{n}^{\prime\prime}$ in (106) to a double sum over binary strings that
have equal digits but, possibly, for the sites $r$ and $s$. We shall then
denote by $\mathcal{N}^{(rs)}_{p}(n^{\prime}_{r},n^{\prime}_{s})$ and
$\mathcal{N}^{(rs)}_{p}(n^{\prime\prime}_{r},n^{\prime\prime}_{s})$ the
combinadic inidices (95) of the binary strings with $p$ ones that have the
same entries $n_{j}$ everywhere but, possibly, for the sites $r$ and $s$.

With this notation, the two-spin density matrix formally reads

$\displaystyle\rho_{(r,s)}$ $\displaystyle:=$ $\displaystyle{\rm
Tr}_{(r,s)}\big{(}\rho_{\infty}\big{)}$ (108) $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}n^{\prime}_{r},n^{\prime}_{s}\\\
n^{\prime\prime}_{r},n^{\prime\prime}_{s}\end{subarray}}\,\sum_{p=0}^{N}\sum_{\begin{subarray}{c}\mathcal{N}^{(rs)}_{p}(n^{\prime}_{r},n^{\prime}_{s})\\\
\mathcal{N}^{(rs)}_{p}(n^{\prime\prime}_{r},n^{\prime\prime}_{s})\end{subarray}}\,\mathcal{S}^{(p)}_{\mathcal{N}^{(rs)}_{p}(n^{\prime}_{r},n^{\prime}_{s})\mathcal{N}^{(rs)}_{p}(n^{\prime\prime}_{r},n^{\prime\prime}_{s})}\times$
$\displaystyle\times$
$\displaystyle\left(X_{n^{\prime}_{r}}\,X^{\dagger}_{n^{\prime\prime}_{r}}\right)\,\otimes\left(X_{n^{\prime}_{s}}\,X^{\dagger}_{n^{\prime\prime}_{s}}\right)\
.$

Since the $N$-tuples indexed by
$\mathcal{N}^{(rs)}_{p}(n^{\prime\prime}_{r},n^{\prime\prime}_{s})$ have the
same entries but, possibly, for the sites $r$ and $s$, it follows that the
allowed values for $n^{\prime}_{r},n^{\prime}_{s}$ and
$n^{\prime\prime}_{r},n^{\prime\prime}_{s}$ must satisfy
$n^{\prime}_{r}\,+\,n^{\prime}_{s}=n^{\prime\prime}_{r}\,+\,n^{\prime\prime}_{s}$.
These latter ones and the corresponding two-spin operators are as follows:

$\displaystyle\left\\{\begin{matrix}n^{\prime}_{r}=0\
,&n^{\prime\prime}_{r}=0\cr n^{\prime}_{s}=0\
,&n^{\prime\prime}_{s}=0\end{matrix}\right.:\quad\frac{1-\sigma_{z}}{2}\otimes\frac{1-\sigma_{z}}{2}\
,$ (109) $\displaystyle\left\\{\begin{matrix}n^{\prime}_{r}=0\
,&n^{\prime\prime}_{r}=0\cr n^{\prime}_{s}=1\
,&n^{\prime\prime}_{s}=1\end{matrix}\right.:\quad\frac{1-\sigma_{z}}{2}\otimes\frac{1+\sigma_{z}}{2}\
,$ (110) $\displaystyle\left\\{\begin{matrix}n^{\prime}_{r}=0\
,&n^{\prime\prime}_{r}=1\cr n^{\prime}_{s}=1\
,&n^{\prime\prime}_{s}=0\end{matrix}\right.:\quad\sigma_{-}\otimes\sigma_{+}\
,$ (111) $\displaystyle\left\\{\begin{matrix}n^{\prime}_{r}=1\
,&n^{\prime\prime}_{r}=1\cr n^{\prime}_{s}=0\
,&n^{\prime\prime}_{s}=0\end{matrix}\right.:\quad\frac{1+\sigma_{z}}{2}\otimes\frac{1-\sigma_{z}}{2}\
,$ (112)

and

$\displaystyle\left\\{\begin{matrix}n^{\prime}_{r}=1\
,&n^{\prime\prime}_{r}=0\cr n^{\prime}_{s}=0\
,&n^{\prime\prime}_{s}=1\end{matrix}\right.:\quad\sigma_{+}\otimes\sigma_{-}\
,$ (113) $\displaystyle\left\\{\begin{matrix}n^{\prime}_{r}=1\
,&n^{\prime\prime}_{r}=1\cr n^{\prime}_{s}=1\
,&n^{\prime\prime}_{s}=1\end{matrix}\right.:\quad\frac{1+\sigma_{z}}{2}\otimes\frac{1+\sigma_{z}}{2}\
.$ (114)

Therefore, for all sites $1\leq r<s\leq N$, the reduced two-spin density
matrix is a $X$-state for the case of a two-spin chain (see (227) in Appendix
H):

$\displaystyle\hskip-22.76228pt\rho_{(r,s)}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}a&0&0&0\cr 0&b&c&0\cr 0&c&d&0\cr
0&0&0&e\end{pmatrix}$ (115) $\displaystyle=$ $\displaystyle
a\,\frac{1+\sigma_{z}}{2}\otimes\frac{1+\sigma_{z}}{2}\,+\,b\,\frac{1+\sigma_{z}}{2}\otimes\frac{1-\sigma_{z}}{2}$
(117) $\displaystyle+$ $\displaystyle
c\,\sigma_{+}\otimes\sigma_{-}\,+\,c\,\sigma_{-}\otimes\sigma_{+}$ (118)
$\displaystyle+$ $\displaystyle
d\,\frac{1-\sigma_{z}}{2}\otimes\frac{1+\sigma_{z}}{2}\,+\,e\,\frac{1-\sigma_{z}}{2}\otimes\frac{1-\sigma_{z}}{2}\
,$ (119)

with off-diagonal entry

$c=\sum_{p=0}^{N}\sum_{\begin{subarray}{c}\mathcal{N}^{(rs)}_{p}(1,0)\\\
\mathcal{N}^{(rs)}_{p}(0,1)\end{subarray}}\,\mathcal{S}^{(p)}_{\mathcal{N}^{(rs)}_{p}(1,0)\mathcal{N}^{(rs)}_{p}(0,1)}\
.$ (120)

Here the combinadic indices of the entries of $\mathcal{S}^{(p)}$ contributing
to the only off-diagonal term $c$ involve different sites $r$ and $s$. The
combinadic indices are instead the same for the entries of $\mathcal{S}^{(p)}$
contributing to the diagonal entries:

$\displaystyle a$ $\displaystyle=$
$\displaystyle\sum_{p=0}^{N}\sum_{\mathcal{N}^{(rs)}_{p}(1,1)}\,\mathcal{S}^{(p)}_{\mathcal{N}^{(rs)}_{p}(1,1)\mathcal{N}^{(rs)}_{p}(1,1)}\
,$ (121) $\displaystyle b$ $\displaystyle=$
$\displaystyle\sum_{p=0}^{N}\sum_{\mathcal{N}^{(rs)}_{p}(1,0)}\,\mathcal{S}^{(p)}_{\mathcal{N}^{(rs)}_{p}(1,0)\mathcal{N}^{(rs)}_{p}(1,0)\
,}$ (122) $\displaystyle d$ $\displaystyle=$
$\displaystyle\sum_{p=0}^{N}\sum_{\mathcal{N}^{(rs)}_{p}(0,1)}\,\mathcal{S}^{(p)}_{\mathcal{N}^{(rs)}_{p}(0,1)\mathcal{N}^{(rs)}_{p}(0,1)}\
,$ (123) $\displaystyle e$ $\displaystyle=$
$\displaystyle\sum_{p=0}^{N}\sum_{\mathcal{N}^{(rs)}_{p}(0,0)}\,\mathcal{S}^{(p)}_{\mathcal{N}^{(rs)}_{p}(0,0)\mathcal{N}^{(rs)}_{p}(0,0)}\
.$ (124)

For such states the concurrence takes the following analytic expression

$C(r,s)=2\,\max\left\\{0\,,\,\Big{(}|c|-\sqrt{a\,e}\Big{)}\right\\}\ .$ (125)

whence the stationary bipartite entanglement corresponding to a non-vanishing
positive $C(r,s)$, can be evaluated as a function of the sites $r$ and $s$ and
their distance $s-r$. The concurrences $C(1,2)$, $C(2,3)$ and $C(1,3)$ for a
three spin chain are studied in Appendix I.

###  VI.3 Two-spin concurrence

In this section we study the stationary two-spin entanglement in a $N$-spin
chain. In doing so, we use Appendix J which shows how the coefficients
$a,b,c,d$ and $e$ appearing in the concurrence $C(r,s)$ in (125) can be
algorithmically reconstructed. The quantity $C(r,s)$ depends on the parameters
$\Delta$ and $g$ of the chain Hamiltonian, on the temperatures $T_{L,R}$, on
the number of spins, $N$, and on the spin sites $0\leq r\leq s\leq N$.

Firstly, although the algorithm developed in Appendix J works for all $N$, its
algorithmic implementation rapidly becomes time-consuming so that, in the
following figures, we shall focus upon a chain consisting of $N=8$ spins. In
full generality, we observe that, similarly to the sink and source terms in
(86), and (87), the bipartite entanglement between any pair of sites scales as
$1/N$; this follows from the fact that, for large $N$, such is the leading
order of the matrix elements
$\mathcal{S}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N}^{\prime\prime}_{p}}$
in (104). In turn, such a behaviour is due to the fact that the transition
frequencies $\omega_{\ell}$ in (47), and thus the eigenvalues (67), are of
order $1$ with respect to $N$, while the quantities
$\mathcal{D}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N}^{\prime\prime}_{p}}$
introduced in Remark 5 are of order $1/(\sqrt{N})^{p}$ and, in each of the
expressions (120)–(124), there appear sums from $p=1$ to $p=N$ of products of
pairs of such terms.

Secondly, as much as in the case of source and sink terms and of heat flows,
we set $T_{L}=0$ and then inspect the dependence on the right temperature
$T_{R}$ only. What one expects by letting $T_{L}>0$ is that when
$T_{R}=T_{L}>0$ one reaches the Gibbs state in (75). This thermal equilibrium
state can not provide transport effects, for
$n_{L}(\omega_{\ell})=n_{R}(\omega_{\ell})$, but may however support bipartite
entanglement at finite non-vanishing temperatures. On the other hand, for
$T_{L}=T_{R}=0$ the state becomes the vacuum state $|vac\rangle$ in (25) which
is clearly separable. close to the maximum value (55) that ensure the
positivity of all transition frequencies $\omega_{\ell}$ in (47))) and plot
various concurrences versus $T_{R}$.

Figure 3: Maximum achievable entanglement between sites 1 and $s=2,3,\cdots 8$
by varying $T_{R}$. Here $N=8,\lambda=1,T_{L}=0,\Delta=15$ and $g=7.8$ is
close to its upper bound in Eq. (55).

Expected features of the concurrence $C(r,s)$ are that by increasing the
distance $s-r$ between the spins with fixed $r$, the maximum achievable
bipartite entanglement $C_{max}(r,s)$ diminishes, as shown in Figure 3 for
$C_{max}(1,s)$ in a chain of size $N=8$, with $T_{L}=0$, $\lambda=1$,
$\Delta=50$ and $g=7.8$ is very close to upper bound in Eq. (55) while the
concurrence itself vanishes at lower temperatures, in agreement with the fact
that distance and temperature play against correlations. Furthermore, the lack
of translational invariance makes $C(r,s)$ depend not only on $s-r$, but also
on the position $r$ of the first spin.

Figure 4: Bipartite entanglement between spins $3$ and $4$, as measured by the
concurrence $C(3,4)$, versus $T_{R}$ for $N=8$, $T_{L}=0$, with
$\Delta=15,30,50$ and $g$ close to the saturation value relative to
$\Delta=15$. Figure 5: Bipartite entanglement between spins $3$ and $4$, as
measured by the concurrence $C(3,4)$, versus $T_{R}$ for $N=8$, $T_{L}=0$,
with $\Delta=15,30,50$ and $g$ close to the corresponding saturation values.

As regards the dependence of the concurrence on the parameters $\Delta$ and
$g$, Figure 4 first shows that, with temperature $T_{L}=0$, and $g$ fixed,
close to the saturation the bound (55) for $\Delta=15$, the entanglement as a
function of $T_{R}$ diminishes while increasing $\Delta$. This behaviour
agrees with the fact that augmenting the transverse external field the spins
tend to become all parallel and thus the stationary state separable. On the
other hand, by increasing $g$ the spins interact more strongly thus favouring
the generation of quantum correlations that may persist asymptotically against
temperature. In fact, the farther is $g$ away from the saturation value at
given $\Delta$, the smaller is the achieved entanglement, the larger is the
temperature at which it appears and the smaller the one at which it
disappears. Specifically, Indeed, the chosen value of $g$ is sufficient to
generate entanglement for $\Delta=15,30$, but not for $\Delta=50$, while
increasing $g$ beyond the saturation value for $\Delta=15$ would violate the
condition assumed throughout the manuscript that the transition frequencies
$\omega_{\ell}$ in (47) be positive.

Instead, in Figure 5, under the same conditions as in Figure 4, the values of
the interaction strength $g$ are chosen close to the saturation bound (55) for
each $\Delta$. The graph shows that in this case the highest possible $g$,
despite the higher values of $\Delta$, contributes to the creation of
entanglement as soon as $T_{R}>0$; moreover, it also makes it last up to
higher values of $T_{R}$.

##  VII Discussion

Spin chains coupled to external baths at their endpoints represent
paradigmatic models for the study of transport properties in quantum many-body
systems, as they allow the precise analysis of the behaviour of spin and heat
flows along the chain. So far, analytic treatments of the asymptotic transport
properties in these systems have been obtained assuming only ad hoc couplings
between the system and the external baths, those that allow expressing the
chain steady state in terms of the so-called “matrix product states” and the
like Prosen5 .

Here instead, taking an arbitrary, energy preserving coupling to the end
baths, we have been able to derive an exact analytic expression for the unique
steady state of a generic $N$-sites spin-1/2 chain, with $XX$-type inter-spin
interaction, in a transverse constant magnetic field. This has allowed
discussing in detail the open transport properties of the model, treated in
the so-called global approach, revealing the presence of sink and source terms
in the spin flow continuity equation, never pointed out before (except in the
$N=3$ case BFM ).

In addition, having the explicit form of the system asymptotic stationary
state allowed analyzing the entanglement properties of the chain. In
particular, a procedure has been devised able to algoritmically provide the
explicit expression of the reduced two-spin density matrix for any two sites
along the chain. The behaviour of the corresponding entanglement content of
the reduced state, as measured by concurrence, has been discussed in some
relevant cases in terms of the parameters of the system Hamiltonian and of the
bath temperatures. While increasing bath temperatures and magnitude of the
external magnetic field counteracts entanglement, sufficiently high values of
the inter-spin interaction coupling constant would always allow the presence
of asymptotic entanglement among any couple of close enough sites.

In addition, these results show that, for generic $N$, there is no apparent
relation between the behaviour of heat flow and two-site entanglement as a
function of the bath temperatures, as claimed in the literature for the
special case $N=2$ Khandelwal : entanglement in the chain is generated
independently from the heat flow and even in absence of it, as in the case of
isothermal baths. Furthermore, the two quantities behave rather differently
with respect to the length of the chain: while the concurrence vanishes as
$1/N$, the heat flow does not.

We are confident that these findings will stimulate further research on the
use of many-body systems, and spin-chains in particular, for modelling quantum
transport processes, in view of possible applications in quantum techonology.

##  Appendix A Diagonalization of the spin-chain Hamiltonian

Let us consider the Hamiltonian

$\widetilde{H}=\gamma\,\sum_{j=1}^{N}a^{\dagger}_{j}a_{j}\,+\,\sum_{j=1}^{N-1}\Big{(}a^{\dagger}_{j}\,a_{j+1}\,+\,a^{\dagger}_{j+1}a_{j}\Big{)}$

and recast it as
$\widetilde{H}=\sum_{j,k=1}^{N}h_{jk}\,a^{\dagger}_{j}\,a_{k}$ where

$h:=\Big{[}h_{jk}\Big{]}=\begin{pmatrix}\gamma&1&0&0&\ldots&0&0\cr
1&\gamma&1&0&\ldots&0&0\cr
0&1&\gamma&1&\ldots&0&0\cr\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\ldots&\dot{\hfil}&\dot{\hfil}\cr\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\ldots&\dot{\hfil}&\dot{\hfil}\cr\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\ldots&\dot{\hfil}&\dot{\hfil}\cr
0&0&0&0&\ldots&\gamma&1\cr
0&0&0&0&\ldots&1&\gamma\end{pmatrix}=\gamma\,\mathbb{I}\,+\,B\ .$ (126)

The rank $N$ tridiagonal matrix

$B=\begin{pmatrix}0&1&0&0&\ldots&0&0\cr 1&0&1&0&\ldots&0&0\cr
0&1&0&1&\ldots&0&0\cr\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\ldots&\dot{\hfil}&\dot{\hfil}\cr\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\ldots&\dot{\hfil}&\dot{\hfil}\cr\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\dot{\hfil}&\ldots&\dot{\hfil}&\dot{\hfil}\cr
0&0&0&0&\ldots&0&1\cr 0&0&0&0&\ldots&1&0\end{pmatrix}\ ,$ (127)

can be diagonalized as follows. We shall emphasize the rank by writing $B_{N}$
instead of $B$. Then, one notices that

$\hbox{det}\Big{(}x-B_{N}\Big{)}=x\,\hbox{det}\Big{(}x-B_{N-1}\Big{)}-\hbox{det}\Big{(}x-B_{N-2}\Big{)}\
,$ (128)

whence the same equation is satisfied by the associated characteristic
polynomial $p_{N}(x)$:

$p_{N+2}(x)\,-\,x\,p_{N+1}(x)\,+\,p_{N}(x)=0\ .$ (129)

Setting $x=2\cos\theta$, one finds the solution

$p_{N}(x)=C_{1}\,{\rm e}^{i\,N\,\theta}\,+\,C_{2}\,{\rm e}^{-i\,N\,\theta}\ .$
(130)

From $p_{1}(x)=x$ and $p_{2}(x)=x^{2}-1$, one fixes the coefficients

$C_{1}=\frac{{\rm e}^{i\theta}}{{\rm e}^{i\theta}\,-\,{\rm e}^{-i\theta}}\
,\qquad C_{2}=-\frac{{\rm e}^{-i\theta}}{{\rm e}^{i\theta}\,-\,{\rm
e}^{-i\theta}}\ ,$ (131)

whence $\displaystyle p_{N}(x)=\frac{\sin((N+1)\theta)}{\sin\theta}$.

Since for $\theta=0$, respectively $\theta=\pi$, $p_{N}(2)=N+1$, respectively
$p_{N}(-2)=(-)^{N}(N+1)$, the only zeroes of $p_{N}(x)$ are at
$\displaystyle\theta_{k}=\frac{k\,\pi}{N+1}$, $k=1,2,\ldots,N$. It thus
follows that the eigenvalues of $h=\gamma\,\mathbb{I}+B$ are

$\lambda_{k}=\gamma\,+\,2\,\cos\left(\frac{k\,\pi}{N+1}\right)\ .$ (132)

Finally, one can check that the symmetric matrix

$U=\Big{[}u_{jk}\Big{]}_{j,k=1}^{N}\ ,\
u_{jk}=\sqrt{\frac{2}{N+1}}\,\sin\left(\frac{j\,k\,\pi}{N+1}\right)\ .$ (133)

is also orthogonal. It can be directly checked that

$U\,B\,U={\rm
diag}\left[2\,\cos\left(\frac{\ell\pi}{N+1}\right)\right]_{\ell=1}^{N}\ .$
(134)

##  Appendix B Energies and eigenvectors of two and three spin chains

###  B.1 Two-spin chain

For a chain consisting of two spins only as in Levy – Haack2 , $N=2$ and (35)
yield the energies,

$E_{00}=-2\,\Delta,\ E_{10}=2\,g,\ E_{01}=-2\,g,\ E_{11}=2\,\Delta\ .$ (135)

Furthermore, from (133), the unitary matrix $U$ results

$U=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1\cr 1&-1\end{pmatrix}\ .$ (136)

Therefore, applying (37) one obtains

$\displaystyle b_{1}$
$\displaystyle=\frac{\sigma^{(1)}_{-}-\sigma^{(1)}_{z}\sigma^{(2)}_{-}}{\sqrt{2}}\
,\ b_{2}=\frac{\sigma^{(1)}_{-}+\sigma^{(1)}_{z}\sigma^{(2)}_{-}}{\sqrt{2}}\
,$ (137) $\displaystyle b_{1}b_{2}$
$\displaystyle=-\sigma^{(1)}_{-}\sigma^{(2)}_{-}\ ,$ (138)

so that using (31) and (25), one can recast the eigenvectors
$|n_{1}n_{2}\rangle$ relative to the eigenvalues $E_{n_{1}n_{2}}$ using the
standard basis, $|\uparrow\rangle$, $|\downarrow\rangle$. Indeed,
$|vac\rangle=|00\rangle=|\downarrow\downarrow\rangle$, whence

$\displaystyle|00\rangle$ $\displaystyle=|\downarrow\downarrow\rangle\ ,$
$\displaystyle|10\rangle$
$\displaystyle=\frac{|\uparrow\downarrow\rangle+|\downarrow\uparrow\rangle}{\sqrt{2}}$
(139) $\displaystyle|01\rangle$
$\displaystyle=\frac{|\uparrow\downarrow\rangle-|\downarrow\uparrow\rangle}{\sqrt{2}}\
,$ $\displaystyle|11\rangle$ $\displaystyle=|\uparrow\uparrow\rangle\ .$ (140)

###  B.2 Three spin chain

In the case of a three-spin chains BFM , setting $N=3$ the eigenvalues of the
Hamiltonian (1) are

$E_{n_{1}n_{2}n_{3}}=\Delta\,\Big{(}2\,\sum_{\ell=1}^{3}\,n_{\ell}-3\Big{)}+4g\sum_{\ell=1}^{3}n_{\ell}\,\cos\left(\frac{\ell\pi}{4}\right)\
.$ (141)

Explicitly they and their corresponding eigenvectors read

$\displaystyle|000\rangle$ $\displaystyle=|vac\rangle\ ,\ $ $\displaystyle
E_{000}$ $\displaystyle=-3\Delta\ ,$ (142) $\displaystyle|100\rangle$
$\displaystyle=b^{\dagger}_{1}|000\rangle\ ,\ $ $\displaystyle E_{100}$
$\displaystyle=-\,\Delta\,+\,2\sqrt{2}g\ ,$ (143) $\displaystyle|010\rangle$
$\displaystyle=b^{\dagger}_{2}|000\rangle\ ,\ $ $\displaystyle E_{010}$
$\displaystyle=-\,\Delta\ ,$ (144) $\displaystyle|001\rangle$
$\displaystyle=b^{\dagger}_{3}|000\rangle\ ,$ $\displaystyle E_{001}$
$\displaystyle=-\,\Delta\,-\,2\sqrt{2}g\ ,$ (145) $\displaystyle|110\rangle$
$\displaystyle=b^{\dagger}_{1}b^{\dagger}_{2}|000\rangle\ ,$ $\displaystyle
E_{110}$ $\displaystyle=\Delta\,+\,2\sqrt{2}g\ ,$ (146)
$\displaystyle|101\rangle$
$\displaystyle=b^{\dagger}_{1}b^{\dagger}_{3}|000\rangle\ ,$ $\displaystyle
E_{101}$ $\displaystyle=\Delta\ ,$ (147) $\displaystyle|011\rangle$
$\displaystyle=b^{\dagger}_{2}b^{\dagger}_{3}|000\rangle\ ,$ $\displaystyle
E_{011}$ $\displaystyle=\Delta\,-\,2\sqrt{2}g\,$ (148)
$\displaystyle|111\rangle$
$\displaystyle=b^{\dagger}_{1}b^{\dagger}_{2}b^{\dagger}_{3}|000\rangle\ ,$
$\displaystyle E_{111}$ $\displaystyle=3\,\Delta\ .$ (149)

The correspondence with the eigenvalues $E_{j}$ in BFM  is as follows

$\displaystyle E_{000}=E_{2}\ ,\ E_{100}=E_{7}\ ,\ E_{010}=E_{4}\ ,\
E_{001}=E_{6}\ ,$ $\displaystyle E_{110}=E_{5}\ ,\ E_{101}=E_{3}\ ,\
E_{011}=E_{8}\ ,\ E_{111}=E_{1}\ .$

Since the $3\times 3$ matrix $U$ in (133) reads

$U=\frac{1}{2}\begin{pmatrix}1&\sqrt{2}&1\cr\sqrt{2}&0&-\sqrt{2}\cr
1&-\sqrt{2}&1\end{pmatrix}\ .$ (150)

Then, $\sigma_{\pm}\sigma_{z}=\mp\sigma_{\pm}$ and (38) yield

$\displaystyle b^{\dagger}_{1}$
$\displaystyle=\frac{1}{2}\Big{(}\sigma^{(1)}_{+}\,-\,\sqrt{2}\,\sigma^{(1)}_{z}\,\sigma^{(2)}_{+}\,+\,\sigma^{(1)}_{z}\,\sigma^{(2)}_{z}\,\sigma^{(3)}_{+}\Big{)}$
(151) $\displaystyle b^{\dagger}_{2}$
$\displaystyle=\frac{1}{\sqrt{2}}\Big{(}\sigma^{(1)}_{+}\,-\,\sigma^{(1)}_{z}\,\sigma^{(2)}_{z}\,\sigma^{(3)}_{+}\Big{)}$
(152) $\displaystyle b^{\dagger}_{3}$
$\displaystyle=\frac{1}{2}\Big{(}\sigma^{(1)}_{+}\,+\,\sqrt{2}\,\sigma^{(1)}_{z}\,\sigma^{(2)}_{+}\,+\,\sigma^{(1)}_{z}\,\sigma^{(2)}_{z}\,\sigma^{(3)}_{+}\Big{)}$
(153) $\displaystyle b^{\dagger}_{1}b^{\dagger}_{2}$
$\displaystyle=-\frac{1}{2}\Big{(}\sigma^{(1)}_{+}\,\sigma^{(2)}_{+}\,-\,\sqrt{2}\,\sigma^{(1)}_{+}\,\sigma^{(2)}_{z}\,\sigma^{(3)}_{+}\,+\,\sigma^{(2)}_{+}\,\sigma^{(3)}_{+}\Big{)}$
(154) $\displaystyle b^{\dagger}_{2}b^{\dagger}_{3}$
$\displaystyle=-\frac{1}{2}\Big{(}\sigma^{(1)}_{+}\,\sigma^{(2)}_{+}\,+\,\sqrt{2}\,\sigma^{(1)}_{+}\,\sigma^{(2)}_{z}\,\sigma^{(3)}_{+}\,+\,\sigma^{(2)}_{+}\,\sigma^{(3)}_{+}\Big{)}$
(155) $\displaystyle b^{\dagger}_{1}b^{\dagger}_{3}$
$\displaystyle=-\frac{1}{\sqrt{2}}\Big{(}\sigma^{(1)}_{+}\,\sigma^{(2)}_{+}\,-\,\sigma^{(2)}_{+}\,\sigma^{(3)}_{+}\Big{)}$
(156) $\displaystyle b^{\dagger}_{1}b^{\dagger}_{2}b^{\dagger}_{3}$
$\displaystyle=-\sigma^{(1)}_{+}\,\sigma^{(2)}_{+}\,\sigma^{(3)}_{+}\ .$ (157)

The expressions of the eigenvectors $|n_{1}n_{2}n_{3}\rangle$ in the spin
standard basis and their correspondence with the eigenvectors obtained in BFM
are reported in Appendix C.2.

##  Appendix C Lindblad operators for two and three spin chains

###  C.1 Two-spin chain

In the case of $N=2$, from (47), one computes the following transition
frequencies

$\omega_{1}=2(\Delta+g)\ ,\ \omega_{2}=2(\Delta-g)\ .$ (158)

Using (48) and (52), the following Lindblad operators ensue for the open two-
spin chain:

$\displaystyle A_{L}^{\dagger}(\omega_{1})$
$\displaystyle=\frac{|10\rangle]\langle 00|+|11\rangle\langle 01|}{\sqrt{2}}\
,$ (159) $\displaystyle A_{L}^{\dagger}(\omega_{2})$
$\displaystyle=\frac{|01\rangle]\langle 00|-|11\rangle\langle 10|}{\sqrt{2}}\
,$ (160) $\displaystyle A_{R}^{\dagger}(\omega_{1})$
$\displaystyle=\frac{|10\rangle]\langle 00|-|11\rangle\langle 0|}{\sqrt{2}}\
,$ (161) $\displaystyle A_{R}^{\dagger}(\omega_{2})$
$\displaystyle=-\frac{|01\rangle]\langle 00|+|11\rangle\langle 10|}{\sqrt{2}}\
.$ (162)

According to the discussion before Remark 3, in order to have all of them
contribute to dissipation, one must set $g\leq\Delta$.

###  C.2 Three-spin chain

Setting $N=3$, from (47) we get the three frequencies

$\omega_{1}=2\Big{(}\Delta+\sqrt{2}g\Big{)},\ \omega_{2}=2\Delta,\
\omega_{3}=2\Big{(}\Delta-\sqrt{2}g\Big{)}\ .$ (163)

They correspond to the three frequencies $\omega_{1}$, $\omega_{0}$ and
$\omega_{2}$ in BFM . Furthermore, (48) and (52) yield the left-bath Lindblad
operators

$\displaystyle A^{\dagger}_{L}(\omega_{1})$
$\displaystyle=\frac{1}{2}\Big{(}|100\rangle\langle
000|\,+\,|110\rangle\langle 010|\,+\,|101\rangle\langle 001|$
$\displaystyle\hskip 28.45274pt+\,|111\rangle\langle 011|\Big{)}$ (164)
$\displaystyle A^{\dagger}_{L}(\omega_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\Big{(}|010\rangle\langle
000|\,-\,|110\rangle\langle 100|\,-\,|111\rangle\langle 101|$
$\displaystyle\hskip 28.45274pt+\,|011\rangle\langle 001|\Big{)}$ (165)
$\displaystyle A^{\dagger}_{L}(\omega_{3})$
$\displaystyle=\frac{1}{2}\Big{(}|001\rangle\langle
000|\,-\,|101\rangle\langle 100|\,-\,|011\rangle\langle 010|$
$\displaystyle\hskip 28.45274pt+\,|111\rangle\langle 110|\Big{)}\ ,$ (166)

and the right bath Lindblad operators

$\displaystyle A^{\dagger}_{R}(\omega_{1})$
$\displaystyle=\frac{1}{2}\Big{(}|011\rangle\langle
001|\,+\,|111\rangle\langle 101|\,-\,|010\rangle\langle 000|$
$\displaystyle\hskip 28.45274pt-\,|110\rangle\langle 100|\Big{)}$ (167)
$\displaystyle A^{\dagger}_{R}(\omega_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\Big{(}|011\rangle\langle
001|\,+\,|111\rangle\langle 101|\,-\,|010\rangle\langle 000|$
$\displaystyle\hskip 28.45274pt-\,|110\rangle\langle 100|\Big{)}$ (168)
$\displaystyle A^{\dagger}_{R}(\omega_{3})$
$\displaystyle=\frac{1}{2}\Big{(}|001\rangle\langle
000|\,+\,|101\rangle\langle 100|\,+\,|011\rangle\langle 010|$
$\displaystyle\hskip 28.45274pt+\,|111\rangle\langle 110|\Big{)}\ .$ (169)

According to the discussion before Remark 3, in order to have all of them
contribute to dissipation, one must set $\displaystyle
g\leq\frac{\Delta}{\sqrt{2}}$.

In terms of the eigenstates $|\uparrow\rangle$ and $|\downarrow\rangle$ of
$\sigma_{z}$, using that
$|vac\rangle=|000\rangle=|\downarrow\downarrow\downarrow\rangle$, one can
reexpress the eigenvectors $|n_{1}n_{2}n_{3}\rangle$ in the spin standard
basis:

$|000\rangle=|\downarrow\downarrow\downarrow\rangle=|E_{2}\rangle\ ,$ (170)

in the case of zero excitations, while for one excitation,

$\displaystyle|100\rangle$
$\displaystyle=b^{\dagger}_{1}|000\rangle=\frac{1}{2}\Big{(}|\uparrow\downarrow\downarrow\rangle\,+\,\sqrt{2}\,|\downarrow\uparrow\downarrow\rangle\,+\,|\downarrow\downarrow\uparrow\rangle\Big{)}\
,$ $\displaystyle=|E_{7}\rangle$ (171) $\displaystyle|010\rangle$
$\displaystyle=b^{\dagger}_{2}|000\rangle=\frac{1}{\sqrt{2}}\Big{(}|\uparrow\downarrow\downarrow\rangle\,-\,|\downarrow\downarrow\uparrow\rangle\Big{)}=|E_{4}\rangle\
,$ (172) $\displaystyle|001\rangle$
$\displaystyle=b^{\dagger}_{3}|000\rangle=\frac{1}{2}\Big{(}|\uparrow\downarrow\downarrow\rangle\,-\,\sqrt{2}\,|\downarrow\uparrow\downarrow\rangle\,+\,|\downarrow\downarrow\uparrow\rangle\Big{)}$
$\displaystyle=|E_{6}\rangle\ ,$ (173)

and for two excitations

$\displaystyle|110\rangle$
$\displaystyle=b^{\dagger}_{1}b^{\dagger}_{2}|000\rangle=-\frac{1}{2}\Big{(}|\uparrow\uparrow\downarrow\rangle\,+\,\sqrt{2}\,|\uparrow\downarrow\uparrow\rangle\,+\,|\downarrow\uparrow\uparrow\rangle\Big{)}\
,$ $\displaystyle=-|E_{5}\rangle\ ,$ (174) $\displaystyle|101\rangle$
$\displaystyle=b^{\dagger}_{1}b^{\dagger}_{3}|000\rangle=-\frac{1}{\sqrt{2}}\Big{(}|\uparrow\uparrow\downarrow\rangle\,-\,|\downarrow\uparrow\uparrow\rangle\Big{)}\
,$ $\displaystyle=-|E_{3}\rangle$ (175) $\displaystyle|011\rangle$
$\displaystyle=b^{\dagger}_{2}b^{\dagger}_{3}|000\rangle=-\frac{1}{2}\Big{(}|\uparrow\uparrow\downarrow\rangle\,-\,\sqrt{2}\,|\uparrow\downarrow\uparrow\rangle\,+\,|\downarrow\uparrow\uparrow\rangle\Big{)}$
$\displaystyle=-|E_{8}\rangle\ .$ (176)

Finally, in the case of three excitations one finds

$|111\rangle=b^{\dagger}_{1}b^{\dagger}_{2}b^{\dagger}_{3}|000\rangle=-|\uparrow\uparrow\uparrow\rangle=-|E_{1}\rangle\
,$ (177)

the difference in the overall sign depending on the chosen ordering of the
creation operators $b^{\dagger}_{\ell}$. The correspondence with the
$|E_{k}\rangle$ obtained in BFM  is confirmed after noticing that there
$|0\rangle=|\uparrow\rangle$ and $|1\rangle=|\downarrow\rangle$.

In the same vain, using (170)– (177), one recasts the Lindblad operators in
(164)– (169) as

$\displaystyle A^{\dagger}_{L}(\omega_{1})$
$\displaystyle=\frac{1}{2}\Big{(}|E_{7}\rangle\langle
E_{2}|\,-\,|E_{5}\rangle\langle E_{4}|\,-\,|E_{3}\rangle\langle E_{6}|$
$\displaystyle\hskip 28.45274pt+\,|E_{1}\rangle\langle E_{8}|\Big{)}\ ,$ (178)
$\displaystyle A^{\dagger}_{L}(\omega_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\Big{(}|E_{4}\rangle\langle
E_{2}|\,+\,|E_{5}\rangle\langle E_{7}|\,-\,|E_{1}\rangle\langle E_{3}|$
$\displaystyle\hskip 28.45274pt-\,|E_{8}\rangle\langle E_{6}|\Big{)}\ ,$ (179)
$\displaystyle A^{\dagger}_{L}(\omega_{3})$
$\displaystyle=\frac{1}{2}\Big{(}|E_{6}\rangle\langle
E_{2}|\,+\,|E_{3}\rangle\langle E_{7}|\,+\,|E_{8}\rangle\langle E_{4}|$
$\displaystyle\hskip 28.45274pt+\,|E_{1}\rangle\langle E_{5}|\Big{)}\ ,$ (180)

and

$\displaystyle A^{\dagger}_{R}(\omega_{1})$
$\displaystyle=\frac{1}{2}\Big{(}|E_{7}\rangle\langle
E_{2}|\,+\,|E_{5}\rangle\langle E_{4}|\,+\,|E_{3}\rangle\langle E_{6}|$
$\displaystyle\hskip 28.45274pt+\,|E_{1}\rangle\langle E_{8}|\Big{)}$ (181)
$\displaystyle A^{\dagger}_{R}(\omega_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\Big{(}|E_{5}\rangle\langle
E_{7}|\,-\,|E_{4}\rangle\langle E_{2}|\,-\,|E_{8}\rangle\langle E_{6}|$
$\displaystyle\hskip 28.45274pt+\,|E_{1}\rangle\langle E_{3}|\Big{)}$ (182)
$\displaystyle A^{\dagger}_{R}(\omega_{3})$
$\displaystyle=\frac{1}{2}\Big{(}|E_{6}\rangle\langle
E_{2}|\,+\,|E_{1}\rangle\langle E_{5}|\,-\,|E_{3}\rangle\langle E_{7}|$
$\displaystyle\hskip 28.45274pt-\,|E_{8}\rangle\langle E_{4}|\Big{)}\ .$ (183)

These expressions coincide with those found in BFM .

##  Appendix D Uniqueness of the stationary state

Consider the commutator of a generic chain operator
$X=\sum_{\mathbf{p},\mathbf{q}}\,X_{\mathbf{p}\mathbf{q}}\,|\mathbf{p}\rangle\langle\mathbf{q}|$,
in the energy eigenbasis (34), with the Lindblad operator in (48) and set it
equal to zero:

$\displaystyle\langle\mathbf{r}|\left[X\,,\,A^{\dagger}_{L}(\omega_{\ell})\right]|\mathbf{s}\rangle=u_{1\ell}\,\sum_{\mathbf{p},\mathbf{q},\widehat{\mathbf{n}}_{\ell}}(-1)^{\sum_{j=1}^{\ell-1}n_{j}}\,X_{\mathbf{p}\mathbf{q}}\,\times$
$\displaystyle\hskip
28.45274pt\times\langle\mathbf{r}|\Big{[}|\mathbf{p}\rangle\langle\mathbf{q}|\,,\,|\mathbf{n}_{1_{\ell}}\rangle\langle\mathbf{n}_{0_{\ell}}|\Big{]}|\mathbf{s}\rangle=0\
.$

This yields

$X_{\mathbf{r}\mathbf{s}_{1_{\ell}}}(-1)^{\sum_{j=1}^{\ell-1}s_{j}}\,\delta_{s_{\ell}0}\,=\,X_{\mathbf{r}_{0_{\ell}}\mathbf{s}}(-1)^{\sum_{j=1}^{\ell-1}r_{j}}\,\delta_{r_{\ell}1}\
,$ (184)

whence, choosing $r_{\ell}=0$ and $s_{\ell}=0$ gives
$X_{\mathbf{r}_{0_{\ell}}\mathbf{s}_{1_{\ell}}}=0$. Changing
$A^{\dagger}_{L}(\omega_{\ell})$ into $A^{\dagger}_{R}(\omega_{\ell})$ yields
$X_{\mathbf{r}_{1_{\ell}}\mathbf{s}_{0_{\ell}}}=0$; thus, the only non-
vanishing $X$ commuting with all Lindblad operators must be diagonal in the
energy eigenbasis: namely,
$X=\sum_{\mathbf{n}}\,X_{\mathbf{n}\mathbf{n}}|\mathbf{n}\rangle\langle\mathbf{n}|$.
On the other hand, choosing $r_{\ell}=1$ and $s_{\ell}=0$, (184) yields
$X_{\mathbf{n}_{1_{\ell}}\mathbf{n}_{1_{\ell}}}=X_{\mathbf{n}_{0_{\ell}}\mathbf{n}_{0_{\ell}}}$
for all $\ell=1,2,\ldots,N$, whence $X$ must be a multiple of the identity.

##  Appendix E Sink and source contributions to the stationary transport
properties

Given the stationary state in (67), in order to compute

$\mathfrak{Q}^{(k)}_{L}=\lambda^{2}\sum_{\ell=1}^{N}\text{Tr}\big{(}\rho_{\infty}\,\widetilde{\mathbb{D}}^{(L)}_{\omega_{\ell}}[\sigma^{(k)}_{z}]\big{)}\
,$ (185)

in (86), we need evaluate mean-values of the form

$\displaystyle\langle\mathbf{n}|\widetilde{\mathbb{D}}^{(L)}_{\omega_{\ell}}[\sigma^{(k)}_{z}]|\mathbf{n}\rangle=C^{(L)}_{\omega_{\ell}}\,\left(\langle\mathbf{n}|A^{\dagger}_{L}(\omega_{\ell})\sigma^{(k)}_{z}A_{L}(\omega_{\ell})|\mathbf{n}\rangle\right.$
$\displaystyle\hskip
28.45274pt\left.-\frac{1}{2}\langle\mathbf{n}|\left\\{A^{\dagger}_{L}(\omega_{\ell})A_{L}(\omega_{\ell}),\sigma^{(k)}_{z}\right\\}|\mathbf{n}\rangle\right)$
(186) $\displaystyle\hskip
71.13188pt+\widetilde{C}^{(L)}_{\omega_{\ell}}\,\left(\langle\mathbf{n}|A_{L}(\omega_{\ell})\sigma^{(k)}_{z}A^{\dagger}_{L}(\omega_{\ell})|\mathbf{n}\rangle\right.$
$\displaystyle\hskip
28.45274pt\left.-\frac{1}{2}\langle\mathbf{n}|\left\\{A_{L}(\omega_{\ell})A^{\dagger}_{L}(\omega_{\ell}),\sigma^{(k)}_{z}\right\\}|\mathbf{n}\rangle\right)\
.$ (187)

Using the expressions (48), (53) and (50) one gets

$\displaystyle\hskip-17.07182pt\langle\mathbf{n}|A^{\dagger}_{L}(\omega_{\ell})\sigma^{(k)}_{z}A_{L}(\omega_{\ell})|\mathbf{n}\rangle=u_{1\ell}^{2}\,\delta_{1n_{\ell}}\,\langle\mathbf{n}_{0_{\ell}}|\sigma^{(k)}_{z}|\mathbf{n}_{0_{\ell}}\rangle,$
(188)
$\displaystyle\hskip-17.07182pt\langle\mathbf{n}|\left\\{A^{\dagger}_{L}(\omega_{\ell})A_{L}(\omega_{\ell}),\sigma^{(k)}_{z}\right\\}|\mathbf{n}\rangle=u_{1\ell}^{2}\,\delta_{1n_{\ell}}\,\langle\mathbf{n}|\sigma^{(k)}_{z}|\mathbf{n}\rangle,$
(189)
$\displaystyle\hskip-17.07182pt\langle\mathbf{n}|A_{L}(\omega_{\ell})\sigma^{(k)}_{z}A^{\dagger}_{L}(\omega_{\ell})|\mathbf{n}\rangle=u_{1\ell}^{2}\,\delta_{0n_{\ell}}\,\langle\mathbf{n}_{1_{\ell}}|\sigma^{(k)}_{z}|\mathbf{n}_{1_{\ell}}\rangle,$
(190)
$\displaystyle\hskip-17.07182pt\langle\mathbf{n}|\left\\{A_{L}(\omega_{\ell})A^{\dagger}_{L}(\omega_{\ell}),\sigma^{(k)}_{z}\right\\}|\mathbf{n}\rangle=u_{1\ell}^{2}\,\delta_{0n_{\ell}}\,\langle\mathbf{n}|\sigma^{(k)}_{z}|\mathbf{n}\rangle.$
(191)

Then, from (26) and (36) it follows that

$\sigma^{(k)}_{z}=-1\,+\,2\,\sum_{r,s=1}^{N}u_{kr}\,u_{ks}\,b^{\dagger}_{r}\,b_{s}\
,$ (192)

whence, using (32)– (34),

$\displaystyle\hskip-14.22636pt\langle\mathbf{n}|\sigma^{(k)}_{z}|\mathbf{n}\rangle$
$\displaystyle=$ $\displaystyle 2\,\sum_{r=1}^{N}\,n_{r}\,u^{2}_{kr}\,-\,1\ ,$
(193)
$\displaystyle\hskip-14.22636pt\langle\mathbf{n}_{0_{\ell}}|\sigma^{(k)}_{z}|\mathbf{n}_{0_{\ell}}\rangle$
$\displaystyle-$
$\displaystyle\langle\mathbf{n}|\sigma^{(k)}_{z}|\mathbf{n}\rangle=-2\,n_{\ell}\,u^{2}_{k\ell}\
,$ (194)
$\displaystyle\hskip-14.22636pt\langle\mathbf{n}_{1_{\ell}}|\sigma^{(k)}_{z}|\mathbf{n}_{1_{\ell}}\rangle$
$\displaystyle-$
$\displaystyle\langle\mathbf{n}|\sigma^{(k)}_{z}|\mathbf{n}\rangle=2(1-n_{\ell})\,u^{2}_{k\ell}\
.$ (195)

Inserting the last three expressions into (186)– (191) yields

$\hskip-8.5359pt\langle\mathbf{n}|\widetilde{\mathbb{D}}^{(L)}_{\omega_{\ell}}[\sigma^{(k)}_{z}]|\mathbf{n}\rangle=u^{2}_{k\ell}u^{2}_{1\ell}\Big{(}\widetilde{C}^{(L)}_{\omega_{\ell}}\delta_{0n_{\ell}}-C^{(L)}_{\omega_{\ell}}\delta_{1n_{\ell}}\Big{)},$
(196)

whence, finally using the stationary state eigenvalues in (67)– (70) and the
explicit form of the constants $C^{(L)}_{\omega_{\ell}}$ and
$\widetilde{C}^{(L)}_{\omega_{\ell}}$ in (15), (16), the sink/source
contribution $\mathfrak{Q}^{(k)}_{L}$ in (86) ensues. Similar arguments lead
to $\mathfrak{Q}^{(k)}_{R}$ in (87).

##  Appendix F Stationary state eigenvalues for two and three spin chains

###  F.1 Two-spin chain

For two spins and $N=2$, the four $2$-digit strings are, in anti-lexicographic
order, $(00),(10),(01)$ and $(11)$. Moreover, using (95) their combinadic
ordering shows to be the same; namely,

$(00)\leftrightarrow\mathcal{N}_{0}=1\ ,\
\left\\{\begin{matrix}(10)\leftrightarrow{\mathcal{N}}_{1}=1\cr(01)\leftrightarrow\mathcal{N}_{1}=2\end{matrix}\right.\
,\ (11)\leftrightarrow\mathcal{N}_{2}=1\ .$ (197)

Then, the combinadic list of the eigenvalues $\Lambda_{n_{1}n_{2}}$ of the
stationary state $\rho_{\infty}$ in (96) are

$\displaystyle\mathcal{L}^{(0)}_{1}$
$\displaystyle=\Lambda_{00}=\frac{R^{(1)}_{0}R^{(2)}_{0}}{R_{1}R_{2}}\ ,\
\hbox{for $p=0$}\ ,$ (198) $\displaystyle\mathcal{L}^{(1)}_{1}$
$\displaystyle=\Lambda_{10}=\frac{R^{(1)}_{1}R^{(2)}_{0}}{R_{1}R_{2}}\ ,\
\mathcal{L}^{(1)}_{2}=\Lambda_{01}=\frac{R^{(1)}_{0}R^{(2)}_{1}}{R_{1}R_{2}}$
(199)

for $p=1$ and, for $p=2$,

$\mathcal{L}^{(2)}_{1}=\Lambda_{11}=\frac{R^{(1)}_{1}R^{(2)}_{1}}{R_{1}R_{2}}\
.$ (200)

###  F.2 Three-spin chains

For $N=3$, let us consider the anti-lexicograhic ordering of the eight
$3$-digit strings:

$(000),(100),(010),(110),(001),(101),(011),(111)\ .$

Application of (95) shows that the previous one and the combinadic ordering
coincide, in the sense that

$\displaystyle(000)\leftrightarrow\mathcal{N}_{0}=1\ ,\
\left\\{\begin{matrix}(100)\leftrightarrow{\mathcal{N}}_{1}=1\cr(010)\leftrightarrow\mathcal{N}_{1}=2\cr(001)\leftrightarrow\mathcal{N}_{1}=3\end{matrix}\right.\
,$ (201)
$\displaystyle\left\\{\begin{matrix}(110)\leftrightarrow{\mathcal{N}}_{2}=1\cr(101)\leftrightarrow\mathcal{N}_{2}=2\cr(011)\leftrightarrow\mathcal{N}_{2}=3\end{matrix}\right.\
,\ (111)\leftrightarrow\mathcal{N}_{3}=1\ .$ (202)

Then, for a three spin chain, the eigenvalues $\Lambda_{n_{1}n_{2}n_{3}}$ of
the stationary state $\rho_{\infty}$ in (96) are

$\mathcal{L}^{(0)}_{1}=\Lambda_{000}=\frac{R^{(1)}_{0}R^{(2)}_{0}R^{(3)}_{0}}{R^{(1)}R^{(2)}R^{(3)}}\
,$ (203)

for $p=0$, while for $p=1$

$\displaystyle\mathcal{L}^{(1)}_{1}$ $\displaystyle=$
$\displaystyle\Lambda_{100}=\frac{R^{(1)}_{1}R^{(2)}_{0}R^{(3)}_{0}}{R_{1}R_{2}R_{3}}\
,$ (204) $\displaystyle\mathcal{L}^{(1)}_{2}$ $\displaystyle=$
$\displaystyle\Lambda_{010}=\frac{R^{(1)}_{0}R^{(2)}_{1}R^{(3)}_{0}}{R_{1}R_{2}R_{3}}\
,$ (205) $\displaystyle\mathcal{L}^{(1)}_{3}$ $\displaystyle=$
$\displaystyle\Lambda_{001}=\frac{R^{(1)}_{0}R^{(2)}_{0}R^{(3)}_{1}}{R_{1}R_{2}R_{3}}\
,$ (206)

for $p=2$,

$\displaystyle\mathcal{L}^{(2)}_{1}$ $\displaystyle=$
$\displaystyle\Lambda_{110}=\frac{R^{(1)}_{1}R^{(2)}_{1}R^{(3)}_{0}}{R_{1}R_{2}R_{3}}\
,$ (207) $\displaystyle\mathcal{L}^{(2)}_{2}$ $\displaystyle=$
$\displaystyle\Lambda_{101}=\frac{R^{(1)}_{1}R^{(2)}_{0}R^{(3)}_{1}}{R_{1}R_{2}R_{3}}\
,$ (208) $\displaystyle\mathcal{L}^{(2)}_{3}$ $\displaystyle=$
$\displaystyle\Lambda_{011}=\frac{R^{(1)}_{0}R^{(2)}_{1}R^{(3)}_{1}}{R_{1}R_{2}R_{3}}\
,$ (209)

and, finally, for $p=3$,

$\mathcal{L}^{(3)}_{1}=\Lambda_{111}=\frac{R^{(1)}_{1}R^{(2)}_{1}R^{(3)}_{1}}{R_{1}R_{2}R_{3}}\
.$ (210)

The quantities

$\displaystyle R^{(\ell)}_{0}$ $\displaystyle=$
$\displaystyle[h_{L}(\omega_{\ell})]^{2}\Big{(}1+n_{L}(\omega_{\ell})\Big{)}$
(211) $\displaystyle\hskip
14.22636pt+[h_{R}(\omega_{\ell})]^{2}\Big{(}1+n_{R}(\omega_{\ell})\Big{)}\qquad\hbox{and}$
$\displaystyle R^{(\ell)}_{1}$ $\displaystyle=$
$\displaystyle[h_{L}(\omega_{\ell})]^{2}\,n_{L}(\omega_{\ell})+[h_{R}(\omega_{\ell})]^{2}\,n_{R}(\omega_{\ell})$
(212)

in (69) coincide with the quantities $s_{\ell}$, respectively $\tau_{\ell}$ in
BFM . According to the correspondence of the frequencies $\omega_{1,2,3}$ in
(163) with the frequencies $\omega_{1,0,2}$ in BFM , the following
correspondence arises among the stationary state eigenvalues
$\Lambda_{\mathbf{n}}$ and the eigenvalues $\mu_{k}$ in BFM :

$\displaystyle\hskip-14.22636pt\mathcal{L}^{(0)}_{1}=\mu_{2}\ ,\
\mathcal{L}^{(1)}_{1}=\mu_{7}\ ,\ \mathcal{L}^{(1)}_{2}=\mu_{4}\ ,\
\mathcal{L}^{(1)}_{3}=\mu_{6}\ ,$ (213)
$\displaystyle\hskip-14.22636pt\mathcal{L}^{(2)}_{1}=\mu_{5}\ ,\
\mathcal{L}^{(2)}_{2}=\mu_{3}\ ,\ \mathcal{L}^{(2)}_{3}=\mu_{8}\ ,\
\mathcal{L}^{(3)}_{1}=\mu_{1}\ .$ (214)

##  Appendix G Stationary state in the spin representation

Given the combinadic labeling (97) of the energy eigenstates in the Fermionic
representation, from (23) one gets

$\displaystyle|\mathcal{N}_{p}\rangle=\sum_{j_{1}<j_{2}<\ldots<j_{p}}D^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}\
\left(\prod_{k=1}^{j_{1}-1}(-\sigma^{(k)}_{z})^{p}\right)\,\times$
$\displaystyle\hskip
5.69046pt\times\left(-\sigma^{(j_{1})}_{+}\sigma^{j_{1}}_{z})\,\right)\,\left(\prod_{k=1}^{j_{2}-1}(-\sigma^{(k)}_{z})^{p-1}\right)\,\left(-\sigma^{(j_{2})}_{+}\sigma^{j_{2}}_{z})\,\right)\cdots\,\times$
$\displaystyle\hskip
5.69046pt\times\,\cdots\left(-\sigma_{+}^{(j_{p-1}}\sigma^{(j_{p-1})}_{z}\right)\,\left(\prod_{k=1}^{j_{1}-1}(-\sigma^{(k)}_{z})\right)\,\sigma^{(j_{p})}_{+}|vac\rangle\
.$ (215)

Since $\sigma_{+}\sigma_{z}=-\sigma_{+}$, one rewrites

$|\mathcal{N}_{p}\rangle_{=}\sum_{j_{1}<\ldots<j_{p}}D^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}\
S_{j_{1}\cdots j_{p}}|vac\rangle\ ,$ (216)

where

$\displaystyle\hskip-28.45274ptS_{j_{1}\cdots
j_{p}}|vac\rangle=\left(\prod_{k=1}^{j_{1}-1}(-\sigma^{(k)}_{z})^{p}\right)\,\sigma^{(j_{1})}_{+}\,\times$
$\displaystyle\hskip
28.45274pt\times\,\left(\prod_{k=j_{1}+1}^{j_{2}-1}(-\sigma^{(k)}_{z})^{p-1}\right)\,\sigma^{(j_{2})}_{+}\,\cdots\times$
(217)
$\displaystyle\hskip-14.22636pt\times\cdots\sigma_{+}^{(j_{p-1})}\,\left(\prod_{k=j_{p-1}+1}^{j_{p}-1}(-\sigma^{(k)}_{z})\right)\,\sigma^{(j_{p})}_{+}|vac\rangle\
.$ (218)

Because of (25), with $|\uparrow\rangle=|1\rangle_{S}$ and
$|\downarrow\rangle=|0\rangle_{S}$, one writes

$S_{j_{1}\cdots j_{p}}|vac\rangle=|\mathcal{N}^{\prime}_{p}\rangle_{S}\ ,$
(219)

where, according to (95), $\mathcal{N}^{\prime}_{p}$ uniquely identifies the
$N$-tuple with $n^{\prime}_{j_{1}}=\cdots=n^{\prime}_{j_{p}}=1$ and
$n^{\prime}_{k}=0$ otherwise: such an $N$-tuple corresponds in its turn to a
spin state vector with spin up at the sites $j_{1}<\cdots<j_{p}$, and down at
all other ones. Thus with respect to the standard spin basis, the stationary
states can be recast as in (103) and (104).

In order to rewrite $\rho_{\infty}$ as a tensor product of on-site spin
operators, one starts from the projectors

$\displaystyle|\mathcal{N}_{p}\rangle\langle\mathcal{N}_{p}|$ $\displaystyle=$
$\displaystyle\sum_{j_{1}<\ldots<j_{p}}\sum_{k_{1}<\cdots<k_{p}}D^{i_{1}<\cdots<i_{p}}_{j_{1}<\cdots<j_{p}}\,D^{i_{1}<\cdots<i_{p}}_{k_{1}<\cdots<k_{p}}\,\times$
(220) $\displaystyle\times$ $\displaystyle S_{j_{1}\cdots
j_{p}}\,|vac\rangle\langle vac|\,S^{\dagger}_{k_{1}\cdots k_{p}}\ .$

Then, writing $\displaystyle|vac\rangle\langle
vac|=\prod_{\ell=1}^{N}\frac{1-\sigma^{(\ell)}_{z}}{2}$ and using that

$\sigma_{+}\frac{1-\sigma_{z}}{2}=\sigma_{+}\
,\quad(-\sigma_{z})^{k}\frac{1-\sigma_{z}}{2}=\frac{1-\sigma_{z}}{2}\quad\forall
k\ ,$

one recast $S_{j_{1}\cdots j_{p}}|vac\rangle\langle vac|$ as

$\displaystyle\left(\prod_{k=1}^{j_{1}-1}\frac{1-\sigma^{(k)}_{z}}{2}\right)\,\sigma^{(j_{1})}_{+}\,\left(\prod_{k=j_{1}+1}^{j_{2}-1}\frac{1-\sigma^{(k)}_{z}}{2}\right)\,\sigma^{(j_{2})}_{+}\,\cdots\times$
$\displaystyle\hskip
14.22636pt\times\cdots\,\sigma_{+}^{(j_{p-1})}\,\left(\prod_{k=j_{p-1}+1}^{j_{p}-1}\frac{1-\sigma^{(k)}_{z}}{2}\right)\,\sigma^{(j_{p})}_{+}\
.$ (221)

As done before, using (95), we identify any given set of indices
$j_{1}<\cdots<j_{p}$ and the corresponding $N$-tuple with
$n_{j_{1}}=\cdots=n_{j_{p}}$ by the unique combinadic integer
$\mathcal{N}^{\prime}_{p}$, whence

$\displaystyle\hskip-14.22636ptS_{j_{1}\cdots j_{p}}|vac\rangle\langle
vac|=\mathcal{Z}^{(p)}_{\mathcal{N}^{\prime}_{p}}\ ,$ (222)
$\displaystyle\hskip-14.22636pt|\mathcal{N}_{p}\rangle\langle\mathcal{N}_{p}|=\sum_{\mathcal{N}^{\prime}_{p},\mathcal{N}^{\prime\prime}_{p}}\,\mathcal{D}^{(p)}_{\mathcal{N}_{p}\mathcal{N}^{\prime}_{p}}\,\mathcal{D}^{(p)}_{\mathcal{N}_{p}\mathcal{N}^{\prime\prime}_{p}}\,\mathcal{Z}^{(p)}_{\mathcal{N}^{\prime}_{p}}\,\left(\mathcal{Z}^{(p)}_{\mathcal{N}^{\prime\prime}_{p}}\right)^{\dagger},$
(223)

where, setting $\displaystyle X_{0}^{(\ell)}=\frac{1-\sigma^{(\ell)}_{z}}{2}$
and $X^{(\ell)}_{1}=\sigma_{+}^{(\ell)}$,

$\mathcal{Z}^{(p)}_{\mathcal{N}^{\prime}_{p}}\,\left(\mathcal{Z}^{(p)}_{\mathcal{N}^{\prime\prime}_{p}}\right)^{\dagger}=\prod_{\ell=1}^{N}\left(X_{n^{\prime}_{\ell}}^{(\ell)}\left(X_{n^{\prime\prime}_{\ell}}^{(\ell)}\right)^{\dagger}\right)\
.$ (224)

It thus follows that in spin operatorial form, the stationary state reads as
in (106).

##  Appendix H Spin representation of the stationary state for two and three
spin chains

###  H.1 Two-spin chain

The case of a two-spin chain is the simplest: as there can be at most two
spins up, the values of $p$ are $0,1,2$. Then, $\mathcal{D}^{(0)}=1$,
$\mathcal{D}^{(2)}=-1$ and $\mathcal{D}^{(1)}=U$, whence, using (198)– (200),
from (105) one gets $\mathcal{S}^{(0)}=\Lambda_{00}$,
$\mathcal{S}^{(2)}=\Lambda_{11}$ and

$\displaystyle\mathcal{S}^{(1)}$ $\displaystyle=$ $\displaystyle
U\begin{pmatrix}\Lambda_{10}&0\cr 0&\Lambda_{01}\end{pmatrix}U$ (225)
$\displaystyle=$
$\displaystyle\frac{1}{2}\begin{pmatrix}\Lambda_{10}+\Lambda_{01}&\Lambda_{10}-\Lambda_{01}\cr\Lambda_{10}-\Lambda_{01}&\Lambda_{10}+\Lambda_{01}\end{pmatrix}\
.$

Finally, using (198)– (200), (106) yields

$\displaystyle\rho_{\infty}$ $\displaystyle=\Lambda_{00}P_{-}\otimes
P_{+}+\Lambda_{11}P_{+}\otimes P_{+}$
$\displaystyle+\frac{\Lambda_{10}+\Lambda_{01}}{2}\Big{(}P_{+}\otimes
P_{-}\,+\,P_{-}\otimes P_{+}\Big{)}$
$\displaystyle+\frac{\Lambda_{10}-\Lambda_{01}}{2}\Big{(}\sigma_{+}\otimes\sigma_{+}\,+\,\sigma_{-}\otimes\sigma_{+}\Big{)}\
.$ (226)

The stationary state in the standard spin basis $|\uparrow\uparrow\rangle$,
$|\uparrow\downarrow\rangle$, $|\downarrow\uparrow\rangle$ and
$|\downarrow\downarrow\rangle$, explicitly reads

$\rho_{\infty}=\frac{1}{2}\begin{pmatrix}2\Lambda_{11}&0&0&0\cr\cr
0&\Lambda_{10}+\Lambda_{01}&\Lambda_{10}-\Lambda_{01}&0\cr\cr
0&\Lambda_{10}-\Lambda_{01}&\Lambda_{10}+\Lambda_{01}&0\cr\cr
0&0&0&2\Lambda_{00}\end{pmatrix}$ (227)

and is thus a so-called $X$-state. Its entanglement content is measured by the
concurrence which has, in this special case, the analytic expression ,

$C(\rho_{\infty})=\max\Big{(}0\,,\,\Big{|}\Lambda_{10}-\Lambda_{01}\Big{|}\,-\,2\,\sqrt{\Lambda_{00}\Lambda_{11}}\Big{)}\
.$ (228)

In the case where $h_{L,R}(\omega_{1,2})=h$, from (71) one gets

$\displaystyle\Lambda_{00}$
$\displaystyle=\frac{1}{4}\frac{2+N_{LR}(\omega_{1})}{1+N_{LR}(\omega_{1})}\frac{2+N_{LR}(\omega_{2})}{1+N_{LR}(\omega_{2})}$
$\displaystyle\Lambda_{10}$
$\displaystyle=\frac{1}{4}\frac{N_{LR}(\omega_{1})}{1+N_{LR}(\omega_{1})}\frac{2+N_{LR}(\omega_{2})}{1+N_{LR}(\omega_{2})}$
$\displaystyle\Lambda_{01}$
$\displaystyle=\frac{1}{4}\frac{2+N_{LR}(\omega_{1})}{1+N_{LR}(\omega_{1})}\frac{N_{LR}(\omega_{2})}{1+N_{LR}(\omega_{2})}$
$\displaystyle\Lambda_{11}$
$\displaystyle=\frac{1}{4}\frac{N_{LR}(\omega_{1})}{1+N_{LR}(\omega_{1})}\frac{N_{LR}(\omega_{2})}{1+N_{LR}(\omega_{2})}\
,$

where we set $N_{LR}(\omega):=n_{L}(\omega)+n_{R}(\omega)$. Then, one finds
entanglement in the stationary state whenever
$\left|N_{LR}(\omega_{1})-N_{LR}(\omega_{2})\right|$ is larger than

$2\,\sqrt{(2+N_{LR}(\omega_{1}))(2+N_{LR}(\omega_{2}))N_{LR}(\omega_{1})N_{LR}(\omega_{2})}\
.$

###  H.2 Three-spin chain

For $N=3$, the number of ones in the binary digits of length $3$ is
$p=0,1,2,3$, the integers $1\leq i_{1}<\cdots<i_{p}\leq N$ denoting the sites
at which the ones occur. If there no ones as when $p=0$ we shall set
$i_{0}=0$: the following ones are then the possible strings:

$\displaystyle(000)\,:\,i_{0}=0\
,\quad\left\\{\begin{matrix}(100)\,:\,i_{1}=1\cr(010)\,:\,i_{1}=2\cr(001)\,:\,i_{1}=3\end{matrix}\right.\
,$
$\displaystyle\left\\{\begin{matrix}(110)\,:\,i_{1}=1,i_{2}=2\cr(101)\,:\,i_{1}=1,i_{2}=3\cr(011)\,:\,i_{1}=2,i_{2}=3\end{matrix}\right.\
,\ (111)\,:\,i_{1}=1,i_{2}=2,i_{3}=3\ .$

The binary strings above are listed anti-lexicographically: it turns out that
their combinadic indices according to (95) provide the same ordering:

$\displaystyle(000)\leftrightarrow\mathcal{N}_{0}=1$
$\displaystyle(100)\leftrightarrow{\mathcal{N}}_{1}=1\ ,\
(010)\leftrightarrow\mathcal{N}_{1}=2\ ,\
(001)\leftrightarrow\mathcal{N}_{1}=3$
$\displaystyle(110)\leftrightarrow{\mathcal{N}}_{2}=1\ ,\
(101)\leftrightarrow\mathcal{N}_{2}=2\ ,\
(011)\leftrightarrow\mathcal{N}_{2}=3$
$\displaystyle(111)\leftrightarrow\mathcal{N}_{3}=1\ .$

The $p\times p$ matrices $\mathcal{D}^{(p)}$ appearing in (105) have entries
$\mathcal{D}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N}^{\prime\prime}_{p}}$
that are the determinants of the sub-matrices $U^{i_{1}<\cdots
i_{p}}_{j_{1}<\cdots j_{p}}$ of rank $p$ that are obtained from the unitary
matrix

$U=[u_{ij}]=\frac{1}{2}\begin{pmatrix}1&\sqrt{2}&1\cr\sqrt{2}&0&-\sqrt{2}\cr
1&-\sqrt{2}&1\end{pmatrix}$

by choosing the rows indexed by $i_{1}<\cdots<i_{p}$ and the columns indexed
by $j_{1}<\cdots j_{p}$. Here follows some instances of the various entries:

$\displaystyle\mathcal{D}^{(0)}_{1}=\hbox{det}\left(U^{0}_{0}\right)=1\ ,$
(229)
$\displaystyle\mathcal{D}^{(3)}_{1}=\hbox{det}\left(U^{123}_{123}\right)=-1\
,$ (230)
$\displaystyle\mathcal{D}^{(1)}_{11}=\hbox{det}\left(U^{1}_{1}\right)=u_{11}=\frac{1}{2}\
,$ (231)
$\displaystyle\mathcal{D}^{(2)}_{11}=\hbox{det}\left(U^{12}_{12}\right)=\left|\begin{matrix}u_{11}&u_{12}\cr
u_{21}&u_{22}\end{matrix}\right|=-\frac{1}{2}\ .$ (232)

It follows that $\mathcal{D}^{(1)}=U$ and $\mathcal{D}^{(2)}=-U$.

Given the diagonal symmetric matrix (96) consisting of the eigenvalues in (67)
of the stationary state $\rho_{\infty}$ , the entries of the blocks
$\mathcal{S}^{(p)}$ of the matrix $\mathcal{S}$ which represents the three-
spin stationary state with respect to the standard basis (103) are (notice
that, due to (102),
$\mathcal{S}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N}"_{p}}=\mathcal{S}^{(p)}_{\mathcal{N}^{\prime\prime}_{p}\mathcal{N}^{\prime}_{p}}$):

$\mathcal{S}^{(0)}_{11}=\mathcal{L}^{(0)}_{1}\
,\quad\mathcal{S}^{(3)}_{11}=\mathcal{L}^{(3)}_{1}\ ,$ (233)

for $p=0$ and $p=3$;

$\displaystyle\mathcal{S}^{(1)}_{11}$
$\displaystyle=\mathcal{S}^{(1)}_{33}=\frac{\mathcal{L}^{(1)}_{1}+2\mathcal{L}^{(1)}_{2}+\mathcal{L}^{(1)}_{3}}{4}\
,\hskip 28.45274pt$ (234) $\displaystyle\mathcal{S}^{(1)}_{12}$
$\displaystyle=\mathcal{S}^{(1)}_{23}=\frac{\mathcal{L}^{(1)}_{1}-\mathcal{L}^{(1)}_{3}}{2\sqrt{2}}\
,$ (235) $\displaystyle\mathcal{S}^{(1)}_{13}$
$\displaystyle=\frac{\mathcal{L}^{(1)}_{1}-2\mathcal{L}^{(1)}_{2}+\mathcal{L}^{(1)}_{3}}{4}\
,$ (236) $\displaystyle\mathcal{S}^{(1)}_{22}$
$\displaystyle=\frac{\mathcal{L}^{(1)}_{1}+\mathcal{L}^{(1)}_{3}}{2\sqrt{2}}\
,$ (237)

for $p=1$ and, for $p=2$,

$\displaystyle\mathcal{S}^{(2)}_{11}$
$\displaystyle=\mathcal{S}^{(2)}_{33}=\frac{\mathcal{L}^{(2)}_{1}+2\mathcal{L}^{(2)}_{2}+\mathcal{L}^{(2)}_{3}}{4}\
,$ (238) $\displaystyle\mathcal{S}^{(2)}_{12}$
$\displaystyle=\mathcal{S}^{(2)}_{23}=\frac{\mathcal{L}^{(2)}_{1}-\mathcal{L}^{(2)}_{3}}{2\sqrt{2}}\
,$ (239) $\displaystyle\mathcal{S}^{(2)}_{13}$
$\displaystyle=\frac{\mathcal{L}^{(2)}_{1}-2\mathcal{L}^{(2)}_{2}+\mathcal{L}^{(2)}_{3}}{4}\
,$ (240) $\displaystyle\mathcal{S}^{(2)}_{22}$
$\displaystyle=\frac{\mathcal{L}^{(2)}_{1}+\mathcal{L}^{(2)}_{3}}{2}\ .$ (241)

Using (106) and the above expressions for
$\mathcal{S}^{(p)}_{\mathcal{N}^{\prime}_{p}\mathcal{N}^{\prime\prime}_{p}}$,
one recovers the following algebraic form for the diagonal contributions to
the stationary state $\rho_{\infty}$ in the spin-operator representation:

$\displaystyle\hskip-28.45274pt\rho^{diag}_{\infty}$ $\displaystyle=$
$\displaystyle\mathcal{L}^{(0)}_{1}\,P_{---}\,+\,\mathcal{L}^{(3)}_{1}\,P_{+++}$
(242) $\displaystyle+$
$\displaystyle\frac{\mathcal{L}^{(1)}_{1}+\mathcal{L}^{(1)}_{3}}{2}\,P_{-+-}\,+\,\frac{\mathcal{L}^{(2)}_{1}+\mathcal{L}^{(2)}_{3}}{2}P_{+-+}$
$\displaystyle+$
$\displaystyle\frac{\mathcal{L}^{(1)}_{1}+2\mathcal{L}^{(1)}_{2}+\mathcal{L}^{(1)}_{3}}{4}\,\Big{(}P_{+--}\,+\,P_{--+}\Big{)}$
$\displaystyle+$
$\displaystyle\frac{\mathcal{L}^{(2)}_{1}+2\mathcal{L}^{(2)}_{2}+\mathcal{L}^{(2)}_{3}}{4}\,\Big{(}P_{++-}\,+\,P_{-++}\Big{)}\
,$

where $P_{ijk}=P_{i}\otimes P_{j}\otimes P_{k}$, $i,j,k=\pm$ and
$\displaystyle P_{\pm}=\frac{1\pm\sigma_{z}}{2}$. The off-diagonal
contributions instead read

$\displaystyle\hskip-14.22636pt\rho_{\infty}^{off}=$
$\displaystyle\hskip-14.22636pt=\frac{\mathcal{L}^{(1)}_{1}-\mathcal{L}^{(1)}_{3}}{2\sqrt{2}}\,\Big{(}\sigma_{+}\otimes\sigma_{-}\otimes
P_{-}+P_{-}\otimes\sigma_{+}\otimes\sigma_{-}\Big{)}$
$\displaystyle\hskip-14.22636pt+\frac{\mathcal{L}^{(2)}_{1}-\mathcal{L}^{(2)}_{3}}{2\sqrt{2}}\,\Big{(}P_{+}\otimes\sigma_{+}\otimes\sigma_{-}+\sigma_{+}\otimes\sigma_{-}\otimes
P_{+}\Big{)}$
$\displaystyle\hskip-14.22636pt+\frac{\mathcal{L}^{(1)}_{1}-2\mathcal{L}^{(1)}_{2}+\mathcal{L}^{(1)}_{3}}{4}\,\sigma_{+}\otimes
P_{-}\otimes\sigma_{-}$
$\displaystyle\hskip-14.22636pt+\frac{\mathcal{L}^{(2)}_{1}-2\mathcal{L}^{(2)}_{2}+\mathcal{L}^{(2)}_{3}}{4}\,\sigma_{+}\otimes
P_{+}\otimes\sigma_{-}\,+\,\hbox{h.c.}$ (243)

Due to the correspondence in (213) and (214) among the eigenvalues of
$\rho_{\infty}$ in (204)–(210) with those in BFM , one checks that the spin-
operator expressions above coincide with those obtained there.

##  Appendix I Two-spin concurrence in a three-spin chain

Aided by the fact that, in the case $N=3$, the lexicographic and combinadic
indices coincide as expressed by (201) and (202), in order to find the indices
$\mathcal{N}^{(rs)}_{p}(n^{\prime},n^{\prime\prime})$, one proceeds as
follows. Since for $p=0$ the only possible string is $(000)$,

$\displaystyle\mathcal{N}^{(12)}_{0}(0,0)=1\ ,\ \mathcal{N}^{(12)}_{0}(0,1)=0\
,$ (244) $\displaystyle\mathcal{N}^{(12)}_{0}(1,0)=0\ ,\
\mathcal{N}^{(12)}_{0}(1,1)=0\ .$ (245)

For $p=1$, the possible strings are $(100)$, $(010)$ and $(001)$, whence

$\displaystyle\mathcal{N}^{(12)}_{1}(0,0)=3\ ,\ \mathcal{N}^{(12)}_{1}(0,1)=2\
,$ (246) $\displaystyle\mathcal{N}^{(12)}_{1}(1,0)=1\ ,\
\mathcal{N}^{(12)}_{1}(1,1)=0\ .$ (247)

For $p=2$, the possible strings are $(110)$, $(101)$ and $(011)$, whence

$\displaystyle\mathcal{N}^{(12)}_{2}(0,0)=0\ ,\ \mathcal{N}^{(12)}_{2}(0,1)=3\
,$ (248) $\displaystyle\mathcal{N}^{(12)}_{2}(1,0)=2\ ,\
\mathcal{N}^{(12)}_{2}(1,1)=1\ .$ (249)

Finally, for $p=3$, the only possible string is $(111)$ so that

$\displaystyle\mathcal{N}^{(12)}_{3}(0,0)=0\ ,\ \mathcal{N}^{(12)}_{3}(0,1)=0\
,$ (250) $\displaystyle\mathcal{N}^{(12)}_{3}(1,0)=0\ ,\
\mathcal{N}^{(12)}_{3}(1,1)=1\ .$ (251)

Similarly, for $r=1$ and $s=3$, that is for computing the concurrence of the
first and third spin, one finds

$\displaystyle\mathcal{N}^{(13)}_{0}(0,0)=1\ ,\ \mathcal{N}^{(13)}_{0}(0,1)=0\
,$ (252) $\displaystyle\mathcal{N}^{(13)}_{0}(1,0)=0\ ,\
\mathcal{N}^{(13)}_{0}(1,1)=0\ ,$ (253)
$\displaystyle\mathcal{N}^{(13)}_{1}(0,0)=2\ ,\ \mathcal{N}^{(13)}_{1}(0,1)=3\
,$ (254) $\displaystyle\mathcal{N}^{(13)}_{1}(1,0)=1\ ,\
\mathcal{N}^{(13)}_{1}(1,1)=0\ ,$ (255)
$\displaystyle\mathcal{N}^{(13)}_{2}(0,0)=0\ ,\ \mathcal{N}^{(13)}_{2}(0,1)=3\
,$ (256) $\displaystyle\mathcal{N}^{(13)}_{2}(1,0)=1\ ,\
\mathcal{N}^{(13)}_{2}(1,1)=2\ ,$ (257)
$\displaystyle\mathcal{N}^{(13)}_{3}(0,0)=0\ ,\ \mathcal{N}^{(13)}_{3}(0,1)=0\
,$ (258) $\displaystyle\mathcal{N}^{(12)}_{3}(1,0)=0\ ,\
\mathcal{N}^{(12)}_{3}(1,1)=1\ .$ (259)

Finally, for the concurrence of the second and third spin, setting $r=2$ and
$s=3$ one finds

$\displaystyle\mathcal{N}^{(23)}_{0}(0,0)=1\ ,\ \mathcal{N}^{(23)}_{0}(0,1)=0\
,$ (260) $\displaystyle\mathcal{N}^{(23)}_{0}(1,0)=0\ ,\
\mathcal{N}^{(23)}_{0}(1,1)=0\ ,$ (261)
$\displaystyle\mathcal{N}^{(23)}_{1}(0,0)=1\ ,\ \mathcal{N}^{(23)}_{1}(0,1)=3\
,$ (262) $\displaystyle\mathcal{N}^{(23)}_{1}(1,0)=2\ ,\
\mathcal{N}^{(23)}_{1}(1,1)=0\ ,$ (263)
$\displaystyle\mathcal{N}^{(23)}_{2}(0,0)=0\ ,\ \mathcal{N}^{(23)}_{2}(0,1)=2\
,$ (264) $\displaystyle\mathcal{N}^{(23)}_{2}(1,0)=1\ ,\
\mathcal{N}^{(23)}_{2}(1,1)=3\ ,$ (265)
$\displaystyle\mathcal{N}^{(23)}_{3}(0,0)=0\ ,\ \mathcal{N}^{(23)}_{3}(0,1)=0\
,$ (266) $\displaystyle\mathcal{N}^{(23)}_{3}(1,0)=0\ ,\
\mathcal{N}^{(23)}_{3}(1,1)=1\ .$ (267)

Given the indices of the entries of the matrices $\mathcal{S}^{(p)}$ in (105)
which are necessary to compute the quantities $a$, $b$ and $c$ in (121), (120)
and (124), in the case of $r=1$, $s=2$ one finds

$\displaystyle a$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(2)}_{11}\,+\,\mathcal{S}^{(3)}_{11}\ ,$ (268)
$\displaystyle c$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(1)}_{12}\,+\,\mathcal{S}^{(2)}_{23}\ ,$ (269)
$\displaystyle e$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(0)}_{11}\,+\,\mathcal{S}^{(1)}_{33}\ ,$ (270)

while, in the case of $r=1$, $s=3$,

$\displaystyle a$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(2)}_{22}\,+\,\mathcal{S}^{(3)}_{11}\ ,$ (271)
$\displaystyle c$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(1)}_{13}\,+\,\mathcal{S}^{(2)}_{13}\ ,$ (272)
$\displaystyle e$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(0)}_{11}\,+\,\mathcal{S}^{(1)}_{22}\ ,$ (273)

and

$\displaystyle a$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(2)}_{33}\,+\,\mathcal{S}^{(3)}_{11}\ ,$ (274)
$\displaystyle c$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(1)}_{23}\,+\,\mathcal{S}^{(2)}_{12}\ ,$ (275)
$\displaystyle e$ $\displaystyle=$
$\displaystyle\mathcal{S}^{(0)}_{11}\,+\,\mathcal{S}^{(1)}_{11}\ ,$ (276)

in the case of $r=2$ and $s=3$. Insertion of (233)– (241) into the previous
expressions finally yields (277)– (279) for both $r=1$, $s=2$ and $r=2$,
$s=3$, while (280)– (282) result for $r=1$, $s=3$.

In order to inspect the stationary entanglement of the first two spins, we set
$r=1$ and $s=2$ and seek the combinadic indices
$\mathcal{N}^{(12)}_{p}(n^{\prime},n^{\prime\prime})$ that select the entries
of $\mathcal{S}^{(p)}$ to be used in (120), (121) and (124). As shown above,
the coefficients contributing to the concurrence $C_{1,2}$ in (125), are

$\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{\mathcal{L}^{(2)}_{1}+2\mathcal{L}^{(2)}_{2}+\mathcal{L}^{(2)}_{3}+4\mathcal{L}^{(3)}_{1}}{4}\
,$ (277) $\displaystyle c$ $\displaystyle=$
$\displaystyle\frac{\mathcal{L}^{(1)}_{1}+\mathcal{L}^{(2)}_{1}-\mathcal{L}^{(1)}_{3}-\mathcal{L}^{(2)}_{3}}{2\sqrt{2}}\
,$ (278) $\displaystyle e$ $\displaystyle=$
$\displaystyle\frac{\mathcal{L}^{(1)}_{1}+2\mathcal{L}^{(1)}_{2}+\mathcal{L}^{(1)}_{3}+4\mathcal{L}^{(0)}_{1}}{4}\
.$ (279)

In a similar fashion, again as shown in Appendix I, in the case of $C_{1,3}$,
one finds

$\displaystyle\hskip-14.22636pta$ $\displaystyle=$
$\displaystyle\frac{\mathcal{L}^{(2)}_{1}+2\mathcal{L}^{(3)}_{1}+\mathcal{L}^{(2)}_{3}}{2}\
,$ (280) $\displaystyle\hskip-14.22636ptc$ $\displaystyle=$
$\displaystyle\frac{\mathcal{L}^{(1)}_{1}-2\mathcal{L}^{(1)}_{2}+\mathcal{L}^{(1)}_{3}+\mathcal{L}^{(2)}_{1}-2\mathcal{L}^{(2)}_{2}+\mathcal{L}^{(2)}_{3}}{4}\
,$ (281) $\displaystyle\hskip-14.22636pte$ $\displaystyle=$
$\displaystyle\frac{\mathcal{L}^{(1)}_{1}+2\mathcal{L}^{(0)}_{2}+\mathcal{L}^{(1)}_{3}}{4}\
.$ (282)

while $C_{2,3}=C_{1,2}$ as the coefficients $a,c,e$ coincide with those for
$\rho_{(1,2)}$. The explicit expressions of the concurrences $C_{r,s}$ are not
particularly suggestive and their dependence on $r$, $s$ and on the bath
temperatures must be addressed numerically: this will be done in the next
section for arbitrarily large chains. Here we shall focus upon the coefficient
$c$ as its vanishing gives zero concurrence and thus excludes the existence of
entanglement: in particular, we look at it under the simplifying assumption
$h_{L,R}(\omega_{\ell})=h$ for all $\ell=1,2,3$. Then, the eigenvalues
$\mathcal{L}^{(p)}_{\mathcal{N}_{p}}$ in (203)– (210) together with (71) yield

$c=\frac{N_{LR}(\omega_{1})-N_{LR}(\omega_{3})}{4\sqrt{2}(1+N_{LR}(\omega_{1}))(1+N_{LR}(\omega_{3}))}\
,$ (283)

where $N_{LR}(\omega_{\ell}):=n_{L}(\omega_{\ell})+n_{R}(\omega_{\ell})$ for
$\rho_{(12)}$ and $\rho_{(23)}$, while

$c=\frac{N_{LR}(\omega_{1})+N_{LR}(\omega_{3})-2N_{LR}(\omega_{2})}{8(1+N_{LR}(\omega_{1}))(1+N_{LR}(\omega_{2})(1+N_{LR}(\omega_{3}))}$
(284)

for $\rho_{(13)}$. Unlike for the spin (88), (89) and heat (94) flows, that
depend on the differences $n_{L}(\omega_{\ell})-n_{R}(\omega_{\ell})$, $c$
need not vanish even for identical baths such that
$n_{L}(\omega_{\ell})=n_{R}(\omega_{\ell})=n(\omega_{\ell})$ for $\ell=1,2,3$.
Indeed, in the latter case the stationary state is the Gibbs thermal state
(75) which can carry two-spin quantum correlations because of the inter-spin
interactions.

##  Appendix J Computing the concurrence for $N$ spin chains

In order to compute the concurrence in the genral case of an $N$-spin chain,
one needs the coefficients in (121)– (123); for that purpose, one has to
select from each $p\times p$ matrix $\mathcal{S}^{(p)}$ the entries specified
by the indices $\mathcal{N}^{(rs)}_{p}(0,0)$, $\mathcal{N}^{(rs)}_{p}(0,1)$,
$\mathcal{N}^{(rs)}_{p}(1,0)$ and $\mathcal{N}^{(rs)}_{p}(1,1)$. These
combinadic indices correspond to a total number of ones equal to $p$: among
them the indices $\mathcal{N}^{(rs)}_{p}(0,0)$ number binary strings with
zeroes at sites $r$ and $s$ and $p$ ones over the remaining $N-2$ sites,
$\mathcal{N}^{(rs)}_{p}(0,1)$ and $\mathcal{N}^{(rs)}_{p}(1,0)$ those binary
strings with $1$ at site $s$, respectively $r$ and $p-1$ ones over the
remaining $N-2$ sites and, finally, the combinadic indices
$\mathcal{N}^{(rs)}_{p}(1,1)$ list the binary strings with two ones at sites
$r$ and $s$ and $p-2$ ones distributed over the remaining $N-2$ sites. We will
reconstruct such combinadic indices by means of the choices $i_{1}<\cdots
i_{p}$, $i_{1}<\cdots<i_{p-1}$ and $i_{1}<\cdots<i_{p-2}$ among the $N-2$
sites whereat the ones not already allocated at $r$ and/or $s$ can be
assigned.

In order to do this, let us first consider $\mathcal{N}^{(rs)}_{p}(0,0)$:
according to (95), $\mathcal{N}^{(rs)}_{p}(0,0)$ labels all $N$-tuple with $p$
ones distributed over all sites but the site $r$ and the site $s$. There are
at most $N-2\choose p$ of such sites if $0\leq p\leq N-2$; let
$i_{1}<\cdots<i_{p}$ be the sites with $i_{\ell}\neq r,s$ chosen among

$1,2,\ldots,r-1,r+1,\ldots,s-1,s+1,\ldots,N\ .$ (285)

Then, the required indices will be of the form

$\mathcal{N}_{p}^{(rs)}(0,0)=1\,+\,\sum_{\ell=1}^{p}{i_{\ell}-1\choose\ell}\
.$ (286)

In the case of the combinadic indices $\mathcal{N}^{(rs)}_{p}(0,1)$, there are
$p-1$ ones to be distributed over the $N-2$ sites in (285) and $N-2\choose
p-1$ such choices, the binomial vanishing if $p<N-2$.
$\mathcal{N}^{(rs)}_{p}(0,1)$ signals a $1$ at site $s$.
Let $i_{1}<\cdots<i_{p-1}$ be the positions of the other $1$’s and let
$i_{s^{*}}$ denote the largest $0\leq i_{\ell}<s$, corresponding to the
following distribution of ones:

$0<i_{1}<\cdots<i_{s^{*}}<s<i_{s^{*}+1}<\cdots<i_{p-1}\ .$

Then, with the proviso that $s^{*}=0$, if $i_{s*}=0$, that is when all other
ones occurs at sites beyond $s$ and that the sums are set to zero if the first
summation index is smaller than the last one, the required combinadic indices
are retrieved as

$\displaystyle\mathcal{N}_{p}^{(rs)}(0,1)$ $\displaystyle=$ $\displaystyle
1\,+\,\sum_{\ell=1}^{s^{*}}{i_{\ell}-1\choose\ell}\,+\,{s-1\choose s^{*}+1}$
(287) $\displaystyle+$
$\displaystyle\,\sum_{\ell=s^{*}+1}^{p-1}{i_{\ell}-1\choose\ell+1}\ .$

Analogously, when a $1$ occurs at site $r$, then

$\displaystyle\mathcal{N}_{p}^{(rs)}(1,0)$ $\displaystyle=$ $\displaystyle
1\,+\,\sum_{\ell=1}^{r^{*}}{i_{\ell}-1\choose\ell}\,+\,{r-1\choose r^{*}+1}$
(288) $\displaystyle+$
$\displaystyle\sum_{\ell=r^{*}+1}^{p-1}{i_{\ell}-1\choose\ell+1}\ .$

Finally, in the case of two $1$’s at $r$ and $s$ there remain other $p-2$ ones
to be distributed among the $N-2$ sites in (285), thus making for $N-2\choose
p-2$ indices $\mathcal{N}_{p}^{(rs)}(1,1)$. Then, the combinadic index

$\displaystyle\mathcal{N}_{p}^{(rs)}(1,1)$ $\displaystyle=$ $\displaystyle
1\,+\,\sum_{\ell=1}^{r^{*}}{i_{\ell}-1\choose\ell}\,+\,{r-1\choose r^{*}+1}$
(289) $\displaystyle+$
$\displaystyle\,\sum_{\ell=r^{*}+1}^{s^{*}}{i_{\ell}-1\choose\ell+1}+{s-1\choose
s^{*}+2}$ $\displaystyle+$
$\displaystyle\,\sum_{\ell=s^{*}+1}^{p-2}{i_{\ell}-1\choose\ell+2}\ .$

corresponds to a choice of $p-2$ ones of the form

$\displaystyle i_{1}<\cdots<i_{r^{*}}<r<i_{r^{*}+1}<\cdots$
$\displaystyle\hskip 28.45274pt<\cdots<i_{s^{*}}<s<i_{s^{*}+1}<\cdots<i_{p-2}\
,$

with the positions $i_{1}<\cdots<i_{p-2}$ of the $1$’s in an $N$-tuple where
two $1$’s are already present at sites $r<s$ and $r^{*}$ and $s^{*}$ denoting
the largest integers such that $i_{r^{*}}<r$ and $i_{s^{*}}<s$.

Notice that by setting $r^{*}=0$ and $s^{*}=0$ if $i_{r*}=0$ and $i_{s*}=0$
and with the convention about the sums introduced before (287), the above
expression accounts also for the cases when $i_{r^{*}}=i_{s^{*}}=0$ which
corresponds to having all ones at sites beyond $s$,

$0<r<s<i_{1}<\cdots<i_{p-2}\ ,$

the cases when there are ones before $r$, but no ones in between $r$ and $s$,
that is when $0<i_{r^{*}}$, $i_{s*}=0$,

$0<i_{1}<\cdots<i_{r^{*}}<r<s<i_{r^{*}+1}<\cdots<i_{p-2}\ ,$

and the cases when there are no ones before $r$, but $s^{*}$ ones before $s$,
that is when $i_{r^{*}}=0$ and $i_{s^{*}}>0$,

$0<r<i_{1}<\cdots<i_{s^{*}}<s<i_{s^{*}+1}<\cdots<i_{p-2}\ .$

## References

  * (1) S. Datta, Quantum Transport: Atom to Transistor, (Cambridge University Press, Cambridge, 2005)
  * (2) J. Gemmer, M. Michel and G. Mahler, Quantum Thermodynamics Lect. Notes Phys. 784, (Springer, Berlin, 2009)
  * (3) V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 3d Ed., (Wiley, Weinheim, 2011)
  * (4) G. Benenti et al., Phys. Rep. 694 (2017) 1
  * (5) S. Lepri, R. Livi and A. Politi, Phys. Rep. 377 (2003) 1
  * (6) L.-A. Wu and D. Segal, J. Phys. A 42 (2009) 025302
  * (7) F.Caruso et al., J. Chem. Phys. 131 (2009) 105106
  * (8) J.T. Barreiro et al., Nature 470 (2011) 486
  * (9) J. Wu and M. Berciu, Phys. Rev. B 83 (2011) 214416
  * (10) F. Giazotto and M.J. Martinez-Perez, Nature 492 (2012) 401
  * (11) J.-P. Brantut et al., Science 342 (2013) 2013
  * (12) R. Labouvie et al., Phys. Rev. Lett. 115 (2015) 050601
  * (13) R. Labouvie et al., Phys. Rev. Lett. 116 (2016) 235302
  * (14) F. Schlawin et al., Nat.Commun. 4 (2013) 1782
  * (15) A. Bermudez, M. Bruderer and M.B. Plenio, Phys. Rev. Lett. 111 (2013) 040601
  * (16) B. Leggio, R. Messina and M. Antezza, Europhys. Lett. 110 (2015) 40002
  * (17) N. Freitas, E.A. Martinez and J.P. Paz, Phys. Scr. 91 (2016) 013007
  * (18) B. Dutta et al., Phys. Rev. Lett. 119 (2017) 077701
  * (19) P. Doyeux, R. Messina, B. Leggio and M. Antezza, Phys.Rev.A 95 (2017) 012138
  * (20) R. Biele et al., npj Quantum Materials 2 (2017) 38
  * (21) B. Bertini et al., Finite-temperature transport in one-dimensional quantum lattice models, arXiv:2003.0334
  * (22) R. Alicki and K. Lendi, Quantum dynamical semigroups and applications, Lect. Notes Phys. 717, (Springer Verlag, Berlin, 2007)
  * (23) A. Rivas and S.F. Huelga, Open Quantum Systems (Springer Verlag, Berlin, 2012)
  * (24) F. Benatti and R. Floreanini, Int. J. Mod. Phys. B 19 (2005) 3063
  * (25) F. Benatti, Dynamics, Information and Complexity in Quantum Systems, (Springer, Berlin, 2009)
  * (26) R. Alicki, Invitation to quantum dynamical semigroups, in: Lect. Notes Phys. 597, P. Garbaczewski and R. Olkiewicz, Eds., (Springer-Verlag, Berlin, 2002), p.239
  * (27) Dissipative Quantum Dynamics, F. Benatti and R. Floreanini, Eds., Lect. Notes Phys. 622, (Springer-Verlag, Berlin, 2003)
  * (28) A. Kossakowski, Bull. Acad. Pol. Sc. 12 (1972) 1021
  * (29) E.B. Davies, Comm. Math. Phys. 39 (1974) 91
  * (30) E.B. Davies, Math. Ann. 219 (1976) 147
  * (31) E.B. Davies, Quantum Theory of Open Systems, (Academic Press, New York, 1976)
  * (32) V. Gorini, A. Kossakowski and E.C.G. Sudarshan, J. Math. Phys. 17 (1976) 821
  * (33) G. Lindblad, Comm. Math. Phys. 48 (1976) 119
  * (34) V. Gorini, A. Frigerio, M. Verri, A. Kossakowski and E.G.C. Sudarshan, Rep. Math. Phys. 13 (1976) 149
  * (35) R. Dümcke, H. Spohn, Z. Phys. B34 (1979) 419
  * (36) H. Spohn, Rev. Mod. Phys. 52 (1980) 569
  * (37) M. Merkli, Quantum markovian master equations: Resonance theory overcomes the weak coupling regime, arXiv:1908.01984
  * (38) E.B. Davies, J. Stat. Phys. 18 (1978) 161
  * (39) R. Alicki, J. Phys. A: Math. Gen. 12 (1979) L103
  * (40) H. Spohn and J.L. Lebowitz, Adv. Chem. Phys. 38 (1978) 109
  * (41) H. Zoubi, M. Orenstien and A. Ron, Phys. Rev. A 67 (2003) 063813
  * (42) H. Wichterich, M.J. Henrich, H.-P. Breuer, J. Gemmer and M. Michel, Phys. Rev. E 76 (2007) 031115
  * (43) A. Rivas, A.Plato, S.F. Huelga and M. Plenio, New J. Phys. 12 (2010) 113032
  * (44) R. Migliore et al., J. Phys. B 44 (2011) 075503
  * (45) J.P. Santos and F.L. Semiao, Phys. Rev. A 89 (2014) 022128
  * (46) T. Werlang and D. Valente, Phys. Rev. E 91 (2015) 012143
  * (47) J.P. Santos and G.T. Landi, Phys. Rev. E 94 (2016) 062143
  * (48) A. Rivas and M.A. Martin-Delgado, Scient. Rep. 7 (2017) 6350
  * (49) M. Michel and O. Hess, Phys. Rev. B 77 (2008) 104303
  * (50) D. Karevski and T. Platini, Phys. Rev. Lett. 102 (2009) 207207
  * (51) M. Znidaric, Phys. Rev. Lett. 106 (2011) 220601
  * (52) T. Prosen Phys. Rev. Lett. 107 (2011) 137201
  * (53) T. Prosen and M. Znidaric, Phys. Rev. B 86 (2012) 125118
  * (54) T. Prosen, Phys. Scr. 86 (2012) 058511
  * (55) T. Prosen, J. Phys. A 48 (2015) 373001
  * (56) V. Popkov, J. Stat. Phys. 2012 (2012) P12015
  * (57) J.J. Mendoza-Arenas, S. Al-Assam, S.R. Clark and D. Jaksch, J. Stat. Phys. 2013 (2013) P07007
  * (58) D. Karevski, V. Popkov and G.M. Schütz, Phys. Rev. Lett. 110 (2013) 047201
  * (59) L.A. Correa, J.P. Palao, G. Adesso and D. Alonso, Phys. Rev. E 87 (2013) 042131
  * (60) V. Popkov and M. Salerno J. Stat. Phys. 2013 (2013) P02040
  * (61) A. Asadian, D. Manzano, M. Tiersch and H. J. Briegel, Phys. Rev. E 87 (2013) 012109
  * (62) V. Popkov, M. Salerno and R. Livi, New J. Phys. 15 (2013) 023030
  * (63) G.T. Landi, E. Novais, M.J. de Oliveira and D. Karevski, Phys. Rev. E 90 (2014) 042142
  * (64) D. Manzano and P.I. Hurtado, Phys. Rev. B 90 (2014) 125138
  * (65) J. Cui, J.I Cirac and M.C. Banuls, Phys. Rev. Lett. 114 (2015) 220601
  * (66) V. Popkov, M. Salerno and R. Livi, New J. Phys. 17 (2015) 023066
  * (67) F. Nicacio, A. Ferraro, A. Imparato, M. Paternostro and F.L. Semiao, Phys. Rev. E 91 (2015) 042116
  * (68) G.T. Landi and D. Karevski, Phys. Rev. B 91 (2015) 174422
  * (69) L. Schuab, E. Pereira and G.T. Landi, Phys. Rev. E 94 (2016) 042122
  * (70) P.H. Guimaraes, G.T. Landi and M.J. de Oliveira, Phys. Rev. E 94 (2016) 03213
  * (71) S. Campbell, G. De Chiara M. Paternostro, Scient. Rep. 6 (2016) 19730
  * (72) D. Manzano, C. Chuang and J. Cao, New J. Phys. 18 (2016) 043044
  * (73) C. Monthus, J. Stat. Phys. 2017 (2017) 043303
  * (74) G. De Chiara et al., New J. Phys. 20 (2018) 113024
  * (75) F. Carollo, J.P. Garrahan, I. Lesanovsky and C. Perez-Espigares, Phys. Rev. E 96 (2017) 052118
  * (76) F. Carollo, J.P. Garrahan and I. Lesanovsky, Phys. Rev. B 98 (2018) 094301
  * (77) T. Chanda et al., Phys. Rev. A 97 (2018) 062324
  * (78) E. Pereira, Phys. Rev. E 97 (2018) 022115
  * (79) M. Brenes et al. Phys. Rev. B 98 (2018) 235128
  * (80) K.V. Hovhannisyan and A. Imparato, New J. Phys. 21 (2019) 052001
  * (81) A. Levy and R. Kosloff, Europh. Lett. 107 (2014) 20004
  * (82) A.S. Trushechkin and I.V. Volovich, Europh. Lett. 113 (2016) 30005
  * (83) J. B. Brask et al., New J. Phys. 17 (2015) 113029
  * (84) S. Khandelwal et al., New J. Phys. 22 (2020) 073039
  * (85) G.L. Decordi and A. Vidiella-Barranco, Opt. Commun. 387 (2017) 366
  * (86) J.T. Stockburger and T. Motz, Fortschr. Phys. 65 (2017) 1600067
  * (87) J.O. Gonzalez et al., Open Syst. Inf. Dyn. 24 (2017) 1740010
  * (88) G.G. Giusteri et al., Phys.Rev. E 96 (2017) 012113
  * (89) P.P. Hofer et al., New J. Phys. 19 (2017) 123037
  * (90) M. Tahir Naseem, A. Xuereb and O.E. Mustecaplioglu, Phys. Rev. A 98 (2018) 052123
  * (91) N. Shammah et al., Phys. Rev. A 98 (2018) 063815
  * (92) J. Kolodynski et al., Phys.Rev. A 97 (2018) 062124
  * (93) M.T. Mitchison and M. Plenio,New J. Phys. 20 (2018) 033005
  * (94) E. Mascarenhas et al., Phys.Rev. B 99 (2019) 245134
  * (95) M. Cattaneo et al., New J. Phys. 21 (2019) 113045
  * (96) D. Farina et al., Phys. Rev. A 102 (2020) 052208
  * (97) F. Benatti, R. Floreanini, L. Memarzadeh, Phys. Rev. A 102 (2020) 042219
  * (98) P. Coleman, Introduction to Many-Body Physics, (Cambridge University Press, Cambridge, 2015)
  * (99) A. B. Siddique, S. Farid, and M.Tahir, Proof of bijection for combinatorial number system, arXiv:1601.05794 (2016)
  * (100) S.Hill and W.K. Wootters, Phys. Rev. Lett. 78 (1997) 5022
  * (101) H. Spohn, Rep. Math. Phys. 10 (1976) 189
  * (102) A. Frigerio, Lett. Math. Phys. 2 (1977) 79
  * (103) A. Frigerio, Comm. Math. Phys. 63, 269 (1977)
  * (104) D.E. Evans, Commun. Math. Phys. 54 (1977) 293
  * (105) F. Fagnola and R. Rebolledo, J. Math. Phys. 42 (2001) 1296
  * (106) F. Fagnola and R. Rebolledo, Infin. Dimens. Anal. Qu. 11 (2008) 467
  * (107) B. Baumgartner, H. Narnhofer and W. Thirring, J. Phys. A 41 (2008) 065201
  * (108) B. Baumgartner and H. Narnhofer, J. Phys. A 41 (2008) 395303
  * (109) F. Benatti, A. Nagy and H. Narnhofer, J. Phys. A 44 (2011) 155303
  * (110) B. Baumgartner and H. Narnhofer, Rev. Math. Phys. 24 (2012) 1250001
  * (111) S. Khandelwal et al., New J. Phys. 22 (2020) 073039