Short-scale quantum kinetic theory including spin-orbit interactions

P ^HYSICAL J ^OURNAL D

Regular Article - Plasma Physics

Short-scale quantum kinetic theory including spin–orbit interactions

R. Ekman ^1,2,a , H. Al-Naseri ^1,b , J. Zamanian ^1,c , and G. Brodin ^1,d

1 Department of Physics, Ume˚ a University, 901 87 Ume˚ a, Sweden

2 Center of Mathematical Sciences, University of Plymouth, Plymouth PL4 8AA, UK

Received 14 August 2020 / Accepted 2 November 2020

© The Author(s) 2021

Abstract. We present a quantum kinetic theory for spin-1/2 particles, including the spin–orbit interaction, retaining particle dispersive effects to all orders in , 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.

1 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 hydro- dynamic [7] and kinetic equations [8], where, in the present paper, the focus is on the latter class of the- ories. 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, sev- eral quantum kinetic models for plasmas have been put forward, within the Hartree approximation [12–15] or the Hartree–Fock approximation [16–18]. In the latter case, exchange effects are also included, which might be important in various applications [1,2, 8]. Depending on the scope of the theory, the density matrix can be based on the Schr¨ odinger equation [8], the Pauli equation [13], or the Dirac equation, where the latter equation can be applied in the weakly [19] or fully [12] relativistic approximation.

In order to produce a quantum kinetic theory of elec- trons from the Dirac equation, certain restrictions need to be applied. Firstly, pair-production must be negli- gible, 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 ² c ³ / |q|, and the characteristic

a e-mail: robin.ekman@plymouth.ac.uk

b e-mail: haidar.al-naseri@umu.se

c e-mail: jens.zamanian@umu.se (corresponding author)

d e-mail: gert.brodin@umu.se

spatial scales of the fields, which should be much longer than the Compton length L _c = /m c. 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 differ- ent 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, par- ticle 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 pre- cession 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-o-rbit interaction. In

Sect. 2, 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. Our main result is the kinetic equation, Eq. (58)

in Sect. 2.5.

In Sect. 3, in order to demonstrate the usefulness of the present theory, we have calculated two exam- ples from linearized theory: electromagnetic waves in a non-magnetized plasma and electrostatic waves prop- agating parallel to an external magnetic field. Classi- cally, an external magnetic field does not affect par- allel 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 damp- ing rate.

Finally, in Sect. 4, we compare our model with some theories from the recent research literature.

2 The gauge invariant Wigner function and the spin transformation

We will consider the semi-relativistic Pauli Hamilto- nian

H = ˆ π ˆ ²

2m + q ˆ φ − μ e σ · ˆ B + μ e

2mc ²



π × ˆE − ˆE × ˆπ ˆ 

· σ, (1) including the spin–orbit interaction, applying the mean field approximation. This means we are dropping inter- actions due to exchange effects and due to correlations (see e.g. [16–18, 24]) associated with the two-particle correlation function, and possibly also higher order cor- relations in the BBGKY-hierarchy [25], limiting our- selves to the Hartree approximation. Here ˆ π = ˆp − q ˆ A is the gauge invariant momentum operator, ˆ φ, ˆ A are the scalar and vector potentials and μ _e is the electron mag- netic moment. Note that the spin–orbit interaction is written in a symmetric form so as to make the Hamil- tonian 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,

W (x, p, t) αβ = tr

 W (x, p, t)ρ ˆ αβ



=



d ³ r r| ˆ W ˆ ρ αβ |r, (2) 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 ˆ W (x, p, t) is time-dependent, but the density matrix ρ is not. The operator ˆ W (p, x, t) can be expressed as a Fourier trans- form

W (x, p, t) = ˆ  F( ˆ T )

 (x, p)

=

 d ³ u (2π) ³

d ³ v

(2π) ³ e ⁻

^i

^(u·p+v·x) T (u, v), ˆ (3)

where the operator ˆ T (u, v, t) is given by T (u, v, t) = exp ˆ

 i

 (u · ˆπ + v · ˆx)



. (4)

To obtain a scalar function from the matrix-valued W ˆ _αβ , one can use the spin transform

f (x, p, s) = 1

4π tr[(1 + s · σ)W ]

= 1 4π

α,β

[δ αβ + s · σ αβ ]W βα , (5)

which is a Husimi Q-function for the spin [26]. This transform was used in Ref. [13] for a theory including only the magnetic dipole interaction.

If we let

S = ˆ 1

4π (1 + s · σ), (6) then

f = Tr[ ˆ S ˆ W ρ], (7) i.e. f is the expectation value of ˆ S ˆ W . Note that ˆ S and W are both Hermitian and commute, which implies ˆ that f is real. The Heisenberg equation of motion then gives

