Spin-Hall effect in topological materials: evaluating the proper spin current in systems with arbitrary degeneracies

Since its discovery two decades ago^1,2 the spin-Hall effect (SHE) has become one of the most actively studied topics in modern physics. Its uses range from the inverse SHE used as a detection and characterization tool^3,4,5,6 to the generation of spin torques that flip magnetic bits^7,8. The spin-Hall torque induces magnetization dynamics in spintronic memory devices^9,10,11 and is strong in many topological materials^{12,13,14,15,16}, including, recently, van der Waals heterostructures coupled to WTe₂¹⁷, Mn₃Sn¹⁵, and heavy metals¹⁸. Despite its manifold uses, the underlying mechanisms of the SHE remain somewhat mysterious and poorly understood. It remains difficult to distinguish experimentally between intrinsic (band structure) and extrinsic (disorder) contributions, as well as between the SHE and other mechanisms such as the orbital Hall and Edelstein effects^11,19. Hence, realistic calculations of the SHE for real materials are urgently needed.

The absence of an experimental blueprint for measuring the SHE accurately is related to the inherent difficulty in calculating the spin-Hall current. The conventional spin current is physically meaningless, since the spin precesses as it is transported^{20,21,22,23,24,25}. The proper spin current, based on the equation of continuity and the Onsager relations, is conserved, but contains a torque dipole term that involves the position operator^{21,24,25,26,27,28,29,30,31}. The intra-band elements of the position operator, which are associated with Fermi surface, group velocity, and dipolar effects, are challenging to evaluate in extended systems, since Bloch electrons are delocalized^21,29,30,32. Additional major challenges include the presence of degeneracies in realistic band structures^30,31, and the difficulty of bridging the gap between analytical approaches and density functional theory (DFT), while minimizing the computational cost.

In this work, we exploit recent theoretical breakthroughs in evaluating the proper spin-Hall current in extended systems to overcome these challenges. We develop a new tight-binding approximation (TBA) methodology to determine the intrinsic proper spin current (IPSC), using the physical definition of ref. ²⁴, for different classes of topological materials as well as for strongly spin-orbit coupled metals. This paper breaks new ground: we develop a method for treating arbitrary matrix elements of the position operator, including its diagonal elements, for Bloch electrons with arbitrary degeneracies, and merge the analytical formalism with the TBA approach to yield a blueprint for evaluating the proper SHE for arbitrary band structures. We determine the size and structure of the spin current and its implications for the spin-Hall torque in topological materials. Of the materials we studied, Pt had the largest intrinsic spin current when using both the conventional and conserved definitions. Remarkably, in all classes of topological materials studied, as well as in metals, the proper SHE generally shares the same features. The torque dipole reduces the spin current and causes its energy-dependence to be relatively flat and featureless. This is in sharp contrast to the conventional spin current, which exhibits sharp peaks and dips as a function of energy. Using the proper definition significantly alters the calculated spin current, this cements the need to use the proper definition when making theoretical predictions.

Our work enables us to connect equilibrium density functional theory with non-equilibrium quantum mechanics based on the density matrix, which provides the most complete description of a quantum mechanical system. The results are presented in a form that is directly comparable to the experiment, where the SHE can be inferred from the spin-Hall angle. Whereas the focus here is on the intrinsic case, the method can be extended to disordered systems, incorporating skew scattering and side jump along the lines of ref. ³¹. The extension to disordered and inhomogeneous systems will enable the study of realistic devices and architectures in the most accurate and least computationally expensive way possible.

Generating a spin current typically requires spin-orbit coupling, which causes spin precession and hence non-conservation. The conventional definition of the spin current is the product of the spin and velocity operators (J={hat{v},hat{s}})³³ or a redefined velocity operator³⁴. However, the non-conservation of spin as it is transported makes the conventional spin current physically meaningless: its relationship to spin accumulation is not obvious^35,36, it does not satisfy an equation of continuity or an Onsager relation, and is nonzero even in thermodynamic equilibrium^24,37.

