Dephasing-induced diffusive transport in anisotropic Heisenberg model

We study transport properties of anisotropic Heisenberg model in a disordered magnetic field experiencing dephasing due to external degrees of freedom. In the absence of dephasing the model can display, depending on parameter values, the whole range of possible transport regimes: ideal ballistic conduction, diffusive, or ideal insulating behavior. We show that the presence of dephasing induces normal diffusive transport in a wide range of parameters. We also analyze the dependence of spin conductivity on the dephasing strength. In addition, by analyzing the decay of spin-spin correlation function we discover a presence of long-range order for finite chain sizes. All our results for a one-dimensional spin chain at infinite temperature can be equivalently rephrased for strongly-interacting disordered spinless fermions.


Introduction
Theory of free fermions has been very successful in explaining properties of many condensed matter systems. Sometimes though, in particular in low dimensional systems, the picture of isolated noninteracting fermions is an oversimplifications and one has to include additional interactions. Integrable model of free fermions can be upgraded by taking into account several effects: (i) interaction -fermions can interact with each other, (ii) disorder -fermions can experience different local eigenenergies and, (iii) coupling -fermions can interact with external degrees of freedom. A non-perturbative inclusion of any of these effects greatly complicates the analysis forcing one to use various successful phenomenological theories [1]. There is nevertheless a desire to understand properties of such strongly-interacting many-body disordered systems from the first principles, that is starting from the governing equations of motion.
In the present work we shall study transport properties of a one-dimensional (1D) antiferromagnetic Heisenberg model at infinite temperature in the presence of all the above mentioned effects: interaction, disorder and coupling to the environment. To our knowledge this is a first study of a large many-body system in such setting by using microscopic equations of motion. In terms of Pauli matrices the 1D antiferromagnetic spin-1/2 Heisenberg model reads Magnetic field strength h j at site j is chosen according to a uniform distribution in the interval [−h, h]. Even though a clean Heisenberg model without disordered magnetic field is explicitly solvable by the Bethe ansatz and has been studied for several decades, its spin transport properties are far from being understood, for an overview see [2]. By using the Jordan-Wigner transformation it can be mapped to a system of strongly interacting spinless fermions, with ∆ being the interaction strength. Therefore, the problem of spin transport in the Heisenberg model is equivalent to that of particle (charge) transport in a system of interacting fermions. In the absence of coupling to the environment, rigorous statements are available only for special cases without anisotropy (interaction), ∆ = 0, or in the case of no disorder, h = 0, see Fig. 1. For ∆ = 0 the model is equivalent to a system of non-interacting fermions and can be described in terms of single-particle states that exhibit exponential localization in the presence of nonzero disorder, h = 0, resulting in a perfect insulating behavior due to Anderson localization [3]. If disorder is zero the Heisenberg model displays ballistic spin transport in a gapless phase for ∆ < 1 [4], while it is likely a diffusive spin conductor for ∆ > 1 [5,6,7,8,9]. This last result about diffusive transport in an integrable model, based mainly on numerical investigations, is somewhat unexpected and still lacks a deeper understanding, see though [10,11] for theoretical arguments supporting diffusive transport. As soon as interaction ∆ and disorder h are both nonzero things get more complicated because one has to deal with a genuine many-body physics, for single-particle treatment see, e.g., [12,13]. There are in fact opposing claims in the literature on whether the localization is destroyed or preserved in the presence of strong interaction; some predict a phase transition from an insulating to a conducting regime at high temperature [14], other numerical simulations, albeit on small systems, suggest localization even at an infinite temperature [15,16] or even ideal conducting behavior [17].
Other non-equilibrium properties, like dynamics in various quantum quenches, has also been studied intensely in recent years, see [18] and references therein. For a recent study in 3D see [19]. Besides its theoretical importance, Heisenberg model is also experimentally realized in spin-chain materials [20]. It is argued that an anomalously high spin conductivity of the Heisenberg model is required to explain large measured heat conductivity. While a detailed engineering of interactions in condensed-matter systems is probably out of question, experiments with cold gases in optical lattices could in near future lead to controlled experiments with strongly correlated many-body systems [21]. It is therefore important to get a better understanding of the influence of the interplay between interaction, disorder and external coupling on the dynamics in such systems.

