Kubo formula for dc conductivity: generalization to systems with spin-orbit coupling

We revise the Kubo formula for the electric dc conductivity in the presence of spin-orbit coupling (SOC). We discover that each velocity operator that enters this formula differs from $\partial H/\partial \boldsymbol p$, where $H$ is the Hamiltonian and $\boldsymbol p$ is the canonical momentum. Moreover, we find an additional contribution to the Hall dc conductivity from noncommuting coordinates that is missing in the conventional Kubo-Streda formula. This contribution originates from the"electron-positron"matrix elements of the velocity and position operators. We argue that the widely used Rashba model does in fact provide a finite anomalous Hall dc conductivity in the metallic regime (in the noncrossing approximation) if SOC-corrections to the velocity and position operators are properly taken into account. While we focus on the response of the charge current to the electric field, linear response theories of other SOC-related effects should be modified similarly.


I. INTRODUCTION
Approximately 25 years ago, out-of-equilibrium phenomena caused by spin-orbit coupling (SOC) started gaining strong attention.Among such phenomena are the anomalous Hall effect (AHE) [1], the spin Hall effect (SHE) [2], the spin-orbit torques (SOT) [3], and their thermoelectric and spin caloritronics counterparts [4].These effects continue to be the subject of very active and fruitful research today.Many ideas for technological applications based on them have been suggested and developed, such as SOT magnetoresistive randomaccess memory [5,6], magnetic racetrack memory [7], and neuromorphic computing [8].AHE is often considered as a hallmark of the nontrivial topology of the system [9].It also serves as an important probe of magnetic order [10][11][12][13].
From the computational point of view, these effects can be analyzed using kinetic-type equations or different versions of the linear response Kubo formula [14].The latter approach normally allows one to compute the corresponding quantities in the most direct and reliable way.In this paper, however, we show that, for systems with SOC, the existing Kubo formulas miss important contributions and should be revised.One of the reasons for this lies in the fact that the unitary transformation that diagonalizes a relativistic Hamiltonian does not, in general, diagonalize other operators of physical observables.While the spectrum of the system and its ground state properties are not affected by this, linear response theories are.Another important issue is that the operator ∂/∂p is not invariant under p-dependent unitary transformations.Below, using the electrical dc conductivity as an example, we demonstrate that a consistent treatment of observables and perturbations reveals previously overlooked contributions to out-ofequilibrium phenomena originating from SOC.
In the nonrelativistic (nr) theory, the Kubo formula for the full dc conductivity tensor of an electron gas was derived in 1971 by Bastin et al., using the Green's function formal- * iv.a.ado@gmail.comism [15].It can be cast in the form where G ± = (ε ± i0 − H) −1 denotes the retarded (+) and advanced (-) Green's functions, H is the Hamiltonian, f ε is the Fermi-Dirac distribution function, Ω represents the system volume, v = [r, H]/i = ∂H/∂p is the velocity operator, p is the canonical momentum, r represents the position operator, and prime denotes a derivative with respect to ε.The system is assumed to be noninteracting.In 1982, Středa published a modified version of the Kubo-Bastin formula, in which he split the nonrelativistic conductivity tensor into a sum of two distinct contributions [16], σ nr αβ = σ nr-I αβ +σ nr-II αβ .By applying a partial integration to Eq. ( 1) and then taking a half-sum of the obtained result and the original formula, one can express the Středa's splitting as At the end of the 1990s, when the interest in SOC-related out-of-equilibrium phenomena increased, Středa's formula [Eqs.(2)] and its modifications became standard tools for computing tensors that describe such effects.For example, to obtain the dc spin Hall conductivity, one would replace the velocity operator v α in Eqs.(2) with the spin velocity operator.For SOT, v α would change to the spin operator.For AHE, the Středa's formula was used directly.
In the vast majority of all such computations, the velocity operators are defined as v = ∂H/∂p.This is because the commutation of the Hamiltonian with r corresponds to differentiating the former with respect to p.For AHE, this is also in line with the conventional definition of the electric current density operator, which is considered proportional to the functional derivative δH/δA.Here A is the vector potential and, according to the Peierls substitution, it enters the Hamiltonian through the combination π = p − (e/c)A [17].
At the same time, it is known that the position operator r acquires a relativistic (or effectively relativistic) correction ∆r both in vacuum [18] and in crystals [19,20].This correction is sometimes called the Yafet term [21,22].Upon commuting the modified coordinate operator with the Hamiltonian, one finds that the velocity operator should be corrected as well, where v an = [∆r, H]/i is what Adams and Blount called the "anomalous velocity" in 1959 [23][24][25].Sometimes it is confused with the spin-dependent part of the velocity operator.Now we can formulate the problem.On the one hand, we want to use Kubo formulas for effects like AHE, SHE, and SOT, which are relativistic (as they originate from SOC).On the other hand, Eqs.(2) were derived for a nonrelativistic system.A priori, it is unclear how the anomalous velocity v an and the Yafet term ∆r enter Kubo formulas in the presence of SOC.Thus, a separate derivation of the Kubo formulas is required in this case to include all effects of SOC consistently.
For the conductivity tensor σ αβ , this task was approached by Crépieux and Bruno in 2001.They found an elegant expression for the Green's function of the Dirac Hamiltonian in terms of the nonrelativistic Green's function [26].In Ref. [27], they used this expression to project the Středa's formula derived in the Dirac theory onto the electronic branch.This was done by expanding the latter formula up to the order 1/c 2 , where c is the speed of light.According to the result of Crépieux and Bruno, Eqs.(2) with v = ∂H/∂p still hold when the leading order relativistic corrections, including SOC, are taken into account.In other words, the anomalous velocity v an is absent in their formulas.The origin of this is that the Yafet term is missing in the definition of the velocity operator (see, e.g., Eqs. ( 28) and (30) of Ref. [26]).In 2015, another work [28] extended some of the results of Crépieux and Bruno, but also disregarded the Yafet term.
Below, we derive the Středa's formula for electrons in the presence of SOC, starting from the Dirac Hamiltonian.We use a different method of projecting the relativistic Kubo formula onto the electronic branch.Our approach can be used to reach any required accuracy with respect to the parameter 1/c in a unified fashion.We find, in particular, that each velocity operator that enters the Středa's formula for electrons contains the anomalous part.This applies not only to the velocity operator corresponding to the computed observable but also to the velocity operator representing the Hamiltonian's perturbation by the electric field.
The important principle that we use is that the unitary transformation applied to the Dirac Hamiltonian to remove positrons should also be applied to all other operators.In particular, SOC corrections must be computed not only for effective Hamiltonians, but for observables and perturbations too.Moreover, we show that a proper linear response theory for a system with SOC should also take into account the "electron-positron" matrix elements of the involved operators.While the former principle was already formulated, e.g., in Ref. [29], the latter one, to the best of our knowledge, was not discussed before.We note that while our analysis focuses on the conductivity tensor, our approach can be applied to other SOC-related effects.In what follows, we also discuss several implications of the obtained results.