The proper spin current was introduced in Phys. Rev. Lett. 93, 046602 (2004) as a conserved spin current, which satisfies the equation of continuity. This led to the distinction between the conventional spin current, which is convenient but physically meaningless, and the proper spin current, which is conserved, despite being difficult to measure. The proper spin current operator includes the conventional spin current but with an extra contribution, the torque dipole (I={hat{r},dhat{s}/dt}) which arises from spin precession^{24,25,26,27,28,29,38,39}. Its calculation is subtle, involving matrix elements of the position operator between Bloch states. Operators containing the position operator are often difficult to deal with, as in a crystal in which the electron states are Bloch states the density matrix is diagonal in the crystal momentum. However, the position operator is not diagonal in the crystal momentum and couples wave vectors that are infinitesimally spaced.

The complications involved with properly defining the spin current can be avoided by calculating the spin response directly without resorting to the spin current^40,41, which is appropriate when the quantity of interest is the spin accumulation. For spin-Hall related phenomena knowledge of the spin accumulation is generally preferable as it is observable. Furthermore, the relationship between spin current and spin accumulation is unclear⁴¹. Spin accumulation can be calculated via a drift-diffusion, Boltzmann, or multi-band density matrix approach^42,43,44. However, for phenomena that involve spin currents in which the spin accumulation is ambiguous or absent such an approach is not appropriate and knowledge of the spin current is needed to gain physical insight. This is the case for spin torques driven by the spin-Hall effect or for the inverse spin-Hall effect. In spin-Hall torques, spins generated by the spin current flow into a magnetic material and generate a torque on the magnetization. Whereas the spin-orbit torque is generally given by the spin accumulation, in such systems the fate of the spins as they travel across the interface and into the magnet is incompletely understood and it is not clear whether a spin accumulation can also be present at the interface. It is in this context that the calculation of spin currents is indispensable in interpreting spin torque experiments: it can reveal whether they are zero or finite, whether they change sign under certain circumstances, as well as their variation in different materials and with system parameters. The inverse spin-Hall effect is the Onsager inverse of the spin-Hall effect, where a spin current gives rise to a transverse voltage^3,45. In this context, the spin-Hall conductivity, which is directly related to the spin current, is the quantity of interest. Hence, calculations of the spin current are crucial for the study of this effect.

Spin precession is incorporated into the proper spin current such that only the flow of the conserved part of the spin is considered. We believe this definition still allows for the incorporation of extrinsic spin relaxation mechanisms, such as the Dyakonov-Perel mechanism, when scattering is added to the kinetic equation. It has previously been shown that adding spin-orbit coupling to the Born-approximation scattering term correctly captures the interplay between spin precession and scattering in both the weak and the strong scattering regimes⁴⁶. Such a scattering term will still be present in the kinetic equation used to calculate the spin current and spin accumulation, and will affect the spin distribution, regardless of which current operator is traced with the density matrix.

Until now, there have been a number of works on the topic of calculating the proper spin current, each work has a distinct approach to the problem. In this work, which is an extension of refs. ^30,31, we take an approach based on a quantum kinetic theory. Ref. ²⁹ used a semiclassical approach based on wave packet dynamics. Lastly, Ref. ²¹ employed a Keldysh approach. There are important differences between each approach, and comparisons between results are not straightforward.

Our approach is fully quantum mechanical and is equivalent to the Kubo linear response formalism, we evaluate the expectation value of the proper spin current operator operator. Whereas, the formalism in ref. ²⁹ constructs the proper spin current by adding a number of distinct contributions. These are formulated using classical physics and Maxwell’s equations while referencing the center of mass of a wave packet. Although the two approaches are quite different the final expressions for the proper spin current in a non-degenerate system are very similar, only differing by the position of spin and Berry connection matrix elements. However, it should be noted this work extends the results of ref. ³⁰ by extending the formula to systems with degenerate states. In this work, we separate the torque dipole into two contributions I₁ and I₂, which represent the contributions from the degenerate/band-diagonal and non-degenerate matrix elements of the spin operator respectively.

The Keldysh approach of ref. ²¹ calculates the torque dipole by employing a fictitious electromagnetic field. The calculation is semiclassical in the sense that it mixes both position and momentum coordinates. Their results find that the spin current vanishes in spin-1/2 systems regardless of the model. This differs from our findings and those of ref. ²⁹. A detailed comparison is hampered by the lack of generic formulas in ref. ²¹. Furthermore, it is unclear what the analogs to the fictitious electromagnetic fields and explicit gradient expansion are in our approach.

In this work, we extend the result of ref. ³⁰ for the intrinsic proper spin current to systems with arbitrary degeneracies, whereas the formula presented in this previous work is only valid in fully non-degenerate systems. Our general analytical expression for the IPSC spin-Hall conductivity in systems with arbitrary degeneracies is

