Semiclassical equations of motion for disordered conductors: extrinsic interband velocity, corrected collision integral and spin-orbit torques

The semiclassical equations of motion are widely used to describe carrier transport in conducting materials. Nevertheless, the substantial challenge of incorporating disorder systematically into the semiclassical model persists, leading to quantitative inaccuracies and occasionally erroneous predictions for the expectation values of physical observables. In the present work we provide a general prescription for reformulating the semiclassical equations of motion for carriers in disordered conductors by taking the quantum mechanical density matrix as the starting point. We focus on external electric fields, without magnetic fields, and spin-independent disorder. The density matrix approach allows averaging over impurity configurations, and the trace of the velocity operator with the disorder-averaged density matrix can be reinterpreted as the semiclassical velocity weighted by the Boltzmann distribution function. Through this rationale the well-known intrinsic group and anomalous velocities are trivially recovered, while we demonstrate the existence of an extrinsic interband velocity, namely a disorder correction to the semiclassical velocity of Bloch electrons, mediated by the interband matrix elements of the Berry connection. A similar correction is present in the non-equilibrium expectation value of the spin operator, contributing to spin-orbit torques. To obtain agreement with diagrammatic approaches the scattering term in the Boltzmann equation is corrected to first order in the electric field, and the Boltzmann equation is solved up to sub-leading order in the disorder potential. Our prescription ensures all vertex corrections present in diagrammatic treatments are taken into account, and to illustrate this we discuss model cases in topological insulators, including the anomalous Hall effect as well as spin-orbit torques.