II. UNITARY TRANSFORMATION OF THE DIRAC PICTURE OPERATORS
Let us outline how, by a unitary transformation, operators of the Dirac theory are projected onto the electronic branch.The Dirac (D) Hamiltonian reads where σ denotes Pauli matrices, π = −i ∇ − (e/c)A, is the diagonal (d) component of the Hamiltonian, while is the off-diagonal (od) one.We want to diagonalize H (D) up to a desired order in 1/c by applying a unitary transformation.Such a procedure is sometimes called the Löwdin partitioning [30,31].We define a unitary operator U = e iW with the help of where w † is a Hermitian conjugate of w.Diagonal and offdiagonal components of the unitary transformed Dirac Hamiltonian H (D) U = U † H (D) U are then represented by [30] where the notation n times (9) has been introduced.We need to solve the equation H (D) U,od = 0 for w.We begin by considering only the first two terms on the right hand side of Eq. (8b).This gives H (D)  od + [H (D) d , iW ] = 0, which implies The difference H (D) + − H (D) − = 2mc 2 is a large parameter and we can solve Eq. ( 10) iteratively.The second iteration gives However, this is not yet the full result for w up to the order 1/c 3 .We need to consider two more terms of the expansion (8b).By doing so, we find All other contributions to w are at least of the order 1/c 5 .Substituting Eq. ( 12) into Eqs.( 8), we deduce that the unitary transformed Dirac Hamiltonian becomes where H (el) and H (p) are determined by Eq. (8a).We see that the off-diagonal elements of the Hamiltonian H (D) have been eliminated up to the order 1/c.Consequently, H (el) and H (p) , which are the Hamiltonians of electrons (el) and positrons (p), respectively, are decoupled from each other with the corresponding accuracy.Such a decoupling procedure can be performed up to any order of 1/c.Importantly, not just the Hamiltonian H (D) , but all other operators transform similarly.However, the essential difference is that, by eliminating the off-diagonal elements of the Hamiltonian H (D) , we produce (or modify) off-diagonal (electronpositron) terms in other operators.For position and velocity operators, r (D) U , we will denote the off-diagonal terms as δr and δv, It is straightforward to relate the elements of r (D) U and v (D) U to each other.Since, in the Dirac theory, velocity is a commutator v (D) = [r (D) , H (D) ]/i (note that r (D) = r) and since unitary transformations preserve commutators, we observe that Hence, for electrons, a familiar relation between velocity and position, v (el) = [r (el) , H (el) ]/i , holds.In addition, we have derived an important result: δv = (δrH (p) − H (el) δr)/i .We will use it below.
Let us also give the explicit results for the velocity and position operators up to the order 1/c 2 , We clearly see that the velocity operator v (el) differs from the p-derivative (or π-derivative) of the electron Hamiltonian [32] The difference between the velocity operator of Eq. (16b) and ∂H (el) /∂p equals (− /4m 2 c 2 )[∇V × σ].This is precisely the expression for the anomalous velocity operator stemming from the leading order relativistic corrections in vacuum.It was conjectured in Appendix A of Ref. [24].Several later works [9,20,29,33,34] have also mentioned this term.