where m, n and ({n}^{{prime} }) are band indices, E is the external electric field, f_mk ≡ f(ε_mk) the equilibrium Fermi–Dirac distribution, ε_mk the band dispersion, k the wave vector, s the spin operator, and ({{bf{{mathcal{R}}}}}_{mn{boldsymbol{k}}}=langle {u}_{m{boldsymbol{k}}}| {rm{i}}frac{partial {u}_{n{boldsymbol{k}}}}{partial {boldsymbol{k}}}rangle) the Berry connection. The check over the spin term indicates band diagonal matrix elements and elements between degenerate states, and the tilde over the Berry connection indicates matrix elements between non-degenerate states. This is the only intrinsic contribution to the proper spin current, which is shown to flow perpendicular to the applied electric field. Since only the band off-diagonal matrix elements of the Berry connection are included in the calculation the proper spin current is gauge invariant. It can also be finite in the insulating gap, with important implications for TI/FM devices, as we show below.

We consider four different types of materials in our calculations: a topological insulator (Bi₂Se₃), 2D quantum spin-Hall insulator (WTe₂), antiferromagnet (Mn₃Sn) and a heavy metal (Pt). We chose these materials due to their prevalence in spintronic research⁴⁷. We employ density functional theory (DFT) to calculate the electronic and spintronic properties of Bi₂Se₃, WTe₂, Mn₃Sn, and Pt.

Our computational approach involves the utilization of the projector augmented-wave (PAW) method implemented within the Vienna ab-initio Simulation Package (VASP)⁴⁸. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional⁴⁹ was employed, accompanied by standard scalar-relativistic PAW pseudopotentials^50,51. The electronic configurations for the valence electrons of Bi, Se, W, Te, Mn, Sn, and Pt were specified as 5d¹⁰6s²6p³, 6s²6p⁴, 5s²5p⁶5d⁴6s², 5s²5p⁴, 3p⁶3d⁵4s², 4d¹⁰5s²5p², and 5d⁹6s¹, respectively. The experimental crystal parameters for Bi₂Se₃⁵², WTe₂⁵³, Mn₃Sn⁵⁴, and Pt⁵⁵ were adopted, with dimensions set as a = b = 4.14 and c = 28.64 Å; a = 6.28, b = 3.50, c = 19.15 Å; a = b = 5.67 and c = 4.53 Å; and a = b = c = 3.92 Å, respectively. For WTe₂, a Hubbard-like correction was applied within the Dudarev scheme⁵⁶, with a U value of 5.1 eV chosen to accurately reproduce the correct band structure. Furthermore, all computations incorporated spin-orbit coupling (SOC), with a plane-wave basis set cutoff energy established at 350 eV. Convergence criteria for energy optimization were set to 10⁻⁶ eV. The Monkhorst-Pack k-point grid was applied, with the k-mesh tailored for different materials.

The determination of spin-Hall conductivity (SHC) entailed the use of a tight-binding Hamiltonian, derived from localized Wannier functions (WFs)^57,58 projected from the DFT Bloch wavefunctions. Specifically, SHC calculations were performed through Qiao’s method⁵⁸. Subsequent postprocessing was conducted using Wannier90⁵⁹ and WannierTools⁶⁰ codes. We used a k-point grid of 101 × 101 × 101 for calculating SHC. The proper spin current formula was evaluated using the Wannier90 and WannierTools by first decomposing it into a combination of spin and velocity operators, the detailed equations used for this calculation can be found in the supplementary material. Notably, atomic-orbital-like WFs encompassing Bi-p/Se-p, W-d/Te-p, Mn-p/Sn-p, and Pt-spd orbitals were selectively chosen for this analysis. The calculated band structures are in good agreement with ab-initio results, plots of the band structures can be found in the supplementary material.

For our numerical calculations, we separate the torque dipole I into two contributions. The first is

and the second is

where ({S}_{{boldsymbol{k}}}^{E}) is the band off-diagonal part of the k diagonal non-equilibrium density matrix. Further details on the derivation of these expressions can be found in section V and the supplementary material.

Bi₂Se₃