The central idea of the semiclassical model is the separation between the dynamics of individual carriers and the carrier distribution.Carrier dynamics between collisions are described by the semiclassical equations of motion, which do not incorporate disorder, while collisions are taken into account through the Boltzmann equation, and affect solely the distribution function, in-ducing changes in the occupation of quantum states 2 .The semiclassical velocity, originally assumed to be simply the band group velocity, is now known to incorporate a transversal anomalous component linear in the driving electric field and proportional to the Berry curvature Ω m of a given band m. 1,52 Although written in terms of the curvature for a single band this anomalous velocity includes inter-band coherence effects, and is associated with band mixing by an electric field. 49The anomalous velocity lies at the heart of the quantum Hall effect and of the intrinsic contribution to the anomalous Hall effect, together with its quantized counterpart.In recent years, however, it has been realised that disorder itself leads to band-mixing effects which are not captured by the Boltzmann equation, and are challenging to include in the wave packet description, since averaging over disorder configurations cannot be done at the level of the wave function.It is well established that a naive application of the semiclassical model to the anomalous and spin-Hall effects in disordered systems makes inaccurate predictions 53,54 .Indeed, the role of disorder in the anomalous Hall effect, [55][56][57][58][59] and its relationship to semiclassical dynamics, remains an intensely researched topic 34,38,[60][61][62][63][64][65][66][67] .
Nevertheless, since transport is fundamentally semiclassical, all transport-related quantities must be expressible in semiclassical terms.The assumptions behind wave packet dynamics and diagrammatic approaches are the arXiv:2109.06214v1[cond-mat.mes-hall]13 Sep 2021 same: external fields are treated classically and are assumed to be slowly-varying in space, a separation is made between scattering processes and the dynamics between scattering events, and the calculation is performed in the regime F τ / 1, where F is the Fermi energy and τ the momentum relaxation time.Recent work has investigated strategies for incorporating the findings of diagrammatic linear response theory into semiclassical dynamics.Xiao and Niu showed that agreement is obtained with diagrammatic approaches if all semiclassical quantities are dressed by disorder, the cost being that the introduction of a disorder-dressed Berry curvature 68 .0][71] These quantities, however, are difficult to work with, since a spin-dependent position operator introduces complications of its own.
In light of the above, in this paper we formulate the semiclassical equations of motion for Bloch electrons so as to include disorder using the density matrix formalism, 49 in the process making a connection to Green's functions approaches.4][35][36][37][38][39][40][41][42][43][44][45][46][47] Our primary aim is to provide a straightforward method to incorporate disorder into such computational strategies once a model of disorder is chosen.Taking the density matrix as the starting point 72 allows one to average over disorder configurations, something that cannot be done using a wave function.We demonstrate that disorder affects not only the state occupation but also the semiclassical equations of motion, and that it generates a correction to the velocity that accounts for band mixing mediated by the Berry connection and disorder.This approach enables one to distinguish disorder effects on the distribution function from disorder effects on carrier dynamics, yet entails a change in one's point of view so as to regard the semiclassical equations as describing carrier propagation averaged over many disorder scattering events.The carrier undergoes transitions between bands as it scatters, and its trajectory can be determined by averaging over impurity configurations.Whereas the equation of motion for the wave vector k follows trivially from the time derivative of the momentum operator, our central result is the revised semiclassical equation of motion for the position of a carrier in band m, with dispersion m , propagating under the action of an electric field E in the presence of disorder We identify a new contribution to the velocity, which we term the extrinsic interband velocity β m k , defined as: where R represents the inter-band Berry connection with matrix elements R mn The extrinsic interband velocity β m k is proportional to the disorder strength, which is typically quantified by the impurity density n i , scattering potential strength u 2 0 , or alternatively 1/τ , where τ is the characteristic scattering time.Formally, β m k is similar to the customary scattering term in the Born approximation, except the distribution function is replaced by the band off-diagonal elements of the Berry connection.Since only the band off-diagonal elements of the Berry connection appear β m k is by construction gauge covariant.In general β m k , being independent of applied fields, can be thought of as a disorderdependent correction to the semiclassical band velocity, or a random inter-band walk on the Fermi surface.We find that β m k is nonzero in systems in which time reversal symmetry is broken by e.g. a magnetization.Whereas β m k is similar to the side jump as defined in Ref. [73], unlike Ref. 73, the present formalism does not employ coordinate shifts, so that the formal and physical position operators coincide.Furthermore, unlike Ref. [68], the Berry curvature is the same as in the clean system, rather than being dressed by disorder.Importantly, we show that the scattering term in the Boltzmann equation, needed to determine the effective distribution function, acquires a correction to first order in the electric field, which is equivalent to a gradient expansion in the electrostatic potential.In addition, the Boltzmann equation needs to be solved up to the sub-leading order in the impurity strength, in order to incorporate processes customarily termed skew scattering and side jump. 57,58,74he method we present here also enables us to calculate spin densities using the semiclassical model and obtain accurate results for spin-orbit torques.In order to accomplish this the bare spin expectation value needs to be supplemented with an electric field contribution, and we find an analogous quantity to β m k in the spin expectation value.
More generally, we present a prescription for mapping steady-state expectation values onto the semiclassical model by expressing traces purely in terms of the band diagonal elements of the density matrix.Since in linear response theory all expectation values are traced back to the equilibrium density matrix, which is banddiagonal, they can all be recast in terms of semiclassical quantities.The band-diagonal elements of the density matrix represent the Boltzmann distribution, which can be evaluated from the much simpler Boltzmann equation.In fact we will argue briefly in the latter part of this work that linear response theories can be thought of as a family tree with its roots in the quantum Liouville equation: the Kubo approach is the integral formulation, the quantum kinetic, or quantum Boltzmann, approach is the integro-differential formulation, while the semiclassical model is an offshoot of the latter, which arises as a result of an additional separation between the carrier dynamics and distribution.Sharing a common origin, these methods yield equivalent results, and, in particular, vertex corrections present in diagrammatic approaches have straightforward equivalents in the semiclassical language.The blueprint presented in this work can be used in the future to incorporate electron-electron interactions into computational approaches in a mean field picture.
The outline of this paper is as follows.In Sec.(II) we introduce the Hamiltonian and the model of disorder.In Sec.(III) we review linear response theory based on the density matrix and introduce the electric-field correction to the collision term.Next in Sec.(IV) we outline the general methodology for deriving the semiclassical equations of motion from the quantum kinetic equation, and discuss also disorder effects on spin expectation values.In Sec.(V) we analysed the relation among different linear response methodologies commonly used to calculate transport coefficients.In Sec.(VI) we discuss at length two model examples, the anomalous Hall effect and spinorbit torques in magnetic topological insulators.We end with a summary and conclusions.

