Short-scale quantum kinetic theory including spin-orbit interactions

We present a quantum kinetic theory for spin-$1/2$ particles, including the spin-orbit interaction, retaining particle dispersive effects to all orders in $\hbar$, based on a gauge-invariant Wigner transformation. Compared to previous works, the spin-orbit interaction leads to a new term in the kinetic equation, containing both the electric and magnetic fields. Like other models with spin-orbit interactions, our model features"hidden momentum". As an example application, we calculate the dispersion relation for linear electrostatic waves in a magnetized plasma, and electromagnetic waves in a unmagnetized plasma. In the former case, we compare the Landau damping due to spin-orbit interactions to that due to the free current. We also discuss our model in relation to previously published works.


I. INTRODUCTION
Dense plasmas, where quantum effects are important, can, for example, be found in different types of solid-state plasmas [1,2], various astrophysical environments [3,4], and certain forms of laser-plasma interactions [5,6]. Such plasmas have been modelled using hydrodynamic [7] and kinetic equations [8], where, in the present paper, the focus is on the latter class of theories. Quantum kinetic equations can be derived, for example, from the Green's function formalism [9,10] or from the density matrix approach [11,12]. Recently, several kinetic models for plasmas have been put forward, within the Hartree approximation [12][13][14][15] or the Hartree-Fock approximation [16][17][18]. Depending on the scope of the theory, the density matrix can be based on the Schrödinger equation [8], the Pauli equation [13], or the Dirac equation, where the latter can be applied in the weakly [19] or fully [12] relativistic approximation.
In order to produce a quantum kinetic theory of electrons from the Dirac equation, certain restrictions need to be applied. Firstly, pair-production must be negligible, such that a Foldy-Wouthuysen transformation [20,21] can be applied to separate electron states from positron states. This puts restrictions on the maximal electric field strength, which should be well below the critical field E cr = m 2 c 3 /|q| , and the characteristic spatial scales of the fields, which should be much longer than the Compton length L c = /mc. Here m and q are the electron mass and charge, respectively, c is the speed of light in vacuum, and is the reduced Planck's constant. Given these restrictions, there are two different regimes that can be studied. Firstly, there is the fully relativistic regime, where the relativistic factor γ fulfills γ − 1 ∼ 1, in which case the de Broglie length λ dB = /p ch is of the same order as L c . Here p ch is the characteristic momentum of the electrons -determined by the Fermi temperature T F or the thermodynamic temperature T whichever is larger. Since, according to the applicability conditions, this means that the spatial scales must be longer than the de Broglie length, particle dispersive effects cannot be included in the fully relativistic regime [12]. However, for the second case, the weakly relativistic regime with γ − 1 1, we have λ dB L c , and as a consequence, a theory including particle dispersive effects based on the Dirac equation can be formulated. This is the main goal of the present paper. We thus generalize the work presented in Ref. [13] based on the Pauli equation, by including several new effects, such as spin-orbit interaction, Thomas precession and the polarization currents associated with the spin. We also generalize the results of Ref. [19], by also covering the short scale physics down to spatial scales of the order of the de Broglie length.
Our approach is as follows. We start from the Pauli Hamiltonian including the spin-orbit interaction. In Section II, we apply a combined Wigner-and Q-transform (for the spin) using the gauge-independent approach of Stratonovich [22,23], to formulate our scalar kinetic theory.
In Section III, in order to demonstrate the usefulness of the present theory, we have calculated two examples from linearized theory: electromagnetic waves in a nonmagnetized plasma and electrostatic waves propagating parallel to an external magnetic field. Classically, an external magnetic field does not affect parallel propagating electrostatic waves. By contrast, the present model predicts a condition, depending on the magnetic field, for resonant wave-particle interaction, generalizing those of previous works. We use this result to calculate the magnetic field dependence of the damping rate.
Finally, in Section IV, we compare our model with some theories from the recent research literature.