The model
We are going to study transport by numerically solving dynamical equations of motion, finding a non-equilibrium stationary state (NESS) to which the system converges after long time. Once we calculate NESS we can evaluate various expectation values in that stationary non-equilibrium state. The dynamics of the system will be described by the Lindblad master equation [22], Bath superoperator L bath takes into account the coupling to unequal heat baths at both ends, inducing a non-equilibrium situation, while L deph represents the effect of dephasing caused by the coupling to some external degrees of freedom. Even though description in terms of Lindblad equation, originating in quantum optics, has not been employed very often in condensed matter physics, it has been recently used in several works studying spin transport in the Heisenberg model [23,24,7,8,25]. While a number of assumptions is involved in the derivation of the Lindblad master equation from unitary evolution of a combined system [26], Markovian approximation can be justified because we are interested in the stationary state reached after long time, e.g., much longer than the bath correlation time. Apart from coupling to the bath the only non-unitary effect we take into account is dephasing. At high temperatures, which is the regime studied in the present work, dephasing is the dominant environmental effect. We are going to solve master equation using time-dependent density matrix renormalization group (tDMRG) method enabling us to study by an order of magnitude larger chains [7,8] than is possible with other methods. Dephasing present in our simulations in fact acts stabilizing on the numerical method enabling us to simulate relaxation in chains of up-to 512 spins. Let us describe various terms in the Lindblad equation (2). The Lindblad bath superoperator L bath affects only two border spins at the left and right end of the chain and involves four Lindblad operators L L k=1,2,3,4 at the left end and four L R k=1,2,3,4 at the right end, Lindblad operators L L,R k are chosen in such way that if they would be the only term present in the master equation, the NESS, i.e., the state ρ(t → ∞), would be equal to the grand canonical state on two border spins. They are obtained by targeting the grand canonical state at an infinite temperature, ρ grand. ∼ exp (−µ L,R Σ z ), Σ z = k σ z k , on two boundary spins, see [7,27] for more details. By choosing different potential µ L and µ R at the left and the right end, a NESS with a nonzero spin current is obtained. The bath operators are such that for µ L = µ R they induce the correct equilibrium grand canonical state in the bulk of an interacting chaotic spin chain [27]. Dephasing part L deph , due to, for instance, coupling with phonons, involves only one Lindblad operator where the sum over j runs over spin sites, σ ± = σ x ± i σ y , and γ is a dephasing strength.
Dephasing (also called phase damping or phase-flip in quantum information) causes an exponential decay of the off-diagonal elements in the diagonal basis of σ z (in fermionic language this corresponds to the off-diagonal decay in the number basis). It might be interesting to note [28] that the dynamics of operators in the Heisenberg model in the presence of dephasing, i.e., evolution by the Lindblad equation (2), can be exactly mapped to the evolution of pure states in a one-dimensional non-hermitian Hubbard model with complex on-site repulsion (iU)n j↑ n j↓ . ‡ The same L deph is obtained for Recently, the influence of dephasing on transport has been studied in short disordered networks modelling the light-harvesting complex [29,30,31,32], where the presence of dephasing can enhance network's ability to transmit excitations. First we are going to show expectation values of some simple observables in NESS obtained as an infinite-time solution of the Lindblad equation (2). We choose bath potentials on left and right ends as µ L,R = ±0.1 (same throughout the paper) and numerically solve the Lindblad equation, starting with some arbitrary initial condition. Details of our numerical implementation of tDMRG can be found in the Appendix. After sufficiently long time the state ρ(t) converges to a time-independent NESS. This "relaxation" time after which we converge to the NESS depends on the bath details, i.e., it is a combined property of the central system and of the coupling. It increases with the increasing chain size n. For each run we carefully check that the simulation time has been long enough so that the convergence to the NESS is indeed reached. We show in Fig. 3 the expectation value of local magnetization σ z k and of local spin current j k = 2(σ x k σ y k+1 − σ y k σ x k+1 ) in a NESS of a disordered Heisenberg model with n = 128 spins. One obtains a linear spin (magnetization) profile, as it is expected for a diffusive conductor, as well as a homogeneous spin current j k throughout the chain. Our numerical experience tells us that the homogeneity of the spin current can be used as a rather sensitive indicator of the convergence to NESS as well as of truncation errors.

