Transport enhancement from incoherent coupling between one-dimensional quantum conductors

We study the non-equilibrium transport properties of a highly anisotropic two-dimensional lattice of spin-1/2 particles governed by a Heisenberg XXZ Hamiltonian. The anisotropy of the lattice allows us to approximate the system at finite temperature as an array of incoherently coupled one-dimensional chains. We show that in the regime of strong intrachain interactions, the weak interchain coupling considerably boosts spin transport in the driven system. Interestingly, we show that this enhancement increases with the length of the chains, which is related to superdiffusive spin transport. We describe the mechanism behind this effect, compare it to a similar phenomenon in single chains induced by dephasing, and explain why the former is much stronger.


Introduction
Since the experimental discovery that quantum coherent dynamics is present in excitation energy transport in biological light harvesting complexes [1] and theoretical work demonstrating that environmental fluctuations can be used to optimise transport efficiency [2,3], a great deal of interest has focused on the beneficial consequences that the unavoidable coupling of a quantum system to its environment can have [4,5,6,7,8,9,10,11]. In particular, it has been found that the optimal regime for excitation transport through various systems consists of a balance between coherent and incoherent phenomena, or, more specifically, the interplay between coherent electronic dynamics and the vibrational environment [2,3,12,13]. Most of this effort to characterise and understand environment-assisted transport has been restricted to single-particle effects [2,3,12,14,15], given that light-harvesting complexes under physiological conditions usually contain very few excitations at the same time due to the low intensity of ambient sunlight [16,17]. Nevertheless, since the interplay between coherent and incoherent phenomena is relevant beyond a biological context, it is important to consider its impact on the transport properties of many-body systems.
It has been known for many years that coherent and incoherent particle transport processes take place in various condensed matter systems. These include cuprates [18,19,20] and several organic conductors such as conjugated polymers [21,22], layered organic metals [23,24] as well as Bechgaard and Fabre salts [25,26,27]. As illustrated in figures 1(a),(b), these systems share a common feature, which is a highly anisotropic structure consisting of lattice sites strongly coupled in one or two directions and weakly coupled in other directions. At not-too-low temperatures, particle transport along the strongly coupled directions is predominantly coherent, while the transport along directions with weak coupling occurs via incoherent hopping processes. Due to the many-body nature of most of these systems, it is expected that they show a rich variety of phenomena emerging from the interplay between the different transport mechanisms and particle-particle interactions.
Several interesting effects resulting from the competition between coherent and incoherent processes in the presence of strong interactions have recently been found [28,29,30]. In particular, previous research by some of us has uncovered a novel mechanism of dephasing-enhanced transport in linear homogeneous strongly interacting systems [31,32]. Nevertheless, the physics resulting from this interplay still represents relatively uncharted territory, and anisotropic condensed matter systems displaying both coherent and incoherent transport processes arise as attractive and natural scenarios to take the exploration even further.
To unravel the competition between coherent and incoherent dynamics in these systems, we propose a concrete minimal model that mimics their structure, as shown in figure 1(c). Specifically, we consider excitation transport through arrays of incoherently coupled one-dimensional quantum spin chains, including the possibility of strong interactions between excitations on the same chain. We study DC transport properties by imposing a net current flowing in one direction through the system. We find that the effective environment furnished by nearby chains significantly enhances intrachain transport for sufficiently strong interactions between excitations. In addition, such a transport enhancement increases with the size of the chains, indicating the relevance of the mechanism for bulk materials. Furthermore, the incoherent hopping of spin excitations between chains results in a much more pronounced effect than that produced by pure dephasing due to, for example, lattice vibrations [31,32].
The paper is organised as follows. In Section 2 we describe the model to be studied and the approximations considered. In Section 3 we study weakly interacting systems, where the incoherent coupling only degrades the transport. In Section 4 we analyse the case of strong interactions, where current enhancement due to incoherent coupling is observed. The origin of this effect is explained, and compared to that resulting from Illustration of systems of incoherently coupled chains. The red solid and green dashed lines represent interchain and intrachain hopping, respectively. (a) Excitation transport in a conjugated polymer system. Interchain hopping of excitons, observed to be incoherent [21,22], occurs where polymeric chains are close to each other. Our lattice model (in the non-interacting case) can be relevant for describing exciton dynamics in conducting polymer chains arranged in a highly ordered configuration, similar to those recently achieved by chemical doping [33,34]. (b) Excitation hopping in organic salts. Coherent hopping occurs along one direction of the system, while incoherent hopping takes place along another direction, mediated by a scaffolding of different molecules. (c) System of three incoherently coupled spin chains, with intrachain hopping τ , interaction strength ∆ and interchain coupling γ. The blue arrows represent the right-to-left driving of the system, while the red arrows correspond to the left-to-right driving.