II. THE GAUGE INVARIANT WIGNER FUNCTION AND THE SPIN TRANSFORMATION
We will consider the semi-relativistic Pauli Hamiltonian +qφ−µ e σ ·B+ µ e 2mc 2 π ×Ê −Ê ×π ·σ (1) including the spin-orbit interaction. Hereπ =p − qÂ is the gauge invariant momentum operator,φ,Â are the scalar and vector potentials and µ e is the electron magnetic moment. Note that the spin-orbit interaction is written in symmetric form so as to make the Hamiltonian Hermitian -which is necessary to have unitary time evolution.
Our goal is to formulate a scalar quantum kinetic theory based on the gauge invariant Wigner function [22,23], that is, where α, β are the Pauli spin indices, and this trace is only over the spatial degrees of freedom. We work in the Heisenberg picture so that the operatorŴ (x, p, t) is time-dependent, but the density matrix ρ is not. The op-eratorŴ (p, x, t) can be expressed as a Fourier transform where the operatorT (u, v, t) is given bŷ To obtain a scalar function from the matrix-valued W αβ , one can use the spin transform which is a Husimi Q-function for the spin [24]. This transform was used in Ref. [13] for a theory including only the magnetic dipole interaction. The remainder of this section presents the calculations involved in deriving an evolution equation for f , with the conclusion presented in Section II C. The calculations in this section are similar to those in Ref. [23], where the spin-independent terms in the kinetic equation are found.
Our notation is as follows. Throughout, quantities with hats are operators and quantities without hats are c-numbers, except the Pauli spin operator σ which is always an operator (a finite dimensional one, i.e., a matrix). We will also use the summation convention for repeated indices, and use commas to denote derivatives, e.g., Various other forms of the operatorT (Baker-Campbell-Haussdorff formulas) can be found [23] and will be useful in the following: One can note that these occur in pairs related by operator orderings.
The time evolution ofŴ is given by the Heisenberg equation of motion, Here the only effect of the partial time derivative is to give the ∂ t A contribution to the electric field in the Lorentz force. In the Heisenberg picture, the operator σ is timedependent, so we must remember to include dσ/dt when deriving the evolution equation for f . Hence, if we let i.e. f is the expectation value ofŜŴ . Note thatŜ and W are both Hermitian and commute, which implies that f is real. The Heisenberg equation of motion then gives where in this and the previous equation, the trace is over both the spin indices and the spatial degrees of freedom.
The commutator is linear inĤ, and the terms in the evolution equation forŴ corresponding to the lowest order ("non-relativistic") part of the Hamiltonian have already been computed in [13]. Therefore we will here be concerned only with the spin-orbit contribution We will the denote the contribution from this Hamiltonian by (∂ t f ) SO . The spin-orbit interaction Hamiltonian has the form σ ·V, and hence we find Then using σ i σ j = δ ij + i ijk , we get For the second commutator, we writeVŴ as the sum of its symmetric and anti-symmetric parts and clearly, the second part cancels with cross product part ofŜ[Ŵ ,Ĥ SO ] in Eq. (12). Now we just need to use which implies and thus The remaining, laborious, task is to evaluate and take the Fourier transform ofT (π ×Ê) ± (π ×Ê)T , then express the result in terms of differential operators -functions of ∂ p and ∂ x -acting onŴ . These operators can be taken outside the trace, so that they act on f , and this will give the kinetic equation. Since we need to evaluate the commutators withπ andÊ. This is sufficient, because for any operatorÔ sinceŴ is Hermitian, and we get the commutator with the Hermitian ordering (which appears in Eq. (9)) by reading off the "imaginary part" of [Ŵ ,π ×Ê].

