Conservation laws in kinetic theory for spin-1/2 particles

(1)

Conservation laws in kinetic theory for spin-1/2 particles

Linn´ ea Gr¨ ans Samuelsson

Abstract

In this thesis a kinetic theory for spin-1/2 particles is given a brief overview, focusing on the derivation of an evolution equation for the quasiprobability distribution function used in the theory to describe certain types of quantum plasma. The current theory is expanded upon by exploring conservation laws. A local conservation law for momentum is derived using two different expressions for electromagnetic momentum, given by Abraham and Minkowski respectively. There has been some controversy over which of these expressions should be used;

in the case considered here the expression given by Minkowski seems to be more suitable.

Based on the conservation law for momentum, a conservation law for angular momentum is also derived.

Bachelor’s thesis in physics Supervisor: Jens Zamanian

(2)

(3)

1 Introduction

Plasmas can be described as ionized gasses dominated by collective effects due to the charged particles’ interaction via electromagnetic forces. They have traditionally been treated classically, but there is an increasing interest in the quantum mechanical treatment needed to understand properties of so called quantum plasmas. These are plasmas where quantum properties are of importance, for which a common criterion is that the Fermi temperature should be at least as high as the thermodynamic temperature. The increased interest can be related to recent experimental development, such as the possibility to create quantum effects in plasmas with the help of lasers [1].

One essentially quantum mechanical aspect of a plasma made up by spin-1/2 particles, such as electrons, is the intrinsic spin possessed by said particles. In some experimentally accessible regimes, the energies and forces related to the dipole moments due to this spin are of a magnitude comparable to the other energies and forces present. To see which regimes are relevant, there are some parameters that are of interest. The criterion that the Fermi temperature T_F should be at least as high as the thermodynamic temperature T can be met by demanding that the parameter TF/T should be large, which relates the temperature of the plasma to the density. In this regime there are quantum effects due to Fermi statistics replacing the classical Boltzmann statistics. In the case considered in this thesis we will instead consider quantum effects due to spin, and then there is no need to demand that T TF. In fact we will assume the opposite to simplify some calculations, although the results do not depend on that assumption either. (See comments in Appendix A.1.) A relevant criterion for having significant quantum effects due to spin can instead be found by comparing the Zeeman energy µB/kBT , i.e. the potential energy of the magnetic dipoles, to the thermal energy kBT . If the parameter µB/kBT is large enough the dipole energy becomes important. Here µ is the magnetic moment of a particle, B is the magnetic field strength, k_B is Boltzmann’s constant and T is the temperature in Kelvin. For further examples of useful parameters, see Ref. [2].

The dynamics of these new regimes can be described by a quantum spin kinetic theory. In this thesis such a theory will be explored, in which a quasiprobability distribution function describes the system. We consider a plasma where the particles of interest (usually electrons) are of spin 1/2 and do not interact directly, but are affected by the average electromagnetic fields created by the other particles – this is the so called mean field approximation of a plasma. Classically, such a system can be described by the Vlasov equation, which together with Maxwell’s equations gives a closed system for the electromagnetic fields, E(r, t) and B(r, t), and the probability distribution function f (r, v, t). We will begin by a short derivation of the Vlasov equation.

The function f (r, v, t) gives the distribution of particles in phase space, so that f (r, v, t)d³rd³v gives the probability of finding a particle in a volume element d³r around position r and d³v around the velocity v. If we imagine a volume Vr in position space and a volume Vv in velocity space, then

N (t) = ˆ

V_r

ˆ

V_v

f (r, v, t)dVrdVv (1.1)

is the number of particles in the volume Vr× Vv in phase space. If we assume that particles can neither be created nor destroyed, and that there are no collisions between the particles, the number N (t) can only change in two ways: by an influx of particles with correct velocities into V_r that we call ∆N_r, or by an influx of particles with correct positions into V_v that we call ∆N_v. In other words,

dN dt =

ˆ

Vr

ˆ

Vv

∂f

∂tdV_rdV_v= ∆N_r+ ∆N_v, (1.2) where ∆N_r and ∆N_v can be described by surface integrals of the flux over the surfaces S_r and S_v that surround the volumes V_r and V_v. These surface integrals can be turned into volume integrals

(6)

by applying Gauss theorem, and we can also use that ^∂r_∂t = v and ^∂v_∂t = a = (q/m)(E + v × B).

The change of particles is then described by

∆N_r+ ∆N_v= − ˆ

V_v

˛

S_r

f∂r

∂t · dSrdV_v− ˆ

V_r

˛

S_v

f∂v

∂t · dSvdV_r

= − ˆ

Vr

ˆ

Vv

∇r· (f v)dVrdVv− ˆ

Vr

ˆ

Vv

∇v·h f q

m(E + v × B)i

dVrdVv. (1.3) In the first integral, ∇r will not act on v and we can write ∇r· (f v) = v · ∇rf . In the last integral

∇v will act on v, but we still get ∇v· [f (q/m)(E + v × B)] = [(q/m)(E + v × B)] · ∇vf since the fields do not depend on velocity and the divergence of v × B is zero. Since the volumes Vr

and Vv are arbitrary, we can conclude that the integrands must be equal in eq. (1.2) and eq. (1.3).

Collecting all integrands on the left hand side we get the Vlasov equation:

∂

∂tf (r, v, t) + v · ∇_rf (r, v, t) + q

m[E(r, t) + v × B(r, t)] · ∇_vf (r, v, t) = 0. (1.4) For other approaches to the derivation of the Vlasov equation see Ref. [3], in which the classical theory of plasma physics is described.

While the intrinsic spin of the electrons is a quantum mechanical feature, it can be introduced semiclassically into eq. (1.4) by considering a volume Vs in spin space as well. This amounts to adding a term ^∂s_∂t· ∇s, where ^∂s_∂t= ^2µ

~s × B. We also need to include the dipole force _m^µ∇r(s · B) in ^∂v_∂t as seen in Ref. [4]:

∂

∂tf + v · ∇_rf + hq

m(E + v × B) + µ

m∇r(s · B)i

· ∇vf + 2µ

~