Model of non-equilibrium incoherently coupled spin chains
In this work we consider a N × Λ rectangular two-dimensional (2D) lattice, consisting of Λ chains with N sites each (see figure 1(c)). The model for coherent intrachain transport should describe conserved excitations that can hop between lattice sites and Transport enhancement from incoherent coupling between one-dimensional quantum conductors4 interact with each other. We therefore choose the simple spin-1 2 XXZ Hamiltonian to govern the dynamics of each chain [35,36] where the super-index (λ) refers to any operator of chain λ, σ k(λ) i (k = x, y, z) are Pauli matrices at lattice site i of chain λ, τ is the exchange coupling between nearest neighbours (in the following we take units of energy and time such that τ = 1 and = 1), and ∆ is the anisotropy (we consider ∆ > 0 only), where both parameters are assumed to be the same for every chain. The presence of an excitation at a certain site corresponds to a spin pointing up, while the absence of an excitation corresponds to a spin pointing down. The hopping is encapsulated by the first two terms of equation (1), while the final term corresponds to an energy penalty for nearest-neighbour lattice sites in the same spin state, creating an interaction between spin excitations. In this sense, we will refer to a strongly interacting model if ∆ > 1, which is an energy-gapped regime.
We assume that hopping also occurs between nearest-neighbour sites of neighbouring chains, with a hopping rate η. Let t φ denote the time taken for a spin excitation to lose its phase coherence, either from collisions with other spin excitations on the same chain or due to dephasing induced by an external bath, e.g. phonons. If η t −1 φ , it is reasonable to neglect quantum correlations between sites of neighbouring chains, and treat the interchain coupling as a purely incoherent hopping process. This is expected to be a good approximation for temperatures intermediate between the two hopping energy scales, i.e. η k B T τ . This is easily satisfied in several systems. For example, η ∼ 100K and τ ∼ 1000K for typical Bechgaard salts [37].
By means of a Jordan-Wigner transformation, the present model can be mapped onto an interacting spinless fermion system [38], as we demonstrate in Appendix A. The parameter τ then corresponds to nearest-neighbour hopping, and τ ∆ to nearest-neighbour Coulomb repulsion, up to factors of order unity. Normally such a transformation is not feasible in a 2D system, due to the appearance of non-local Jordan-Wigner string operators (see equation (A.8)) which enforce the correct exchange phase between fermions at different sites. However, due to the purely incoherent nature of the interchain coupling, there is no need to maintain a definite phase relation between fermion states localised on different chains. This equivalence between fermion and spin representations makes our model relevant for describing not only spin transport, but also particle transport of hard-core bosons or fermions.
We describe the combination of coherent and incoherent dynamics by a quantum master equation of Lindblad form [39]: where ρ is the density matrix of the total system, H = λ H (λ) is the total Hamiltonian, and L(ρ) is the dissipator describing the interaction of the spin chains with the environment and each other. The dissipator is given by with L k the jump operators describing each incoherent process and {., .} the anticommutator of two operators. The incoherent coupling between two spin chains is modelled by the jump operators representing the transfer of a spin excitation from site i of chain λ to site i of chain µ = λ ± 1, with rate γ. Simple Golden Rule arguments [25] indicate that the incoherent hopping rate is of order γ ∼ η 2 /t φ . Due to the large number of factors that can contribute to this hopping rate (e.g. temperature, collision rate, interchain distance etc.), we treat γ as a free parameter that can be varied independently.
To analyse the transport properties of this system, we drive it into a non-equilibrium configuration by coupling its boundaries to unequal reservoirs, as depicted in figure 1(c). This driving scheme imposes a magnetisation imbalance on each chain, and thus induces a spin current. We assume that the correlation time of the reservoirs is negligibly small, so that the energy dependence of the incoherent transition rates may be neglected. We also assume that the points of contact between the bath and any pair of neighbouring chains are further apart than the correlation length ‡, leading to independent driving reservoirs for each chain. Under these conditions, the action of the reservoirs driving the system out of equilibrium is represented by the following Lindblad operators [40,41] L +(λ) L,R = Γ(1 ± f )/2σ Here, σ ), Γ is the strength of the coupling to the reservoirs (we choose Γ = 1 in all our numerical calculations), and f is the driving parameter. The driving operators are such that when applied to the boundary spins in isolation, they induce a state with magnetisation σ z(λ) 1 = f and σ z(λ) N = −f . At f = 0 there is no magnetisation imbalance between the boundaries of the chains, so there is no net spin transport. As f increases, so does the imbalance between the boundaries, forcing a spin current to flow from the left to the right boundary of each chain. Our simulations are performed with the weak driving f = 0.1, thus staying in the linear response regime [31,40].
Due to the finite temperature, we would also expect local dephasing processes to exist, described by jump operators of the form ‡ The correlation length is defined by λ c = c/ω c , where c is a characteristic velocity of bath excitations and ω c is the frequency of the most energetic bath mode that the system interacts with appreciably. For example, for a fermionic bath consisting of a macroscopic metal lead the correlation length is essentially the Fermi wavelength, which gives λ c < 10 −9 m for typical carrier concentrations greater than 10 28 m −3 .
with γ d the dephasing rate. However, in order to simplify the analysis, in most of the paper we will assume that, apart from driving, the only effect of the environment is to generate incoherent hopping between chains. We show in Section 4.2 that the large current enhancement induced by incoherent coupling cannot result solely from dephasing processes. We also show in Appendix C that our qualitative conclusions about the enhancement due to incoherent coupling remain valid in the simultaneous presence of dephasing.