DFT first-principles calculations were performed to analyze the band structure and spin-Hall conductivity (SHC) of Bi₂Se₃, a well-known topological insulator⁶¹. The crystal structure of Bi₂Se₃ is defined by the space group (Roverline{3}m) and associated with the Laue group (overline{3}m1{1}^{{prime} }). We defined an energy range for the Wannier functions (WFs), spanning an inner frozen window from −5.8 to 1.8 eV and an outer disentanglement window from −5.8 to 14.2 eV relative to the Fermi level. This approach yielded 30 spinor WFs, exhibiting p-like characteristics, which were used to construct a tight-binding Hamiltonian that replicates the ab-initio band structure with high fidelity, as illustrated in Fig. S1 and corroborated by the work of ref. ⁶².

In Fig. 1, we plot the SHC ({sigma }_{zx}^{y}) vs energy. The conventional SHC value at the Fermi level is calculated to be 93(ℏ/e)S/cm, decreasing to − 749(ℏ/e)S/cm at 2.82 eV below the Fermi level. The similarity between Fig. 1a and previously presented results⁶³ validates the accuracy of our calculation. The introduction of the torque dipole corrections, as shown in Fig. 1b, significantly alters the SHC spectrum, mainly flattening it due to I₂ canceling with the conventional spin current. We find the proper SHC ({sigma }_{zx}^{y}) to be − 22(ℏ/e)S/cm at the Fermi level. Additionally, the introduction of the torque dipole has shown a similar effect on another spin conductivity tensor component, ({sigma }_{xy}^{z}), where the flattening of the spin conductivities energy dependence is also seen; however, the proper SHC value at the Fermi level increases to 165(ℏ/e)S/cm, contrasting with the 80(ℏ/e)S/cm determined through the conventional spin current.

Panel (a) shows the SHC calculated using the conventional spin current formula, while panel (b) shows the SHC calculated using the proper spin current formula. Panels (c) and (d) contain the two torque dipole contributions I₁ and I₂, respectively.

Full size image

WTe₂

The second material that we investigated was (1{{rm{T}}}^{{prime} })-WTe₂, a monolayer 2D quantum spin-Hall insulator distinguished by its stability and electronic properties among (1{{rm{T}}}^{{prime} }) phase monolayer transition metal dichalcogenides⁶⁴. The choice of WTe₂ was motivated by its promising potential in spintronic and orbitronic devices. WTe₂’s crystal structure is defined by the space group Pm and classified within the Laue group ({2}^{{prime} }/{m}^{{prime} }). Through a meticulously chosen energy window we extracted 44 spinor WFs, we used an inner frozen window spanning from −7.2 to 2.5 eV and an outer disentanglement window ranging from −7.2 to 9.3 eV relative to the Fermi level. Our spinor WFs have W-d/Te-p-like Gaussian forms and were used to construct a tight-binding Hamiltonian that precisely mirrors the ab-initio band structure, as shown in Fig. S1.

Our calculated conventional SHC for WTe₂, ({sigma }_{xy}^{z}), displayed in Fig. 2, shows strong agreement with previously reported values⁶⁵, with the conventional SHC at the Fermi level calculated to be − 174(ℏ/e)S/cm, increasing in magnitude to − 364(ℏ/e)S/cm at 1.92 eV below the Fermi level. The inclusion of the torque dipole corrections, as depicted in Fig. 2b, substantially modifies the SHC spectrum. This is again primarily due to I₂ canceling with contribution from the conventional spin current. Using the proper spin current definition we find ({sigma }_{xy}^{z}) to be − 257(ℏ/e)S/cm at the Fermi level. Additionally, we calculated ({sigma }_{zx}^{y}) which exhibited a similar cancellation due to the introduction of I₂. However, for this tensor component, the proper SHC value at the Fermi level exhibited negligible variation from the conventional definition adjusting to 10(ℏ/e)S/cm from 9(ℏ/e)S/cm.

Panel (a) shows the SHC calculated using the conventional spin current formula, while panel (b) shows the SHC calculated using the proper spin current formula. Panels (c) and (d) contain the two torque dipole contributions I₁ and I₂, respectively.

Full size image

Mn₃Sn

