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

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 [1]- [21].
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 [22]- [37].
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 [38]- [48]). These gives rise to new, global effects in transport phenomena that can not be captured using other, simplified approaches to the chain open dynamics. 1 * memarzadeh@sharif.edu The authors contributed equally to this study and are listed in alphabetical order. 1 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 equation (e.g. see [49]- [80]). Although the two approaches, named global and local, are regularly adopted in applications and com-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.
pared [41,[43][44][45][46]70], [81]- [96], it turns out that the local approach might not be able to capture all the correct system transport properties [97]. arXiv:2102.10036v1 [quant-ph] 19 Feb 2021 2 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, . . . , 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: with free boundary ocnditions, where ∆ > 0 is the intensity of the constant transverse magnetic field, σ ( ) x,y,z are the Pauli matrices at site , 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, = κ 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 where b α (ν), b † α (ν) are Bosonic operators satisfying the canonical commutation relations The coupling of the baths to the left and right spins are described by the interaction Hamiltonian where λ << 1 is a dimensionless coupling constant, are spin ladder operators at site , whence σ (L) are bath operators, with * meaning complex conjugation and h L,R (ν) are suitable smearing functions. Referring to [97] 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 ≡ 1/β L and T R ≡ 1/β R , the state of the environment is then given by It is invariant under the bath dynamics generated by H env = α=L,R H α and exhibits thermal expectations of the form Tr B ρ env b α (ν)b † α (ν ) = δ αα δ(ν − ν ) (1 + n α (ν)), (8) with thermal mean occupation numbers Finally, choosing the initial state of the compound system chain plus baths of the form ρ tot (0) = ρ(0) ⊗ ρ env , with ρ(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 [22]- [37]. In practice, the initial state of the compound system spin-chain plus baths evolves into ρ tot (t) = e −itHtot ρ tot (0) e itHtot , where H tot = H +H env +H is the total system Hamiltonian. The state of the open chain at time t, ρ(t), is then retrieved by tracing over the baths degrees of freedom, ρ(t) = Tr env ρ tot (t) . Then, one rescales the physical time variable t to τ = t λ 2 and takes the limit λ → 0 in ∂ t Tr env ρ tot (t) . In doing so, too fast oscillations with respect to the chain transition frequencies ω = E i − E j are suppressed, where E i and |E i solve the spin Hamiltonian eigenvalue equation H|E i = E i |E i . This procedure corresponds to the so-called rotating wave approximation leading to a master equation of the GKSL form On the right hand side of the above time-evolution equation, one distinguishes a Hamiltonian term λ 2 H LS which provides a Lamb-shift correction to the spin-chain Hamiltonian H and a purely dissipative term D[ρ(t)]. As we shall see, in the specific physical context here considered, ω [ρ(t)] resulting from positive transition frequencies ω ≥ 0, only: Their explicit form reads whose coefficients with n α (ω) 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 where, unlike the dissipative term, the sum runs now over all positive and negative transition frequencies and whose coefficients read with P denoting the principal value. In all the previous expressions there appear Lindblad operators of the form together with their Hermitean conjugates 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 weakcoupling 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 weakcoupling procedure.
Remark 1. The fact that the dissipative contribution to the generator, D[ρ(t)], involves only transition frequencies ω ≥ 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 † α (ω)b α (ν) contribute with time oscillations exp(±i(ν − ω)t). On the time scale τ = t/λ 2 and in the weak-coupling limit when λ → 0, fast oscillations select contributions with ω = ν ≥ 0. Negative transitions frequencies, ω ≤ 0, would also be selected if in (3) there were interaction terms of the form σ (α) which, together with their Hermitian conjugates, would correspond to the presence of terms of the form A † α (ω)b † α (ν), and Hermitian conjugates, contributing with time oscillations exp(±i(ν + ω)t) .

A. 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 -th spin ladder operators in (4) one rewrites By means of the Jordan-Wigner transformation [98], one introduces Fermionic annihilation and creation operators with the convention that z ) = 1 for j = 1, satisfying the anti-commutation relations Let | ↑ and | ↓ be the eigenvectors of σ z , σ z | ↑ = | ↑ , σ z | ↓ = −| ↓ . Since σ − | ↓ = 0, the vacuum vector such that a j |vac = 0, for all j = 1, 2, . . . , N , amounts to Using that a † j a j = σ , one inverts the transformation (23): finally turning the spin Hamiltonian into a Fermionic one, As shown in Appendix A, H can then be diagonalized, where the operators are also Fermionic, b j , b † k = δ jk with the same vacuum as the operators a j : b |vac = 0 for all = 1, 2, . . . , N , while the coefficients form an orthogonal and symmetric matrix U = [u k ].
In the following, we shall denote by n the N -tuple n 1 , n 2 , . . . , 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 . The eigenvectors of the Hamiltonian (1) have thus the form Indeed, according to (31), where, n ± denote the N -tuples n 1 , . . . , n ± 1, . . . , n N . Then, one verifies that H |n = E n |n , where Remark 2. Since the matrix U = [u k ] with entries u k as in (30) is real and symmetric, from (29) one obtains Then, by using (26), one expresses the Fermionic operators of type b in terms of spin operators: 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 Their explicit form can be derived by expressing the spin operators σ (α) + first in terms of the Fermionic operators a j , a † j , σ (L) and then in terms of the operators b , b † . Using (36), one immediately derives while the presence of + requires some preliminary manipulation. Firstly, using that 1 − 2 a † j a j = exp iπ a † j a j and that the relations (36) one gets whence, finally, By means of (41) and (33), one then computes the transition amplitudes m|σ (L) Let n 0 respectively n 1 denote the N -tuples with fixed digits n = 0, respectively n = 1, at site . Then, the only contributing transition amplitudes are n 1 |σ (L) with = 1, 2, . . . , N . Also, from (35) and (21), the transition frequencies associated with such amplitudes are while the corresponding Lindblad operators in (39) read where the symbol n means that the summation is performed over all binary 2 N −1 -tuples of indices n j = 0, 1 with j = . It thus follows thats In a similar way, from (44), one obtains that the only contributing transition amplitudes associated to the right bath are with Lindblad operators whence Notice that, in the spin representation, all Lindblad operators A † α (ω ) 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 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 † α (ω ) in (48) and (52) need be such. This means that, the Lindblad operators A † α (ω ) with ω < 0 can only contribute to the Lamb-shift Hamiltonian and not to the dissipative part of the generator. The sign of ω depends on the strength of the inter-spin coupling constant g; indeed, Correspondingly ω < 0 for g > ∆ 2 cos The second observation is that, should any of the transition frequencies ω in the list (47) be negative, the opposite one, −ω = E n0 − E n1 , not being in the list, would not give rise to a dissipative contribution of the form In the following we shall assume so that ω ≥ 0, = 1, 2, . . . , 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 L[ρ ∞ ] = 0. This follows from the fact that, as shown in Appendix D, the only operator commuting with all Lindblad operators A L,R (ω ) and A † L,R (ω ) and with the Hamiltonian must be multiples of the identity [101]- [106].
Because of the diagonal form of H + λ 2 H LS in the energy eigenbasis, H + λ 2 H LS , P k = 0, for all energy eigenprojections P k := |k k|. On the other hand, by inserting ρ(t) = P k into (13) and (14), one obtains From (12), using the two previous expressions one finds where, using (15) and (16), Consider the diagonal expression X diag = n x n P n ; then, the dissipator maps it into with whence Therefore, D[X diag ] = 0 is obtained by the factorized expressions All x ( ) ≥ 0 and, after normalization, the uniqueness of the stationary state together with the expressions (59) and (60) yield With the simplifying assumption h L,R (ω ) = h, for each = 1, 2, . . . , N , one retrieves If we further restrict to identical baths, by imposing equal temperatures and thus β L = β R = β, one computes so that On the other hand, using (35), Then, (47) implies that the open chain stationary state ρ ∞ is the Gibbs state at inverse temperature β with Hamiltonian H as given in (1):

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.