(s × B) · ∇_sf = 0. (1.5) Here the unit vector s gives the direction of the spin and the magnetic moment µ is for electrons given by the signed Bohr magneton as µ = −µB = e~/2m^e. However, it turns out that treating the system quantum mechanically instead of semiclassically gives rise to correction terms to eq. (1.5) and to new effects such as the Stern-Gerlach effect [5]. Such effects can become important in cases characterized by different parameters as described above. Treating the system in terms of quantum mechanical states and operators has the drawback that direct comparison with the classical theory can be difficult. To deal with this difficulty there are different ways of bridging the gap between quantum and classical, one of which is to use transforms on the density matrix used in quantum mechanics (the density matrix will be introduced in Section 2) to get a quasiprobability distribution function similar to the probability distribution function seen in the equations above.

Apart from easy comparision, there is also the advantage of being able to use methods developed for the classical case, e.g. numerical schemes for simulation of the system.

In order to simplify the calculations we here deal with plasmas in the so called long scale length limit, in which the scale lengths of the physics investigated are much longer than the de Broglie wavelength of the electrons. The scale length L of a physical quantity such as the electric potential Φ is defined from Φ⁰/Φ ∼ 1/L, with Φ⁰ ≡^∂Φ_∂r, and is a measure of at which length scale potentials, field strengths or other relevant values change significantly. Meanwhile, the de Broglie wavelength is defined as λdB≡ h/mvT, where vT is defined from v²_T/2 ≡ kBT , and gives an indication of the length scale needed to notice the wave nature of the electrons. The concepts scale length and de Broglie wavelength are illustrated in Fig. 1.1. If the scale length is much longer than the de Broglie wavelength, the potentials et cetera are approximately constant over the average wavelength in the wave packet describing the electron. From the definitions above, the long scale length limit can be written formally as

h/mv_T L. (1.6)

(7)

Potential or field strength Wave packet

Scale length

de Broglie wavelength

Figure 1.1: Illustration of the concepts scale length and de Broglie wavelength.

This limit will motivate taking corrections only to first order in the evolution equation derived in Section 3. We also assume that the velocities are low enough to be non-relativistic, and that our system is close to thermal equilibrium.

The structure of this thesis is as follows: After some comments on notation and necessary theory in Section 2, a modified evolution equation similar to eq. (1.5) will be derived in Section 3 with first order corrections due to quantum mechanical effects. So far, the results are known from before, see e.g. Ref. [6]. In Section 4 the known theory is expanded upon by investigating what form the local conservation laws for momentum and angular momentum will take with first order corrections to the evolution equation. Here, two different expressions for electromagnetic momentum are considered, given by Abraham and Minkowski respectively. There is some uncertainty over which of these two expressions is the most natural, see for example Ref. [7], Ref. [8]. While this matter cannot be resolved by one example it is still of interest to determine which expression is best suited for achieving total conservation of momentum in the specific case considered in this thesis. Finally, the results are briefly summarized and discussed in Section 5. For further details on some of the calculations there is also an appendix.

(8)

(9)

2 Preliminaries

2.1 Notation

In several of the calculations index notation will be used. Unlike in the covariant formulation used in general relativity, only down indices will be used, the difference between up and down indices having no significance in 3-space. The exception to the usage of down indices is the vector of Pauli matrices (see below), where up indices are used to differentiate between the three Pauli matrices and down indices are reserved for specifying the components (row and column) of the matrices. We also differentiate between different physical quantities by writing for instance P^P for momentum associated with particles compared to P^EMfor momentum associated with electromagnetic fields.

For the indices, Latin letters i, j, k, ... are used for coordinates in phase space as in for example ri and vi, i = x, y, z, while Greek letters α, β, ... are used for describing spin up or down in the z-direction as in for example ψα, α = 1, 2.

Together with index notation we use Einstein’s summation convention. The summation convention means that an index that is repeated twice in a product implies a summation over said index, even if no summation sign is written out. For example, a_ib_i meansPz

i=xa_ib_i. When doing calculations in index notation we will also frequently use the Kronecker delta δ_ij and the Levi-Civita symbol _ijk which are defined as:

δij =