II. MODEL HAMILTONIAN
We consider a Hamiltonian of the form: with H 0 the low energy effective band Hamiltonian, in principle assumed to include the Zeeman interaction with an external magnetic field, V (r) is the electrostatic potential, and U = U (r) represents the disorder scattering potential.We emphasize the Hamiltonian is always Hermitian.Non-Hermitian systems were considered in Ref. [75].
We work in the crystal momentum representation |m, k = e ik•r |u m k .The matrix elements of a scalar disorder potential U (r) are given by the equation where we have defined the Fourier transform of the spatial function in d dimensions with q = k − k .The impurity average is defined by with where • • • refers to an average over impurity configurations.For concreteness we will use a model of disorder whose spatial correlations function is defined as: Then, it follows that where u 2 0 is a parameter that takes into account the strength of the disorder potential.

III. QUANTUM KINETIC EQUATION
In this section we give a brief presentation of the quantum kinetic equation in a somewhat different language than that used in Ref. [49].We note that similar densitymatrix based approaches have been used recently to describe carrier dynamics in the semiclassical regime 50,51 .The starting point is the quantum Liouville equation for the single particle density operator ρ, namely: For the sake of convenience, we introduce at this stage the free retarded Green function In the frequency domain where we introduced the factor e −ηt to ensure convergence.The advanced Green function follows by Hermitian conjugation.

A. Kinetic equation in equilibrium
For the sake of simplicity, let us for now ignore the effect of the driving electric field in the kinetic equation.Using a decomposition of the density matrix as ρ = ρ + g 0 in the quantum Liouville equation 49 , we get for the disorder averaged part in equilibrium: while for g 0 we get the equation: In order to solve the kinetic equation for ρ , we first solve Eq.( 15) for g 0 and then we use it in Eq.( 14).In the first Born approximation 49 we neglect the last two terms on the left hand side of eq.( 15).We are left with Solving for g 0 In terms of Green's functions g 0 can be expressed as This solution is substituted into Eq.( 14).We arrive at the equation with the collision integral J( ρ ) defined as:

B. Adding an electric field
Let us now consider the effect of the driving electrostatic potential up to linear order.For simplicity we take this potential to have the form V (r) = eE • r, implying a uniform electric field, which corresponds to the overwhelming majority of experimental setups.The case of inhomogeneous systems, including systems in inhomogeneous electric fields, entails additional subtleties which we postpone for later consideration [76][77][78][79] .Adding an electric field to the Hamiltonian implies a correction to the function g, which can then be written as g = g 0 + g E , where g 0 was found in the previous section, and The notation g E reflects the fact that eventually it is the electric field that appears in the final expressions, rather than the electrostatic potential.For g E we find explicitly The function g E is off-diagonal in the momentum as well as in the band index.We solve Eq.( 22) by introducing Markovian approximation which reads g 0 (t − t ) ≈ g 0 (t), whereupon in the frequency domain we obtain and in the commutator, we should use Eq.( 18) as a functional of the equilibrium distribution function f 0 ( ) to fulfill linear response.The kinetic equation for the disorder averaged density matrix ρ is now modified to with the collision integral J E ( ρ ) defined as: As we show below, the electric field correction to the collision integral in Eq. ( 25) with the off-diagonal density function as given in Eq. ( 23) will provide results in agreement with previous calculations based on diagrammatic perturbation theory 62,80 .We consider such an agreement as a positive test of the Markovian approximation.We note also that Eq. ( 23) was used in a different but equivalent form in a previous paper 81 in order to calculate side jump effects in a system with extrinsic spin-orbit coupling.In this paper, we will focus on systems with intrinsic spin-obit coupling.