A. Stationary spin flow: sinks and sources
The spin flow at site k = 1, 2, . . . , N along the spin chain corresponds to the rate of change in time of the average of σ (k) z given by the quantity In the first equality L is the generator at the right hand side of (11), while the second equality follows from the cyclicity of the trace, Tr(XY ) = Tr(Y X) and defines the so-called dual generator L The Hamiltonian contribution to the rate of change in time of the average of σ z can be expressed in terms of the dimensionless spin currents: as where the Lamb-shift contribution is characterized by a constant 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 ρ(t) = ρ ∞ in the right hand side of (76) and find J (k,k+1) ∞ := Tr ρ ∞ J (k,k+1) = 0. Indeed, passing from spin to Fermionic operators, by (26) and (36), one finds Hence, all their averages with respect to the energy eigenstates vanish, while the columns |u k of the orthogonal and symmetric matrix U (see Remark 2) are orthogonal. Thus the stationary left, J (k−1,k) ∞ , and right, J (k,k+1) ∞ , spin currents through site k both vanish.
Clearly, being ρ ∞ time-independent, the right hand side of (76) then yields Tr ρ ∞ D[σ (k) ] = 0. However, the left and right purely dissipative contributions, do not separately vanish; indeed, as shown in Appendix E, with R as in (70). Furthermore, since from (30) one finds that . (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 (ω ) = n R (ω ) and are thus interpretable as spin flow source or sink contributions, depending on whether they are positive or negative.
contribution to the rate of change in of z can be expressed in terms of pin currents: ⇣ ift contribution is characterized by a ences in (81) thus contribute to the (76) as current divergence terms ension of energy. Since we are intionary transport properties, we set right hand side of (76) and find ⇢ 1 J (k,k+1) = 0. Indeed, passing onic operators, by (26) and (36), one ages with respect to the energy eigen- n j ) and , while the columns |u k i of the oretric matrix U (see Remark 2) are the stationary left, hJ (k 1,k) i 1 , and spin currents through site k both vantime-independent, the right hand ields Tr ⇢ 1 e D[ (k) ] = 0. However, urely dissipative contributions, with R`as in (70). Furthermore, since from (30) one finds that . (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 di↵erent 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 (!`) = n R (!`) and are thus interpretable as spin flow source or sink contributions, depending on whether they are positive or negative. Note that, due to the scaling as 1/N 2 of the products u 2 k`u 2 N`( see (30)) and the presence of N of them in (88) and (89), the source and sink terms scale as 1/N with ribution to the rate of change in z can be expressed in terms of urrents: g+) s in (81) thus contribute to the ) as current divergence terms on of energy. Since we are inry transport properties, we set ht hand side of (76) and find (k,k+1) = 0. Indeed, passing operators, by (26) and (36), one with respect to the energy eigen- n j ) and ile the columns |u k i of the orc matrix U (see Remark 2) are tationary left, hJ (k 1,k) i 1 , and currents through site k both vane-independent, the right hand with R`as in (70). Furthermore, since from (30) one finds that u 1`= ( )`u N`f or all`= 1, 2, . . . , N, it follows that Q (k) .
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 di↵erent 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 (!`) = n R (!`) and are thus interpretable as spin flow source or sink contributions, depending on whether they are positive or negative. Note that, due to the scaling as 1/N 2 of the products u 2 k`u 2 N`( see (30)) and the presence of N of them in (88) and (89), the source and sink terms scale as 1/N with  Note that, due to the scaling as 1/N 2 of the products u 2 k u 2 N (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 8 chain with N = 8 spins, set T L = 0 so that n L (ω ) = 0 and show the dependence of the source term Q (4) in the middle of the chain as a function of the right temperature T R and various values of the transverse magnetic field ∆ and the inter-spin coupling strength g. The values of g associated with ∆ 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 = 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 [81]. Notice that in the global approach developed before sinks and sources are present even in the limit where the inter-spin coupling g → 0; indeed, g appears in the thermal factors n L,R (ω ) through the transition frequencies ω (see (9)). These terms remain different and non zero whenever β L = β R , even for g = 0.

B. 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 [39,40], the heat flow through an open quantum system due to its weak coupling to a thermal bath, is measured by where ρ → ρ(t) is the dissipative evolution due to the bath, generated by L, while H is the open system timeindependent 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 Certainly, D[ρ ∞ ] = 0 implies H st L + 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 Notice that the heat flow is positive, namely it flows from the left bath into the chain if n L (ω) > n R (ω), that is (see (9)) if the left bath is at higher temperature than the right one. Furthermore, the simplifying assumption Furthermore, the transition frequencies ω in (47) are of order 1 with respect to increasing N , whence each of the N contributions ω u 2 1 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 ∆ and g as in Figure 1, the behaviour of the heat flow H st R as a function of T R is reported in Figure 2.
contribution to the rate of change in e of z can be expressed in terms of pin currents: ⇣ ift contribution is characterized by a rences in (81) thus contribute to the n (76) as current divergence terms ension of energy. Since we are intionary transport properties, we set e right hand side of (76) and find ⇢ 1 J (k,k+1) = 0. Indeed, passing onic operators, by (26) and (36), one rages with respect to the energy eigen- j , while the columns |u k i of the ormetric matrix U (see Remark 2) are the stationary left, hJ (k 1,k) i 1 , and spin currents through site k both van-1 time-independent, the right hand ields Tr ⇢ 1 e D[ (k) ] = 0. However, urely dissipative contributions, do not separately vanish; indeed, as shown in Appendix E, with R`as in (70). Furthermore, since from (30) one finds that . (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 di↵erent 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 (!`) = n R (!`) and are thus interpretable as spin flow source or sink contributions, depending on whether they are positive or negative. Note that, due to the scaling as 1/N 2 of the products u 2 k`u 2 N`( see (30)) and the presence of N of them in (88) and (89), the source and sink terms scale as 1/N with ry entanglement.
nk and source contributly related to the global rs in (48) and (52) Notice that the heat flow is positive, namely it flows from the left bath into the chain if n L (!) > n R (!), that is (see (9)) if the left bath is at higher temperature than the right one. Furthermore, the simplifying assumption .
Notice that the transition frequencies !`in (47) are of order 1 with respect to increasing N , whence each of the N contributions !`u 2 1`t o 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 and g as in Figure 1, the behaviour of the heat flow H st R as a function of T R is reported in Figure 2.  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 ρ 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.

A. Stationary state: spin representation
In this respect, instead of the standard lexicograhic ordering, it proves convenient to regroup the binary strings n in terms of the number of ones they contain. We then introduce the enumeration of the 2 N binary N -tuples n known as combinatorial numbering of degree p [99], that we shall refer to as combinadic ordering for sake of shortness. For fixed p = 0, 1, 2, . . . , N , one bijiectively associates to each of the N p N -tuples with n i = 1 at sites i 1 < · · · < i p the integers where the binomial coefficients i −1 are set to vanish if i − 1 < . According to such a numbering, we identify n with a unique N p for some p = 0, 1, . . . , N ; then, the stationary state ρ ∞ may be written as where L

(p)
Np denotes the eigenvalue Λ n in (67) corresponding to the binary N -tuple n with p 1's, indexed by the combinadic integer 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, · · · , N , the integers N p in (95) correspond to the Fermionic excitations of the modes i 1 < · · · < i p of type b. Indeed, N p identifies the binary N -tuple n, where n i1 = · · · = n ip = 1 while the remaining n j vanish. We can thus consistently label: Then, using (29) one writes (98) Notice that, unlike the indices i 1 , . . . , i p , the indices j 1 , j 2 , . . . , j p are in general not ordered; their reordering such that j 1 < j 2 < · · · < j p yields where is the determinant of the p × p sub-matrix U i1<···<ip j1<···<jp of the orthogonal and symmetric matrix U in Remark 2 with p rows indexed by i 1 < · · · < i p and p columns by j 1 < · · · < j p . Its entries are thus given by It is convenient to introduce the N p × N p matrices D (p) , where D (0) = 1 and D (N ) = −1 are scalars, otherwise D (p) has entries D (p) N p N p corresponding to the determinants D i 1 <···<i p i 1 <···<i p , where N p identifies an N -tuple with 1's at sites i 1 < · · · < i p , N p identifies an N -tuple with 1's at sites i 1 < · · · < i p . Because of (102), the matrices D (p) are symmetric for all p = 0, 1, . . . , N .
In the spin representation |vac = | ↓ ⊗N ; therefore, as shown in Appendix G, setting | ↑ = |1 S and | ↓ = |0 S so that σ + |0 S = |1 S , one can express the Fermionic states |N p with p excitations at sites i 1 < · · · < i p as linear combinations of the spin states |N p S with p spins flipped up at the sites i 1 < · · · < i p identified by the combinadic index N p . It follows that, with respect to the standard spin basis, the stationary state ρ ∞ in (67) can be recast as where L Np are the eigenvalues of ρ ∞ 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 where L (p) is the diagonal matrix whose entries are the eigenvalues in (67) labelled by the combinadic integers N p while the matrices D (p) are as in the previous remark.
Finally, again in Appendix G it is shown that the stationary state ρ ∞ has the following structure in terms of spin operators where X The above expression of the stationary state ρ ∞ 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.

B. Two-spin entanglement
The spin-operator expression of the stationary state ρ ∞ is useful to investigate the entanglement content of any pair of spins along the chain and its dependence on their positions N ≥ s > r ≥ 1. We shall quantify the two-spin entanglement by means of the concurrence [100] 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 Tr (r,s) . Considering the expression (106) one readily computes 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 n and n 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 N (rs) p (n r , n s ) and N (rs) p (n r , n 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 ρ (r,s) := Tr (r,s) ρ ∞ = n r ,n s n r ,n s N p=0 N (rs) p (n r ,n s ) Since the N -tuples indexed by N (rs) p (n r , n s ) have the same entries but, possibly, for the sites r and s, it follows that the allowed values for n r , n s and n r , n s must satisfy n r + n s = n r + n s . These latter ones and the corresponding two-spin operators are as follows: n r = 0 , n r = 0 n s = 1 , n s = 1 : n r = 0 , n r = 1 n s = 1 , n s = 0 : n r = 1 , n r = 1 n s = 0 , n s = 0 : and n r = 1 , n r = 0 n s = 0 , n s = 1 : n r = 1 , n r = 1 n s = 1 , n s = 1 : Therefore, for all sites 1 ≤ r < s ≤ N , the reduced twospin density matrix is a X-state for the case of a two-spin chain (see (H3) in Appendix H): .
Here the combinadic indices of the entries of 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 S (p) contributing to the diagonal entries: For such states the concurrence takes the following analytic expression 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.

C. 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 ∆ and g of the chain Hamiltonian, on the temperatures T L,R , on the number of spins, N , and on the spin sites 0 ≤ r ≤ s ≤ 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 S (p) N p N p in (104). In turn, such a behaviour is due to the fact that the transition frequencies ω in (47), and thus the eigenvalues (67), are of order 1 with respect to N , while the quantities D 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 (ω ) = n R (ω ), 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 in (25) which is clearly separable. close to the maximum value (55) that ensure the positivity of all transition frequencies ω in (47))) and plot various concurrences versus T R . .
For such states the concurrence takes the following analytic expression C(r, s) = 2 max 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.

C. 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 coe cients 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 and g of the chain Hamiltonian, on the temperatures T L,R , on the number of spins, N , and on the spin sites 0  r  s  N .
Firstly, although the algorithm developed in Appendix J works for all N , its symbolic numerical 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 S (104). In turn, such a behaviour is due to the fact that the transition frequencies !`in (47), and thus the eigenvalues (67), are of order 1 with respect to N , while the quantities D  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, = 1, = 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.
As regards the dependence of the concurrence on the parameters and g, Figure 4 first shows that, with temperature T L = 0, and g fixed, close to the saturation the bound (55) for = 15, the entanglement as a function of T R diminishes while increasing . This behaviour agrees with the fact that augmenting the transverse ex- .

(124)
For such states the concurrence takes the following analytic expression 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.

C. 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 coe cients 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 and g of the chain Hamiltonian, on the temperatures T L,R , on the number of spins, N , and on the spin sites 0  r  s  N .
Firstly, although the algorithm developed in Appendix J works for all N , its symbolic numerical 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 S (p) N 0 p N 00 p in (104). In turn, such a behaviour is due to the fact that the transition frequencies !`in (47), and thus the eigenvalues (67), are of order 1 with respect to N , while the quantities D separable. close to the maximum value (55) that ensure the positivity of all transition frequencies !`in (47))) and plot various concurrences versus T R . 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, = 1, = 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.
As regards the dependence of the concurrence on the parameters and g, Figure 4 first shows that, with temperature T L = 0, and g fixed, close to the saturation the bound (55) for = 15, the entanglement as a function of T R diminishes while increasing . This behaviour agrees with the fact that augmenting the transverse ex- 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, λ = 1, ∆ = 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.
As regards the dependence of the concurrence on the parameters ∆ and g, Figure 4 first shows that, with temperature T L = 0, and g fixed, close to the saturation the bound (55) for ∆ = 15, the entanglement as a function of T R diminishes while increasing ∆. 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 ∆, 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 cho- hat, due to the scaling as 1/N 2 of the products see (30)) and the presence of N of them in (88) , the source and sink terms scale as 1/N with R as it should physically be. suming h L (!) = h R (!) = h, one gets . (89) sees that, while the continuity equation (76) tionary case does not contain any current dit site k, it does however contain terms of a rigin that are due to the presence of the two ese terms vanish only if the temperatures are so that n L (!`) = n R (!`) and are thus inters spin flow source or sink contributions, den whether they are positive or negative. at, due to the scaling as 1/N 2 of the products ee (30)) and the presence of N of them in (88) the source and sink terms scale as 1/N with rium state can not provide transport ef-`) = n R (!`), but may however support glement at finite non-vanishing temperaother hand, for T L = T R = 0 the state cuum state |vaci in (25) which is clearly ected features of the concurrence C(r, s) easing the distance s r between the spins e maximum achievable bipartite entangles vanishing at decreasing temperatures, in the fact that distance and temperature relations. Furthermore, the lack of transce makes C(r, s) depend not only on s r, position r of the first spin. sen value of g is sufficient to generate entanglement for ∆ = 15, 30, but not for ∆ = 50, while increasing g beyond the saturation value for ∆ = 15 would violate the condition assumed throughout the manuscript that the transition frequencies ω 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 ∆. The graph shows that in this case the highest possible g, despite the higher values of ∆, 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 [21].
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 interspin 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 [97]).
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 [111]: 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 spinchains 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 can be diagonalized as follows. We shall emphasize the rank by writing B N instead of B. Then, one notices that (A3) whence the same equation is satisfied by the associated characteristic polynomial p N (x): Setting x = 2 cos θ, one finds the solution From p 1 (x) = x and p 2 (x) = x 2 − 1, one fixes the coefficients whence p N (x) = sin((N + 1)θ) sin θ .
Since for θ = 0, respectively θ = π, p N (2) = N + 1, respectively p N (−2) = (−) N (N + 1), the only zeroes of p N (x) are at θ k = k π N + 1 , k = 1, 2, . . . , N . It thus follows that the eigenvalues of h = γ I + B are Finally, one can check that the symmetric matrix is also orthogonal. It can be directly checked that (A9) Appendix B: Energies and eigenvectors of two and three spin chains

Three spin chain
In the case of a three-spin chains [97], setting N = 3 the eigenvalues of the Hamiltonian (1) are (B7) Explicitly they and their corresponding eigenvectors read The correspondence with the eigenvalues E j in [97] is as follows Since the 3 × 3 matrix U in (A8) reads Then, σ ± σ z = ∓σ ± and (38) yield The expressions of the eigenvectors |n 1 n 2 n 3 in the spin standard basis and their correspondence with the eigenvectors obtained in [97] are reported in Appendix C 2.
Appendix C: Lindblad operators for two and three spin chains

Two-spin chain
In the case of N = 2, from (47), one computes the following transition frequencies ω 1 = 2(∆ + g) , ω 2 = 2(∆ − g) . (C1) Using (48) and (52), the following Lindblad operators ensue for the open two-spin chain: According to the discussion before Remark 3, in order to have all of them contribute to dissipation, one must set g ≤ ∆.

Three-spin chain
Setting N = 3, from (47) we get the three frequencies They correspond to the three frequencies ω 1 , ω 0 and ω 2 in [97]. Furthermore, (48) and (52) yield the left-bath Lindblad operators and the right bath Lindblad operators According to the discussion before Remark 3, in order to have all of them contribute to dissipation, one must set In terms of the eigenstates | ↑ and | ↓ of σ z , using that |vac = |000 = | ↓↓↓ , one can reexpress the eigenvectors |n 1 n 2 n 3 in the spin standard basis: in the case of zero excitations, while for one excitation, and for two excitations Finally, in the case of three excitations one finds the difference in the overall sign depending on the chosen ordering of the creation operators b † . The correspondence with the |E k obtained in [97] is confirmed after noticing that there |0 = | ↑ and |1 = | ↓ .
Appendix D: Uniqueness of the stationary state Consider the commutator of a generic chain operator X = p,q X pq |p q|, in the energy eigenbasis (34), with the Lindblad operator in (48) and set it equal to zero: This yields whence, choosing r = 0 and s = 0 gives X r0 s1 = 0. Changing A † L (ω ) into A † R (ω ) yields X r1 s0 = 0; thus, the only non-vanishing X commuting with all Lindblad operators must be diagonal in the energy eigenbasis: namely, X = n X nn |n n|. On the other hand, choosing r = 1 and s = 0, (D1) yields X n1 n1 = X n0 n0 for all = 1, 2, . . . , 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 in (86), we need evaluate mean-values of the form Using the expressions (48), (53) and (50) one gets Then, from (26) and (36) it follows that whence, using (32)- (34), Inserting the last three expressions into (E2)-(E7) yields whence, finally using the stationary state eigenvalues in (67) 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, (F1) Then, the combinadic list of the eigenvalues Λ n1n2 of the stationary state ρ ∞ in (96) are for p = 1 and, for p = 2, L

Three-spin chains
For N = 3, let us consider the anti-lexicograhic ordering of the eight 3-digit strings: Application of (95) shows that the previous one and the combinadic ordering coincide, in the sense that Then, for a three spin chain, the eigenvalues Λ n1n2n3 of the stationary state ρ ∞ in (96) are for p = 0, while for p = 1 and, finally, for p = 3, L The quantities in (69) coincide with the quantities s , respectively τ in [97]. According to the correspondence of the frequencies ω 1,2,3 in (C6) with the frequencies ω 1,0,2 in [97], the following correspondence arises among the stationary state eigenvalues Λ n and the eigenvalues µ k in [97]: 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 Since σ + σ z = −σ + , one rewrites where Because of (25), with | ↑ = |1 S and | ↓ = |0 S , one writes where, according to (95), N p uniquely identifies the Ntuple with n j1 = · · · = n jp = 1 and n k = 0 otherwise: such an N -tuple corresponds in its turn to a spin state vector with spin up at the sites j 1 < · · · < 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 ρ ∞ as a tensor product of on-site spin operators, one starts from the projectors Then, writing |vac vac| = N =1 1 − σ ( ) z 2 and using that one recast S j1···jp |vac vac| as As done before, using (95), we identify any given set of indices j 1 < · · · < j p and the corresponding N -tuple with n j1 = · · · = n jp by the unique combinadic integer N p , whence where, setting X 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

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 ≤ i 1 < · · · < i p ≤ 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: by choosing the rows indexed by i 1 < · · · < i p and the columns indexed by j 1 < · · · j p . Here follows some instances of the various entries: where P ijk = P i ⊗ P j ⊗ P k , i, j, k = ± and P ± = 1 ± σ z 2 . The off-diagonal contributions instead read in the case of r = 2 and s = 3. Insertion of (H9)-(H17) into the previous expressions finally yields (I34)-(I36) for both r = 1, s = 2 and r = 2, s = 3, while (I37)-(I39) 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 N (12) p (n , n ) that select the entries of 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 a = L 1 + 2L (2) 2 + L 3 + 4L while C 2,3 = C 1,2 as the coefficients a, c, e coincide with those for ρ (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 (ω ) = h for all = 1, 2, 3. Then, the eigenvalues L where N LR (ω ) := n L (ω ) + n R (ω ) for ρ (12) and ρ (23) , while for ρ (13) . Unlike for the spin (88), (89) and heat (94) flows, that depend on the differences n L (ω ) − n R (ω ), c need not vanish even for identical baths such that n L (ω ) = n R (ω ) = n(ω ) for = 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.