∂ _t f = Tr



S∂ ˆ _t W ˆ − i

 S[ ˆ ˆ W , ˆ H] − i

 [ ˆ S, ˆ H] ˆ W 

ρ

 ,

(8) where in this and the previous equation, the trace is over both the spin indices and the spatial degrees of freedom. The only effect of ∂ t W is to give the ∂ ˆ t A con- tribution to the electric field in the Lorentz force.

The remainder of this section presents the calcula- tions involved in expressing the right-hand side of this equation in terms of f , with the conclusion presented in Sect. 2.5. The steps are similar to those in Ref. [23], where the spin-independent terms in the kinetic equa- tion are found.

2.1 Gross structure of the calculation

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., E _i,j = ∂E _i /∂x _j .

The main thrust of the calculation is to take the matrix elements x ^ | · · · |x between position eigenstates of the operators on the right-hand side of (8) then to move as much as possible outside the inner product using commutation relations and the Baker–Campbell–

Haussdorff formulas for ˆ T , Eq. (4), below, to put every-

thing in terms of phase-space differential operators on

f . For a simplest-cases demonstration, we recommend

the treatment in Ref. [23] of the scalar potential.

Baker–Campbell–Haussdorff formulas for ˆ T can be found [23] and will be useful in the following, as we shall often want to put ˆ π:s or ˆ x:s on the left or the right, as needed. They are:

T (u, v) = exp ˆ

 i 2  u · v



× exp

 i

 v ·  ˆ

x − qu · ₁

0 dτ ˆ A(ˆ x + τ u) 

× exp

 i

 u · ˆp



= exp



− i 2 u · v

 exp

 i

 u · ˆp



× exp

 i

 v ·  ˆ

x − qu · ₁

0 dτ ˆ A(ˆ x + τ u) 

= exp

 i 2  u · v

 exp

 i

 v · ˆx

 exp

 i

 u · ˆπ



= exp



− i 2  u · v

 exp

 i

 u · ˆπ

 exp

 i

 v · ˆx

 .

One can note that these occur in pairs related by oper- ator orderings.

The commutator is linear in ˆ H, and the terms in the evolution equation for ˆ W 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

H ˆ _SO = μ e

2mc ² ( ˆ π × ˆE − ˆE × ˆπ) · σ. (9) We will denote the contribution from this Hamiltonian by (∂ _t f ) _SO .

The spin–orbit interaction Hamiltonian has the form σ · ˆ V, and hence we find

[ ˆ W , ˆ H _SO ] = σ i [ ˆ W , ˆ V i ], (10) [ ˆ S, ˆ H _SO ] = 2i s i ijk V ˆ j σ k . (11) Then using σ i σ j = δ ij + i ijk σ k , we get

S[ ˆ ˆ W , ˆ H _SO ] = (σ j + s j + i ijk s i σ k )[ ˆ W , ˆ V j ]. (12)

For the second commutator, we write ˆ V ˆ W as the sum of its symmetric and anti-symmetric parts

2i _ijk s _i V ˆ _j σ _k W ˆ

= i ijk s i σ k

 V ˆ j W + ˆ ˆ W ˆ V j

 + [ ˆ V j , ˆ W ]

 (13) and clearly, the second part cancels with cross product part of ˆ S[ ˆ W , ˆ H _SO ] in Eq. (12). Now we just need to use

∂ s

_i

S = ˆ 1

4π (σ i − s i s j σ j ) , (14)

(see footnote ¹ ) which implies

(∂ _s

_i

+ s _i ) ˆ S = 1

4π (σ _i + s _i ) , (15) _ijk ∂ _s

_i

(s _j S) = ˆ 1

4π _ijk σ _i s _j , (16) and thus

(∂ _t f ) _SO = − i

 Tr 

(∂ _s

_i

+ s _i ) ˆ S[ ˆ W , ˆ V _i ] + ijk ∂ s

_i

s j S[ ˆ ˆ W , ˆ V k ] ₊ 

ρ



. (17)

The remaining, laborious, task is to evaluate and take the Fourier transform of ˆ T ( ˆ π × ˆE) ± (ˆπ × ˆE) ˆ T , then express the result in terms of differential operators—

functions of ∂ _p and ∂ _x —acting on ˆ W . These operators can be taken outside the trace, so that they act on f , and this will give the kinetic equation. Since

[ ˆ T , ˆ π × ˆE] = [ ˆ T , ˆ π] × ˆE + ˆπ × [ ˆ T , ˆ E], (18)

we need to evaluate the commutators with ˆ π and ˆE.

This is sufficient, because for any operator ˆ O

[ ˆ W , ˆ O ^† ] = −[ ˆ W , ˆ O] ^† , (19)