A. The commutator [Ŵ ,Ê]
In Sections II A to II C, we will typically omit boldface for vector quantities; the type of each quantity should be apparent from context. Using Eq. (4b), we can write the operatorT witĥ Θ(u) = exp( i u ·p) on the right. Sincep is the generator of translations, the operator Θ(u) acts likê and thus, using Eq. (4b) to put the translation operator on the left, with the other factors commuting withÊ, Now we need an expression for the productπ × [T ,Ê] in terms ofT . By putting the commutator in the form Eq. (22), we can use any expression forT as is most appropriate. The one that is most appropriate here is Eq. (4e) since we can lift the result from Ref. [23] their (4.30b), and its Hermitian conjugate: We thus have that the operator π i [T , E j ] can be expressed as Now let |x , |x be eigenstates ofx (at the appropriate time). Use Eq. (4b) so that thex operators are on the left At this point we see that the Fourier transform in v can be carried out and will give an expression proportional to δ(x − x + u/2). Hence the matrix elements ofQ B ij (x, p) remain unchanged if we change x to x − u/2 in the argu-ments of the fields. In the integral over τ , we also change variables to τ − 1/2 for a somewhat simpler expression.
We can then restore the operators making upT , resulting in Since the Fourier transform sends u → i ∂ p , we can now express this operator in terms ofŴ and its derivatives. First, use that and note that this quantity is real. In the integral with B, 1/2 is even and −τ is odd, so even and odd powers of iτ ∂ p in the expansion of B survive the integration, respectively. Since we want the imaginary part, it is the term with −τ that we should keep, and in conclusion We now turn toÔ D ij . Note that from the expression forT in Eq. (4e) we have and henceÔ The first term can be handled using the argument above, except that the Fourier transform in v will send v i to i ∂ xi δ(x). The result is that but this quantity is real, and so makes no contribution to the evolution equation. For the term with the u-derivative, we integrate by parts in the Fourier transform, i.e., In the first term, p i can be taken outside the Fourier transform, and the remaining operator is treated like above resulting in We can write this in a somewhat simpler form as where we have introduced the notationẼ To summarize, we obtain The quantitiesẼ and ∆p are the same as in Refs. [13,23]. Since Eq. (38) is to be contracted with the Pauli matrices, in the kinetic equation it represents the spin-orbit force.

B. The commutator [Ŵ ,π]
We use Eq. (4e) forT , so that also using (4.25) from Ref. [23]. We can now look at the matrix elements of F Ô A ij , but it is of the form in Eq. (33) except we have only oneÊ operator,Ê(x). Therefore, We see that the imaginary part of E i,j term cancels the imaginary part of the corresponding term in Eq. (35), since when dividing by i, even powers of i ∂ p should be kept. For the second term, we can integrate by parts to show that Again keeping even powers of i ∂ p since we are dividing by i, This is the "hidden momentum" velocity correction term ∝ E × s [25] (to be discussed more below), but with E =Ẽ + ∆Ẽ where ∆Ẽ is higher order in .
Following steps like those above, and using Eq. (42) again, one readily expresses the matrix elements ofÔ C ij in terms of those ofŴ . The result is that which, in the kinetic equation, will give i.e., it is the correction to the magnetic force due to the relation between momentum and velocity. HereB is defined in the same way asẼ, Eq. (37).

