Block-Lanczos density-matrix renormalization-group approach to spin transport in Heisenberg chains coupled to leads

We adapt the block-Lanczos density-matrix renormalization-group technique to study the spin transport in a spin chain coupled to two non-interacting fermionic leads. As an example, we consider leads described by two-dimensional tight-binding models on a square lattice. Although the simulations are carried out using a chain representation of the leads, observables in the original two-dimensional lattice can be calculated by reversing the block-Lanczos transformation. This is demonstrated for leads with Rashba spin-orbit coupling.


Introduction
Magnetic insulators are potentially useful for spintronics applications because of the greater decay length of spin currents compared to conductors. In this light, there have been several experiments on the spin-current generation in insulating materials. [1][2][3][4] For example, Kajiwara et al. achieved the injection of spin currents into the ferrimagnetic insulator yttrium iron garnet (YIG) from an attached platinum (Pt) film, where a spin current was generated electrically using the spin Hall effect. 1) Spin currents have also been experimentally investigated for other types of magnets, including various antiferromagnets, 5) the paramagnet Gd 3 Ga 5 O 12 6) and the spin nematic liquid LiCuVO 4 . 7) Relevant to the present study, a spin current was induced in the spin-1/2-chain material Sr 2 CuO 3 by using the longitudinal spin Seebeck effect 8) at an interface with a Pt film. 4) In the same work, the spin transport for this setup has been studied theoretically using a Green's function technique. 9) See also Refs. 7, 10 for similar analyses of spin transport in other magnetic systems.
Motivated by these developments, we numerically study the spin transport in a one-dimensional (1D) antiferromagnetic spin-1/2 chain coupled to two non-interacting fermionic leads. In numerical simulations, the leads are typically represented by a finite number of sites which causes a discretization of the Hamiltonian in energy space. To reach a larger number of sites and a finer discretization of energy, it is advantageous to map the non-interacting leads to a chain representation that permits the use of matrix-product state (MPS) techniques, e.g., the density-matrix renormalizationgroup (DMRG) or the numerical renormalization group method. 11,12) There are different ways to obtain such a chain mapping. Here, we use a Lanczos technique, 13) which is convenient when starting from a model defined on a real-space lattice.
Previously, we investigated the spin transport for a junction in which the leads are modelled by tight-binding chains with uniform hopping parameters. 14,15) The present study differs in that we start from a higher-dimensional tight-binding model and obtain the chain representation numerically. This serves two purposes: (i) we can evaluate how important the assumption of homogeneous 1D leads is for the conclusions of Refs. 14, 15, (ii) we outline a general strategy for the numerical study of spin transport in 1D systems coupled to non-interacting leads of arbitrary dimension. Here, we specifically consider a junction with two-dimensional (2D) leads defined on a square lattice. Within the 1D representation of this model, we calculate the spin conductance with the DMRG and the time-evolving block decimation methods, 16) finding qualitatively similar results to the setup with uniform 1D leads. Observables in the original 2D lattice can also be calculated by reversing the block-Lanczos transformation. We demonstrate this for leads with Rashba spin-orbit coupling, which is important for semiconductor heterostructures. 17) Note that although we focus on a specific type of spin chain and 2D leads, the presented numerical approach could be applied to other models, e.g., different types of spin chains or other lattices for the leads.
The rest of this paper is organized as follows. In Sec. 2, we introduce the model. Section 3 describes the block-Lanczos transformation applied to facilitate the use of MPS methods. The results of the MPS calculations are presented in Sec. 4. Section 5 contains a summary and a discussion of the results. The green spheres represent spin-chain sites, the blue ones lead sites.

Model
We consider a junction of a 1D interacting region and two non-interacting leads (see Fig. 1). The interacting part is modeled by a spin-1/2 Heisenberg chain of length N S , with antiferromagnetic exchange coupling (J > 0). For the lead HamiltoniansĤ L 1 andĤ L 2 , we assume 2D tight-binding models, possibly with additional Rashba spin-orbit coupling, i.e.,Ĥ whereĉ ( †) jσ is the annihilation (creation) operator of an electron with spin σ ∈ {↑, ↓} at site j [located at (x j , y j )] on an infinite square lattice with an open edge and lattice constant a = 1, and i, j indicates all pairs of nearest-neighbor sites i and j. The spin-orbit coupling is parametrized by λ, σ x and σ y are Pauli matrices, and µ is the chemical potential. At each end, the spin chain is coupled via an exchange interaction to a site at the open edge of one of the leads. Denoting the indices of these sites by j 0 , the coupling terms are of the form where an additional index a ∈ {L 1 , L 2 } was introduced to distinguish between the two leads. With Eq. (3), the complete Hamiltonian takes the form