since ˆ W is Hermitian, and we get the commutator with the Hermitian ordering [which appears in Eq. (9)] by reading off the “imaginary part” of [ ˆ W , ˆ π × ˆE]. Sec- tions 2.2 to 2.4 present the details.

2.2 The commutator [ ˆ W , ˆ E]

In Sects. 2.2 to 2.4, we will typically omit boldface for vector quantities; the type of each quantity should be apparent from context.

Using Eq. (4), we can write the operator ˆ T with Θ(u) = exp( ˆ _ ⁱ u · ˆp) on the right. Since ˆp is the gen- erator of translations, the operator Θ(u) acts like

Θ(u)|x = |x − u ˆ (20)

Θ(u)f (ˆ ˆ x) = f (ˆ x + u) ˆ Θ(u) (21)

1 Equations (14)–(16) can be derived using the definition

Eq. (6), together with the expression for the gradient in

spherical coordinates, expressed in terms of the Carte-

sian basis. For example ∇ s

x

= ˆ x · ∇ s = cos θ cos φ∂ θ −

(sin φ/ sin θ)∂ φ . Alternatively, we may temporarily relax

the normalization requirement on the spin s and define

ˆ s = s/|s|. Taking the derivative results in ∂ s

i

ˆ s j = δ ij /|s| −

s i s j /|s| ³ . Now, applying the constraint |s| = 1 we arrive at

Eq. (14).

and thus, using Eq. (4) to put the translation operator on the left, with the other factors commuting with ˆ E,

[ ˆ T , ˆ E i ] = ˆ T (u, v)

 E ˆ i (ˆ x) − ˆ E i (ˆ x − u) 

(22)

Now we need an expression for the product ˆ π×[ ˆ T , ˆ E]

in terms of ˆ T . By putting the commutator in the form Eq. (22), we can use any expression for ˆ T as is most appropriate. The one that is most appropriate here is Eq. (4) since we can lift the result from Ref. [23]

their (4.30b), and its Hermitian conjugate:

ˆ π i exp

i

 u · ˆπ 

= 

i

∂

∂u i − q

 ₁

0

dτ (1 − τ)

×[u × ˆ B(ˆ r + τ u)] i

 exp

i

 u · ˆπ 

, (23)

× exp

− i

 u · ˆπ 

ˆ π _i

= exp

− i

 u · ˆπ  ⎛ ⎝−

i

← ∂

∂u _i

−q

 ₁

0

dτ (1−τ)[u × ˆ B(ˆ r+τ u)] i

 . (24) We thus have that the operator π _i [ ˆ T , E _j ] can be expressed as π _i [ ˆ T , E _j ] = ˆ O ^B _ij + ˆ O _ij ^D where

O ˆ ^B _ij = −q

 ₁

0

dτ (1 − τ)[u × ˆ B(ˆ x + τ u)] i T (u, v) ˆ

× 

E ˆ j (ˆ x) − ˆ E j (ˆ x − u) 

, (25)

O ˆ _ij ^D = e ⁻

^2i

^{u ·v}

  i

∂

∂u _i exp i

 u · ˆπ 

× exp i

 v · ˆx  

E ˆ j (ˆ x) − ˆ E j (ˆ x − u)  . (26)

Now let |x ^ , |x ^  be eigenstates of ˆx (at the appro- priate time). Use Eq. (4) so that the ˆ x operators are on the left in ˆ T .

Let ˆ Q ^B _ij be the Fourier transform in u, v of ˆ O ^B _ij . Then, letting e ^{i ˆ Φ} = exp

 − ^iq _ u · ₁

0 dτ ˆ A(ˆ x + τ u)



its matrix elements are

x ^ | _Q ˆ ^B _{ij |x} ^ _ = F 

x ^ | _O ˆ _ij ^B _|x ^ _ 

= −qF



e

^i

^v·(x

^

^+u/2)

 ₁

0 dτ (1 − τ )

× [ u × B ( x ^ + τ u )] _i x ^ |e ^iΦ e

^i

^{u·ˆ p}

 _E ˆ _j (ˆ x ) − _E ˆ _j (ˆ x − u ) 

|x ^ 



= −qF



e

^i

^v·(x

^

^+u/2)

 ₁

0 dτ (1 − τ )

× [ u × B ( x ^ + τ u )] i x ^ |e ^iΦ 

E ˆ j (ˆ x + u ) − _E ˆ _j (ˆ x ) 

|x ^ −u 

= −qF



e

^i

^v·(x

^

^+u/2)

 ₁

0 dτ (1 − τ )

× [ u × B ( x ^ + τ u )] i

E j ( x ^ + u ) −E j ( x ^ )

x ^ |e ^iΦ |x ^ −u 

.

(27)

