We study a non-interacting, one-dimensional lattice of L sites with Hamiltonian given by the X X model H = J ∑ j = 1 L − 1   σ j x σ j + 1 x + σ j y σ j + 1 y + ∑ j = 1 L ε j σ j z . (1)

In Equation 1 { σ j x , σ j y , σ j z } j = 1 L are Pauli matrices, J is the hopping term and ε j is a local field which will be specified shortly. In order to induce a non-equilibrium state in the system, we couple the boundary sites to high-temperature reservoirs with different chemical potentials; in addition, each site is also exposed to its own local reservoir which induces dephasing. All the reservoirs are assumed to be ideal and couple weakly to the system. The dynamics of the system’s density operator ρ ( t ) at time t is given by the Lindblad master equation [2, 3, 30]. d ρ t d t = − i ℏ H , ρ t + ∑ α = l , r , d L α ρ t . (2)

Here the superoperator L α describes dissipation induced by the reservoirs α = l , r , d corresponding to left boundary, right boundary and dephasing, respectively. It is a sum of local Lindblad jump operators of the form L α   ⋅ = ∑ j , β L j α β ⋅ L j α β † − 1 2 L j α β † L j α β , ⋅ , (3) where j in Equation 3 labels the site acted on by the reservoir α and β labels different jump operators on that same site. The effect of each boundary reservoir ( α = l , r ) is represented by two jump operators L 1 l ± = Γ 1 ± f 2 σ 1 ± , L L r ± = Γ 1 ∓ f 2 σ L ± , (4) where Γ in Equation 4 is the coupling rate at the boundaries, f is the chemical potential bias and σ j ± = ( σ j x ± i σ j y ) / 2 . For forward bias ( 0 ≤ f ≤ 1 ) excitations are mostly created on the first site and annihilated on the last site; the opposite reasoning applies for reverse bias ( − 1 ≤ f ≤ 0 ) . In order to describe the effects of dephasing, we consider L jump operators acting on individual sites as if each of them were coupled to its own reservoir L j d = γ 2 σ j z , j = 1,2 , … , L , (5) where γ in Equation 5 is the dephasing rate.

The key observables in the system are the magnetization and its associated current. We denote by ⟨ A ⟩ ( t ) ≡ T r [ A ρ ( t ) ] the expectation value of an observable A at time t . Differentiating with respect to time and using Equation 2 yields a continuity equation (see Supplementary Material for more details) d ⟨ σ j z ⟩ t d t = ⟨ I j − 1 ⟩ t − ⟨ I j ⟩ t + δ j , 1 Γ f − ⟨ σ j z ⟩ t − δ j , L Γ f + ⟨ σ j z ⟩ t , (6) involving the expectation value of the magnetization current ⟨ I j ⟩ ( t ) flowing from j to j + 1 , where Equation 7 I j = 2 J ℏ σ j x σ j + 1 y − σ j y σ j + 1 x (7) defines the corresponding current operator. Note that the dephasing reservoir does not change the average magnetization and therefore does not contribute explicitly to Equation 6, although it still affects the average magnetization and currents implicitly. Moreover, the boundary reservoirs do not contribute to the expression of the current in the bulk of the chain.

When the steady state is reached, the left-hand side of Equations 2, 6 vanishes. We denote by ⟨ A ⟩ ∞ ≡ T r [ A ρ ( ∞ ) ] the expectation value of A at the steady state ρ ( ∞ ) . From Equation 6 we obtain Equations 8–10 for the steady-state current: 〈 I 1 〉 ∞ = Γ f − 〈 σ 1 z 〉 ∞ , (8) 〈 I j − 1 〉 ∞ = 〈 I j 〉 ∞ , j = 2,3 , … , L − 1 (9) 〈 I L − 1 〉 ∞ = − Γ f + 〈 σ L z 〉 ∞ . (10) This current is uniform across the chain, with its value dictated by the average magnetization at the boundaries; from now on, we remove the subscript and denote it simply by ⟨ I ⟩ ∞ . In the absence of any bias ( f = 0 ) , the boundary reservoirs induce an infinite temperature (maximally-mixed) steady state ρ ( ∞ ) = I / 2 L , in which the values of magnetization and current all vanish.

In order to study the transport properties of the model at the steady state, both numerically and analytically, we make use of a perturbative ansatz characterized by an expansion in terms of potential bias f as [31, 32, 36]. ρ = 1 2 L I + f H + B + O f 2 , (11) H = ∑ r = 1 L ∑ j = 1 L + 1 − r h j r H j r , (12) B = ∑ r = 2 L ∑ j = 1 L + 1 − r b j r B j r , (13) where h j ( r ) and b j ( r ) where in Equations 11–13 are expansion coefficients, and the corresponding operators are given by Equation 14 H j 1 = − σ j z , H j r = σ j x Z j + 1 r − 2 σ j + r − 1 x + σ j y Z j + 1 r − 2 σ j + r − 1 y , r ≥ 2 B j r = σ j x Z j + 1 r − 2 σ j + r − 1 y − σ j y Z j + 1 r − 2 σ j + r − 1 x , r ≥ 2 . (14)