C. Kinetic equation and linear response
When Eq. ( 24) is expressed in the crystal momentum representation we obtain the quantum kinetic equation 82 We have written the matrix elements of ρ in this representation as f k .We refer to f k henceforth as the density matrix, noting that it has matrix elements connecting different bands, although the band index n has not been written explicitly.The covariant derivative To solve Eq. 26, the density matrix is separated into a band diagonal and a band off-diagonal part, namely, we write f k = n k + S k .The band diagonal term n k represents the fraction of carriers in a specific band and is essentially the solution of the ordinary Boltzmann equation, while S k contains the effect of inter-band coherence, or band mixing.All our effort in recovering the semiclassical theory consists of eliminating S k .The effective Boltzmann equation that we shall derive is simply what is obtained for n k once all references to S k have been eliminated.Fortunately, as we recapitulate below, the solution for S k in an electric field is relatively simple, making it straightforward to express expectation values in terms of n k alone.
The equilibrium density matrix is band diagonal, its elements represented by the Fermi-Dirac distribution for each band n F D ( m k ).In an electric field one may expand to linear order f mn k = n F D ( m k )δ mn + f mn Ek , with corresponding expressions for n Ek and S Ek .The kinetic equation is split into two coupled equations for n Ek and S Ek , whose solution, based on an expansion in the small parameter /( F τ ), is explained in detail in Ref. [49].It was shown that n Ek starts at order −1 in this small parameter, since it is proportional to the scattering time τ , while S Ek starts at order 0. Consequently, the subleading correction to n Ek , referred to as n (0) Ek , is also required.
To leading order in /( F τ ), the diagonal part reads: where the Born approximation collision integral is The solution of Eq. ( 27) is in general rather complicated 83 .For a system with isotropic dispersion it reduces to the simple form where the transport time τ m p is defined as: (30) The solution for S (0) Ek takes the simple form 49 with the intrinsic and anomalous driving terms, 49 Since S Ek starts at zeroth order in the parameter /( F τ ), we also require the sub-leading term n (0) Ek , which is found from the equation where the right hand side acts as the driving term, whose constituents will be explained shortly.Solving this equation will yield two different contributions to the subleading diagonal density matrix n Ek , which we write as n Ek are of zeroth order in /( F τ ), they are parametrically different with respect to magnetisation and Fermi energy, as we will see later on.We can solve for these two terms separately as follows.

The contribution n (sk)
Ek stems from D and is associated with skew scattering in the semiclassical theory.It is solved in an analogous manner to Eq. ( 27), namely, the driving term is found by substituting Eq. ( 31) into a collision integral of the form of Eq. ( 28), obtaining which can be solved for n (sk) Ek using the standard techniques of Boltzmann theory 83 .The driving term in this equation can be written explicitly as a function of the leading-order density matrix n and we recall that n Ek was found in Eq. ( 29).The second contribution to the driving term in Eq. ( 34) is due to the electric field correction of the collision inte-gral J E acting on the equilibrium distribution function.Since this contribution is associated with side jump scattering in the semiclassical theory, we will refer to it as Ek we need to solve the equation In the crystal momentum representation the electric field correction to the collision integral takes the form where the derivatives act only on U mm kk .This equation should be compared with the side jump velocity calculated from a coordinate shift introduced in Ref. [73].The balance between two collision integrals, as stated in Eq. ( 37) provides the necessary information to calculate a new subleading density function n (sj) E that in the semiclassical language 84 is interpreted as an anomalous distribution due to coordinate shift of the scattered particle after many collisions.