Block-Lanczos Transformation
The Lanczos algorithm is a way to obtain a unitary transformation that tridiagonalizes a given Hermitian matrix H 0 . 18) One starts with a single unit vector v 1 which is the first column of the transformation matrix P = (v 1 , v 2 , ...). All remaining v j are then obtained by setting v j ← H 0 v j−1 ( j = 2, 3, 4, . . . ) and orthogonalizing against previous vectors. We use a block version of the Lanczos method, 19) in which one chooses the first M orthonormal vectors v 1 , v 2 , . . . , v M and sets v n ← H 0 v n−M , again followed by an orthogonalization. From this construction and the Hermicity of H 0 , it follows that P † H 0 P is a band matrix with bandwidth 2M + 1: where E j and T j are Hermitian and lower-triangular M × M matrices, respectively. If the block-Lanczos transformation is applied to the matrix H 0 describing a single-particle HamiltonianĤ 0 of a sys- is a vector of creation operators (here the spin index is implicitly assumed), the banded structure ofH 0 means that, in terms of the new operatorsâ † = (â † 1 ,â † 2 , ...,â † N ) =ĉ † P,Ĥ 0 describes an open chain with short-ranged hopping and a site-dependent potential. Explicitly stated: where N is the order of the matrix H 0 , corresponding to the number of single-particle states. The first M operatorŝ a † 1 ,â † 2 , ...,â † M can be fixed through the choice of initial vectors. Increasing M thus allows for a greater flexibility in the transformation but it also increases the maximum hopping range.
For an interacting HamiltonianĤ =Ĥ 0 +V, in which the interaction is a function of only M electron operators, i.e., , we can use a block-Lanczos transformation with block size M and setâ † n =ĉ † j n for n = 1, 2, ..., M. The interaction is then restricted to the first M sites after the transformation and the chain representation of H contains only short-ranged terms. An exact solution of the problem will still be impossible in general but the 1D form makes the model suitable for a numerical treatment with MPS techniques. The Lanczos method has been used in this manner mostly in the context of impurity problems. 13,20) Our application here to the spin-chain junction is similar to it, with the sites coupled to the spin chain taking the role of the impurities.

Infinite boundary conditions
In the following discussion, we assume that the Hamiltonian is originally defined on a lattice in real space and has only short-ranged hopping terms. The interactionV shall act on one site whose corresponding fermion operators (with internal degrees of freedom, such as spin, and thus M being larger than 1) will be invariant under the transformation. While we are ultimately interested in the thermodynamic limit, the matrix H 0 and thus the system size need to be finite in a numerical calculation. However, for a given j max one can always choose the original system large enough so that the transformed operatorsâ † j with j ≤ j max are not affected by its finite size. This follows from for j ≤ rM, which implies thatâ † j is supported only on sites connected to the interacting site through at most r − 1 hopping operations. Of course, finite-size effects will eventually appear if the transformation is carried out to completion. We can, however, stop the Lanczos recursion before that happens and work with a truncated transformation by ignoring the remaining sites in the chain representation. For a fixed number of chain sites these infinite boundary conditions will be closer to the thermodynamic limit with regard to the physics at the interacting site. 20) On the downside of this procedure, P is no longer unitary which complicates the measurements for the original lattice. The one-body expectation values in both representations are for the full transformation. Typical quantities in the original system thus have to be reconstructed from a large number of correlation functions in the effective 1D model. If the Lanczos recursion is stopped prematurely, the above relation is not fulfilled because states are missing on the right-hand side. However, it is still possible to calculate with reasonable accuracy the change in the expectation values that is induced by a perturbation at the interacting site, e.g., an injected current. The reason is that, for sufficiently many sites in the effective 1D model, the only missing terms in Eq. (8) are then between sites that are both outside the range of the perturbation, or whose distance so large that the contribution from the correlation function can be neglected. We thus expect the following relation to hold during the time evolution: where τ denotes the time at which the expectation value is calculated and indicates that the sum over m and n is truncated.
Note that another more accurate way to calculate expectation values away from the interacting site is to increase the block size M and include additional sites in the initial block-Lanczos basis set for measurements. This has been used in Ref. 13 to calculate two-point correlation functions between impurity and conduction sites. For our problem this scheme is less suitable, since we intend to carry out measurements at a lot of different positions, which would require the simulations to be repeated many times.