These operators are thus strings of Pauli matrices starting at site j and having length r , where Z k + 1 ( r − 2 ) = σ k + 1 z σ k + 2 z   …   σ k + r − 3 z σ k + r − 2 z are strings of Pauli z matrices of length r − 2 . The defining characteristic of this ansatz is that all operators in the expansion are orthogonal according to the Hilbert-Schmidt inner product ⟨ A , B ⟩ ≡ T r [ A † B ] . This implies that the expectation values of H j ( r ) and B j ( r ) are exactly determined within first order in f . In particular, H j ( 1 ) and B j ( 2 ) are respectively proportional to the magnetization and current 〈 σ j z 〉 ∞ = − f h j 1 , 〈 I j 〉 ∞ = f   4 J b j 2 ℏ . (15)

Once we determine the first order expansion coefficients h j ( r ) and b j ( r ) , we have access to the exact current and magnetization profile, as dictated by Equation 15 [32, 36, 38, 39]. In order to obtain them, we insert Equation 11 into Equation 2 and set the left-hand side to zero as demanded at stationarity (see Refs. Žnidarič [31]; Žnidarič and Horvat [32] and Supplementary Material for more details). This gives rise to a set of equations for the coefficients. ℏ Γ 1 + h 1 1 − 4 J b 1 2 = 0 , (16) ℏ Γ 1 − h L 1 − 4 J b L − 1 2 = 0 , (17) b j 2 − b j − 1 2 = 0 , j = 2 , … , L − 1 , (18)

for r = 1 . Note that Equations 16–18 just express the fact that the current is uniform in the steady state, with its value dictated by the boundary reservoirs as discussed before. For r ≥ 2 we have. J h j r − 1 − h j + 1 r − 1 + h j r + 1 − h j − 1 r + 1 + ε j − ε j + r − 1 h j r + ℏ γ + Γ 4 δ j , 1 + δ j + r − 1 , L b j r = 0 , (19) J b j r − 1 − b j + 1 r − 1 + b j r + 1 − b j − 1 r + 1 + ε j − ε j + r − 1 b j r − ℏ γ + Γ 4 δ j , 1 + δ j + r − 1 , L h j r = 0 . (20)

In total, Equations 16–20 are a closed set of L 2 coupled equations and will be the starting point of our study. This set of equations is equivalent to a Lyapunov-type equation for the correlation matrix in the presence of dephasing [32]. Instead of the exponential scaling with the system size, the quadratic scaling makes the problem computationally tractable. Moreover, as we show below, it also allow us to find analytical expressions for the current under some mild assumptions.

Having shown the degrading effect of the tilt and dephasing separately on the spin transport across the system, we now study their interplay which gives rise to dephasing-enhanced transport. In this section, we focus on transport for small system sizes L = 2,3 . An exact and compact solution for the current can be found from Equations 16–20 in these cases. For L = 2 we obtain ⟨ I ⟩ ∞ = f   16 J 2 Γ 2 γ + Γ 16 J 2 2 γ + Γ + 4 Γ U 2 + ℏ 2 Γ 2 γ + Γ 2 . (22)

Dephasing-assisted transport can be characterised through the maximum of Equation 22 with respect to the dephasing rate γ , yielding an optimal dephasing rate γ max = U ℏ − Γ 2 ≥ 0 . (23)

The behaviour of Equation 22 is shown in Figure 4. Note that since the dephasing is always positive, the maximum only occurs if U / ℏ ≥ Γ / 2 .

For L = 3 we have a similar expression, namely, ⟨ I ⟩ ∞ = f   16 J 2 Γ 4 γ + Γ 16 J 2 4 γ + Γ + 16 Γ U 2 + ℏ 2 Γ 4 γ + Γ 2 , (24) where the maximum is given by Equation 25 γ max = U ℏ − Γ 4 ≥ 0 . (25)

The complexity of these expressions increases rapidly for L > 3 , so the exact evaluation of the current and its maximum is too involved to be displayed here for larger systems. In the following section, we resort to an approximate solution for the current.