We then calculated the spin conductivity of Mn₃Sn, a non-collinear antiferromagnetic material at room temperature^66,67. This material’s crystal structure is classified by the space group P6₃/mmc and the Laue group 6/mmm. The derivation of 72 spinor WFs, with Mn-d/Sn-p-like Gaussian forms, was done using inner frozen and outer disentanglement windows spanning from −7.1 to 0.7 eV and −7.1 to 15.1 eV relative to the Fermi level, respectively. We again used the WFs to construct a tight binding model that accurately replicated the ab-initio band structure (Fig. S1)⁶⁸, demonstrating the accuracy of our tight-binding model.

As depicted in Fig. 3, conventional SHC at the Fermi level exhibits sharp variation across different energy levels. At the Fermi level, the conventional SHC is calculated to be 61(ℏ/e)S/cm changing to − 1044(ℏ/e)S/cm at 1.26 eV below the Fermi level. Our results for the conventional SHC are consistent with previously published data⁶⁹. The inclusion of the torque dipole significantly alters the SHC spectrum, as illustrated in Fig. 3b, flattening out the sharp energy dependence. There is a plateau in the spin conductivity around the Fermi energy, this is likely due to the Fermi energy being in the bulk band gap. However, there may be other contributing factors, as the plateau extends past the band gap. Among the materials studied, Mn₃Sn seems to have the most profound correction from I₂, it almost completely cancels with the conventional spin current. Hence, I₁ is the primary factor in determining the proper SHC spectrum, this is clear when comparing Fig. 3b, c. We find the proper SHC ({sigma }_{xy}^{z}) at the Fermi level to be 53(ℏ/e)S/cm. Furthermore, another component of the spin conductivity tensor, ({sigma }_{zx}^{y}), has similar cancellation between I₂ and the conventional spin current. We find the value of ({sigma }_{zx}^{y}) to be 47(ℏ/e)S/cm using the proper definition as opposed to − 91(ℏ/e)S/cm using the conventional definition.

Panel (a) shows the SHC calculated using the conventional spin current formula, while panel (b) shows the SHC calculated using the proper spin current formula. Panels (c) and (d) contain the two torque dipole contributions I₁ and I₂, respectively.

Full size image

Pt

The last material we considered was Pt, a heavy metal with a significant role in spintronics due to its strong spin-orbit coupling^7,70. The crystal structure of Pt is classified by the space group (Fmoverline{3}m) and the Laue group m(overline{3})m. We derived 18 spinor WFs, which exhibit Pt-spd-like Gaussian forms, using an inner frozen window and an outer disentanglement window, spanning from −11 to 6.7 eV and −11 to 33.7 eV relative to the Fermi energy, respectively. We again used the WFs to create a tight binding model that replicates the ab-initio band structure accurately, as depicted in Fig. S1.

Figure 4 plots the SHC, ({sigma }_{xy}^{z}), at different energy levels. At the Fermi level, the conventional SHC reaches 2180(ℏ/e)S/cm declining and changing the sign to a value of − 2177(ℏ/e)S/cm at 4.28 eV below the Fermi level. Our calculation of the conventional SHC aligns with previous calculations⁵⁸. The inclusion of the torque dipole results in a proper SHC ({sigma }_{xy}^{z}) of 1305(ℏ/e)S/cm at the Fermi level. The energy dependence of the proper SHC, depicted in Fig. 4b, shows that the torque dipole introduces a small but non-negligible correction to the SHC in Pt, differing from the more pronounced effects observed in materials like Mn₃Sn. This highlights the nuanced role these corrections play in the overall spectral behavior.

Panel (a) shows the SHC calculated using the conventional spin current formula, while panel (b) shows the SHC calculated using the proper spin current formula. Panels (c) and (d) contain the two torque dipole contributions I₁ and I₂, respectively.

Full size image

In general, our spin conductivity results show that the magnitude of the proper spin current will generally be smaller than the conventional spin current. This is important as it highlights that previous theoretical works employing the conventional definition may overestimate the size of the spin-Hall effect. Our figures also show that the energy dependence of the proper spin current is flatter than the conventional spin current, which seems to have sharp fluctuations that the proper spin current does not. The reduction in magnitude and flatter energy dependence is due to the torque dipole canceling exactly with terms in the conventional spin current. This is prominently seen in Figs. 1 to 4 when comparing the second part of the torque dipole with the conventional spin current.