IV. RECOVERING THE SEMICLASSICAL THEORY
In this section we decompose the kinetic equation into a part representing carrier dynamics and a part representing the distribution, which is found from a modified Boltzmann equation.Since the equation of motion for the carrier wave vector, yielding k = −eE, follows immediately from the operator commutator [p, V (r)], the bulk of our effort is devoted to finding the disorderaveraged velocity, which will yield the time evolution of the carrier position ṙn .The prescription for recovering the semiclassical theory from the quantum kinetic equation proceeds as follows: • Determine the velocity expectation value as the op-erator trace Tr ( ṙf ), where ṙ = i [H, r] represents the matrix elements of the velocity operator.In the crystal momentum representation these are given by the covariant derivative ṙ = 1 DH Dk .
• Reduce the trace to a form in which only banddiagonal elements of the density matrix appear.These will contain either the equilibrium Fermi-Dirac distribution n 0k , or the correction to the band-diagonal part n Ek , which we recall has three constituents: n Ek = n Ek .
• The result follows a natural separation into a contribution associated with the equation of motion ṙn and one associated with the Boltzmann equation.
• For the spin density, we follow similar steps, namely, we take the trace of the spin operator in the Bloch basis with the averaged density matrix.We will also find an extrinsic spin matrix element that accounts for spin rotations during scattering events.
The Hamiltonian is H = H 0 + V (r) + U (r), and since the last two terms commute with the position operator they do not contribute to the velocity operator.The band Hamiltonian yields The Berry connection R mm The first term gives the usual group velocity v m k = ∇ k m k / which is diagonal, while the second term gives a contribution due to band mixing and is purely off diagonal.We will concentrate on the second factor or band mixing velocity.The off-diagonal density matrix is composed of two terms: an intrinsic one and an extrinsic one.Let us first consider the intrinsic one.It is After replacing the driving term D m m Ek by exchanging m → m in the first term, and summing over intermediate states, the intrinsic contribution can be written as the average of the transverse velocity with the Berry curvature The extrinsic contribution reads After some algebra it can be written as with the function β m k formally defined as where the prime in R indicates that only the band offdiagonal matrix elements of the Berry connection enter.In Eq. ( 45) we have written directly the electric-field dependent correction to the distribution function, since Eq. ( 44) makes it obvious that this contribution vanishes when f is replaced by the equilibrium distribution n F D .This is because, for scalar scattering as studied in this work, the equilibrium distribution causes the entire collision integral to vanish.For computational evaluations it will be useful to list the explicit equation for β m k : Note that β m k is proportional to the disorder strength quantified here by u 2 0 , making it first order in /( F τ ).It represents a disorder-dependent correction to the semiclassical band velocity, which is independent of the applied electric field.Physically, β m k represents the average value of the random changes in the carrier velocity that occur every time the carrier is scattered between bands.Since β m k has units of velocity and depends on the disorder potential we will refer to it as the extrinsic interband velocity.Given that β m k is formally of first order in /( F τ ), we are only interested in its product with the leading term in the distribution function, n Ek , so that its overall contribution to the current is formally zeroth order in disorder.Moreover, with n Ek representing a Fermi surface contribution, the net effect of β m k can be thought of as a random inter-band walk on the Fermi surface.Interestingly, β n has the same mathematical form as the Born approximation scattering term J 0 , except the band off-diagonal elements of R appear instead of n k .The presence of only the band off-diagonal matrix elements of R ensures β m k is gauge covariant.In the examples we study below we find that β m k is nonzero in systems in which time reversal symmetry is broken by e.g. a magnetization.It is similar to the side jump appearing in Ref. [73], although we stress that our approach makes no reference to any coordinate shifts, and the formal position operator is identical to the physical position operator.
Since all contributions to the current density are now expressed in terms of the distribution function (the equilibrium as well as the leading and sub-leading terms in an electric field), we are able to write the semiclassical equations of motion as The distribution function is found from the Boltzmann equation, with the caveat that we require both the leading and subleading order terms in the disorder strength.The procedure is as follows.First the leading-order term in the distribution function n while the sub-leading correction n (0) E is given by where the left hand side is the quantity to be found, and the right hand side plays the role of a driving term.
Finally, we are able to write the full expectation value of the current in terms of semiclassical quantities Let us consider the expectation value of the spin operator in the presence of an electric field where s nm k represent the matrix elements of the spin operator.Writing explicitly the off diagonal terms of the density matrix in the average of the spin operator we can separate the intrinsic and extrinsic contributions as If we define the quantity the intrinsic contribution can be rewritten as where we only have to take off-diagonal components inside the commutator.This enables us to define an intrinsic spin expectation expectation value as The extrinsic part can be written as This is mathematically analogous to the expression for the extrinsic interband velocity Eq. ( 44), and using a similar manipulation we can re-express it as where n Ek is the leading order distribution function and we introduce the new extrinsic spin γ m k given by In exact analogy with β m k , since the Born approximation scattering term vanishes when the distribution function is replaced by the Fermi-Dirac distribution, γ m k also vanishes in equilibrium in the presence of scalar scattering.This quantity represents spin rotations during scattering events.Again, the quantity T only has inter-band matrix elements.Explicitly γ m k is evaluated as In analogy to the above, we can write a modified intraband spin matrix element as The first factor is the bare matrix element of the spin operator in band m, while the remaining two represent non-equilibrium corrections, the first being intrinsic and the second extrinsic.We note that γ k and β k are mathematically very similar and both contribute at the Fermi energy.The intrinsic contribution to the spin density from the Fermi sea appears in the second term in Eq. ( 62).The full semiclassical expression for the spin density is given by V. LINEAR RESPONSE FAMILY TREE In closing the methodological discussion we remark briefly on the relationships between the various linear response theories.The most common strategy for solving Eq. ( 11) in an electric field is via Kubo linear response theory, which is discussed in detail in many textbooks 85 , hence we only dwell upon its fundamental aspects.Briefly, the Hamiltonian is decomposed as H = (H 0 + U ) + V , and the non-equilibrium part of the density matrix is likewise singled out as ρ = ρ 0 + ρ E .Then in linear response one can write This equation is solved immediately to yield where G R , G A are the retarded and advanced Green's functions, respectively, for the disordered system: We refrain from writing out the energy dependence in full.Note that the Green's functions have not been averaged over disorder configurations at this stage.To obtain the customary Kubo formula one must trace over the velocity operator and average over impurity configurations The procedure is standard, so we do not cover it here in detail, the purpose of this description is illustrative.The important point to notice is that the Kubo formula is the integral approach to solving the Liouville equation, whereas the kinetic equation we follow in this work represents the differential approach, or integro-differential in view of the complex scattering term.The Keldysh theory follows a similar path, and although it takes as its starting point a series of Green's functions, its ultimate origin lies in the quantum Liouville equation.The Keldysh theory is formally non-local in time, although in the vast majority of practical applications the non-locality is removed and the Keldysh Green's function, which is analogous to the density matrix employed in this work, depends only on the difference in time variables.Thus,  for the purposes of the present comparison, the quantum Boltzmann equation derived in the Keldysh theory is indistinguishable from the quantum kinetic equation derived from the density matrix.The Kubo formula, Keldysh theory and quantum kinetic equation may be regarded as holistic approaches, in which both the carrier dynamics and the carrier distribution are accounted for in the density matrix (or Keldysh Green's function), and the net result is the expectation value of a physical observable.In contrast, the semiclassical theory involves a separation between the carrier dynamics and carrier distribution, which can help to build an intuitive picture of the underlying physics.The relationship between the different approaches is summarized in the family tree of Fig. 1.The most important observation in this context is the common origin of all linear response theories, which reinforces the expectation that they should all lead to the same results.

VI. APPLICATIONS
We now turn to applications of the theory, which are intended to illustrate the way the extrinsic interband velocity, extrinsic spin terms, and additional scattering terms in the Boltzmann equation appear in the explicit evaluations of physical observables for model systems.
In particular, we emphasize the relationship between these various contributions and the analogous quantities appearing in diagrammatic theories, which enables us to reconcile the semiclassical and diagrammatic results.
Our focus will be on topological insulators 86 , where we discuss the anomalous Hall effect as well as spin-orbit torques, and compare the semiclassical results with previous work.

A. Anomalous Hall effect in topological insulators
In this section we calculate the anomalous Hall conductivity in topological insulators.The anomalous Hall conductivity is basically expressed in terms of four contribution: intrinsic contribution that takes into account the whole Fermi sea of the system, extrinsic contribution due to the extrinsic velocity β m k at the Fermi surface, side jump like contribution at the Fermi energy due to an electric field correction to the collision integral and a skew scattering contribution.
The Hamiltonian that describes low energy excitations in the surface of 3D topological insulators reads: where v F is the effective Fermi velocity and σ i are Pauli matrices, while M is the magnetization.The eigenvalues are , where ± labels conduction/valence and the eigenstates are with a parameter From the eigenstates we determine the Berry connection vector R mm We decompose this vector into a diagonal and an off-diagonal contribution, namely, We have defined a unit vector parallel to momentum k = (cos θ k , sin θ k ) and a unit vector perpendicular to momentum θ = (− sin θ k , cos θ k ).Also, these vectors are related as θ = ẑ × k.
Let us start by calculating the intrinsic contribution to the Hall conductivity.It reads: Explicit evaluation for a topological insulator gives the Berry curvature After a straightforward integration we find the intrinsic contribution to the Hall conductivity to be The extrinsic inter-band velocity makes the following contribution to the anomalous Hall conductivity The leading order correction to the distribution function is n with the transport time given by the expression with the scattering time defined as 1/τ = πu 2 0 ρ( k )/ and the density of states ρ The extrinsic inter-band velocity is a transverse velocity since it is proportional to the unit vector θ.After explicit integration we arrive at the expression Notice that n (−1)m Ek is inversely proportional to the impurity density while the extrinsic velocity β is proportional to the impurity density.As a result, the overall effect in Eq. ( 75) is independent of disorder.The extrinsic inter-band velocity β m k comprises the effect of disorder on carrier dynamics, namely, it can be interpreted as an effective velocity of the electron after many collisions, in contrast to the group velocity, which is a velocity between collisions.In this sense our extrinsic velocity can be associated to the side jump velocity encountered in previous semiclassical results 84 but with the difference that β m k is entirely due to band mixing mediated by the off-diagonal components of the Berry connection vector and that it is constructed from a collision integral without introducing any quantity related to coordinate shift.
As we discussed earlier, there are also two contributions to the anomalous Hall conductivity related to two different diagonal subleading density matrix functions.Let us first calculate the anomalous Hall conductivity related to the term n (sj) Ek in the distribution function.It reads The diagonal velocity reads while the correction to the distribution function is Replacing all elements we find for the conductivity This term doubles the contribution in Eq.( 78) due to the extrinsic velocity although in this case the effect of disorder is completely captured by n (sj) Ek .In previous semiclassical studies an anomalous distribution function was introduced as a result of a coordinate shift 84 .In contrast, we have derived n (sj) Ek from Eq. ( 37) without any need for introducing a coordinate shift.
The contribution to the anomalous Hall conductivity related to the skew scattering correction to the distribution function, n (sk) Ek , reads The diagonal velocities are After integration we obtain Adding all the contributions to the Hall conductivity we get the final expression in exact agreement with previous results using the noncrossing approximation and diagrammatic perturbation theory 62,70,71,74,84 .

B. Spin density and spin-orbit torques in topological insulators
In this section we determine the spin density and spinorbit torques in topological insulators with an out-ofplane magnetization described by the effective Hamiltonian Eq. (68).As for the conductivity, the spin density has five contributions: a dominant contribution from the Fermi surface leading to the Edelstein effect 87 , an intrinsic contribution from the Fermi sea, a contribution due to the extrinsic spin expectation value γ m k at the Fermi surface, a side-jump like contribution at the Fermi energy due to the electric field correction to the collision integral and a skew scattering contribution.
The leading order contribution to the spin density is Using n the Edelstein effect contribution to the spin density is The fraction of the spin density weighted by the intrinsic driving term reads: The dot product is between the electric field and the Berry connection vector.With the Berry connection and the off-diagonal spin expectation values we get the intrinsic correction to the spin density The extrinsic correction is defined as yielding This extrinsic contribution to the spin density is the counterpart of the extrinsic inter-band velocity contribution in the anomalous Hall effect.The extrinsic spin γ m k as defined in Eq. ( 60) is an interband coherence effect mediated by an effective off-diagonal spin operator defined in Eq. ( 56).
The side-jump contribution is given by Using the spin expectations values above and the density function given by Eq.( 81) we find that the side-jump contribution to the spin density reads This term doubles the contribution of Eq. (99).This is the counterpart of the anomalous distribution function introduced in the semiclassical theory 84 and also calculated in Eq. ( 82) for the anomalous Hall conductivity.The skew scattering contribution takes the form Using the diagonal spin expectations values and the distribution function found in Eq.( 84) we find for the skew scattering contribution to the spin density This is the counterpart of the skew scattering term in the anomalous Hall conductivity calculated in Eq. (85).Finally, the total spin density of the system then reads The spin-orbit torque is defined as τ = (2M/ )m × s , where m is a unit vector in the direction of magnetisation that we took here as m = ẑ.Restoring the factor of /2 in the spin matrix elements and substituting the density of states ρ( F ) = F /2π 2 v 2 F , the spin-orbit torque is finally given by in agreement with previous results. 80,88A peculiarity of topological insulators is that the velocity operator is directly related to the spin operator as v = −v F ẑ × σ.
Then the current density is j = −e v what implies that the anomalous Hall conductivity and the spin density are related by σ yx = J y /E x = ev F s x .This fact can indeed be verified from Eq. ( 86) and Eq.(104).

VII. CONCLUSIONS AND OUTLOOK
We have demonstrated that the semiclassical dynamics of electrons in disordered solids can be determined using linear response theory by taking the density matrix and quantum Liouville equation as the starting point.This results in a disorder-dependent correction to the semiclassical equation of motion for the carrier position, which we have termed the extrinsic inter-band velocity β m k , and which accounts for the effect of disorder on the carrier velocity after many collisions.In analogy to the extrinsic inter-band velocity, an extrinsic correction to the spin expectation value γ m k is also present, which accounts for spin rotations during scattering events.This is accompanied by an intrinsic, electric field dependent non-equilibrium correction to the spin expectation value, which is analogous to the anomalous velocity present in the semiclassical equations of motion.At the same time, the Boltzmann equation must be solved up to sub-leading order in the disorder strength.The collision integral in the Boltzmann equation includes an electric field dependent correction which is analogous to the side-jump scattering term, as well as an additional correction analogous to the skew scattering term in systems with intrinsic spin-orbit interactions.We have applied this theory to describe the anomalous Hall effect and spin-orbit torques in topological insulators, obtaining exact agreement with quantum mechanical results using the Kubo formula.
The general prescription we have formulated in this work can be straightforwardly generalised to include extrinsic spin-orbit and magnetic impurity scattering, as well as magnetic fields, which, however, require more effort since the description relies on the Wigner function.Our prescription paves the way towards a systematic semiclassical picture encompassing a host of effects that conventionally lie beyond the purview of semiclassical theory: disorder effects beyond the Born approximation, such as weak localization, electron-electron interactions in the mean-field approximation, as well as Kondo physics in magnetic systems.
and |u m k is the lattice-periodic part of the Bloch wave function.

Figure 1 :
Figure 1: Family tree of linear response theories.

k 1 / 2
cos θ k and the extrinsic inter-band velocity in the conduction band takes the form