In order to find an exact solution for systems of arbitrary size, we perform the tridiagonal approximation which consists of retaining strings of Pauli matrices up to r = 2 , that is, ignoring higher order coefficients in Equations 16–20 so that h j ( r ) , b j ( r ) = 0 for r ≥ 3 . In the correlation matrix approach, this is equivalent to retaining only the diagonal and first off-diagonal elements [32]. Physically, this implies that spin and current correlations can be ignored within first order in f , which allows the derivation of an exact expression for the current as a function of all parameters as we now show. We only have to consider Equations 19, 20 for r = 2 . Summing Equation 19 from j = 1 to j = L − 1 we obtain J h 1 1 − h L 1 + ∑ j = 1 L − 1 ε j − ε j + 1 h j 2 + ℏ b γ L − 1 + Γ / 4 = 0 , (26) where b ≡ b j ( 2 ) is the coefficient associated with the current, which is independent of j as a consequence of Equation 18. First, we can write h 1 ( 1 ) − h L ( 1 ) as a function of b by using Equations 16, 17, which yields h 1 1 − h L 1 = 8 J b ℏ Γ − 2 . (27)

Second, we can also express h j ( 2 ) as a function of b by using Equation 20, which leads to Equation 28 h j 2 = ε j − ε j + 1 b ℏ γ + δ j , 1 + δ j , L − 1 Γ / 4 . (28)

We can now substitute the previous two equations into Equation 26 and obtain an expression for b , namely, b = 4 J 2 ℏ γ L − 1 + ℏ Γ + 16 J 2 ℏ Γ + 2 ℏ ∑ j = 1 L − 1 ε j − ε j + 1 2 γ + δ j , 1 + δ j , L − 1 Γ / 4 . (29)

The current is now obtained from Equations 15 and 29. For the local tilted potential in Equation 21 it becomes ⟨ I ⟩ ∞ = f   16 J 2 Γ 16 J 2 + ℏ 2 Γ Γ + 2 γ L − 1 + 4 Γ U 2 1 / 2 γ + Γ L = 2 L − 3 / 2 γ + 4 / 4 γ + Γ L ≥ 3 . (30)

We note that this expression is exact for L = 2 and L = 3 , in which case we recover Equations 22, 24, respectively. It is also exact in the absence of tilt U = V = 0 [31, 32, 37], which is verified in Figure 5A. In Figures 5B–D we plot the exact solution and the tridiagonal approximation for non-zero values of total tilt V > 0 . We observe that the approximation works remarkably well, deviating from the exact solution only for very small dephasing rates and system sizes; when γ = 0 , the approximation fails dramatically and only the exact solution is shown. The latter case corresponds to a Wannier-Stark localized system which we already discussed. Moreover, we see that for fixed γ the current can now increase with the system size up to a critical size L max where it is maximal, and only then decreases. In the same way, we see that for fixed L there is a critical dephasing rate γ max at which the current is maximal.

Since the validity of Equation 30 is established, we now use it to determine both L max and γ max . The maxima of Equation 30 as a function of L (for fixed V and γ ) and γ (for fixed L and V ) can be determined exactly. However, as before, the expressions are too lengthy to be displayed. More physically intuitive expressions are obtained by assuming that L ≫ 1 and L ≫ 8 γ / ( 4 γ + Γ ) in the denominator of Equation 30 which leads to the approximate maxima. L max ≃ V ℏ γ , (31) γ max ≃ V ℏ L = U ℏ . (32)

These maxima characterize dephasing-assisted transport as we show in Figure 6A. Below the critical length ( L max ≫ L ≫ 1 ) , dephasing is not yet effective and an increase in L decreases the local tilt U = V / ( L − 1 ) and thus localization, which increases the current according to Equation 30. Above this length scale ( L ≫ L max ≫ 1 ) dephasing becomes effective and diffusive transport clearly emerges ⟨ I ⟩ ∞ ∼ 1 / ( γ L ) . Analogously, below the critical dephasing rate ( γ max ≫ γ ) dephasing jump events happen at a much lower frequency than Bloch oscillations, so the current flowing remains hampered by localization. Increasing the rate of dephasing destroys localization and increases the current. When the dephasing rate passes this threshold ( γ ≫ γ max ) diffusive behaviour emerges. We illustrate this in Figure 6B where we plot the (finite size) diffusion coefficient defined by ⟨ I ⟩ ∞ ≡ f D / L as a function of dephasing rate for the largest system size we reached L = 1000 .

The imprint of dephasing-enhanced transport is also evident in the magnetization profile, as shown in Figure 7. Here, varying the dephasing strength towards its critical value Equation 32 decreases the magnetization gradient, rendering the magnetization profile almost uniform. This is in contrast with the results obtained in the absence of tilt, where an increasing dephasing rate tends to increase the magnetization gradient and decrease uniformity as shown in Figure 3A. Uniformity as a feature of environment-assisted transport has been pointed out for quantum systems subject to on-site disorder [9]; our results confirm that the same phenomenon appears in the presence of a linear potential inducing Wannier-Stark localization.