At this point we see that the Fourier transform in v can be carried out and will give an expression propor- tional to δ(x ^ − x + u/2). Hence the matrix elements of Q ˆ ^B _ij (x, p) remain unchanged if we change x ^ to x − u/2 in the arguments 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 up ˆ T , resulting in

x ^ | _Q ˆ ^B _{ij |x} ^ _ = −qF

 _1/2

−1/2 dτ

 1 2 ^{− τ}



[ u × B ( x + τ u )] i

× 

E j ( x + u/ 2) − E j ( x − u/ 2) 

×x ^ |e

^2i

^u·v e

^i

^v·x e ^{i ˆ Φ} |x ^ 

= −qF

 _1/2

−1/2 dτ

 1 2 ^−τ



[ u × B ( x + τ u )] i  E j ( x + u/ 2) −E j ( x−u/ 2) 

x ^ | _T ˆ ( u, v ) |x ^  

(28)

Since the Fourier transform sends u → i∂ p , we can now express this operator in terms of ˆ W and its deriva- tives. First, use that

E j (x + u/2) − E j (x − u/2) = u k

 _1/2

−1/2 dτ E j,k (x + τ u), (29) 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

Im

 Q ˆ ^B _ij



= −q ²

 _1/2

−1/2 dτ τ [∂ _p × B(x + iτ∂ p )] _i

×

 _1/2

−1/2 dσ E j,k (x + iσ∂ p )∂ p

_k

W . ˆ (30)

We now turn to ˆ O _ij ^D . Note that from the expression for ˆ T in Eq. (4) we have

∂

∂u i

T = ˆ − i

2 v i T + e ˆ ⁻

^2i

^u·v ∂

∂u i

e

^i

^u·ˆπ e

^i

^v·ˆx (31)

and hence

O ˆ _ij ^D = ˆ T

⎡

⎣v ⁱ 2 + 

i

← ∂

∂u i

⎤

⎦  ˆ E _j (ˆ x) − ˆ E _j (ˆ x − u) 

. (32)

The first term can be handled using the argument above, except that the Fourier transform in v will send v i to i∂ x

_i

δ(x). The result is that

F 

x ^ | ˆ T v _i 2

 E ˆ _j (ˆ x) − ˆ E _j (ˆ x − u) 

|x ^  

= −  2i ∂ x

_i

 E j (x + i∂ p /2)

−E j (x − i∂ p /2) 

x ^ | ˆ W |x ^  

(33) but this quantity is real, and so makes no contribu- tion to the evolution equation. For the term with the u-derivative, we integrate by parts in the Fourier trans- form, i.e.,

F

 ∂

∂u _i T ˆ   E ˆ _j (ˆ x) − ˆ E _j (ˆ x − u) 

= F



p i T ˆ  E ˆ j (ˆ x) − ˆ E j (ˆ x − u) 

− 

i T ˆ ˆ E j,i (ˆ x − u)

 . (34) In the first term, p _i can be taken outside the Fourier transform, and the remaining operator is treated like above resulting in

F 

x ^ | ˆ O _ij ^D |x ^  

=  p i

 E j (x + i∂ p /2)

−E j (x − i∂ p /2) 

− 

i E _j,i (x − i∂ p /2)



x ^ | ˆ W |x ^ .

(35) We can write this in a somewhat simpler form as

E _j (x + i ∂ p /2) − E j (x − i∂ p /2) = i (∂ x

k

E ˜ _j )∂ _p

_k

(36) where we have introduced the notation

E ˜ _j =

 _1/2

−1/2

E _j (x + i τ∂ p ) dτ. (37)

To summarize, we obtain

Im( ˆ π × [ ˆ W , ˆ E]) = (p + Δ˜p) × ˜E( ∂ ^← x · ∂ ^→ p ) ˆ W (38) where

Δ˜ p = − q

i ∂ p ×

 _1/2

−1/2 dτ τ B(x + iτ∂ p ). (39)

The quantities ˜ E 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.

2.3 The commutator [ ˆ W , ˆπ]

We use Eq. (4) for ˆ T , so that

[ ˆ T , ˆ π i ]E j = [−ˆπ i , ˆ T ]E j = −e ⁻

^2i

^u·v

 ˆ π i , exp

i

 u · ˆπ 

× exp i

 v · ˆx 

E j

= e ⁻

^2i

^{u ·v}

−v i + q

 ₁

0

dτ [u × ˆ B(ˆ x + τ u)] _i 

× exp i

 u · ˆπ 

exp i

 v · ˆx 

E j

= −v i T (u, v)E ˆ _j + q

 ₁

0

dτ [u × ˆ B(ˆ x + τ u)] i T (u, v)E ˆ j

= ˆ O ^A _ij + ˆ O _ij ^C , (40) also using (4.25) from Ref. [23]. We can now look at the matrix elements of F 