(1 if i = j

0 if i 6= j (2.1)

and

ijk=







1 if (i, j, k) is an even permutation of (x, y, z)

−1 if (i, j, k) is an odd permutation of (x, y, z) 0 if (i, j, k) is not a permutation of (x, y, z).

(2.2)

As an example of their usefulness, the ith component of the curl of a vector ∇r×a(r) is ijk ∂

∂r_jak(r).

When working with spin we will use the vector of Pauli matrices

←→σ =←→σ^(x), ←→σ^(y), ←→σ^(z)

=0 1 1 0

,0 −i i 0

,1 0 0 −1

. (2.3)

In index notation, a component of a Pauli matrix is referred to as σ⁽ⁱ⁾_αβ, e.g. σ₁₂^(y) = −i. We will often use a unit vector s to indicate possible spin directions. In spherical coordinates, s = (sin θ cos φ, sin θ sin φ, cos θ) and integrating over spin directions means taking an integral on the form ´

d²s = ´

dθdφ sin θ. When integrating over both spin and velocity, we use the notation

´dΩ =´

d²sd³v =´

dθdφ sin θ dvxdvydvz.

Finally, we use arrows on operators to indicate what direction they act in. As an example,

A

←−∂

∂xB

−−→∂

∂v_yC = ∂A

∂xB∂C

∂v_y. (2.4)

Operators can appear as arguments of a function such as cos(_∂x^∂ ) which is then defined by its Taylor expansion:

cos ∂

∂x

= 1 − 1 2

∂²

∂x² + .... (2.5)

(10)

2.2 Density matrix, reduced density matrix

Quantum mechanics is usually introduced with examples where the exact state of the system is known. In a some cases, however, the exact state of the system is not known, e.g. when dealing with unpolarized light or thermal distributions. (See further Ref. [9].) Instead, one can consider a statistical ensemble of systems each characterized by a different state |ψii with a corresponding probability pi. The probabilities fulfilP

ipi= 1 and could as an example be given by Boltzmann factors if the system is a thermal distribution.

When dealing with pure ensembles the wave function completely describes the state of the system. For N particles in one spatial dimension it is written ψ(x1, x2, ..., xN, t). It follows the Schr¨odinger equation

i~∂

∂tψ(x₁, x₂, ..., x_N, t) = ˆHψ(x₁, x₂, ..., x_N, t) (2.6) for a given Hamiltonian operator ˆH that may depend on position, spin or other relevant quantities.

In a statistical ensemble one describes the state of the system with the density matrix ρ, which is defined as

ρ(x1, x2, ..., xN, y1, y2, ..., yN, t) ≡X

i

piψi(x1, x2, ..., xN, t)ψ^∗_i(y, y1, y2, ..., yN, t), (2.7)

instead of a wave function. It should be noted that x and y are not coordinates along different axes in a two dimensional system; we are still considering only one spatial dimension where xi and yi correspond to two possible positions for the ith particle. In the simple case of one particle the density matrix is defined as

ρ(x, y, t) =X

i

piψi(x, t)ψ_i^∗(y, t). (2.8) The density matrix can be seen as an infinite continuous matrix with x corresponding to row, y to column. It is quadratic in the wave function and is as such related to the probability to find electrons in a certain region given a specific state |ψ_ii, but it also takes into account the probability to have that state in the first place. In fact, from the definition it is easy to see that the the diagonal elements can be written as

ρ(x, x, t) =X

i

pi|ψi(x, t)|²= n(x, t), (2.9)

where n(x, t) it the probability density.

We can also define the density operator ˆρ =P

ipi|ψiihψi|. With this definition, we can write the density matrix as ρ(x, y, t) = hx| ˆρ|yi =P

ipihx|ψiihψi|yi. For a Hermitian operator ˆO(ˆx, ˆp) the expectation value of the corresponding observable O is given by

hOi = Tr( ˆO ˆρ) = ˆ

dx hx| ˆO ˆρ|xi = ˆ

dx X

i

p_iψ_i^∗(x, t) ˆOψ_i(x, t). (2.10)

For further motivation of the definition of the density matrix along with some useful properties, see Ref. [9].

We now consider the time evolution of the density matrix in eq. (2.8). With only an electric field present, the potential energy is given by V (x, t) = qΦ(x, t) resulting in a Hamiltonian of the form

H(x, t) =ˆ pˆ²

2m+ qΦ(x, t) (2.11)

(11)

where ˆp = −i~∂x^∂ . Since the sum over i in eq. (2.8) is a linear function we can consider it term by term or, for notational simplicity, consider a density matrix where all probabilities pi are zero except for one. With the Hamiltonian in eq. (2.11), together with eq. (2.6) and its conjugate, we get the evolution of the density matrix as

i~∂

∂tρ(x, y, t) =

i~∂

∂tψ(x, t)

ψ^∗(y, t) + ψ(x, t)

i~∂

∂tψ^∗(y, t)

=

−~² 2m

∂²

∂x² + qΦ(x, t)

ψ(x, t)ψ^∗(y, t) + ψ(x, t)

~² 2m

∂²

∂y² − qΦ(y, t)

ψ^∗(y, t)

= −~² 2m

∂²

∂x² − ∂²

∂y²

ψ(x, t)ψ^∗(y, t) + q [Φ(x, t) − Φ(y, t)] ψ(x, t)ψ^∗(y, t)

= −~² 2m

∂²

∂x² − ∂²

∂y²

ρ(x, y, t) + q [Φ(x, t) − Φ(y, t)] ρ(x, y, t)

(2.12) As a side note, this is a special case of the general von Neumann equation:

i~∂

∂tρ = [ ˆˆ H, ˆρ], (2.13)

where the brackets denote a commutator.

The density matrix in eq. (2.8) can be written out in three spatial dimensions by letting x → x,

∂

∂x → ∇xand so forth. One can also consider spin degrees of freedom by letting the wave function go over into a spinor

ψ(x, α) = ψα(x) =ψ1(x) ψ2(x)

(2.14) where the Greek letters denote spin states and can take the values 1, 2 corresponding to spin up/down along the z axis. We then have a 2 × 2 matrix as the density matrix:

←→ρ (x, y) =ρ₁₁(x, y) ρ12(x, y) ρ21(x, y) ρ22(x, y)

=X

i

pi

ψ_i1(x)ψi∗

1(y) ψ_i1(x)ψi∗ 2(y) ψ_i2(x)ψi∗

1(y) ψ_i2(x)ψi∗ 2(y)

. (2.15)

The time evolution for such a density matrix, along with the relevant Hamiltonian for our further purposes, will be discussed in Section 3.

A plasma is an example of a system for which a density matrix may be suitable. We will here assume that the particles of interest are of the same type and do not directly interact with each other - the interaction is taken into account through the fields generated by the charge density and current density caused by the particles. Since the particles in question are usually electrons, the terms particle and electron will be used interchangeably. Due to the sheer number of electrons in a typical plasma, describing the system with an N particle density matrix is not feasible, so one uses a so called reduced density matrix. This approach is similar to what is done for reduced distribution functions in Ref. [3]. Roughly, it is accomplished by using the evolution equation for an N particle density matrix, writing it as a so called BBGKY hierarchy and making certain approximations. The exact details are beyond the scope of this thesis but the interested reader can consult Ref. [10].

Compared to the full density matrix ρ(x₁, x₂, ..., x_N, y₁, y₂, ..., y_N, t) describing all particles, our reduced density matrix ρ(x, y, t) describes ”any particle”. While n(x₁, x₂, ..., x_N, t) is the probability density for having particle 1 at x₁ and particle 2 at x₂ ... and particle N at x_N, the probability density n(x, t) given by the reduced density matrix gives the probability of finding any of the particles at x regardless of where the other particles are.

(12)

The reduced density matrix for a many particle system should not be confused with the density matrix for a one particle system; while in the one particle case we must have´

dxdp n(x, t) = 1, the correct normalization for a system of N particles is´

dxdp n(x, t) = N . From now on we will for convenience use the term density matrix which should be understood to mean the reduced density matrix; it should be clear from the context that we are dealing with a many particle system. The reduced density matrix will follow the same evolution equation as the one particle density matrix, see eq. (2.12), the only difference being that Φ is now the potential created by all particles in the system. This in turn makes eq. (2.12) together with Maxwell’s equations a self consistent system of equations which together with initial conditions contains all information needed for describing the behaviour of the plasma within the approximations made, similarly to the classical case in Ref. [3].

2.3 Wigner transform

To be able to compare with and use methods from classical kinetic theory we would like to transform the density matrix into something similar to a classical probability distribution function. We define the Wigner transform [11]

W (r, p, t) = ˆ dλ

2π~e^−ipλ/~ρ r + λ

2, r −λ 2, t

(2.16) This can be seen as a coordinate change from x, y to r, λ,

x = r +λ

2, y = r −λ

2 ↔ r = x + y

2 , λ = x − y, (2.17)

and then a Fourier transform over λ. From this it is easy to see that it is invertible.

The Wigner function W should, to be a perfect analogy to classical theory, have the usual properties of a probability distribution function such as being real-valued and non-negative. It is, in fact, a quasiprobability distribution function and, while always being real-valued, it can take on negative values. This seems less strange when one takes into account the Heisenberg uncertainty principle forbidding exact simultaneous knowledge of some observables; for integration over large enough regions in phase space, the value will be non-negative. Furthermore, the expectation values will always be given correctly. [12]

An important property that will be used in the calculations is that ˆ

dp W (r, p, t) = ˆ

dp ˆ dλ

2π~e^−ipλ/~ρ r +λ

2, r −λ 2, t

= ˆ dλ

2π~ρ r +λ

2, r −λ 2, tˆ

dp e^−ipλ/~= ˆ dλ

2π~ρ r +λ

2, r −λ 2, t

[2π~δ(λ)]

= ρ(r, r, t) = n(r, t), (2.18) where the delta function δ(λ) picks out the integrand evaluated at λ = 0. That is, integration over all possible values of the momentum gives the spatial density function. In the same way,

n(p, t) = ˆ

dr W (r, p, t) (2.19)

gives the density in momentum space [13].

The expectation value of an operator ˆA(ˆx, ˆp), corresponding to some observable quantity, can be calculated from the Wigner function. It is then in general necessary to write the operator in the so called Weyl ordering, where the operators ˆx and ˆp are arranged symmetrically using the

(13)

commutation relation. With this ordering one can let ˆx → x and ˆp → p to obtain a function A(x, p) which can be integrated together with the Wigner function to obtain the expectation value [14]:

hAi(t) = ˆ

drdp A(r, p)W (r, p, t). (2.20)

In many situations, such as when the the quantities of interest can be computed without using products of ˆx and ˆp, Weyl ordering is not necessary. This is the case in this thesis, where the source terms to Maxwell’s equations will be computed from equations similar to eq. (2.21) and eq. (2.22) below. Furthermore, the integrals will not be taken over d³r; this way we get the expectation values as a function of position and also introduce a factor of density due to the Wigner function.

Classically the Wigner function is often given in terms of r, v, t instead of r, p, t; we will make this change of variables in the end of Section 3. Written in this form we can introduce two important quantities that will appear later on: the charge density ρ_f(r, t) and the current density J_f(r, t).

ρf(r, t) = q ˆ

dv W (r, v, t) (2.21)

and

J_f(r, t) = q ˆ

dv vW (r, v, t). (2.22)

We would like to look at the evolution of the Wigner function in eq. (2.16), when having a Hamiltonian on the form seen in eq. (2.11). To do this, we first need the derivatives of the variables x and y written as the derivatives of the variables r and λ as given in eq. (2.17):

∂

∂x = ∂r

∂x

∂

∂r+∂λ

∂x

∂

∂λ = 1 2

∂

∂r+ ∂

∂λ (2.23)

and, similarly:

∂

∂y = 1 2

∂

∂r− ∂

∂λ. (2.24)

From this, one finds that

∂²

∂x² − ∂²

∂y² =1 4

∂²

∂r² + ∂²

∂r∂λ+ ∂²

∂λ²−1 4

∂²

∂r² + ∂²

∂r∂λ− ∂²

∂λ² = 2 ∂²

∂r∂λ (2.25)

so that eq. (2.12) turns into i~∂

∂tρ r + λ

2, r −λ 2, t

= −~² m

∂²

∂r∂λρ + q

Φ

r +λ λ

− Φ r − λ

λ

ρ (2.26)

and we are ready to write out the evolution of the Wigner function:

i~∂

∂tf (r, p, t) = i~∂

∂t ˆ dλ

2π~e^−ipλ/~ρ = ˆ dλ

2π~e^−ipλ/~i~∂ρ

∂t

= ˆ dλ

2π~e^−ipλ/~n

−~² m

∂²

∂r∂λ

| {z }

(i)

+ q

Φ

r +λ 2

− Φ r −λ

λ

| {z }

(ii)

o

ρ. (2.27)

In (i), understood to be the integral with only the term marked with (i) as integrand, we can just take out the r derivative since it is independent of the integration variable λ:

(i) = − ˆ dλ

2π~e^−ipλ/~~² m

∂²

∂r∂λρ = −~² m

∂

∂r ˆ dλ

2π~e^−ipλ/~ ∂

∂λρ. (2.28)

(14)

This can be integrated by parts, noting that ρ → 0 as λ → ±∞ since the system is finite:

(i) = −~² m

∂

∂r ˆ dλ

2π~e^−ipλ/~ ∂

∂λρ =~² m

∂

∂r ˆ dλ

2π~

∂

∂λe^−ipλ/~ρ

= −ip

~

~² m

∂

∂r ˆ dλ

2π~e^−ipλ/~ρ = −i~p m

∂

∂rW (r, p, t) (2.29) For the (ii) term it is advantageous to rewrite the potential in order to get rid of λ. Since

∂

∂pe^−ipλ/~= −^iλ

~e^−ipλ/~ and nothing apart from the exponential depends on p one can write:

q

Φ

r +λ 2

− Φ r − λ

λ

e^−ipλ/~ρ = q

Φ

r +i~

2

∂

∂p

− Φ r −i~

2

∂

∂p

e^−ipλ/~ρ. (2.30) The integral then becomes

(ii) = ˆ dλ

2π~q

Φ

r +i~

2

∂

∂p

− Φ r −i~

2

∂

∂p

e^−ipλ/~ρ

= q

Φ

r +i~

2

∂

∂p

− Φ r −i~

2

∂

∂p

ˆ dλ

2π~e^−ipλ/~ρ

= q

Φ

r +i~

2

∂

∂p

− Φ r −i~

2

∂

∂p

W (r, p, t). (2.31) Thus, after a division by i~, eq. (2.27) becomes

∂

∂tW (r, p, t) + p m

∂

∂rW (r, p, t) = q i~

Φ

r +i~

2

∂

∂p

− Φ r −i~

2

∂

∂p

W (r, p, t), (2.32) where Φ(r ±^i~₂ _∂p^∂) is defined by the Taylor expansion of the the function Φ around r:

Φ r ±i~

2

∂

∂p

= Φ(r) +∂Φ

∂r

±i~

2

∂

∂p + 1 2!

∂²Φ

∂r²

±i~

2

2 ∂²

∂p²+ 1 3!

∂³Φ

∂r³

±i~

2

3 ∂³

∂p³ + .... (2.33) As seen by this expansion, there will appear terms such as Φ(←−

∂

∂r~

−→

∂

∂p)ⁿf in the evolution equation.

In the long scale length limit the fields and potentials vary significantly on a scale much longer than the thermal de Broglie wavelength. More exactly, h/mv_T L, where the speed v_T is defined by k_BT = ^mv₂^T² and _L¹ ∼ ^Φ_Φ⁰. It can be shown (see Appendix A.1) that assuming we are close to thermal equilibrium, the ratio between two terms in an expansion such as in eq. (2.33) is small:

Φ(←−

∂

∂r_k~

−−→

∂

∂p_k)ⁿ⁺¹W Φ(←−−

∂

∂x_k~

−−→

∂

∂p_k)ⁿW

∼ h

Lmv_T. (2.34)

Assuming that _Lmv^h

T 1 the ratio between consecutive terms in eq. (2.33) is small; we will keep such terms only up to first order. With this approximation, the time evolution of the Wigner function is:

∂W

∂t + p m

∂W

∂r = q∂Φ

∂r

∂W

∂p (2.35)

Substituting E = −^∂Φ_∂r, p = mv we get the classical Vlasov equation, eq. (1.4), in one dimension with no magnetic field:

(15)

∂

∂tW (r, v, t) + v ∂

∂rW (r, v, t) + q mE ∂

∂vW (r, v, t) = 0. (2.36) In three dimensions, the Wigner transform takes the form

W (r, p, t) =

ˆ d³λ

(2π~)³e^−ip·λ/~ρ r +λ

2, r −λ 2, t

(2.37) and the results for the one dimensional case carry over to the three dimensional case in a natural way. Introducing spin will turn the Wigner function into a Wigner matrix that can be transformed into a new scalar function by the spin transform shown below. The dependence of spin, as well as the presence of a magnetic field, will add terms to the three dimensional evolution equation, compared to eq. (2.35), as will be seen in Section 3.

2.4 Spin transform

As we have seen above, the Wigner transform turns the density matrix ρ(x, y, t) seen in eq. (2.15) into a quasiprobability distribution function W (r, p, t) that can be used to find probability of having a particle in a certain region in coordinate space or momentum space. When dealing with spin degrees of freedom we would also like to have a spin transform that turns the 2 × 2 matrix ρ_αβ(x, y, t) into some scalar function Q(x, y, s, t) that can be used to find the probability of having a particle with spin up in a given direction specified by s = (sin θ cos φ, sin θ sin φ, cos θ). More precisely, Q(x, x, s, t) should give the probability density of finding particles with spin up in s direction at position x. Defining the transform as in Ref. [6],

Q(x, y, s, t) = 1

4π(δαβ+ s · σ)ρβα(x, y, t), (2.38) we will show that it does indeed fulfil this criterion. Like the Wigner transform it is invertible, as shown in Ref. [6]. Taking s = ez to be the unit vector in the z direction, we wish to find the probability of spin up in z-direction. The result is general since the z-direction can be chosen freely.

The probability to find spin up in z-direction is proportional to |ψ1(x)|². We want Q(x, x, ez, t) to give this probability, which means we want Q(x, y, ez, t) = Aψ1(x)ψ^∗₁(y) for some constant A.

The factor of 1/4π in eq. (2.38) is for ensuring normalization, as will be seen later. We treat the two terms on the right hand side of eq. (2.38) one by one, starting with the first. For simplicity we will assume that all pi in the density matrix are zero except for one. With s = ezand suppressing t we get

δαβρβα(x, y) = Trs

ψ1(x)ψ₁^∗(y) ψ1(x)ψ₂^∗(y) ψ2(x)ψ₁^∗(y) ψ2(x)ψ₂^∗(y)

= ψ1(x)ψ₁^∗(y) + ψ2(x)ψ^∗₂(y), (2.39) where Trsmeans taking the trace over spin space only. Meanwhile, for the second term, we get

σ^(z)_αβρ_βα(x, y) = Tr_s1 0 0 −1

ψ₁(x)ψ₁^∗(y) ψ₁(x)ψ₂^∗(y) ψ₂(x)ψ₁^∗(y) ψ₂(x)ψ₂^∗(y)

= Tr_s ψ₁(x)ψ₁^∗(y) ψ₁(x)ψ^∗₂(y)

−ψ₂(x)ψ₁^∗(y) −ψ₂(x)ψ^∗₂(y)

= ψ₁(x)ψ^∗₁(y) − ψ₂(x)ψ₂^∗(y). (2.40)

Altogether (δ_αβ+ e_z· σ)ρ_βα(x, x) = 2|ψ₁(x)|² and we only need to check that the normalization is correct. We want´

d²s Q(x, x, s, t) = n(x, t), i.e. the probability to find a particle that has spin

(16)

up in any direction possible is the same as finding any particle at all. As can be seen from the definition of the density matrix, n(x, t) is given by the trace of ρ(x, x, t). When integrating over all possible spin directions we obtain

ˆ

d²s (δαβ+ s · σ)ρβα= 4πδαβρβα= 4πn(x, t) (2.41) since only s depends on the integration variable and will integrate to zero over the whole sphere.

As seen from eq. (2.41) the factor 1/4π in eq. (2.38) ensures proper normalization.

We now want to find an expression for the expectation value of spin in terms of the function Q(x, y, s, t). Similarly to in eq. (2.10) the expectation value for spin is given as a trace; this can also be seen from the definition of the density matrix in eq. (2.15). To get the total average of spin one should take the total trace over both position space (as in 2.10) and spin space (as in 2.42). However, we are here interested in the local value for a given position r, which means taking only a partial trace over spin space. This will introduce a factor of density as well, motivating the definition of the density of spin S as

S(x, t) ≡ Trs(←→ρ (x, x, t)←→σ ) (2.42) We make an ansatz to find the trace in terms of Q by taking

ˆ

d²s s_iQ(x, x, s, t) = ˆ

d²s s_i 1

4π(δ_αβ+ s_jσ_αβ^(j))ρ_βα(x, x, t) = ˆ

d²s s_i 1

4πs_jσ^(j)_αβρ_βα(x, x, t), (2.43) where the last step follows from the fact that when integrating a unit vector s around a sphere, the contributions from opposite sides will cancel. Taking the factors not depending on spin out of the integral, this can be rewritten as

1 4π

σ_αβ^(j)ρβα

ˆ

d²s sisj

= 1 4π

σ^(j)_αβρβα

4πδij

3

=1

3(σ_αβ⁽ⁱ⁾ρβα) = 1

3(σ⁽ⁱ⁾_βαραβ) (2.44) where the relabelling of dummy indices turns our result into ¹₃Trs(←→ρ←→σ ). Comparing with eq. (2.42) the correct expression for the expectation value of spin must be

S(x, t) = 3 ˆ

d²s sQ(x, x, s, t). (2.45)

From this, the magnetization can be obtained by M(x, t) ≡ µS(x, t) = 3µ

ˆ

d²s sQ(x, x, s, t). (2.46)

(17)

3 Spin kinetic theory

When we include spin degrees of freedom, we do not immediately get a scalar function after a Wigner transform of the density matrix as in eq. (2.37). Instead, we get the matrix

Wαβ(r, p, t) =

ˆ d³λ

(2π~)³e^−ip·λ/~ραβ

r +λ

2, r −λ 2, t

. (3.1)

We can combine the spin transform with the Wigner transform to turn the matrix Wαβ(r, p, t) into a scalar function f (r, p, s, t). This scalar function then has the property that integrating over p (or r) gives the probability density to find the particle in a certain position (or with a certain velocity) that is also having spin up in the direction given by s. Such a function is similar to classical distribution functions and s plays a role similar to that of r and p. We get

f (r, p, s, t) = 1

4π(δ_αβ+ s · σ)W_βα (3.2)

which is now a function in the 8-dimensional extended phase space that includes spin degrees of freedom. To make comparison with the classical distribution function easier we will later change variables to get f (r, v, s, t) instead of f (r, p, s, t). The expectation values shown in eq. (2.21), eq. (2.22), eq. (2.45) and eq. (2.46) then take the form

ρ_f(r, t) = q ˆ

dΩ f (r, v, s, t), (3.3)

Jf(r, t) = q ˆ

dΩ vf (r, v, s, t), (3.4)

S(r, t) = 3 ˆ

dΩ sf (r, v, s, t), (3.5)

M(r, t) = 3µ ˆ

dΩ sf (r, v, s, t). (3.6)

(3.7)

3.1 Evolution equation

Similarly to the one dimensional example shown in eq. (2.35) we want the evolution equation for eq. (3.2). The Hamiltonian operator will now include interaction between the magnetic moments of the electrons and the magnetic field. Furthermore, the vector potential must be included in the expression for kinetic energy due to the difference between kinetic momentum mv = p − qA(x, t) and canonical momentum p. Altogether, we have [4]

H =ˆ 1

2m[ˆp − qA(ˆx, t)]²+ qΦ(ˆx, t) − µB(ˆx, t) · σ, (3.8) where A is the vector potential, ˆx = x and ˆp = −i~∇^x. In the Coulomb gauge, ∇x· A = 0 and the Hamiltonian can be rewritten as

H = −ˆ ~²

2m∇²_x+i~q

mA · ∇_x+ q²

2mA²+ qΦ(ˆx, t) − µB(ˆx, t) · σ. (3.9) As in eq. (2.12) we first find the evolution equation for the density matrix by using the Schr¨odinger equation i~∂t^∂ψ(x, α) = ˆHψ(x, α) together with its conjugate. This procedure follows the same

(18)

pattern as in the one dimensional case, albeit being more tedious in three dimensions. For the last term, note that for the Pauli matrices are Hermitian, i.e. σ_αβ^∗ = σβα.

i~∂ραβ

∂t = −~²

2m[∇²_x− ∇²_y]ραβ+i~q

m [A(x, t) · ∇x+ A(y, t) · ∇y]ραβ

+ q²

2m[A²(x, t) − A²(y, t)]ραβ+ q[Φ(x, t) − Φ(y, t)]ραβ− µ[B(x, t) · σαγργβ− B(y, t) · σγβραγ]

| {z }

(i)

(3.10) The next step is to use the Wigner transform and the spin transform. Here, the treatment of the term (i) in eq. (3.10) deserves some elaboration, since it has no counterpart in the one dimensional example. For this purpose, some helpful notes are included in Appendix A.2. To simplify notation, let r^±≡ r ±^i~₂∇p. After applying both transforms we get

∂

∂t + 1 mp · ∇r

f =

(1

i~Φ(r⁺) − Φ(r⁻) +1 2

q

m[A(r⁺) + A(r⁻)] · ∇r

+ 1 i~

q²

2m[A²(r⁺) − A²(r⁻)] − 1 i~

q

mp ·A(r⁺) − A(r⁻)

− 1

i~µB(r⁺) − B(r⁻) · (s + ∇s) −1

~

µ s ×B(r⁺) + B(r⁻) · ∇s

)

f (3.11) where we have also divided through with i~. As in the one dimensional example, the long scale length limit allows us to do a Taylor expansion to first order. The three dimensional equivalent to eq. (2.34) is

Φ(←−

∇r· ~−→

∇p)ⁿ⁺¹f Φ(←−

∇r· ~−→

∇p)ⁿf

∼ h

LmvT

1 (3.12)

and is motivated in Appendix A.1. Ignoring higher order terms we obtain

∂

∂t + 1 mp · ∇r

f (r, p, t) = (

qΦ(r, t)(←−

∇r·−→

∇p) + q

mA(r, t) ·−→

∇r

+ q²

2mA²(r, t)(←−

∇r·−→

∇p) − q

m[p · A(r, t)](←−

∇r·−→

∇p)

− µ [B(r, t) · (s + ∇s)] · (←−

∇r·−→

∇p) −2

~

µ [s × B(r, t)] ·−→

∇s

)

f (r, p, t). (3.13) We now change variables from p to v in order to make the comparison with the classical Vlasov equation easier. The position remains the same, as does the time, while the velocity is related to the canonical momentum through v = [p − qA(r, t)]/m. In terms of the new variables r, v, t the derivatives become

∂

∂r_i =∂rj

∂r_i

∂

∂r_j +∂vj

∂r_i

∂

∂v_j + ∂t

∂r_i

∂

∂t = ∂

∂r_i − q m

∂Aj

∂r_i

∂

∂v_j,

∂

∂p_i =∂rj

∂p_i

∂

∂r_j +∂vj

∂p_i

∂

∂v_j + ∂t

∂p_i

∂

∂t = 1 m

∂

∂v_i,

∂

∂t =∂rj

∂t

∂

∂r_j +∂vj

∂t

∂

∂v_j +∂t

∂t

∂

∂t = −q m

∂Aj

∂t

∂

∂v_j + ∂

∂t.

(3.14)

(19)

Together with the relations between potentials and fields, E = −∇rΦ −^∂A_∂t, and B = ∇r× A these new variables gives an equation very similar to the semiclassical Vlasov equation with spin:

∂f

∂t + v · ∇rf +hq

m(E + v × B) + µ

m∇x(B · (s + ∇s))i

· ∇vf +2µ

~

(s × B) · ∇sf = 0. (3.15) When comparing eq. (3.15) to eq. (1.5), the only difference is an extra term involving ∇s. As seen in Ref. [5] this extra term is related to the Stern-Gerlach effect.

(20)

(21)

4 Conservation laws

In physics, there are often quantities that are conserved, which means they are neither created or destroyed. The simplest example in plasma kinetics is probably the conservation of the probability density n(r, t), which we will show as an introduction to the concept of local conservation laws in which the change in time is related to the flow of the quantity. Due to the three dimensional equivalent of eq. (2.18), to find _∂t^∂ n(r, t) we just need to integrate the evolution equation for the Wigner function, eq. (3.15), over dΩ = d²sd³v. We then get

∂

∂tn(r, t) + ˆ

dΩn

v · ∇rf +hq

m(E + v × B) + µ

m∇x(B · (s + ∇s))i

· ∇vf +2µ

~ (s × B) · ∇sfo

= 0 (4.1) where only the first term in the integral will survive. Apart from the two terms containing cross products, the rest of the terms can easily be integrated over d³v due to the factor ∇v and since f → 0 as v → ±∞ they are zero. The term involving (v × B) · ∇vf can be integrated by parts;

the boundary term disappears and ∇v· (v × B) = B · (∇v× v) = 0. Finally, the last term in the integral will integrate to zero as shown in Appendix A.3. Taking out the spatial divergence we get

∂

∂tn(r, t) + ∇_r· j(r, t) = 0 (4.2)

where j(r, t) =´

dΩ vf (r, t) can be thought of as the flow of particle current density. Equation (4.2) is the so called continuity equation for particle conservation. If multiplied by the charge q it becomes the continuity equation for charge conservation, relating the charge density to the flow of current density. It is is written in the form the conservation laws below will take, and its interpretation is as follows: the only way in which the probability density can change in a volume is if there is a net flow of probability density in to or out of that volume. No probability density is created or destroyed, it is a conserved quantity.

Another example is conservation of energy, taking into account contributions from particles, fields and interaction between the spin of the particles and the fields. To show what form the energy conservation will take when dealing with a spin-1/2 plasma, we take the non-relativistic limit of eqs. (24)-(26) in Ref. [15] rewritten in SI units and in terms of r, v, t instead of r, p, t. To keep with the notation in Ref. [15] we call the energy density W ; this letter is understood to mean the Wigner function outside of the following three equations. We define the energy density W as

W = 1

20(E²+ c²B²) + ˆ

dΩ mv²

2 − 3µs · B

f, (4.3)

where the first term is the electromagnetic energy and the terms inside the integral are the kinetic energy and the spin dipole energy respectively. We then define the energy flux vector K as

K = 1

µ₀E × (B − M) + ˆ

dΩ mv²

2 + 3µs · B

v, (4.4)

where the first term is the Poynting vector in a magnetized field, which accounts for the flow of electromagnetic energy out of the volume, and terms inside the integral account for the flow of kinetic and spin dipole energy out of the volume. The conservation law takes the form

∂W

∂t + ∇r· K = 0, (4.5)

i.e. the decrease in energy density in a volume is equal to the flow of energy out of the volume.

(22)

4.1 Conservation of momentum

We will now attempt to find an equation similar to eq. (4.2) that describes the conservation of momentum, taking into account momentum from the particles as well as electromagnetic momentum. Here, we run into an uncertainty that has plagued physics for a long time, namely which expression one should use for the electromagnetic momentum. In vacuum, it is well established that the Poynting vector_µ¹

0E × B should be used, but for a medium where the permeability or per- mittivity differs from the vacuum case two different expressions have been proposed by Minkowski and Abraham, respectively. For a very brief summary of the controversy, see Ref. [7]. We will try to derive an equation for momentum conservation using first one expression and then the other, to see which one seems best suited for the situation.

We define momentum as

P(r, t) = P^P(r, t) + P^EM(r, t) = ˆ

dΩ mvf (r, v, t) + P^EM(r, t), (4.6) where P^EM is either D × B as given by Minkowski or _c¹2E × H as given by Abraham. The time derivative will then take the shape

∂P

∂t = ˆ

dΩ mv∂f

∂t + ∂

∂t(P^EM) (4.7)

and we wish to rewrite it in the form of a local conservation law:

∂Pi

∂t + ∂

∂rj

Tij= 0, (4.8)

where Tij is a stress tensor.

4.1.1 Particle momentum

In order to define the part of the stress tensor in eq. (4.8) that will be associated with the particle motion, we start by considering a situation where no electromagnetic fields are present. In this case only the first two terms in eq. (3.15) will be relevant. Upon performing integration by parts over d³v we obtain

∂Pi

∂t = ∂P_i^P

∂t = ˆ

dΩ mvi

∂f

∂t = − ˆ

dΩ mvivj

∂

∂rj

f = − ∂

∂rj

ˆ

dΩ mvivjf (4.9) which motivates defining a tensor

T_ij^P= ˆ

dΩ mv_iv_jf (4.10)

such that _∂r^∂

jT_ij^Prepresents the flow of the particle momentum. We can then write

∂P_i^P

∂t + ∂

∂r_jT_ij^P= 0.

We now introduce nonzero fields. For P^P we then need to take into account some additional terms from the evolution equation and we get

∂P^P

∂t = ˆ

dΩ mv∂f

∂t = − ˆ

dΩ mvv · ∇_rf

− ˆ

dΩ mvnhq m( E

|{z}(i)

+ v × B

| {z }

(ii)

) + µ m∇r

B · (s + ∇_s)

| {z }

(iii)

i· ∇v+2µ

~

(s × B) · ∇_so

f (4.11)

(23)

As is shown in Appendix A.3,´

d²s(s × B) · ∇sf = 0 and the last term vanishes. The first term has already been dealt with in the case of no fields, leaving the three numbered terms. Integration by parts over d³v is used on each of these terms, noting that that f = 0 at the boundary ±∞ in velocity space. We also use that ^∂v_∂vⁱ

j = δij and that jklδjl= 0. Then:

(i)i= − ˆ

dΩ viqEj

∂

∂v_jf = ˆ

dΩ ∂vi

∂v_jqEjf = ˆ

dΩ qEif = Eiρf (4.12) and

(ii)i= − ˆ

dΩ viqjklvjBk

∂

∂vl

f = − ˆ

dΩ qvivjjklBk

∂

∂vl

f

= ˆ

dΩ q ∂vi

∂v_lv_j+∂vj

∂v_lv_i

_jklB_kf = ˆ

dΩ qδ_ilv_j_jklB_kf = ˆ

dΩ qv_j_jkiB_kf. (4.13) From these two terms we thus get a contribution of Eρf +´

dΩ q(v × B)f = Eρf + Jf× B that can be recognized as the force exerted by electromagnetic fields on a volume with charge density ρf and current density Jf. The terms in (iii) depend on the spin and we use that

ˆ

dΩ v∇r[B · (s + ∇s)] · ∇vf = 3 ˆ

dΩ v∇r(B · s) · ∇vf (4.14) which is proved in Appendix A.4. With this,

(iii)i= − ˆ

dΩ vi

µ ∂

∂r_j(B · (s + ∇s))

∂

∂v_jf

= 3 ˆ

dΩ

µ ∂

∂r_j(B · s) ∂vi

∂v_jf = 3 ˆ

dΩ µ∂Bj

∂r_i sjf. (4.15) The factor ^∂B_∂r^j

i depends only on position and time. Taking it outside the integral eq. (4.15) can be rewritten as

∂Bj

∂ri

3 ˆ

dΩ µsjf = Mj

∂Bj

∂ri

, (4.16)

where M = 3´

dΩ µsf is the magnetization due to the spin of the electrons. This term can be seen as the dipole force exerted on the magnetic moment by the magnetic field.

Collecting all the terms, we have

∂P_i^P

∂t +∂T_ij^P

∂xj

= Eiρf+ (Jf× B)i+ Mj

∂Bj

∂xi

(4.17) motivating the definition

Fi≡ Eiρf+ (Jf× B)i+ Mj

∂Bj

∂r_i. (4.18)

The equation can then be interpreted as follows: the momentum of the particles changes due to a total force Fi that is exerted on the particles by the fields. This is what would be expected from a physical point of view, it can be recognized as the sum of the Lorentz force and the magnetic dipole force.

To prove that the total momentum conserved, we need the change in momentum from the fields to generate the force F_i so that the final equation takes the form

∂Pi

∂t +∂Tij

∂r_j = 0 (4.19)

(24)

where Tij = T_ij^P+ T_ij^EM and Pi= P_i^P+ P_i^EM. The electromagnetic momentum must satisfy

∂P_i^EM

∂t +∂T_ij^EM

∂rj

= −F_i. (4.20)

We thus search a tensor with a divergence ^∂T

EM ij

∂rj = −_∂PEM i

∂t + F_i

, trying first P˜^EM= 1

c²E × H (4.21)

as given by Abraham and then

P^EM= D × B (4.22)

as given by Minkowski.

It is interesting to note that the difference between the Abraham momentum and the Minkowski momentum is a term (E × M)/c² that is neither purely electromagnetic momentum nor purely particle momentum. This term can be related to so called hidden momentum. For a discussion on hidden momentum, see Ref. [16].

4.1.2 Electromagnetic momentum as given by Abraham The time derivative of eq. (4.21) is:

∂ ˜P_i^EM

∂t = 1 c²

∂

∂t(E × H) = 1 c²

∂E

∂t × H + 1

c²E ×∂H

∂t (4.23)

Due to the spin we have a magnetization. Relativistically we would also have a polarization density, but since we are treating the system non-relativistically we can ignore it’s contribution.

For a further discussion, see Ref. [15]. We thus have H = _µ¹

0B − M and D = 0E, and we can use Maxwell’s equations to rewrite eq. (4.23):

∂E

∂t = 1

0

∂D

∂t = 1

0

(∇r× H − Jf) (4.24)

and ∂H

∂t = 1 µ0

∂B

∂t −∂M

∂t = − 1 µ0

∇r× E −∂M

∂t . (4.25)

Note that from SI convention, 1/0µ0= c². Adding the force Fi yields

∂ ˜P_i^EM

∂t + F_i=

−µ₀H × (∇_r× H − J_f) + E ×

−₀∇_r× E − 1 c²

∂M

∂t

i

+ Eiρf+ (Jf× B)i+ Mj

∂Bj

∂ri

= − (Jf× B)i+

µ0Jf× M − 1

c²E × ∂M

∂t

i

− [0E × (∇r× E) + µ0H × (∇r× H)]_i + E_iρ_f+ (J_f× B)_i+ M_j∂Bj

∂r_i

= − [₀E × (∇_r× E) − Eρ_f+ µ₀H × (∇_r× H)]_i + M_j∂Bj

∂ri

+

µ₀J_f× M − 1

c²E × ∂M

∂t

i

.

(4.26)

Conservation laws in kinetic theory for spin-1/2 particles