C. The spin transform
It remains to evaluate the terms in the kinetic equation arising from the spin transform, proportional to where H. c. denotes the Hermitian conjugate. Using Eq. (23) again we establish that where D k is the operator Hence the product (π ×Ê −Ê ×π) jT is given by where so that Putting this into Eq. (49), all terms are similar to ones already seen, and we will provide only an outline of the remaining calculations. The operator of the typeÊ l (x) ∂T ∂u k we can handle as in Eq. (34), resulting in and following the argument leading to Eq. (43), when taking the real part, we will get For the term proportional to v k in Eq. (49), we use that F : v k → i ∂ x k δ(x) as leading up to Eq. (33). When taking the real part, the term where the derivative acts on E cancels with the derivative term from Eq. (52) and what remains is which can be simplified using Eq. (36). Following steps like those forÔ B (Eq. (25) and following), the final term in Eq. (49), proportional to I k , will give us a contribution of the form When taking the real part of the Fourier transform of this, there will be two terms, which is in line with previous generalisations, and which is of a new type. However, both are order 2 , so their limits could not have been found in previous works that were either or didn't include the spin-orbit interaction.
Thus, the spin torque is modified to add a term proportional to (p + ∆p) × (Ẽ + ∆Ẽ) . The overall sign and coefficient of this term is found by matching its long-scale limit to the model in Ref. [19] D. Conclusion The evolution equation for f follows directly by taking the expectation value of the evolution equation forŴ . It is For reference, the notation used is ∆Ẽ(x) = i withB, ∆B defined in the same way asẼ, ∆Ẽ. Functions of operators are defined in terms of their Taylor series.
There are five new terms in Eq. (58) as compared to Ref. [13]. In order, they represent: the "hidden momentum" (see below) correction to the velocity in the diffusion and magnetic force terms, respectively; the spinorbit force; the spin torque due to spin-orbit interaction, including the Thomas precession; and the last term, which is non-linear in the fields and higher-order in . This term lacks an analog in Ref. [19]; the others are short-scale generalizations of terms found in Ref [19] similar to how Ref. [23] generalizes the Vlasov equation.
We have denoted the electron magnetic moment by µ e to allow for an anomalous magnetic moment, i.e., µ e = g 2 µ B where µ B is the Bohr magneton and g = 2 + α π + . . .. That g = 2 leads to a mismatch between the cyclotron frequency and Larmor frequency, resulting in non-classical resonances [26], as will be seen below.
The system is closed with Maxwell's equations with polarization and magnetization, The free charge ρ f is given by introducing dΩ = d 3 p d 2 s. Finding the continuity equation for ρ f using Eq. (58), the free current J f is given by (65) and while it may appear strange that we do not have p v, this is the function on phase space in Weyl correspondence withv = i [Ĥ,x]. This is an example of "hidden momentum", common to systems with magnetic moments [25,[27][28][29], and is found already in the longscale length limit [19], and also in Ref. [12]. Finally, the polarization and magnetization densities are given by Here, a factor 3 appears due to the normalization of the spin transform; the magnetization is proportional to the expectation value of σ. It should be noted that a model very similar to Eq. (58) was recently presented by Hurst et al. [30], starting from the same Hamiltonian. Their model, however, is formulated in terms of a scalar quantity f 0 and a vector quantity f , the mass and spin densities, respectively. Mathematically, this is writing the Hermitian matrixvalued Wigner function in the basis {I, σ}, and the correspondence between models is that f 0 = f d 2 s and f = 3 sf d 2 s. Hurst et al. find the same charge and current densities as we do [31] including the same "hidden momentum" in the free current or velocity. They also find the new non-linear term in the spin torque.
Hurst et al. [30] derive their model using a completely different method, based on a gauge invariant expression for the Moyal bracket [32]. That our results agree, up to trivial transformations, lends confidence to the results. We return to this in Section IV.

III. LINEAR WAVES
To illustrate the usefulness of Eq. (58) we study linear waves in a homogeneous plasma. We consider two cases: first electrostatic waves in a magnetized plasma, then electromagnetic waves in an unmagnetized plasma. From now on, we will let c = 1 to simplify the notation.
In the sums below it will be implicit that ν takes the values ±1.

A. Electrostatic waves
In this geometry we consider: B 0 = B 0ẑ , k = kẑ and E = Eẑ. The perturbed parameters follow the plane-wave ansatz according to f 1 =f e ik·x−ωt . The left hand side of Eq. (68) can then be written where ω ce is the cyclotron frequency, ω cg = g 2 ω ce and RHS is the right hand side of Eq. (68). Making an expansion off 1 in eigenfunctions of the operators of the right hand side of Eq. (68) [33] we writẽ Substituting Eq. (71) in Eq. (70), we can now express f 1 in terms of f 0 where and ∆ω ce = ω cg − ω ce . Next, noting that the magnetization current J m = ∇ × M = 0, we calculate the total current J = J f + J p , where J f and J p are the free and polarization current respectively, in order to obtain the dispersion relation. Note that the operators contained inẼ and ∆Ẽ act on F 0 as translation of the kinetic momentum, see Appendix A for more details. Introducing the notation where Here χ f L and χ pL are the contributions due to the free and polarization currents, respectively. We stress that in the first factor in χ pL , the sum over ± and ν corresponds to four terms in total. Below we will use the indices f and p with the same meaning also for other quantities.
We note that χ f L is the same susceptibility as for spinless quantum plasmas [34]. However, χ pL generalizes previous results [19] to cover also the short-scale physics.
Next, in order to examine how the spin affects the dynamic of the plasma, the Landau damping will be calculated and the contributions of χ f L and χ pL will be compared. Starting from Eq. (76), we Taylor expand (k, ω r + iγ) around ω r [35], which leads to γ = − Im (k, ω r ) ∂(Re )/∂ω| ω=ωr for |γ| |ω r |.
Analogously to Eq. (76), Im can be divided into Im = Im(χ f L ) + Im(χ pL ). Using the Plemelj formula, we get To see how Im χ pL contributes to γ compared to Im χ f L , we define .
The next step is to specify F 0ν . For simplicity, we will consider a Maxwell-Boltzmann distribution F 0ν = C ν e −p 2 /p 2 th e −νµeB0/k B T where C ν is the normalization constant according to d 3 pd 2 s F 0ν = N 0ν , where N 0ν are the densities of spin-up (ν = +1) and spin-down particles (ν = −1). This condition implies the textbook result that the magnetization is proportional to tanh µeB k B T . We then have Finally, the real part of the dispersion relation ω r (k) has to be specified explicitly in order to obtain Γ. Assuming that χ pL in Eq. (76) is small compared to χ f L [36] we neglect it for the purposes of computing ω r . This gives us Assuming ω kv th /m ± k 2 /2m, we can make a Taylor expansion of the denominators up to first non-vanishing order. The real part of the dispersion relation ω r (k) is then Using this expression of ω in Eq. (81), we can plot Γ, see Fig. 1. In general the spin contribution to the damping dominates (i.e., Γ → 1) for long wavelengths, whereas the free current contribution dominates (i.e., Γ → 0) for shorter wavelengths. The transition between these regimes is dependent on the normalized magnetic field B n = µ B B/m and the quantum parameter H = ω p /mv 2 th . We have performed the computation in Fig. 1 for various other values of v th , however, this reveals that the location of the transition region, expressed in k n = kv th /ω p , is almost independent of v th .