Mean-field approximation
To gain insight into the properties of the system, we calculate its non-equilibrium steady state (NESS), which emerges in the long-time limit of equation (2) from the interplay between coherent and incoherent processes. Computing the NESS for a strongly correlated two-dimensional system represents a formidable challenge, therefore an approximation scheme is necessary. In Appendix B we present an exact solution for two incoherently coupled chains in the non-interacting limit ∆ = 0, demonstrating that the NESS factorises as ρ = ρ 1 ⊗ ρ 2 + O(f 2 ). Since magnetisation and current expectation values are of order O(f ), the lowest order contribution to the NESS is sufficient to compute transport observables accurately. This observation motivates the following mean-field approximation (MFA), according to which the state of the entire system is a direct product of the states of each spin chain thus discarding both quantum and classical correlations between different chains. Using this mean-field ansatz, we can obtain the master equation for each chain separately after tracing out the state of the other chains. This provides a considerable advantage in numerical simulations, which would be very demanding if all correlations between the chains are kept in the description. The resulting MFA master equation for chain λ is where the superscripts d and i refer to driving and incoherent coupling terms, respectively. Within the MFA, the incoherent interchain coupling turns into local gain and loss processes at each chain, with rates depending on the magnetisation of the neighbouring chains. In this form, the density matrix of each chain can be evolved separately from the others by the simulation of its own master equation, with the coupling to neighbouring chains being effectively described by expectation values of local operators. This type of evolution can be performed efficiently by means of a parallel implementation of the mixed-state Time Evolving Block Decimation (TEBD) algorithm [42,43]. Our code is based on the open-source Tensor Network Theory (TNT) library [44].
In order to verify that the MFA gives reasonable results, we have also performed TEBD simulations of two coupled chains without the MFA for comparison. In figure 2 we plot the steady-state currents (defined in Section 2.3) obtained within each approach for a pair of incoherently coupled chains of length N = 20. The two sets of results are clearly in close agreement, however the accuracy of MFA calculations is higher for smaller values of γ and ∆. The maximum error is 3.8%, when ∆ = 2 and γ = 1.2. We are therefore confident that this approximation gives accurate results for greater numbers of coupled chains, when quasi-exact TEBD simulations are not feasible.

Approximation for an infinite number of coupled chains
In the case of an infinite number of chains, the reduced density operators of all the chains are exactly the same at any time. So as observed from equation (8), the problem of simulating the evolution of the entire system is reduced to that of performing the calculation for a single chain coupled twice with itself. The resulting Lindblad master equation of each chain, describing an effective non-linear self-consistent time evolution, is where the index (λ) has been dropped for simplicity, and are the effective gain and decay rates at site m.
Transport enhancement from incoherent coupling between one-dimensional quantum conductors8