Our results for the spin conductivity in Bi₂Se₃ for ({sigma }_{zx}^{y}) predict the spin conductivity to have a different sign to previous results using a k ⋅ p Hamiltonian. These previous results found the spin conductivity to be 8(ℏ/e)S/cm³⁰, whereas we find the spin conductivity to be ({sigma }_{zx}^{y}=-22(hslash /e))S/cm. This difference could be due to previous calculations employing a Hamiltonian that had a magnetization, which was not considered in this work. Nonetheless, our results remain two orders of magnitude lower than the larger spin torque measurements observed in Bi₂Se₃ devices^12,71, hence, effects such as the spin-transfer torque and Rashba–Edelstein effect are likely responsible for these larger torques⁷². The spin conductivity ({sigma }_{zx}^{y}) has opposite sign when calculated using the proper definition. Hence, spins generated via the SHE have opposite sign to spins generated via the REE^73,74,75,76, this is consistent with the experimental results that find the spin torque efficiency increase for thinner TI samples where the SHE is smaller⁷⁷.

A spin-orbit torque experiment recorded an out-of-plane field-like torque conductivity in WTe₂ of 45(ℏ/e)S/cm¹³. This number can be compared with our results for ({sigma }_{zx}^{y}). Comparing the experimental result with our result of 10(ℏ/e)S/cm, our calculated spin conductivity is the same order of magnitude, though the experimental result is substantially larger. However, it should be noted that these experimental results are for bilayer WTe₂, whereas our results are for a monolayer structure. Another study of spin torques in WTe₂ devices found the spin conductivity to be up to two orders of magnitude larger⁷⁸, however, it ascribes this large torque to surface state effects. The spin conductivity calculated using the conventional definition is of a very similar magnitude ({sigma }_{zx}^{y}=9(hslash /e))S/cm.

Our results in Mn₃Sn, shown in Fig. 3, find the magnitude of proper spin conductivity at the Fermi level ({sigma }_{xy}^{z}=53(hslash /e))S/cm to be smaller than the conventional spin conductivity ({sigma }_{xy}^{z}=61(hslash /e))S/cm. Our results for the energy dependence of the spin conductivity using the conventional definition agree with previous calculations^69,79,80. The spin conductivity in Mn₃Sn has been measured to be σ_SH = 46.99 ± 20.63(ℏ/e)S/cm in an ISHE experiment⁸¹. Both our proper and conventional spin current results fall within the range of the measured spin conductivity, and our results have the same sign as the experiment.

The conventional and proper spin conductivies in Mn₃Sn exhibit interesting behavior as a function of the relativistic spin-orbit strength, we have provided figures of the spin conductivity at different spin-orbit strengths in the supplement. Whereas the conventional and proper spin currents would normally reduce to zero as the spin-orbit coupling is reduced to zero, the proper and conventional spin conductivities remain nonzero. This is likely due to an additional source of non-relativistic spin-orbit coupling that arises from the non collinear antiferromagnetic properties of Mn₃Sn, such spin-orbit coupling has been attributed to observations of the anomalous Hall effect in Mn₃Sn and other antiferromagnets⁸². The conventional spin conductivity and second torque dipole term I₂ remain largely unchanged for different values of the relativistic spin-orbit coupling. This implies that the spin texture in Mn₃Sn gives rise to a large spin-orbit coupling. Furthermore, the torque dipole almost exactly cancels with the conventional spin current at all values of the spin-orbit coupling strength, hence, this non-relativistic spin-orbit field must also cause a large amount of spin precession. This result deserves further investigation, however, it lies somewhat beyond the scope of the present paper.

Analysis of spin-Hall effect measurements in Pt finds its spin conductivity to be (0.7–1.7) × 10³(ℏ/e)S/cm⁸³. Our result for the spin conductivity using the proper definition is within this range, whereas the result using the conventional definition falls just above the upper end of the range. As is shown in Fig. 4, at the Fermi energy we find the proper spin conductivity to be ({sigma }_{xy}^{z}=1.31times 1{0}^{3}(hslash /e))S/cm and the conventional spin conductivity to be ({sigma }_{xy}^{z}=2.18times 1{0}^{3}(hslash /e))S/cm.

The spin current formula we present here extends the work in ref. ³⁰ to systems with arbitrary degeneracies. Whereas, the previous formula was only valid in fully non-degenerate systems. Furthermore, in this work, we have demonstrated that this method for calculating the proper spin current is not restricted to k ⋅ p models and can be straightforwardly used with DFT calculations. The use of more accurate models that are accurate beyond the band center is necessary for spin current calculations, since all wavevectors in the Brioullin zone must be summed over and filled bands can generate nonzero spin currents.