B. Electromagnetic waves
In this geometry we look at electromagnetic waves in unmagnetized plasmas with k = kẑ, E = Ex and B = Bŷ. We will allow for a background distribution depending on both p ⊥ and p z , but assume it to be spinindependent, as there is no background magnetic field. Dividing f , E and B into perturbed and unperturbed quantities, Eq. (68) becomes Using an expansion of f 1 like that in Eq. (71) f 1 is computed as, where and C = kµ e B + p z 2mẼ sin θ s ∂f 0 ∂p z .
un-magnetized plasma. The electrostatic dispersion relation possesses two roots. The usual Langmuir root and a root induced by the spin resonances with ω ∆ω c [19].
Here we have focused on the Langmuir root, showing that even when the real part of the frequency is almost unaffected by the spin dynamics, the wave-particle damping can still be significantly modified.
For the electromagnetic wave mode we note that all the new terms can be expressed in terms of the longitudinal susceptibility. Physically it means that the spin-terms couples the longitudinal and the transverse degrees of freedom. However, for scale lengths much longer than the Compton length (k mc/ ) this coupling is relatively weak.
A new feature of the present theory as compared to previous works is seen already for linear wave propagation, and comes from the denominators (ω − kp z /m ± ∆ω c ± k 2 /2m) −1 seen in Eq. (77). We note the importance of not letting µ e = µ B here, as otherwise we would get ∆ω c = 0. The physical origin of the resonances corresponding to (ω − kp z /m ± ∆ω c ± k 2 /2m) = 0 comes from energy momentum conservation written in the form Here (k 1 , ω 1 ) and (k 2 , ω 2 ) are the wavenumbers and frequencies of free particles before and after the interaction with the wave field, fulfilling the free particle dispersion relation ω 1,2 = k 2 1,2 /2m. The generalization from pre-vious results (e.g., Refs. [34,37]) comes from the term ± ∆ω c in (92), corresponding to resonant particles gaining an amount ω c of perpendicular kinetic energy, at the same time as the magnetic dipole energy drops by ω cg , or vice versa. For our particular case of electrostatic waves propagating parallel to B 0 , the kinetic energy and magnetic dipole energy always change in opposite directions, as reflected by the appearance of ∆ω c in (92). However, in a more general setting, the energy relation will be of the form whereas Eq. (93) will be unaffected. Here n is an integer and all sign combinations for the last three terms in (94) are possible. Finally we point out a problem of future interest; to compute the ponderomotive force starting from Eqs. (58) to (67). Previous results [38] indicates that the spin-orbit interaction can give a substantial contribution. Hence it would be of interest to see to what extent the shortscale physics of the present model modify the results of Ref. [38].