Spin conductivity
To determine whether the system is a normal or an anomalous conductor one has to study the scaling of the spin current with the system size. First, we are going to study clean system without any disorder, h = 0. Dephasing strength is always γ = 1 and three different anisotropies ∆ = 0, 0.5, and ∆ = 1.5 are used. Keeping chemical potentials fixed we have calculated the spin current in the NESS for systems between n = 8 and n = 256 spins long. Results are shown in Fig. 4. We can see that at γ = 1, regardless of the value of the anisotropy, the current always scales as ∼ 1/n with the system size, meaning that the system is diffusive (i.e., normal) conductor. The asymptotic value of the product j · n determines the coefficient of spin conductivity κ through j = −κ( σ z 1 − σ z n )/n. For large n we have σ z 1 − σ z n ≈ 0.153 giving κ = 3.79, 3.20, 1.53 for γ = 1 and ∆ = 0, 0.5, 1.5, respectively. With the increasing interaction ∆ the conductivity decreases. For ∆ = 1.5 and without dephasing and disorder, γ = 0 and h = 0, the coefficient of spin conductivity has been calculated in Refs. [7,9] to be κ = 2.30 and is therefore larger than with dephasing.  Figure 4. Scaling of the spin current with system's size in a clean system, h = 0 and γ = 1, at a fixed driving potential. For all anisotropies the presence of dephasing induces normal (diffusive) transport, signaled by the scaling j ∼ 1/n (straight lines in the main plot). In the inset the convergence of the product j · n to its asymptotic value is shown.
Similar results are obtained also in the presence of disorder. In Fig. 5 we show the scaling of the spin current with n when each site, in addition to the dephasing with γ = 1, experiences also a random magnetic field h j with strength h = 2 §. The difference in § Because fluctuations between different realizations of disorder were found to be small for large n all magnetizations between left and right ends is again almost independent of ∆ and equal to σ z 1 − σ z n ≈ 0.153, resulting in κ = 1.63, 1.50, 1.00 for ∆ = 0, 0.5, 1.5, respectively. Compared to the clean case disorder decreases conductivities, however, transport is still diffusive. Dephasing therefore breaks any possible many-body localization present at γ = 0. In the inset of Fig. 5 we also show data for the case of single-particle Anderson localization (γ = 0, ∆ = 0) for which the current exponentially decreases with n. Dephasing therefore induces normal transport irrespective of the anisotropy and disorder strength. Decreasing dephasing strength to 0 one of course has to recover known behavior for γ = 0 which is, depending on ∆ and h (see also Fig. 1), superconducting (for ∆ < 1 and h = 0), diffusive (for ∆ > 1 and h = 0) or insulating (for ∆ = 0 and h = 0). To elucidate this transition we have studied how the coefficient of spin conductivity κ changes as the dephasing strength γ is varied. For each value of parameters ∆ and h we calculated NESS for sizes n = 64, 128 and n = 256, from which we determined κ. Spin current in all cases scaled as ∼ 1/n, the transition to asymptotics though happening at larger n for smaller values of γ. Dependence of κ in the absence of disorder is shown in Fig. 6. One can see that for ∆ = 0.5 spin conductivity increases with decreasing γ. Divergence is scaling as κ ∼ 1/γ which is consistent with the superconducting (i.e., ballistic) limit at γ = 0. For ∆ = 1.5 though, κ converges to a finite value as γ → 0, in accordance with the normal transport in the Heisenberg model at ∆ > 1 without dephasing. data shown are for a single realization of disorder. In the presence of disorder the limit γ → 0 corresponds to Anderson localization if ∆ = 0 while for nonzero ∆ things are not so clear. From our results seen in Fig. 7 we see that for small γ spin conductivity decreases, however, the limit value κ(γ = 0) is rather difficult to predict. Spin conductivity reaches a maximum at some intermediate value of γ and again decreases for large γ. Note that there is no qualitative difference between ∆ = 0.5 and ∆ = 1.5. A nontrivial maximum of conductivity is found also when studying dephasing in photosynthetic complex [29,30,31,32].