Application to leads
We now apply the block-Lanczos transformation to a 2D tight-binding lead with Rashba spin-orbit coupling described in Sec. 2. To keep the exchange interaction with the spin chain local, both up and down spin states of the affected site on the lead are included in the initial block-Lanczos basis set, and the block size is therefore M = 2. Different transformations and the resulting effective 1D models are obtained, depending on the position of the interacting site. Here, the site is assumed to lie on the open edge of the lead (see Fig. 2).
It turns out thatĤ L 1 (L 2 ) expressed in terms of the new fermion operatorsâ † j describes two decoupled tight-binding chains with the same bond-dependent hopping parameterst j . Definingâ 2 j−1 →â j↑ andâ 2 j →â j↓ allows us to writê Therefore, even for the finite Rashba spin-orbit interaction the only difference to a regular tight-binding chain is the positiondependent hoppingt j . The chemical potential term remains the same under the block-Lanczos transformation but the density changes in general unless the system is at half-filling. The conservation of the new spin introduced in Eq. (10) can be exploited in numerical simulations. Note, however, that it corresponds to the physical spin only at the first site. The simple form of Eq. (10) for the Hamiltonian in the block-Lanczos basis implies that the Krylov subspaces generated by the |↑ =â † 1↑ |0 single-particle state at the interacting site are orthogonal to those generated by the |↓ =â † 1↓ |0 state. To prove this, one needs to show that ↑|Ĥ n L 1 (L 2 ) |↓ = 0 for all natural numbers n. This follows from the time-reversal symmetry of the Hamiltonian: where |↓ =T |↓ = −|↑ , |↑ =T |↑ = |↓ , andT is the time-reversal operator, assuming that t is real in Eq. (2). The Lanczos method with block size M = 1 therefore would have been sufficient to obtain the transformation. This is also true if there is both Rashba and Dresselhaus spin-orbit coupling, but not in the presence of a magnetic field. Figure 3 shows the position dependence of the hopping amplitudes in Eq. (10) for the infinite boundary conditions described in the previous section. Already after a few sites ( j ≈ 8), the hopping amplitude approaches a constant value that depends on the spin-orbit interaction λ. This asymptotic value agrees with the hopping with k m = arctan(λ/ √ 2t), which leads to the same bandwidth as in the original Hamiltonian in Eq. (2). Note that if we applied the transformation to a finite and smaller system, the hopping parameter would become non-uniform for sufficiently larger site indices.
From Eq. (10) and Fig. 3, one can see that, regardless of the spin-orbit coupling strength λ, the effective model is quite similar to a regular tight-binding chain with essentially uniform hopping amplitude except for the first few sites. The differences due to the Rashba spin-orbit coupling only become apparent when transforming back to the original representation.