Spin current
We now derive the expression for the spin current through the system. It is obtained from the local magnetisation rate of change, calculated from the master equation directly. For site i of chain λ, we have This expression is equivalent to that of the spin current through a 1D spin chain [40].
Here, in contrast to that case, the longitudinal spin current is site-dependent in the NESS. For example, in the bulk of the system, 1 < i < N , the difference of spin currents through nearest neighbours is Now consider, for example, the left boundary i = 1. Since equation (12) is not defined for i = 0, in the NESS equation (11) gives A similar equation holds for the right boundary i = N . This leads to a natural definition of boundary currents j N , which allow equation (13) to be valid along the entire chain. These currents, which indicate the direct injection and ejection of spin excitations on the chain by the boundary reservoirs, are thus given by We can associate the right hand side of equation (13) to the difference of spin flows between chain λ and its neighbouring chains, by defining a transversal spin current In this form, equation (13) can be rewritten as This balance between the longitudinal and transversal currents is illustrated in figure  3(a). From equation (13) it follows that in the absence of incoherent coupling, the current through each chain is homogeneous in the NESS, j (λ) i = const, i = 0, . . . , N . In addition, from straightforward calculation, it is easily shown that in the presence of incoherent coupling, the total current per site that J i − J i−1 = 0. Also note that due to the symmetry of the system considered, ⊥,i = 0 for all sites i. A concrete example of the longitudinal spin current profiles (equation (12)) for three chains is shown in figure 3(b), including the boundary currents defined in equation (15). Due to the symmetry of the incoherent coupling, the state and thus the currents of chains λ = 1, 3 are equal. The corresponding transversal spin currents, defined in equation (16), are shown in figure 3(c). When moving from the boundary sites i = 1, N towards the centre of the system, the currents through chains λ = 1, 3 significantly increase at the expense of the current in the middle chain. This strong site dependence is reflected in large transversal currents, flowing in opposite directions. In the central sites of the system, the transversal currents are very small since the local magnetisations of neighbouring chains are very similar (see equation (16)). This is expected since the magnetisation of each chain must pass through zero at the same position in the centre due to the symmetric driving.
The NESS spin currents through the system, together with the magnetisation profile, determine the nature of the transport. If it is diffusive, the currents satisfy a diffusion equation: j Transport enhancement from incoherent coupling between one-dimensional quantum conductors10 where κ (λ) is the (N -independent) spin conductivity of chain λ and ∇σ z(λ) i is the magnetisation gradient. On the other hand, if the transport is ballistic κ (λ) diverges, resulting in a size-independent spin current. Ballistic transport has been observed in single dephasing-free chains when |∆| < 1 [45,40], while diffusive transport has been found in the linear response regime for |∆| > 1 and no dephasing [45,40], and for finite dephasing and any interaction strength [41,46,47]. Note that since the transversal current of equation (16) is proportional to the local magnetisation gradient along the transversal direction, it is diffusive by construction.
We now discuss the spin transport properties of the system in both the weakly (|∆| < 1) and strongly (|∆| > 1) interacting regimes, which show a completely different behaviour in the presence of incoherent interchain coupling. For this, instead of observing the spin current through each chain, we consider the total spin current per chain, noted by J, i.e., J = J i /Λ. This quantity, to which we refer in the rest of the paper simply as the spin current, is shown in figure 3(b) for a particular example. Its homogeneity along the system is a good indication of the obtention of the NESS.

Transport in weakly interacting incoherently coupled spin chains
We initially consider the non-interacting case ∆ = 0. The analytical method presented in Refs. [41,46] can be extended to two incoherently coupled non-interacting chains, as explained in Appendix B. This allows us to extract the exact current and magnetisation expectation values. The former is given by while the magnetisation profile is linear in the bulk, see equation (B.11). These results agree with TEBD simulations, as indicated in figure 4. Note that if γ = 0, the current is independent of the size of the spin chains, indicating ballistic transport [41,46]. On the other hand, a finite incoherent coupling induces a decay of the current with the length of the system ∝ N −1 , typical of a diffusive conductor with homogeneous magnetisation gradient. In fact the bulk conductivities are easily shown to be κ (1,2) = 8/γ. So, similarly to dephasing processes on a single non-interacting chain [41,46], incoherent interchain couplings induce a non-equilibrium phase transition between ballistic and diffusive regimes, with a spin current monotonically degraded by the interchain hopping.
In the limit of an infinite number of chains, described by the self-coupled chain (see Section 2.2), the analytical method used for two chains can also be applied (see Appendix B). We thus obtain exact expressions for the current and magnetisation, given by equations (19) and (B.11) respectively, with γ/2 instead of γ/4. The conductivity is thus κ = 4/γ, reduced compared to the case of Λ = 2. This is because each chain has two nearest neighbours rather than one, leading to a stronger degrading effect of the incoherent hopping.
In the intermediate case, i.e. for a finite number of chains Λ > 2, we use the TEBD method to obtain the NESS of the system. Characteristic results for the current and for the magnetisation profiles are shown in figure 4. The same qualitative features of the cases Λ = 2 and Λ → ∞ are found, namely the spin current monotonically decreases with γ and the magnetisation profile is almost linear in the bulk. In addition, for a fixed incoherent coupling, the current decreases with the number of chains, rapidly for small values of Λ and very slowly for large values. An extrapolation of these results to the limit Λ → ∞ agrees with the analytical approach for the self-coupled chain, as shown in figure 4.
Similarly to the cases of two and an infinite number of chains, the transport response induced by incoherent coupling on a system of several chains is characteristic of diffusive conductors. To see this explicitly, we observe that each chain satisfies the local diffusion equation (18), with ∇σ the local magnetisation gradient. Since the spin current through each chain is sitedependent, we have verified that the ratio j (i.e. the conductivity) is homogeneous for each λ, so the diffusion equation (18) holds. In addition, similarly to the analytically-solvable cases, the conductivity of each chain decays monotonically with the incoherent coupling rate, with a behaviour very well described by a decay κ (λ) ∝ 1/γ, as shown in figure 5. We also note that for 2 < λ < Λ − 1 the conductivity is almost indistinguishable from that of λ = 2, Λ − 1, due to the weak effect of the boundary chains. Now we consider weak interactions 0 < ∆ < 1. We find that the effect of incoherent interchain coupling on the system is very similar to that on non-interacting chains. Namely a finite incoherent coupling induces a transition from ballistic (γ = 0) to diffusive (γ > 0) behaviour. The magnetisation profiles become linear, and the spin current and the conductivities of each chain decrease monotonically with γ, the latter following a power law as shown in figure 5. The current also decreases with Λ, approaching a limiting value when Λ → ∞. In figure 5 it is also seen that for fixed values of Λ and γ, the spin transport diminishes as ∆ increases, a known result for single chains in the massless regime [35,47].