III. RELATIVISTIC CORRECTIONS TO THE ST ŘEDA'S FORMULA
The relativistic (r) Kubo-Bastin formula can be derived in the Dirac theory by using standard time-dependent perturbation formalism.No matrix projections are needed for this.Hence no ambiguities appear.We can also assume that the unitary transformation U has already been applied to the derived formula so that it is expressed as Eq. ( 18) has precisely the same form as Eq. ( 1) does, the only difference is the notations.The velocity operators here do not contain the vector potential corresponding to the electric field because the latter is already included in Eq. ( 18) through the diamagnetic term.Integration in Eq. ( 18) is performed over the entire real line.We can split this single integration domain into two disjoint domains: one that corresponds to positrons and another one that corresponds to electrons.The two integrals that we obtain from Eq. ( 18) by performing such a procedure determine Kubo-Bastin formulas for positrons and electrons, respectively.If the electron and positron branches are well separated, each of these two Kubo-Bastin formulas can be formulated exclusively in terms of projections of the involved operators onto the corresponding branch.Below, we derive such a "projected" Kubo-Bastin formula for electrons.For this we also assume that the Fermi level lies in the spectrum of H (el) [35] and that the temperature of the system is sufficiently small as compared to mc 2 .Under these assumptions, the Kubo-Bastin formula for positrons simply provides a constant contribution to the total conductivity tensor.For brevity, we do not write it explicitly anywhere below.
Let us analyze the Green's functions that enter Eq. ( 18), We have already shown that H (D) U can be assumed diagonal with any desired accuracy.Therefore its inverse is diagonal too.Moreover, the positronic component of the latter is the same for the retarded and advanced Green's functions (we consider only the Kubo-Bastin formula for electrons).Hence, from Eq. ( 13) one immediately derives Next we substitute Eqs. ( 14) and ( 19) in Eq. ( 18) and apply the matrix trace.This leads us to where tr represents the operator trace over 2-component spinors.The first term of the integrand above generalizes the nonrelativistic Kubo-Bastin formula of Eq. ( 1), the second term requires further treatment.In Eq. ( 15) we have established a relation between δr and δv.In terms of the Green's functions, it is reformulated as Employing this result, its conjugate, and we obtain for the second term of the integrand in Eq. ( 20): From the commutation relation [r (D)  α , r (D)  β ] = [r (D) U,α , r (D) U,β ] = 0, one can also deduce We substitute Eqs. ( 22) and ( 23) in Eq. ( 20) and apply the same procedure that was used to derive Eqs.(2) from Eq. ( 1).This allows us to finally express the total dc conductivity tensor as σ αβ = σ I αβ + σ II αβ + σ III αβ , where and we refrained from using additional letter superscripts to denote the tensors' elements.Eqs.(24) provide the Středa's formula for electrons that consistently takes into account all effects of SOC.The first two terms, σ I αβ , σ II αβ , correspond to the conventional nonrelativistic Eqs. ( 2), albeit with modified velocities.The third term, σ III αβ , to the best of our knowledge has been missing from the Středa's formula so far.To be more precise, it was never derived as a contribution that should be added to Eqs. (2).
Ref. [36] derived it from Eqs. ( 2) and argued that it was the only intrinsic antisymmetric contribution to the dc conductivity tensor.This is however not the case because another important term was lost during that derivation.Moreover, as we have explained, Eqs. ( 2) should not be used as a starting point for analysis of systems with SOC.Contributions to transport effects from noncommuting coordinates were also derived in Ref. [34].The same work presented a version of the Kubo formula for SHE that takes into account the anomalous velocity.Both results were obtained by adding relativistic corrections to the velocity and position operators in the effective ("projected") model for electrons.However, the electron-positron matrix elements of the same operators were not taken into account in that work.We note that the term σ III αβ defined in Eq. (24c) originates from precisely these matrix elements.
In the context of SHE, the anomalous velocity operator in effective models was also considered, using kinetic-type equations, in Refs.[33,37].
We emphasize that each velocity operator in Eqs. ( 24) contains the anomalous velocity and, in general, is not equal to ∂H/∂p.We also note that the position operator in Eq. (24c) differs from r, as it contains the Yafet term.Expansions of the velocity and position operators up to 1/c 2 (for electrons) are provided by Eqs.(16).The relation between them that is valid in all orders reads Eqs. ( 24) is the first central result of this paper [38].