For the proper spin current to be strictly conserved, the expectation value of the torque, ({rm{Tr}}(hat{t}hat{rho })), needs to cancel so that globally there is no net spin generation in the system^24,26. In other words, the torque density ({rm{Tr}}(hat{t}hat{rho })) vanishes but the torque dipole density ({rm{Tr}}({hat{t},hat{{boldsymbol{r}}}}hat{rho })) is finite. This is true for the models we consider in this paper.

It is well known that the spin accumulation depends crucially on boundary conditions^24,39,43,84, and ref. ⁴¹ demonstrated quantitatively that the spin accumulation can be determined without reference to the spin current. Hence, in systems where spin accumulation is the quantity of interest direct calculation of the spin accumulation without reference to the spin current can be favorable. However, accurately defining the boundary conditions, which are often unknown, is a limiting factor for such calculations.

Our calculation is indispensable in systems in which the spin current does not lead to a spin accumulation. Such systems, which include TI/FM interfaces, are in fact used to infer the presence of a spin current. Since the spin current does not couple to any measurable quantity, its detection is primarily through indirect processes, for example by measuring spin-torque driven magnetization precession^85,86, spin-current induced second-harmonic optical effects^87,88, the inverse spin-Hall effect^3,4,5, and X-ray pump-probe measurements⁸⁹.

The next important step in gaining a complete understanding of the spin-Hall effect is to extend the theory to account for spin currents generated via extrinsic mechanisms. Our formalism can be straightforwardly extended to the case of disorder^42,90,91. A blueprint for the calculation of extrinsic spin currents has been presented in ref. ³¹, and a calculation was done with a k ⋅ p Hamiltonian. In this work, it was shown that spin currents due to spin-orbit scattering effects such as skew and side jump scattering are of a similar order of magnitude to intrinsic spin currents. Furthermore, these contributions will appear at zeroth order in disorder, making them indistinguishable from intrinsic spin currents. This is consistent with previous results for the anomalous Hall effect^92,93, and for the conventional spin current^21,94. Recently it was shown that extrinsic mechanisms in a similar effect, the orbital Hall effect, dominate the intrinsic orbital Hall effect⁹⁵, further highlighting the need to extend our theory to include disorder. The remaining challenge is in extending this formalism to ab-initio and DFT calculations as was done in this work for intrinsic spin current. Including general disorder effects in such calculations is notoriously difficult. We suggest the use of the same simple disorder model as refs. ^31,42,90 and to calculate the disorder contribution iteratively as was done in Ref. ³¹ or by using a simple relaxation time approximation or Bloch lifetime.

We outline in this section the derivation leading to Eq. (1). We first discuss our methodology for dealing with arbitrary matrix elements of the position operator, which is vital in obtaining the correct expression for the torque dipole. Next, we apply this methodology to determine the full expression for the proper spin current in systems with arbitrary degeneracies. Finally, we discuss briefly the procedure for extending this formalism to disordered systems.

The position operator between Bloch states

For a system described by a single-particle density operator (hat{rho }), the expectation value of an arbitrary operator (hat{O}) is given by (langle hat{O}rangle =,{text{Tr}},,hat{O}hat{rho }). Operators containing the position operator can be difficult to deal with, as in a crystal in which the electron states are Bloch states the density matrix (hat{rho }) is diagonal in the crystal momentum. Whereas the matrix elements of the position operator in the crystal momentum representation couple wave vectors that are infinitesimally spaced. Here we outline the derivation of a general expression for the trace of an operator with the position operator. We consider the trace of the position operator with some operator (hat{Lambda }),

Since the wavefunctions are Bloch states we can express them as (leftvert {psi }_{n{boldsymbol{k}}}rightrangle ={e}^{i{boldsymbol{k}}cdot {boldsymbol{r}}}leftvert {u}_{m{boldsymbol{k}}}rightrangle). Therefore the position operator can be expressed as

Evaluating (4) and making the substitution k = k₊ and ({{boldsymbol{k}}}^{{prime} }={{boldsymbol{k}}}_{-}) where k_± = k ± Q/2 we find that

where ({mathcal{D}}) is a covariant derivative defined as³²