Transport enhancement for strong intrachain coupling
We now consider the effect of incoherent interchain coupling on the transport properties of strongly interacting spin chains (∆ > 1).

Environment assisted transport
Due to the strong correlations between spin excitations, the regime ∆ > 1 presents a completely different response to environmental effects to the case of weak interactions ∆ < 1. It has been found [31] that for single 1D chains, dephasing processes can lead to a surprisingly large enhancement of the current even at weak driving. Now we show that the ability of excitations to jump incoherently across different chains leads to an even larger transport enhancement, which constitutes the main result of our work.
As shown in figure 6, the presence of incoherent interchain coupling increases the spin current through the system, compared to that of γ = 0, for a wide range of rates γ. The optimal coupling maximizing the current γ opt , which is obtained by fitting the peak to a polynomial function, strongly depends on the interaction strength, increasing with ∆. Similarly, the current enhancement grows with ∆. For example, the current increases by up to 39% for ∆ = 1.2 and up to 91% for ∆ = 2.
We now consider the effect of the system size on the transport enhancement. The optimal current J opt remains almost constant for all values of Λ considered. In addition, the optimal coupling monotonically decreases with Λ as shown in figure 7, and the range of beneficial couplings narrows. Importantly, extrapolations to Λ → ∞ strongly suggest the existence of a finite optimal incoherent coupling γ ∞ opt in this limit. We have confirmed this result from simulations of a self-coupled chain with ∆ = 2, as shown in the inset of figure 7. The results of both approaches agree very well, giving optimal couplings of γ ∞ opt = 0.146 for the extrapolation and γ opt = 0.144 for the self-coupled chain. Because of this agreement, we henceforth denote by γ ∞ opt the optimal coupling for the self-coupled chain.
To observe how the enhancement effect scales with the length of the system, we have analysed both the optimal current and the optimal coupling for self-coupled chains with different values of N . We found that an exponential decay of the optimal coupling of the form γ ∞ opt = ae −bN + c yields a finite optimal coupling in the thermodynamic limit of γ ∞ opt (N → ∞) = c = 0.057 (9). Nevertheless, a power law decay of the form γ ∞ opt = aN −b also fits well to our results. This means that we are not able to assess whether the optimal coupling is finite for an infinite system. The scaling results indicate, however, that for very large but finite systems, an enhancement effect of the current is still expected for very small incoherent couplings. This does not mean that the increase of the current becomes less important as the system gets larger. In fact, although it is restricted to a narrower range on incoherent coupling rates, the enhancement effect becomes stronger as the size of the system increases. This is shown in figure 8, where the enhancement factor J opt /J(γ = 0) is seen to increase with N . We therefore expect that spin transport can be significantly enhanced by environmental processes even in macroscopic (anisotropic) two-dimensional systems, and thus can be observed experimentally in bulk materials.
To understand the origin of the increase of the enhancement ratio with N , it is important to study the nature of the spin transport in the enhancement regime. To address this point we analyse the scaling with N of the spin current through a strongly interacting self-coupled chain. The results are shown in figure 9. In the presence of diffusive spin transport, the magnetisation profile of the chain is linear in the bulk. The magnetisation gradient is thus homogeneous, and defined as Figure 8. Spin current enhancement factor J opt /J(γ = 0) as a function of the length of the system for a self-coupled chain with ∆ = 2 (•). As N increases, the enhancement factor is well described by a linear scaling, as indicated by the solid line. The optimal current is seen to scale as J opt = aN −b , with a = 0.172(3) and b = 0.830 (5) leading to the increase of the enhancement factor with N .
where N − 5 corresponds to removing two sites from either end to diminish boundary effects. Diffusive spin transport is evidenced if the system satisfies the diffusion equation with κ the spin conductivity and α = 1. As shown in figure 9 for weak incoherent coupling, the results of our TEBD simulations are well described by this equation, but with α < 1. This means that in the regime of transport enhancement, the system presents a spin current which decreases with N slower than normal diffusion, i.e. it shows superdiffusive behaviour. Thus the optimal current also shows a slower decrease than that of a diffusive conductor. Since for the diffusive regime J(γ = 0) ∝ N −1 , a divergence of the enhancement ratio with N is found. Finally, as shown in the upper inset of figure 9, the exponent α gets closer to 1 when increasing γ, the transport thus tending towards being described by normal diffusion when the incoherent effects become stronger. When γ is too large the enhancement effect disappears, since the system is perturbed so frequently that it is prevented from evolving, i.e. the Zeno effect emerges [39].