Spin conductance
With the 2D Rashba system mapped to a chain representation, we are now in the position to examine the spin transport in the junction by means of MPS techniques. Because of the Rashba spin-orbit coupling in the leads, we could induce a spin current by exploiting the spin Hall effect and applying an electric field. However, to simulate an electric field in the Rashba system, we need to either switch on a static potential at the start of the time evolution or add a time-dependent phase factor to the hopping terms. In both cases, the perturbation would be highly non-local in the block-Lanczos basis, rendering MPS simulations inefficient. We therefore neglect the Rashba spin-orbit coupling in this section and assume that the spin current is driven by an effective spin-voltage. Namely, we add a potential termĤ V = (V/2) j∈L 1 (ĉ † j↑ĉ j↑ −ĉ † j↓ĉ j↓ ) to the first lead, which induces a spin current polarized in the z direction.
The spin conductance G S = I/V is defined as the ratio of the spin current I that flows through the junction in the nonequilibrium steady state, and the spin voltage V. For the steady-state spin current, one can write I = lim τ→∞ ĵ z (τ) (1 ≤ < N S ), where τ is the time andĵ z = i(J/2)(Ŝ + Ŝ − +1 − S − Ŝ + +1 ) the spin-current operator between sites and + 1 of the spin chain. To obtain the zero-temperature spin conductance G S of the junction numerically, we first calculate the ground state with the DMRG method and then simulate the time evolution with switched-on spin voltage using the timeevolving block decimation algorithm. Because of the finite lead sizes, a true steady state is not reached in the simulations. However, it is nevertheless possible to estimate I accurately from the time-dependence of the local spin currents ĵ z (τ) .
More details on the numerical method are given in Refs. 14, 15, where a similar setup with uniform 1D leads was studied. In these works, it was found that the spin conductance G S depends sensitively on the model parameters near the interfaces.
Generally, G S is significantly reduced compared with the spin conductance G 0 S = 1/(4π) for a homogeneous tightbinding chain, i.e., a system without spin chain. Through finetuning of the parameters at the interfaces, however, a so-called conducting fixed point may be reached. 21,22) There, the zerotemperature spin conductance is not reduced by interface effects and takes the maximum value G 0 S determined by the leads. 23,24) Using the specific 2D leads (in their 1D representation) should not qualitatively affect this result but the site dependence of the hopping parameters near the interface could move the system away from or towards a conducting fixed point.
We have confirmed this for some values of the model parameters by explicit numerical calculations, as shown in (2) without spin-orbit coupling (λ = 0). In both cases, the overall energy scale is chosen so thatt ∞ = lim j→∞t j = J. Each lead is truncated to a finite length of 400 sites (without including spin degrees of freedom) in the MPS simulations. Fig. 4. At a small but finite spin voltage V S /J = 0.1, the spin conductance G S shows a sharp peak as a function of the interface-coupling strength J . This is observed for both types of leads but the position of the maximum is different. For the block-Lanczos leads, the peak is shifted to smaller J , as may be expected because of the reduced hopping strength near the spin chain. The maximum value of G S is in both cases approximately the ideal value G 0 S = 1/(4π) for the linear conductance. In Refs. 14, 15, only half-filled leads were considered. As shown in Fig. 4, conducting fixed points occur for a finite chemical potential µ as well, although their position (i.e., J /J) is µ dependent.
It should be noted that the coupling strength J c between spin chain and metallic leads corresponding to the conducting fixed point in Fig. 4 is larger than the exchange coupling J inside the spin chain. In a real experiment, on the other hand, J c may be much smaller than J, so that according to our simple model the system would be far away from the conducting fixed point. However, our calculations were done for the zero temperature limit in which the reduction of the conductance away from the conducting fixed point is particularly strong. 21,22) An interesting open question is, how the interface effects change quantitatively when a more realistic description is used that, e.g., also takes finite temperature into account.
Finally, let us briefly comment on the case of ferromagnetic exchange interaction J and J . Numerical calculations indicate that the spin currents are much smaller than for the antiferromagnetic spin-1/2 chains studied here. This may be explained by considering the Kondo model, i.e., a single spin coupled to a fermionic bath. It is known that for ferromagnetic interaction, the Kondo spin effectively decouples from the lead in the low temperature limit. 25) Adding spin sites that interact ferromagnetically with the first spin will not change this behavior. Accordingly, we expect the linear spin conductance of the junction to vanish at zero temperature in the case of ferromagnetic exchange interaction.