Correlation function
We have also calculated the spin-spin correlation function for non-equilibrium steady states from Figs. 4 and 5. Results are similar for all values of ∆ and we show only the case with ∆ = 0.0, for which tDMRG is the most efficient. Spin correlation function Figs.8 and 9, in all cases has a plateau. This indicates the presence of a long-range order. Scaling analysis shows that the plateau in C(i, j) scales with the system size as ∼ 1/n and with the potential difference as ∼ (∆µ) 2 = (µ L − µ R ) 2 . The same scaling is obtained also for ∆ = 0 (data not shown), however, because of less favorable scaling of the entanglement with n we could reliably calculate C(i, j) only for sizes n ∼ 128. Note that decreasing values of the correlation function for increasing size make a precise calculation of C(i, j) rather demanding. Scaling of the correlation function plateau C ∼ (µ L − µ R ) 2 /n ∼ j · (µ L − µ R ) means that the long-range order is a purely non-equilibrium phenomenon. In equilibrium, when the spin current is zero, j = 0, it goes to zero. In addition, correlations go to zero also in the thermodynamic limit n → ∞, as one would expect for a true normal (diffusive) conductor. Similar scaling of the correlation plateau has been recently observed in a non-equilibrium phase transition in the XY spin chain in a homogeneous magnetic field, studied by exactly solving the corresponding master equation [33,34].

Conclusion
We have studied spin transport in the spin-1/2 Heisenberg model at infinite temperature in the presence of disorder and dephasing. By solving master equation for nonequilibrium driving potential we numerically obtained a non-equilibrium steady state and calculated expectation values of different observables. Studying the scaling of the spin current with system's size for systems of up-to 256 spins we found that the dephasing always induces diffusive transport regardless of the anisotropy or disorder strength. The size of spin conductivity coefficient has a non-trivial maximum at intermediate values of the dephasing strength for a disordered model while it monotonically decreases with dephasing in the absence of disorder. Studying the spin-spin correlation function in a non-equilibrium steady state we have found the evidence for a long-range order at all anisotropies. The value of the correlation plateau scales as ∼ (µ L − µ R ) 2 /n, meaning that the correlations disappear in the thermodynamic limit as well as in the equilibrium situation. It would be interesting to study more in detail how the transition from diffusive transport to either superconducting, diffusive or insulating behavior happens as the dephasing strength is decreased to zero. More studies would be also needed to clarify the situation with a non-zero anisotropy and without dephasing, which corresponds to strongly interacting fermions.
Any state ρ from a Hilbert space H which is a tensor product of n local spaces of dimension d, H = H ⊗n loc. , dim(H loc. ) = d, can be written in a matrix product operator (MPO) form, where expansion coefficients are expressed in terms of product of matrices. If we denote basis states of H loc. at j-th site by σ ν j , ν = 0, . . . , d − 1, we can write |ρ = ν i x † M ν 1 1 M ν 2 2 · · · M νn n x|σ ν 1 1 σ ν 2 2 · · · σ νn n . solving the Lindblad equation (2), |ρ represents a density operator which is viewed as a member of a Hilbert space of operators, with local basis being Pauli matrices σ 0 = σ x , σ 1 = σ y , σ 2 = σ z , σ 3 = ½. For a given |ρ the choice of matrices M ν j j is not unique, however, as we will see, there is a preferred choice corresponding to Schmidt decomposition of |ρ .