O ˆ ^A _ij



, but it is of the form in Eq. (33) except we have only one ˆ E operator, ˆ E(ˆ x).

Therefore, F 

x ^ | ˆ O _ij ^A |x ^  

= x ^ |  i ∂ _x

_i



E _j (x + i∂ p /2) ˆ W

 |x ^ 

= x ^ |  i



E j,i (x + i∂ p /2) ˆ W +E _j (x + i ∂ p /2)∂ _x

_i

W ˆ

 |x ^ . (41)

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

 _1/2

−1/2 E j (x + τ u) dτ + u k

 _1/2

−1/2 τ E j,k (x + τ u) dτ

= 1 2

 E j (x + u/2) + E j (x − u/2) 

(42) Again keeping even powers of i∂ p since we are dividing by i,

Im 

i E _j (x + i ∂ p )

= −  2

 E j (x + i∂ p /2) + E j (x − i∂ p /2) 

= − ˜ E j − i ² ∂ p

_k

 _1/2

−1/2 τ E j,k (x + iτ∂ p ) dτ.(43)

This is the “hidden momentum” velocity correction term ∝ E × s [ 27] (to be discussed more below), but with E = ˜E + Δ˜E where Δ˜E is higher order in .

Following steps like those above, and using Eq. (42) again, one readily expresses the matrix elements of ˆ O _ij ^C in terms of those of ˆ W . The result is that

x ^ | ˆ O ^C _ij |x ^  = inm B ˜ _m ( ˜ E _j + Δ ˜ E _j )∂ _p

_n

x ^ | ˆ W |x ^  (44) which, in the kinetic equation, will give

ijk inm B ˜ m ( ˜ E j + Δ ˜ E j )∂ p

_n

σ k W ˆ

∝ 

( ˜ E + Δ˜E) × σ 

× ˜ B 

· ∂ p W , ˆ (45) i.e., it is the correction to the magnetic force due to the relation between momentum and velocity. Here ˜ B is defined in the same way as ˜ E, Eq. ( 37).

2.4 The spin transform

It remains to evaluate the terms in the kinetic equation arising from the spin transform, proportional to

_ijk s _i σ _k 

( ˆ π × ˆE − ˆE × ˆπ) j W + H. c. ˆ 

where H. c. denotes the Hermitian conjugate. Using Eq.

(23) again we establish that ˆ

π k E ˆ l (ˆ x) ˆ T = (D k T )E ˆ l (ˆ x − u), (46) E ˆ l (ˆ x)ˆ π k T = ˆ ˆ E l (ˆ x)D k T , ˆ (47) where D k is the operator

D k =  i

∂

∂u k

+ v k

2 −q

 ₁

0

dτ (1−τ)[u× ˆ B(ˆ x+uτ )] k . (48)

Hence the product ( ˆ π × ˆE − ˆE × ˆπ) j T is given by ˆ _jkl (ˆ π _k E ˆ _l − ˆ E _j π ˆ _k ) ˆ T = _jkl (ˆ π _k E ˆ _l + ˆ E _k ˆ π _l ) ˆ T

= jkl

  i



∂ ˆ T

∂u k

E ˆ l (ˆ x − u) + ˆ E l (ˆ x) ∂ ˆ T

∂u k



+v k E ˆ l (ˆ x) ˆ T − 2I k E l (ˆ x) ˆ T



, (49)

where

I k =

 ₁

0

dτ (1 − τ)[u × ˆ B(ˆ x + uτ )] k . (50)

Differentiating ˆ E _l (ˆ x) ˆ T = ˆ T ˆ E _l (ˆ x−u) with respect to u k

we obtain E ˆ l (ˆ x) ∂

∂u k

T = ˆ ∂ ˆ T

∂u k

E ˆ l (ˆ x − u) − ˆ T ˆ E l,k (ˆ x − u). (51)

so that

∂ ˆ T

∂u _k E ˆ l (ˆ x − u) + ˆ E l (ˆ x) ∂ ˆ T

∂u _k

= 2 ˆ E l (ˆ x) ∂ ˆ T

∂u k

+ ˆ T ˆ E l,k (ˆ x − u)

= 2 ˆ E _l (ˆ x) ∂ ˆ T

∂u k

+ ˆ E _l,k (ˆ x) ˆ T . (52)

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 E ˆ l (ˆ x) _∂u ^{∂ ˆ T}

k

we can handle as in Eq. (34), resulting in

∼ p k E ˜ _j (x − i∂ p /2)W,

and following the argument leading to Eq. (43), when taking the real part, we will get

∼ p k ( ˜ E l + Δ ˜ E l )W. (53) 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

∼ i (E l (x + i ∂ p /2) − E l (x − i∂ p /2)) ∂ _x

_k