The wavevector off-diagonal nature of the position operator is taken care of by the substitution of the infinitesimal wavevector Q. The differential terms in (7) represent the phase of the wave function in more conventional evaluations, while the Berry connection represents the contribution due to the change of the basis states between infinitesimally-separated wave vectors.

In (7) we consider density matrix elements that are infinitesimally off-diagonal in wavevector, whereas we consider the Berry connection to be purely diagonal in wavevector. This is because the Berry connection is defined as

and definitionally cannot have elements off-diagonal in the wave vector. It is straightforward to show that all formulas are gauge invariant.

Proper spin current in degenerate bands

The system is described by a single-particle density operator (hat{rho }), which obeys the quantum Liouville equation

where (hat{H}) is the total Hamiltonian of the system. We will consider an arbitrary Hamiltonian and focus on a clean system, deferring the treatment of disorder to a future publication. The conserved spin current operator ({hat{{mathcal{J}}}}_{j}^{i}=d/dt,({hat{r}}_{j}{hat{s}}_{i})) is the time derivative of the spin dipole operator. Taking trace with the single particle density matrix operator (hat{rho }), we separate the conserved spin current into the conventional spin current and torque dipole contributions ({{mathcal{J}}}_{j}^{i}=frac{1}{2},{rm{Tr}},hat{rho },{{hat{s}}_{i},{hat{v}}_{j}}+frac{1}{2},{rm{Tr}},hat{rho },{{hat{t}}_{i},{hat{r}}_{j}}), with the velocity operator ({hat{v}}_{j}=d{hat{r}}_{j}/dt) and the torque ({hat{t}}_{i}=d{hat{s}}_{i}/dt), both diagonal in wave vector in the crystal momentum representation. The conventional spin current ({J}_{j}^{i}=frac{1}{2}{rm{Tr}}{hat{s}}_{i}{{hat{v}}_{j},hat{rho }}) is straightforward to evaluate. In contrast ({I}_{j}^{i}=frac{1}{2}{rm{Tr}}hat{rho }{{hat{t}}_{i},{hat{r}}_{j}}) stemming from the torque requires some work to deal with the position operator ({hat{r}}_{j}). In order to deal with the position operator in the Bloch representation and evaluate the torque dipole we perform the manipulations outlined in the previous section, taking (hat{Lambda }={hat{t}}_{i}hat{rho }) and evaluating (6) we find

Where, ({{mathbf{Xi }}}_{{boldsymbol{k}}}={mathcal{D}}{{rho }}_{{boldsymbol{k}}})³². It is easy to prove Eq. (10) is gauge invariant. The torque ({t}_{i}=frac{i}{hslash }[{H}_{0},{s}_{i}]) is purely off-diagonal in band index. Hence, when evaluating the trace (10) only band off-diagonal elements of Ξ are required. To find Ξ the covariant derivative is applied to the quantum kinetic equation from refs. ^42,90, the resulting kinetic equation is

Where the k diagonal density matrix ρ_k is found using the original kinetic equation. The new kinetic equation is solved for Ξ_k in an identical manner, the details of this calculation can be found in the supplementary material.

The process of applying the covariant derivative to the kinetic equation and solving for Ξ_k is equivalent to the approaches of refs. ^30,31. In these works the quantum kinetic equation was expanded to linear order in the infinitesimal off-diagonal wavevector Q and then solved for ({rho }_{{{boldsymbol{k}}}_{+}{{boldsymbol{k}}}_{-}}). The quantity Ξ is just the linear order correction to the density matrix in Q, such that ({rho }_{{{boldsymbol{k}}}_{+}{{boldsymbol{k}}}_{-}}={rho }_{{boldsymbol{k}}}+{boldsymbol{Q}}cdot {{mathbf{Xi }}}_{{boldsymbol{k}}}). Furthermore, the expanded kinetic equation used in Refs. ^30,31 is essentially identical to (11). We would like to emphasize that the approach employed for the position operator in this work has the advantage of being more general and can be applied to evaluate any dipole operator in a homogeneous system³².

Solving (11) and taking the trace (10) yields two contributions to the torque dipole. The first contribution is (2) when expanding this expression half of the terms will exactly cancel with the conventional spin current, the other half will add to the spin current and yield the expression in (1). The second contribution to the torque dipole is (3), all of the terms in the expression will cancel exactly with the remaining contributions from the conventional spin current. Hence, the intrinsic contribution to the conventional spin current is contained in (1).