Enhancement mechanism
We now discuss the origin of the transport enhancement in the strongly interacting regime, which is similar to that found in single chains due to dephasing processes [31,32]. For interaction strengths ∆ > 1, the spectrum of the XXZ Hamiltonian (1) consists of several bands whose energetic separation is proportional to ∆ (see figure 10). The highest bands are almost flat, possessing very low conductivity. These bands correspond to bound states of spin excitations, where several spins are clumped together, thus having large potential energy.
When the system is driven out of equilibrium, population is transferred to various eigenstates, depending on the strength of the driving. For example, at f = 1 only the highest bands are populated, leading the system to an insulating NESS [31,40]. Even in the weak-driving regime as considered here, some population is transferred to the highest bands, which then gives a small contribution to the conduction of the system. However, if energy-dissipating processes take place in the system, transitions from these slow bands to lower bands of larger conductivity are induced, leading to an enhancement of the current. In other words, if the energy of spin bound states is dissipated, these break into states of lower potential (and total) energy, but of much higher kinetic energy, thus increasing the conductivity.
The enhancement described in our work emerges from the energy dissipation induced by the incoherent interchain coupling. To clarify this point, consider for simplicity the self-coupled chain configuration described by equation (9) §. A straightforward calculation of the energy dissipation rate corresponding to the § The calculation of the rate of energy dissipation is also easily performed for a finite number of chains. In this case, we obtain a sum over λ of terms like those of equation (23), but instead of the magnetisation at site i + 1 of each chain, the magnetisations of the neighbouring chains at site i + 1 appear. Nevertheless, since the magnetisations of all chains are similar, equation (23) corresponds to a good approximation. incoherent coupling gives where L inc (ρ) is the dissipator describing the incoherent coupling. The first term in the energy dissipation rate appears due to the loss of phase coherence between neighbouring sites [31], and is proportional to the hopping energy The second term is proportional to the sum of nearest-neighbour connected spin-spin correlations and corresponds to a direct dissipation of the interaction energy due to the incoherent hopping, which rips spin excitations away from their nearest neighbours .
Intuitively, if C is positive, the spin excitations of the system are bunched together on average, while if C < 0 the excitations are spread out. Therefore, the sign of C gives a simple indication of how population is distributed between the bound states (bunched) and mobile states (spread out). In the absence of incoherent coupling, we have found As observed during the derivation of equation (23), the terms σ z i σ z i+1 are a direct consequence of the non-conserving nature of the jump operators (i.e. they result from any incoherent process described by jump operators σ + or σ − ). In addition, the terms σ z i σ z i+1 appear because the effective ratesγ + m andγ − m are different (see equation (10)). that C undergoes a marked transition from taking negative to positive values as the interaction strength crosses the critical point ∆ = 1 (see figure 11(a)). This behaviour is a manifestation of the well known non-equilibrium phase transition from ballistic to diffusive conduction at ∆ = 1 [31,40], and demonstrates a tendency of the spin excitations to clump together in the strongly interacting regime of the driven system. Importantly, even when the incoherent coupling is incorporated, our simulations always show that C > 0 when ∆ > 1 (see figures 11(b),(c)). More precisely, C diminishes with γ, indicating the decrease of population in bound (correlated) states with incoherent processes, but remains positive. Similarly, we always found that K > 0 in the strongly Figure 12.
Comparison of the rate of energy dissipation due to bulk dephasing processes (equation (24)), and incoherent self coupling (equation (23)). Both cases correspond to ∆ = 2 and N = 40. For these parameters, the current for a single isolated chain is J ≈ 4.2 × 10 −3 . Inset: Corresponding spin currents.
interacting regime. This indicates that energy is being dissipated from the system due to the incoherent interchain hopping (see equation (23)), transferring population from bound to mobile states and thus leading to transport enhancement.
It is important to note that the transport enhancement induced by incoherent interchain coupling is much larger than that of pure dephasing processes described by equation (6). For example, as shown in figure 12 for ∆ = 2, the spin current is increased up to 37% by dephasing [31], and up to 91% by incoherent coupling. This difference can be explained by looking at the energy dissipation rate due to dephasing, with L d (ρ) the dephasing dissipator. This rate is compared to the dissipation rate from incoherent coupling in figure 12. Since its maximal magnitude is significantly smaller than that of incoherent coupling, more energy is dissipated by the latter, resulting in more population transfer from flat to mobile energy bands and thus to a larger current. This also shows that for ∆ = 0, the effects of incoherent coupling cannot be reproduced just by dephasing processes. In Appendix C we also show that in the simultaneous presence of dephasing and incoherent coupling, the latter dominates the energy dissipation, and the current enhancement is still larger than that of dephasing alone.