IV. ANOMALOUS HALL EFFECT IN THE RASHBA MODEL
We have developed a scheme suitable for derivation of the Středa's formula for Hamiltonians with a block structure where one of the blocks needs to be removed.Our scheme can be applied to crystals, particularly to narrow gap semiconductors, where spin-orbit effects are often much stronger than in vacuum [20,31,39].A widely used description for such systems is provided by the 8-band Kane model [40].It describes, in particular, electrons in the conduction band of semiconductors with zinc blende symmetry (e.g., GaAs, InSb).In such systems, the effective SOC originates from the interband matrix elements between conduction and valence bands [19,31].The latter are analogous to the matrix elements denoted by h in Eq. ( 6).One can eliminate them and then derive the Středa's formula, Eqs.(24), just like we did for the Dirac system.
For a sample confined along one spatial direction by an asymmetric potential, the conduction electron Hamiltonian in the Kane model reduces to the famous two-dimensional Rashba (R) Hamiltonian [41] (if the bulk inversion asymmetry is weak) [21,31,42].In the presence of the internal magnetization in the z-direction, the latter reads where m * represents the effective mass, M describes magnetization (or the exchange field), and α R = λ ∇ z V conf is proportional to the averaged z-derivative of the confining potential (we assume that the confinement is along z).The proportionality coefficient λ is the analog of the relativistic prefactor /4m 2 c 2 in front of the last term in Eq. ( 17).The analog of the π 4 term is disregarded in H (R) , while the confining potential itself and its second derivative correspond to a constant energy shift.Here we ignore "relativistic" corrections associated with the term M σ z .However, they can also be important.
The model of Eq. ( 26) is particularly famous for predicting massive cancellations of SOC-related effects [43][44][45].The absence of dc AHE in the metallic regime of this model (in the noncrossing approximation, for Gaussian white-noise disorder) has been a subject of a very extensive debate [1,[45][46][47][48][49][50][51][52].We argue that the cancellation of the dc AHE in the Rashba model does not occur if correctly defined velocity (and position) operators are used.Instead of the operator v = π/m * − α R [e z × σ], which was used before, one should be using the operator analogous [53] to that of Eq. (16b) [54]: where the last term is the anomalous velocity.
After setting A = 0 in Eq. ( 27) we can go along the lines of, e.g., Refs.[48,55], and from Eqs. (24a), (24b) find for the Rashba model: where δ SO = α 2 R m * /M , µ M = µ/M , and µ is the chemical potential [56].This result assumes the following conditions: the noncrossing approximation, Gaussian white-noise disorder, zero temperature, and metallic regime (both bands partly filled).Detailed explanations of these conditions' meanings can be found in Refs.[1,48,55].We note that, in the derivation of Eq. ( 28), we disregarded relativistic corrections from the disorder potential.We also do not compute constant contributions to conductivity from the fully occupied lower bands of the Kane model.
In addition to the already finite σ I xy in the Rashba model, we find a completely new contribution to the anomalous Hall dc conductivity from Eq. (24c).Indeed, we have to assume that the coordinates are defined as and thus do not commute: [x (R) , y (R) ] = −2iλ σ z , where we have disregarded magnetic fields.Hence, σ III xy is proportional to the equilibrium spin polarization in the z-direction.This quantity was computed in a different context in Appendix E of Ref. [55].Using that result, we find from Eq. (24c), in the metallic regime, In the Kane model, λ/ = δg/4mǫ, where ǫ is of the order of the band gap, δg is the effective g-factor minus the bare g-factor, and m is the free electron mass [19,57].By substituting this relation into Eq.( 30), we obtain an expression for σ III xy with a more transparent meaning: Up to a numerical prefactor, this intrinsic contribution to the anomalous Hall dc conductivity equals the ratio between the magnetization-induced spin splitting and the characteristic band gap energy.This contribution survives even without the Rashba coupling as long as λ is finite.Indeed, in the latter case, relativistic corrections to the Hamiltonian H (R) and to the velocity operator vanish.However the spin-dependent part of the position operator still provides finite σ III xy proportional to the relativistic parameter λ.

V. DISCUSSION
The Středa's formula that we derived, Eqs. ( 24), has two differences in comparison with its nonrelativistic (or fully relativistic [58]) version given by Eqs.(2).The first difference is the deviation of the velocity operators from ∂H/∂p.It roots in the fact that the unitary transformation used to project the total Hamiltonian to a band (or a set of bands) depends on the canonical momentum p.Therefore, in the projected theory, the meaning of the operator ∂/∂p changes.The second difference is the novel term σ III αβ determined by the "electronpositron" matrix elements of the velocity operator.Depending on the system, it may happen that only one of the two differences is present.For example, in the absence of the asymmetric potential in the Kane model, the corresponding anomalous velocity vanishes.The same happens in most of the models describing topological insulators.
If one does not perform any band projections, the differences do not occur at all.Thus, for example, the Středa's formula of Eqs.(2) remains a valid tool for a fully relativistic computation of AHE [59].The conventional Berry-phase theory of AHE [9] is also justified, if the full matrix structure of the Hamiltonian is considered.For effective models, however, this theory has to be modified.
Finally, it should be noted that the importance of the interband components of the velocity operator was understood in several works on AHE, for two-level systems [60].Computation of all Feynman diagrams involving such components represents a particular case of the analysis of σ III .

VI. CONCLUSIONS
Starting from the Dirac theory, we have derived the Kubo formula for the dc conductivity tensor (or the Středa's formula) of electrons in the presence of SOC.Each velocity operator in this formula contains the anomalous part.We also found an additional contribution to the Hall conductivity from noncommuting coordinates.Our scheme can be used to derive the Středa's formula for other Hamiltonians with a block structure.We have shown that AHE does not vanish in the Rashba model.We argue that Kubo formulas for other SOCrelated effects should also be revised.