W, (54) which can be simplified using Eq. (36).

Following steps like those for ˆ O ^B (Eq. (25) and fol- lowing), the final term in Eq. (49), proportional to I k , will give us a contribution of the form

−

 _1/2

−1/2 dτ 1

2 − τ 

[u × B(x + τu)] k E l (x − u/2) ˆ T . (55) When taking the real part of the Fourier transform of this, there will be two terms,

∼

 _1/2

−1/2 τ [i∂ p × B(x + iτ∂ p )] k dτ (E l (x + i∂ p /2) +E _l (x − i∂ p /2))

W ∝ Δ˜p × (˜E + Δ˜E)W, (56)

which is in line with previous generalisations, and

∼

 _1/2

−1/2 [i ∂ p × B(x + iτ∂ p )] _k dτ (E _l (x − i∂ p /2)

−E l (x + i∂ p /2))

W ∝ ( ˜ B × ∂ p ) × (˜E ∂ ^← x · ∂ ^→ p )W, (57)

which is of a new type. However, both are order  ² ,

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 pro- portional to (˜ p+Δ˜p)×(˜E+Δ˜E)+ ^q ₂ ( ˜ B×∂ p )×·(˜E ∂ ^← x ·

∂ → _p ). The overall sign and coefficient of this term is found by matching its long-scale limit to the model in Ref. [19]

2.5 Conclusion

The evolution equation for f follows directly by taking the expectation value of the evolution equation for ˆ W . It is

∂ t f + 1 m

p + Δ˜ p + μ e

2mc ² ( ˜ E + Δ ˜ E) × (s + ∇ s )

· ∇ x f +q 