Appendix A.2. Evolution
The above setting is used to solve master equationρ = Lρ, where L is some linear operator. In our simulation L is just the operator corresponding to the right side of the Lindblad equation (2). Solution can be formally written as ρ(t) = exp (L∆t) ρ(0) in terms of propagators exp (L∆t) for a small time step ∆t; in our simulations we use ∆t = 0.05. Each small-step propagator is then decomposed according to the Trotter-Suzuki formula. Writing L = A + B as a sum of two parts, each of which is a sum of mutually commuting terms, we have We use a fourth order (k = 4) decomposition [35] with . In our Lindblad equation we have only nearest neighbor unitary terms and at most two-qubit nearest neighbor dissipative terms. We put into A all even bonds, into B all odd bonds, while single qubits terms are equally split between the two terms. A and B are therefore sums of commuting two-qubit terms. Basic operation we have to do is then exp (A j,j+1 ǫ), where ǫ ∝ ∆t. For time-independent A j,j+1 this is a fixed linear operator on the Hilbert space of operators. Its matrix representation A j,j+1 in the basis |σ ν j j can be calculated in advance, A j,j+1 is orthogonal for unitary exp (A j,j+1 ǫ). Every bipartite state can be written in the form of the Schmidt decomposition. Splitting a system into the first j sites (1, . . . , j) and the rest, we can write with an explicit form of orthogonal states |w L(j) α and |w R(j) j+1 · · · M νn n ] α,β x β |σ ν j+1 j+1 · · · σ νn n . (A.5) In terms of matrices M ν j j the orthogonality of |w L(j) α is reflected in the condition for all j, while for |w R(j) must hold. Coefficients λ (j) α are called the Schmidt coefficients. How well the tDMRG method works crucially depends on how fast the ordered λ (j) α decrease with α. In general time evolution can not be done in an exact way unless we increase the dimension of matrices M ν j j exponentially with time. To keep dimensions fixed we have to make a truncation. This is done in an optimal way (in terms of probabilities (λ (j) α ) 2 ) if we make sure that we drop states corresponding to the smallest Schmidt coefficients. To do this we first calculate the action of A j,j+1 on the matrix product operator (MPO) form (A.1); it transforms a product of two matrices M ν j j and M ν j+1 j+1 at sites j and j + 1 into a single bigger matrix Θ, The reason to multiply it with λ (j−1) α is to recover Schmidt decomposition after doing a singular value decomposition on Θ, decomposing it as Θ = UdV , with unitary U and V and diagonal d. It turns out that if we prescribe new matricesM and new coefficients λ (j) γ asλ One Trotter-Suzuki step is optimal if all transformations preserve Schmidt decomposition, i.e., if they are unitary. Coupling to the bath and dephasing obviously violate this condition. Because the bath part L bath destroys unitarity only at the two boundaries we included it in A and B terms. Dephasing part though, acting on each spin, would affect unitarity also in the bulk of the chain. As a matter of numerical convenience we have not included it in the Trotter-Suzuki decomposed part (A.2) but have instead applied it only at the end of each Trotter-Suzuki step as exp (L deph ∆t). This corresponds to dephasing acting in a kicked manner and not in a continuous way and has no physical consequences for our results. After each step of length ∆t, we are left with a MPO |ρ whose matrices do not correspond to Schmidt decomposition anymore (due to non-unitary propagators in between). Before proceeding to the next step we reorthogonalize matrices M ν j j in order to recover the optimality of the Schmidt decomposition. This takes of order O(ndD 3 ) steps and is done in the following way [37,38]. First, we recursively construct matrices V L j and V R j as starting with the matrix [V L 0 ] k,l = x k x * l and [V R n ] k,l = x k x * l at the boundaries. If conditions (A.6,A.7) are satisfied the matrices V L j V R j are diagonal. In general this is not the case so we have to rotate each M ν j j in order to recover orthogonality. To do this we calculate a square root of a non-negative V L j and diagonalize the matrix