Lead dynamics
In the previous section, we only considered expectation values in the effective 1D model, which is sufficient to character- ize the spin transport through the spin chain. As described in Sec. 3.1, however, the block-Lanczos transformation can be reversed to obtain expectation values of observables defined on the original real-space lattice. To this end, we calculate the single-particle expectation values â † iσâ jσ in the 1D representation of the lead for all sites i, j ≤ L with some finite number L. This can be done by two nested sweeps in the MPS calculation so that the computational cost scales quadratically with L.
As an example, we investigate the magnetization m i = 1 2 ĉ † i↑ĉi↑ −ĉ † i↓ĉi↓ in the second lead after a spin current polarized in the z direction is injected. Figure 5 shows m i for 2D leads with and without spin-orbit coupling. The spin current entering the lead induces a small position-dependent magnetization that depends strongly on both the chemical potential µ and the strength of the spin-orbit coupling λ. Most strikingly, the Rashba precession for λ 0 leads to oscillations as a function of the position. As can be seen in Fig. 6, the timedependence of the local magnetization m i is relatively small after the wave front has passed, and the Rashba oscillation pattern in particular does not change.
Near half-filling at µ = 0, the induced magnetization is largest along the diagonal directions, while it becomes more uniformly spread when the chemical potential is decreased. This can be understood by looking at the shape of the Fermi surface for λ = 0. At half-filling, it takes a diamond form and thus the energy gradient ∇ k E(k)| |k|=k F , i.e., the group velocity, points in one of the diagonal directions. In the limit of a nearly empty band k F → 0, on the other hand, the Fermi surface becomes a circle and ∇ k E(k)| |k|=k F is proportional to the momentum k. For finite but small λ, this picture remains qualitatively valid and the observed angular dependence of the magnetization is indeed similar to the case for λ = 0.
Above, it was assumed that the spin current in the spin  chain is polarized in the z direction, i.e., orthogonal to the plane of the 2D leads. If we choose a different polarization, the spin current through the chain will have the same magnitude because of the pseudo-spin-rotation symmetry of the 1D representation, but the expectation values in the original lattice will differ. Instead of carrying out a separate MPS calculation, one could calculate these quantities by evaluating formula of Eq. (9) with the same correlation functions â † mân and a different transformation matrix P = PR, where R is a unitary matrix that describes the rotation of the pseudo spins.

Conclusions
We have applied the block-Lanczos DMRG technique to investigate the spin transport in a two-terminal setup consisting of a spin chain and 2D tight-binding leads. As long as the spin chain couples only to a single site of each lead, the Lanczos transformation yields an effective 1D model where the leads are semi-infinite chains with nearest-neighbor hopping. While the hopping amplitudes are not uniform, their site dependence is negligible except in the vicinity of the chain edge. The Lanczos transformation done here can be regarded as a specific case of the chain mappings for non-interacting baths based on orthogonal polynomials. 26) Since it is known that these mappings result for typical environments with finite bandwidth in asymptotically homogeneous chains, 27) the effective Hamiltonian we obtained is not surprising and its explicit calculation mostly amounts to determining the strength of the inhomogeneities near the spin chain. These inhomogeneities can appreciably affect the spin transport in the junction because the parameters at the interface need to be finetuned to achieve a sizeable spin current at low temperatures. Qualitatively, however, the behavior of the spin conductance is the same as the case when the hopping strength is assumed to be uniform. The Lanczos transformation thus does reveal any new phenomena in the setup studied here regarding the spin conductance. One could apply the method also to more complicated interfaces, e.g. involving multiple coupled spin chains. In that case, the lead part would become a ladder model after the transformation, with the number of legs equal to the number of spin chains. As realizations of spin chains in solids typically consist of many weakly coupled chains, the block-Lanczos method could be a way towards a more realistic junction model.
Interestingly, the form of the effective Hamiltonian after the transformation does not change when the Rashba spin-orbit coupling is taken into account. Phenomena characteristic of the Rashba model, such as the spin Hall effect, are hidden in the definition of the Lanczos basis states. As a consequence, a spin current entering the lead appears more or less the same in the 1D representation regardless of the strength of the spinorbit coupling λ. The inverse Hall effect and the Rashba precession in the original real-space lattice, on the other hand, occur only for finite λ. By reversing the Lanczos transformation one can calculate quantities in real space and thereby observe these effects. However, since the tight-binding leads by themselves are non-interacting, the approach used here only makes sense if the interacting region plays an important role. Otherwise, many-body techniques such as the DMRG method are not necessary and more efficient methods, e.g., based on Green's functions, are available.