Summary & Conclusions
We have studied the spin transport in an anisotropic 2D spin− 1 2 lattice driven out of equilibrium by Markovian boundary reservoirs. The assumption of highly anisotropic coupling allowed us to consider the system as an array of incoherently coupled chains. Each chain is described by an XXZ Hamiltonian, which contains the basic elements that constitute many-body lattice systems, namely particle hopping and interactions. Employing a mean-field approximation, we calculated the spin current and magnetisation of the NESS of the system for several parameters. This approximation, found to reproduce the transport properties of two coupled chains, facilitates an accurate and efficient dynamical simulation of the system using a parallel implementation of the TEBD algorithm [44].
We found that in the presence of weak intrachain interactions, the incoherent coupling monotonically degrades the spin conductivity of the chains. However, in the strongly interacting regime we found a significant transport enhancement due to the incoherent coupling. This enhancement can be understood as the result of incoherent transitions from bound states to mobile bands of energy eigenstates, similar to the effect of dephasing on spin and heat transport in 1D systems [31,32]. However, the direct breakdown of bound states by the incoherent hopping between neighbouring chains, which opens more paths for spin excitations to flow, provides a greater improvement than dephasing effects alone.
A self-consistent extension of the mean-field approximation enabled us to perform simulations directly in the limit of an infinite number of chains. In this configuration, we found that the enhancement of the spin current increases with the size of the system, which reveals the importance that the enhancement effect can have in bulk materials. The origin of this scaling was related to the existence of superdiffusive transport in the regime of current enhancement, becoming closer to normal diffusion as the incoherent coupling increases.
We emphasise that the model considered in the present work does not account for the exact processes taking place in a particular material. Instead, it incorporates in simplified form the basic elements identified to exist in several systems such as organic conductors and cuprates, namely coherent and incoherent transport mechanisms in different directions. We thus expect that the observations that we have presented can be relevant to qualitatively explain physical processes in these materials. We start from a fermionic tight-binding model, describing Λ incoherently coupled 1D chains of N sites each. Using the Jordan-Wigner transformation, we can map this system onto the spin model described in the main text. We consider only the incoherent interchain coupling, since the Jordan-Wigner mapping for the boundary driving Lindblad operators is detailed in Ref. [40]. To carry out the proof, it is simplest to work in the Heisenberg picture. The evolution of an operator O is given by The Hamiltonian of each chain is where the ladder operators satisfy {c j . The adjoint dissipator describing the hopping between sites j of chains λ and λ + 1 is Each site of the system is associated to an index pair (j, λ) specifying its position, which we order by the following prescription. If µ < λ then (k, µ) < (j, λ) for all j and k. If µ = λ then (k, µ) < (j, λ) only if k < j. Now we can define the spin representation of the fermion ladder operators which satisfy the required anticommutation relations by construction. Applying this transformation to the Hamiltonian of each chain and discarding a constant energy shift yields where H (λ) is given by equation (1). The second term represents a homogeneous magnetic field, which makes no difference to steady-state magnetisation and current expectation values and can be thrown away [32]. The third term represents magnetic fields acting on the boundary sites, which can be neglected for large N . We have checked numerically that the effect of incorporating these boundary fields disappears as N increases. Now we consider how the transformation acts on the dissipators (A.3). The Lindblad operators transform as where the Jordan-Wigner string is defined by is defined by equation (4). Note also that [L The "sandwich" terms transform, for example, as The only observables we consider are linear combinations of the operators n and must therefore be considered in more detail. The string S (λ) j contains all σ z operators acting on the sites from (j + 1, λ) to (j − 1, λ + 1), inclusive. Therefore, the strings appearing in the sandwich terms of A j+1 unless (k, µ) = (j, λ) or (k, µ) = (j + 1, λ − 1). Nevertheless, in these two special cases the sandwich terms identically vanish and therefore and similar relations hold for the other potentially troublesome sandwich terms. Gathering all the results, we find that the evolution equation for each observable of interest can be written as Transforming back to the Schrödinger picture, we arrive at the master equation described in Section 2. Dephasing terms can also be straightforwardly included, since the dephasing Lindblad operators σ z(λ) j = 2n (λ) j − 1 have a local representation in terms of both spins and fermions.
Appendix B. Analytic approach for incoherently coupled noninteracting spin chains Following the method proposed byŽnidarič [41,46], we derive an analytic approximation for the state of two limit cases of the system of non-interacting (∆ = 0) incoherently coupled spin chains, namely for Λ = 2 and Λ → ∞. We consider first the case of two chains, and propose the following ansatz for the NESS of the entire system: (B.1) Here I is the identity operator of the entire system, and A and B are functions of spin operators in z direction and spin currents, respectively: As seen below, A and B scale with f , so the approximation only gives the NESS up to first order in the driving strength. Nevertheless, this is enough to obtain exact results of the current and magnetisation of the chains for all values of f , since to obtain r-point correlation functions, an expansion up to order r is needed [41]. It is easily shown that the local magnetisation and the spin current are given by where L H (ρ) = −i[H, ρ], and L driv (ρ) and L inc (ρ) are the Lindblad terms corresponding to driving at the boundaries and interchain coupling, respectively, with the jump operators of equations (5) and (4). Now we introduce the ansatz (B.1) in the master equation. Since both chains are equal, they have the same dynamics and steady state, so we set a (λ) j = a j . We then obtain the following results for each process: Transport enhancement from incoherent coupling between one-dimensional quantum conductors24 where we divided the contribution of the driving into its left (L) and right (R) components (L driv (ρ) = L L driv (ρ) + L R driv (ρ)), and To obtain the N + 1 coefficients b and a i we need N + 1 different equations, which result from equating to zero the coefficients in front of each operator. Explicitly, in front of σ z(λ) 1 and σ z(λ) N we have Similarly, the coefficient in front of j The solution of this system of equations gives the magnetisation in the bulk (i > 1) Note that up to O(f ), the total state of the system ρ is a product of the states of each chain (ρ 1 and ρ 2 ). So, if we have it follows that, up to O(f ), ρ = ρ 1 ⊗ ρ 2 (mean-field approximation), with ρ given by the ansatz of equation (B.1). For the case Λ → ∞ with homogeneous incoherent coupling γ, described in Section 2.2 (self-coupled chain), we follow a similar process. Assuming an ansatz for the NESS of the chain like that of equation (B.1), with normalisation to 2 N , we obtain 14) The NESS is then equivalent to that of of two chains, but with γ/2 instead of γ/4 in the factors a i and b. Figure B1. Spin current through a strongly interacting system (∆ = 2 and N = 40) when both dephasing and incoherent self-coupling take place, as a function of γ, for different dephasing rates (γ d = 0.52 is the optimal dephasing in the absence of incoherent coupling. Note that here γ d is defined to be twice that of Ref. [31]).

Appendix C. Simultaneous presence of incoherent coupling and dephasing
We have performed simulations of chains with incoherent self-coupling and dephasing at the same time, and verified that the former tends to be dominant. In figure B1 we show the spin current through the system as a function of the incoherent coupling, for fixed dephasing rates. For all the cases the currents are reduced from those of γ d = 0, but are still larger than those of dephasing alone (compare to the inset of figure  12). This indicates that the current enhancement induced by the incoherent coupling is still the dominant mechanism. This conclusion is reinforced by looking at the rates of energy dissipation corresponding to both the incoherent coupling (equation (23)) and dephasing (equation (24)). We found that for both a large and a small dephasing rate, the amplitude of the energy dissipation rate coming from the incoherent coupling is significantly larger than that of dephasing, except for very small incoherent couplings.