( ˜ E+

p+Δ˜ p+ μ e

2mc ² ( ˜ E + Δ ˜ E) × (s + ∇ s )

× ˜ B 

×∇ p f + μ e



B ˜ − p + Δ˜ p 2mc ² × ˜E



· (s + ∇ s )( ∇ ^← x · ∇ ^→ p )f

+ 2μ e

 

s ×



B ˜ + Δ ˜ B − 1

2mc ² (˜ p + Δ˜ p) × (˜E + Δ˜E)

− q 4mc ² ( ˜ B × ∂ p ) × ·(˜E ∇ ^← x · ∇ ^→ p ) 

· ∇ s f = 0. (58)

For reference, the notation used is

E(x) = ˜

 _1/2

−1/2 E(x + iτ∇ p /2) dτ

=

 E

 _1/2

−1/2 cos(τ ∇ ^← x · ∇ ^→ p ) dτ

 (x), (59)

Δ ˜ E(x) = i

 _1/2

−1/2

τ E(x + iτ∇ p ) dτ ( ∇ ^← x · ∇ ^→ p )

= −  E

 _1/2

−1/2 sin(τ ∇ ^← x · ∇ ^→ p ) dτ ∇ ^← x · ∇ ^→ p

 (x),

(60) Δ˜ p = −iq

 _1/2

−1/2 τ B(x + iτ∇ p ) dτ × ∇ p

= q   B

 _1/2

−1/2

sin( τ ∇ ^← x · ∇ ^→ p ) dτ



( x) × ∇ p , (61) with ˜ B, Δ ˜ B defined in the same way as ˜E, Δ˜E. Func- tions 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 momen- tum” (see below) correction to the velocity in the diffu- sion and magnetic force terms, respectively; the spin–

orbit force; the spin torque due to spin–orbit inter- action, 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 equa- tion.

We have denoted the electron magnetic moment by μ e to allow for an anomalous magnetic moment, i.e., μ e = ^g ₂ μ 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 [28], as will be seen below.

The system is closed with Maxwell’s equations with polarization and magnetization,

∇ · E = 1

₀ (ρ _f − ∇ · P), (62)

∇ × B = μ 0 J f + μ ₀ ∇ × M + 1 c ²

∂E

∂t + μ ₀ ∂P

∂t . (63) The free charge ρ f is given by

ρ _f = q



f dΩ, (64)

introducing dΩ = d ³ p d ² s. Finding the continuity equa- tion for ρ _f using Eq. (58), the free current J f is given by

j f = q



vf dΩ = q

 p

m + 3μ _e E × s 2m ² c ²



f dΩ, (65)

and while it may appear strange that we do not have p  v, this is the function on phase space in Weyl cor- respondence with ˆ v = _ ⁱ [ ˆ H, ˆ x]. This is an example of

“hidden momentum”, common to systems with mag- netic moments [27,29–31], and is found already in the long-scale length limit [19], and also in Ref. [12].

Finally, the polarization and magnetization densities are given by

P = −3μ e

 s × p

2mc ² f dΩ, (66) M = 3μ e



sf dΩ. (67)

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. [32], starting from the same Hamiltonian. Their model, however, is formulated in terms of a scalar quantity f ₀ and a vec- tor quantity f, the mass and spin densities, respectively.

Mathematically, this is writing the Hermitian matrix- valued Wigner function in the basis {I, σ}, and the cor- respondence between models is that f ₀ =

f d ² s and f = 3

sf d ² s. Hurst et al. find the same charge and current densities as we do ² including the same “hidden

2 Taking into account that their convention for the sign of

the charge is different from ours.

momentum” in the free current or velocity. They also find the new non-linear term in the spin torque.

Hurst et al. [32] derived their model using a com- pletely different method, based on a gauge invariant expression for the Moyal bracket [33]. While our results are equivalent to those of Ref. [32], apart from some details to be discussed in the next section, we note that the derivation made here has a value in itself. In case generalizations of the governing equations are sought for, aiming to extend the regime of validity, having more than one route to follow increases the possibility for suc- cess. We will return to related issues in Sects. 3 and 4.

Finally, we note that several models somewhat simi- lar to Eq. (58), but allowing for fully relativistic motion have been published in the literature, see e.g. [34–36].

While these models allow for p ∼ mc, they are subject to other limitations and/or drawbacks. Specifically, this includes not separating between electron and positron states and using a perturbative quantum mechanical treatment [36].

3 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.

The momentum is expressed in cylindrical coordi- nates p = (p ⊥ , ϕ p , p z ) and the spin in spherical coordi- nates s = sin θ s cos ϕ s z + sin θ ˆ s sin ϕ s y + cos θ ˆ s ˆ z. In the linearization of Eq. (58) variables are separated accord- ing to f = f ₀ (p _⊥ , p z , θ s ) + f ₁ (x, p, s, t), E = E 1 and B = B 0 ˆ z + B 1 , where the subscripts 0 and 1 denote unperturbed and perturbed quantities respectively. Eq.

(58) then becomes

∂f 1

∂t + p

m · ∇ x f 1 + p m × B 0

· ∇ p f 1 + 2μ e

h (s × B 0 ) · ∇ s f 1

= −q  E ˜ +

 p + Δ˜p m + μ e

m ( ˜ E + Δ ˜ E ) × (s + ∇ s )



×B 0 + p m × ˜ B ₁

 · ∇ p f 0

+μ e ∇ x

 B ˜ ₁ − p × ˜E 2m

· (s + ∇ s )

· ∇ p f 0

− 2μ e

 s ×



B ˜ + Δ ˜ B − q  ² 4m

 B ₀ × ∇ p 

×  ˜ E ∇ ^← x · ∇ ^→ p 

− (p + Δ˜p) × (˜E + Δ˜E) 2m

· ∇ s f 0 . (68)

In general the background distribution can be divided into its spin-up and spin-down components [37]

f ₀ (p _⊥ , p z , θ s ) = 1 4π

ν=±1

F _0ν (p _⊥ , p z )(1+ν cos θ s ). (69)

In the sums below it will be implicit that ν takes the values ±1.

3.1 Electrostatic waves

In this geometry we consider: B ₀ = B ₀ ˆ z, k = kˆz and E = Eˆz. The perturbed parameters follow the plane- wave ansatz according to f ₁ = ˜ f e ^ik·x−ωt . The left hand side of Eq. (68) can then be written



iω − i kp _z

m + ω _ce ∂

∂ϕ p

+ ω _cg ∂

∂ϕ s



f ₁ = RHS, (70)

where ω ce is the cyclotron frequency, ω cg = ^g ₂ ω ce and RHS is the right hand side of Eq. (68). Making an expansion of ˜ f ₁ in eigenfunctions of the operators of the right hand side of Eq. (68) [37] we write

f ˜ ₁ ( p, s) =

α,β

1

2π g _α,β (p _⊥ , p _z , θ _s )e ^−i(αϕ

^s

^+βϕ

^p

⁾ . (71)

Substituting Eq. (71) in Eq. (70), we can now express f ₁ in terms of f ₀

f ₁ = A

ω − kp z /m + B ₊ e ^i(ϕ

^s

^−ϕ

^p

⁾ ω − kp z /m + Δω ce

+ B ₋ e ^−i(ϕ

^s

^−ϕ

^p

⁾ ω − kp z /m − Δω ce

, (72)

where

A =

ν

−iq ˜ E ∂F _0ν

∂p z

, (73)

B _± = iμ e

ν

sin θ s

− qB ₀

4m ( ˜ E + Δ ˜ E) ∂F _0ν

∂p _⊥

− νp _⊥

2m  ( ˜ E + Δ ˜ E)F _0ν ± kp _⊥ E ˜ 4m

∂F _0ν

∂p _z

∓ νk qB 0 E ˜ 4m

∂ ² F _0ν

∂p _⊥ p z



, (74)

and Δω _ce = ω _cg − ω ce . Next, noting that the magne- tization 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 ˜ E and Δ˜E act on F 0 as translation of the kinetic momentum, see “Appendix A” for more details.

Introducing the notation

p q = k 2 , p ^± _res = mω

k ± p q ,

the dispersion relation is then

(k, ω) = 1 + χ _{f L} + χ _pL = 0, (75) where

χ _{f L} = 4π ² q ² p _q k ² m

 d ² p

1

p ⁺ _res − p z

− 1

p ⁻ _res − p z

 F _0ν ,

(76) χ pL = − π ² μ ² _e

km ²

 ±,ν

d ² p

 p ² _⊥

m  (ν ∓1) − ω ce (1∓2ν)



1

p ^∓ _res −p z +mΔω ce /k − 1

p ^± _res − p z − mΔω ce /k 

F _0ν . (77) 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 sus- ceptibility as for spinless quantum plasmas [38]. How- ever, χ pL generalizes previous results [19] to cover also the short-scale physics. Specifically, this means that the denominators of Eq. (77), responsible for the resonant wave-particle interaction to be studied below, are a new property of the present theory. To obtain such denomi- nators we need to keep all extra features of the present theory; spin–orbit interaction, the polarization current associated with the spin, particle dispersive effects, and also use the QED-corrected expressions for the spin- precession frequency.

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. (75), we Taylor expand (k, ω r + iγ) around ω r [39], which leads to

γ = − Im (k, ω r )

∂(Re )/∂ω| ω=ω

r

for |γ|  |ω r |. (78)

Analogously to Eq. (75), Im can be divided into Im = Im(χ _{f L} ) + Im(χ _pL ). Using the Plemelj formula, we get

Im χ f L = 8π ³ q ² m

k ³



dp _⊥ p _⊥



F _0ν (p _z = p ⁺ _res ) − F _0ν (p _z = p ⁻ _res )

 , (79) Im χ _pL = π ³ μ ² _e

k

 ±,ν

dp _⊥ p _⊥

 p ² _⊥

m  (ν ± 1) − ω ce (1∓2ν) 



F _0ν (p z = p ^∓ _res + mΔω ce /k)

− F 0ν (p z = p ^± _res − mΔω ce /k)



. (80)

To see how Im χ _pL contributes to γ compared to Im χ f L , we define

Γ = γ p

γ _p + γ _f = Im χ pL

Im(χ _{f L} + χ _pL ) . (81) The next step is to specify F _0ν . For simplicity, we will consider a Maxwell–Boltzmann distribution F _0ν = C _ν e ^−p

²

^/p

^2th

e ^−νμ

^e

^B

⁰

^/k

^B

^T where C _ν is the normalization constant according to

d ³ pd ² s F _0ν = N _0ν , where N _0ν are the densities of spin-up (ν = +1) and spin-down particles (ν = −1). This condition implies the text- book result that the magnetization is proportional to tanh _k ^μ

^e

^B

B

T . Note, however, that the above is just a spe- cial case, as the theory presented here is applicable for arbitrary ratios of T /T _F . For the case of completely or partially degenerate electrons, we just need to adopt the thermodynamic background distribution function given by Eq. (60) of Ref. [13]. We then have

Im χ fL =

√ πm ² ω _p ²

 k ³ p th



e ^−(p

^−res

^/p

^th

⁾

²

− e ^−(p

^+res

^/p

^th

⁾

²



, (82)

Im χ pL =  ² ω _p ² √ π 32 m ² kp th



±

−p ² _th m 



tanh ^{μ B B 0} k B T ± 1



− ω ce



1 ± 2 tanh ^{μ B B 0} k B T



×  e ⁻

p

^∓res

+mΔω

ce

/k

2

/p

²th

− e ⁻

p

^±res

−mΔω

ce

/k

2

/p

²th



.

(83) Finally, the real part of the dispersion relation ω _r (k) has to be specified explicitly in order to obtain Γ.

Assuming that χ pL in Eq. (75) is small compared to χ f L (see footnote ³ ) we neglect it for the purposes of computing ω _r . This gives us

_r ≈ 1 − mω _p ² kω √

πp _th

 dp z

 p z /k − 1/2 p ⁺ res − p z

− p z /k + 1/2 p ⁻ res − p z



e ^−p

^2z

^/p

^2th

. (84) Assuming ω  kv th /m ±k ² /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

ω r (k) = ±



ω ² _p + k ² v _th ² +  ² k ⁴

4m ² . (85) Using this expression of ω in Eqs. (82) and (83), we can plot Γ, see Fig. 1. In general the spin contribution to the

3 This will be the case at least if ω p /m  1.