Quantum kinetic theory of plasmas

As is well known, for plasmas of high density and modest temperature, the classical kinetic theory needs to be extended. Such extensions can be based on the Schrödinger Hamiltonian, applying a Wigner transform of the density matrix, in which case the Vlasov equation is replaced by the celebrated Wigner–Moyal equation. Extending the treatment to more complicated models, we investigate aspects such as spin dynamics (based on the Pauli Hamiltonian), exchange effects (using the Hartree–Fock approximation), Landau quantization, and quantum relativistic theory. In the relativistic theory, we first study cases where the field strength is well-beyond Schwinger critical field. Both weakly relativistic theory (gamma factors close to unity) and strongly relativistic theory are investigated, using assumptions that allow for a separation of electron and positron states. Finally, we study the so-called Dirac–Heisenberg–Wigner (DHW) formalism, which is a fully quantum relativistic theory, allowing for field strengths of the order of the Schwinger critical field or even larger. As a result, the quantum kinetic theory is extended to cover phenomena such as Zitterbewegung and electron–positron pair creation. While the focus of this review is on the quantum kinetic models, we illustrate the theories with various applications throughout the manuscript.


Introduction
Historically, the vast majority of plasma physics has been limited to classical (nonquantum) phenomena. With important applications, such as magnetically confined fusion plasmas and space plasmas, where the density is modest and the temperature

The density matrix
To describe an N-particle system, where the initial state is known, we can use a many-particle wave function = ( 1 , … N , t) . Here, i is the position for the i:th particle and | ( 1 , … N ), t| 2 d 3 r 1 … d 3 r N , is the probability of simultaneously measuring the particles at respective positions in the volume d 3 r i around i . Assuming that the particles have mass m and are moving in a potential V = V( 1 , … N ) , the evolution follows the Schrödinger equation: supplemented by the initial state of the system The exact state of a multi-particle system cannot be known in general; instead, the best we can hope for is the knowledge of a statistical distribution of states. To handle this situation, it is then necessary to use the so-called density matrix. It is defined as where p i is the probability of finding the system in state i , and the sum is over all states that the system can be in. From normalization, we require that which can be interpreted that the total probability of finding the system in one of the states i is unity. E.g., in the case where the system is in thermal equilibrium with temperature T, the density matrix is where the diagonal elements yields the density of particles at a given position in space with the normalization (10) ( , t) → ( , t) + ∇Λ( , t) (11) ( , � , t) → e iqΛ( ,t)∕ℏ ( , � , t)e −iqΛ( � ,t)∕ℏ , n( , t) = ( , , t), 4 Page 6 of 58 where N is the total number of particles. The mean-field Hamiltonian formally looks like it is describing a one-particle system, but is modified in that the fields are the self-consistent fields created by all the particles. For example, for an electrostatic plasma, we have where is the self-consistent field, which is relate to the particle density via Poisson's equation Here, we have added a neutralizing and homogeneous charge density −qn 0 . The mean-field approximation is usually applicable in cases where particle correlations can be neglected. E.g., in cases where particle-particle collisions can be neglected. By including further terms in the BBGKY-hierarchy, it is possible to include particle collisions, see e.g., Ref. (Bonitz 2016) In the discussion above, we have also neglected any reference to the particle statistics. Since we are dealing with fermions, one should really take into account the antisymmetry of the wave function. In Sect. 3, we generalize the mean-field approximation to account for this. The resulting approximation is usually called the Hartree-Fock approximation.

The Wigner transformation
A key tool when deriving quantum kinetic equations from the Schrödinger or Pauli Hamiltonian is the Wigner transform, introduced by Eugene Wigner in his famous paper from 1932 (Wigner 1932). Using this transform on the density matrix, a kinetic evolution generalizing the classical Vlasov equation can be derived. This connection with classical theory is helpful for guiding the physical intuition, and the transformed quantities are typically more attractive for practical calculations. Wigner's original approach has been further developed by others, in particular Moyal (1949), and the Wigner equation is sometimes also referred to as the Wigner-Moyal equation. The Wigner transform of the density matrix is defined by We can see that e.g.
gives the local density of particles, and similarly is the momentum distribution of the particles. However, in contrast to the Vlasov distribution, the Wigner distribution can be negative in certain regions. This is attributed to the Heisenberg uncertainty principle, and we should hence be careful in interpreting the Wigner function W( , , t) as a classical phase-space distribution. Irrespective of that, the Wigner function can still be used to calculate macroscopic properties of the system, such as, e.g., the charge and current densities.
If the system is described by the Hamiltonian, (15), we can derive the evolution equation for W using the Schrödinger equation as follows: The evolution equation for the density matrix is given by Using the variable change the equation is cast into We now want to take the Fourier transformation of this and re-identify the Wigner function in each term. To do this, we can use the following trick: where f ( , iℏ∇ ) is defined by its Taylor expansion. Here, the arrow over ⃖�� ∇ indicates that the operator should act to the left on the exponential function only. For the term containing ∇ , we make a partial integration. Doing this, we finally get the equation Page 8 of 58 In the mean-field approximation, the equation above describes a system of particles interacting via the electric potential created by all the particles via where we have assumed that there is a homogeneous, neutralizing background of particles with charge density −qn 0 (note that q = −e < 0 for electrons with our conventions). These two equations together with the appropriate initial conditions yield the dynamics of the system. The Wigner distribution W function can be used to calculate any expectation value in a similar way to how this is done in classical kinetic theory. For example, we have For a general operator Ô = O(̂ ,̂ ) , depending on both the position and momentum variables, we must first put all the position and momentum operators in completely symmetric form using the commutation relations, and then replace all operators with the corresponding phase-space variables. For example, to calculate the expectation value of the operator Ô =xp x , we first write Therefore, the phase-space function corresponding to Ô =xp x is O(x, p x ) = xp x + iℏ . Note that this particular operator is not a physical observable, since it is not Hermitian; we only chose it to illustrate the procedure. It also works in the other direction. Therefore, to find out which quantum mechanical operator corresponds to the phasespace function O(x, p x ) = xp x , we first write the function in symmetric form and then replace the variables with their corresponding operators This correspondence between operators and phase-space functions is called Weyl ordering, see e.g. (Zachos et al. 2005).
As an example of the usefulness of Eq. (24), consider linear waves where the distribution function can be written W( , , t). W( , , t).
(29) W( , , t) = W 0 ( ) + W 1 ( )e ikz−i t , where W 0 ( ) is a spatially homogeneous equilibrium, e.g., a Maxwell-Boltzmann distribution. Similarly, we write the potential as Linearising the evolution equations above and solving for 1 , we get the dispersion relation where the last equation is obtained by a change of variables. In the limit ℏ → 0 , this reduces to the classical Vlasov dispersion relation. This dispersion relation has been investigated in detail in Ref. (Eliasson and Shukla 2009). There, it is shown that for a fully degenerate background distribution ( T = 0 ), the wave particle damping disappears in the Vlasov limit, since the phase-velocity ∕k always exceeds the Fermi velocity. Furthermore, the critical wave number k c , at which the Landau damping sets in, was computed.
Returning to the general case, Eq. (24), it is straightforward to show that this equation reduces to the classical Vlasov equation for long macroscopic scale lengths. In particular, if the potential varies on a scale length L, i.e., ∇ ∼ ∕L , and the characteristic velocity of the system is v, ∇ p W ∼ W∕mv , then in the limit we may keep the first nonvanishing terms in a Taylor expansion of the potential in Eq. (24). We then get i.e., the Vlasov equation, where = −∇ . By keeping further terms in the expansion, we may use this method to derive quantum corrections to an arbitrary order, see, e.g., Ref. (Manfredi 2005).

The gauge-invariant Wigner function
The Hamiltonian describing a charged particle of mass m and charge q ( q = −e for electrons) interacting with a magnetic field is given by where and are the (mean-field) electromagnetic potentials. Note that, here, ̂ = −i�∇ is the canonical momentum operator, and it is related to the velocity via the vector potential in the usual way. Here, we take and to be the mean-field potentials created by all the particles. The density matrix ( , � , t) describing the state of the system is then the reduced density matrix; see Appendix A for a detailed discussion.
Under a gauge transformation, the density matrix changes according to Eq. (11). Due to this, the Wigner function W defined in Eq. (17) will not be gauge-invariant. Furthermore, the momentum variable, will be related to the velocity in a gauge-dependent way. To get something that is more attractive to work with, we will use a modified version of the Wigner transformation, first constructed by Stratonovic (1970). The definition is An important aspect of this transformation is that the momentum variable , is the kinetic momentum related to the velocity via = m . This can be seen by calculating the momentum density where {⋅, ⋅} denotes the anticommutator. The right-hand side can be identified as the (kinetic) momentum density. The evolution equation for the reduced density matrix can be derived in a similar fashion to how Eq. (20) was obtained. We can then obtain the evolution equation for the gauge invariant Wigner function as in the previous subsection. The result is Serimaa et al. (1986) where We note that, since the equation is completely expressed in terms of the electric and magnetic field, it is manifestly gauge-invariant. Also, since the momentum variable involved is the kinetic momentum, we may make the trivial variable change = m and express the Wigner function as W = W( , , t) . In the classical limit, i.e., when Eq. (32) applies, this reduces to the Vlasov equation. By Taylor-expanding the functions Δ ,̃ , and ̃ , it is possible to obtain approximations to arbitrary order in ℏ.
The models we have explored in this section are of importance in their own right, and have been used, for example, to consider quantum dispersive effects, see, e.g., Refs. (Shukla and Eliasson 2011;Vladimirov and Tyshetskiy 2011) and references therein. However, in what follows, we will continue with more elaborate models including additional physical phenomena, such as exchange effects, spin, and relativistic effects.

Exchange effects
Exchange interaction in plasmas follows from electrons being fermions, with a totally antisymmetric wave-function. Since a full many-body wave-function is too difficult to study without approximations, the first step is usually to introduce Slater determinants to construct an antisymmetric many-body wave-function from the single-particle ones. This is done as follows: where 1 , 2 ,... n denote different single-particle wave-functions. However, this construction in itself will not be sufficient to cover exchange effects. Instead, one must go beyond the simplest mean-field description. Covering exchange effects is usually referred to as the Hartee-Fock approximation (MacDonald and Bryant 1987), in contrast to the Hartree approximation, where exchange effects are ignored.
The relative importance of exchange effects in plasmas is proportional to the parameter H 2 = (ℏ p ∕E k ) 2 (Crouseilles et al. 2008;Haas 2021;Zamanian et al. 2013Ekman et al. 2015;, where the characteristic kinetic energy E k is given by E k = k B T F for a degenerate plasma (with T F > T ) and by E k = k B T for the nondegenerate case (with T > T F ). The above suggests that exchange effects are as important as the more basic particle dispersive effects, described by the Wigner-Moyal equation of the previous section. If so, the common approach of including particle dispersive effects through the Wigner equation, but simultaneously neglecting exchange effects would be highly questionable. Fortunately, while the relative importance scales with temperature and density as given by the H-parameter, the overall importance is also proportional to another dimensionless constant that often is much smaller than unity, see, e.g., Ref.  for a more detailed discussion. Hence, the use of the Wigner equation from the preceding section, neglecting exchange effects, can still be a good approximation. What complicates the picture is that the importance of exchange effects is not just dependent on the background plasma parameters (temperature and density), but also of the specific problem under study. Below, we will illustrate this by considering simple examples of high-frequency Langmuir waves and low-frequency ion-acoustic waves.

General electrostatic theory
Since the theory becomes more complicated for the electromagnetic case, we here present a general quantum kinetic approach to exchange effects in the electrostatic limit. The full electromagnetic case will be discussed briefly in the end of Sect. 3.4. Our treatment will follow Ref. (Zamanian et al. 2013), but leaving out some of the technical details. We here consider a completely ionized electron-ion plasma with the particles interacting through a mean-field scalar potential. Quantum effects for the ions will be completely neglected as will effects due to the self-energy and particle correlations (Bonitz 2016). The state of the N-electrons is described by the density operator ̂1 …N (see, for example, Ref. (Bonitz 2016)), and the dynamics is given by the von Neumann equation with the Hamiltonian The last term accounts for the interaction with the electric potential created by the ions. In Appendix A, we derive the mean-field approximation, in the case where we can neglect exchange effects. Here, we generalize that method to take into account the antisymmetry of the wave function for the electrons. We introduce the reduced density operators according to where Tr s+1…N denotes the trace over particles s + 1 to N (i.e., integrating over the position degree of freedom and summing over the spins), and Λ 1…s is the antisymmetrization operator that takes an s-particle state and makes it completely antisymmetric (Boercker and Dufty 1979). We will only need to know that Λ 12 = 1 −P 12 where P 12 interchanges particle 1 and 2, i.e. P 12 ( 1 , 2 ) = ( 2 , 1 ) (see, e.g., Ref. (Bonitz 2016) for further details). The evolution for the one-particle density operator is given by where ĥ 1 =̂ 2 ∕(2m e ) and V 12 = V(̂ 1 −̂ 2 ) = e 2 ∕(4 0 |̂ 1 −̂ 2 |) and ̂1 2 is the twoparticle density operator. The effects of two-particle correlations ĝ 12 can be separated out of the two-particle density operator by writing it in the form 12 =̂1̂2 +ĝ 12 , see, e.g., Ref. (Wang and Cassing 1985). We are interested in the collisionless limit where a mean-field approximation will suffice. This approximation is obtained by neglecting the correlation ĝ 12 (Wang and Cassing 1985). Utilizing this in Eq. (41), we obtain where V 1 = Tr 2V12̂2Λ12 is the Hartree-Fock potential operator. This is a closed system for the one-particle density operator.
To obtain a connection to the classical kinetic theory, we use the Wigner representation (Wigner 1932) of this equation. Using the complete set of states � , ⟩ , where is the position and = 1, 2 is the spin along the axis of quantization, this representation is obtained as where ( , ; � , ) = ⟨ , �̂1� � , ⟩ is the density matrix. This is a slight generalization of Eq. (17), where we trivially include the spin variables. Note, however, that the resulting Wigner function depends on the two spin variables and is hence a 2-by-2 matrix. Writing Eq. (43) first in the position representation and Wigner transforming the result (using, e.g., the method outlined in Sect. 2.2), we obtain where is the total (mean-field and the ionic field) potential and is the Coulomb potential. The left-hand side of Eq. (45) represents the electrostatic limit of the Wigner-Moyal equation, but keeping the spin dependence (as encoded , while the right-hand side is the correction due to exchange effects. This term is nonlocal in phase space and nonlinear in the distribution function. Two steps of our general treatment remain, taking the long-scale limit and averaging over the spin states ( , ) . The first step is straightforward, just expanding the equations in ℏ∇ r ∇ p . The second step takes a little more work, but can be done with the help of the spin transform. The spin transform will be discussed in detail in the next section (Section IV), and the details regarding the procedure for our specific case can be found in Ref. (Zamanian et al. 2013). Thus, we will omit these details here, and proceed directly to the end result, the spin averaged evolution equation (assuming all spin directions equally probable) in the long-scale limit ( �∇ x ∇ p ≪ 1 ). The evolution equation then reads where i r ≡ ∕ r i and analogously for i p , and an arrow above an operator indicates in which direction it acts. We have also used the summation convention, so that a sum over indices occurring twice in a term is understood.

High-frequency Langmuir waves
Equation (48) is derived using a perturbative approach where exchange effects are considered to be small. Thus, there is little reason trying to solve Poisson's equation together with (48) exactly. Instead, we first solve Eq. (48) dropping the right hand side altogether (i.e., we solve the Vlasov equation), and then substitute these solutions into the right-hand side, to evaluate the exchange correction to first order. Even this simplified treatment might be rather difficult, unless the zero-order Vlasov solution is fairly simple. To focus on a case that can be treated analytically to a large degree, we now consider the case of linear Langmuir waves in a homogeneous plasma. The general procedure is as follows: 1. Pick a background distribution function (e.g., a Maxell-Boltzmann or a Fermi-Dirac distribution), linearize the left-hand side Vlasov equation, make a plane wave ansatz, and compute the perturbed distribution function. 2. Linearize the right-hand side exchange term (keep terms where one factor is the linear perturbation, and the other is the background), substitute the distribution functions from the previous step and compute the integrals, and find the exchange correction to the perturbed distribution function. (48) 3. Make a final momentum integration to find the exchange contribution to the charge density, and use this in Poisson's equation to find the exchange correction to the susceptibility.
Treating the ions as immobile, considering a Maxwell-Boltzmann background distribution for the electrons, and performing the steps outlined above, Ref.  was able to derive the following result: Here, cl is the classical electron susceptibility, given by cl = (e 2 ∕ 0 m) ∫ f 0 d 3 v∕ − kv z 2 , and exc is the exchange correction to the susceptibility, given by where the velocity integrations (normalized against the thermal velocity) covers all of velocity space. In general, the last integrals to find exc must be solved numerically. However, the case of most interest is when kv T ∕ ≪ 1 , such that the Landau damping is weak, in which case we can expand the integrals in powers of kv T ∕ . To leading order in kv T ∕ , neglecting the pole contribution associated with Landau damping altogether, Eq. (49) reduces to After correcting a slight numerical error in Eq. (49) of Ref. (Von Roos and Zmuidzinas 1961), we note that the result there is exactly a factor two larger than our result above. The difference is due to that Ref. (Von Roos and Zmuidzinas 1961) does not take into the spin part of the wave-function. We note that it is the full many-body wave-function that should be antisymmetric with respect to particle interchange, not just the spatial part. Ignoring this over-estimates the exchange correction by a factor of two. When noting the negative sign of the exchange term in Eq. (51), it is worth recalling the use of a perturbative treatment in the derivation, which means that the exchange contribution must be a correction. Thus, the second term of the right-hand side cannot change sign. If it did, the Langmuir waves would be unstable, which of course is not physical for a background state in thermodynamic equilibrium. In this context, it is also worth noting that the contribution from electron correlation tends to be of comparable magnitude for Langmuir waves (Crouseilles et al. 2008).
While the scaling with temperature and density is the same as for many other quantum phenomena (proportional to H 2 = ℏ 2 2 p ∕m 2 v 4 T ) , we note that the overall factor also contains the small dimensionless number 1/90 which appears for (49) 1 = cl + exc .
geometrical reasons. Due to this small number, it makes sense to ignore exchange effects at the same time as other quantum effects are kept. However, as we will see below, there is no general principle guaranteeing the relative smallness of exchange effects. Thus, whether or not exchange effects can be ignored for other specific problems, as compared to particle dispersive effects, has to be determined on a case-bycase basis. Moreover, we note that the case of degenerate electrons requires different integrals to be solved. We omit this case here, but point out that the exchange contribution to the Langmuir wave dispersion relation with a fully degenerate Fermi-Dirac background distribution can be found in Ref (Ekman et al. 2015). Before ending this subsection, we point out that although the scaling is not exactly the same, nevertheless, electron correlations often are significant at the same time as exchange effects. A way to cover both strong coupling effects (electron correlations) and exchange effects beyond the mean-field level is the path integral quantum Monte Carlo approach (Dornheim et al. 2018;Hamann et al. 2020;Zhang et al. 2016). In particular, results for the dynamic structure factor, closely related to the Langmuir susceptibility, have been presented in Refs (Dornheim et al. 2018;Hamann et al. 2020).

Low-frequency ion-acoustic waves
The general procedure numbered 1-3 of the previous subsection still applies for ionacoustic waves. However, the concrete integrals that need to be solved depend on whether the electrons are degenerate or nondegenerate (assuming we limit ourselves to a thermodynamic background distribution). Also, the integration will be simplified in the quasi-neutral regime, that applies for frequencies ≪ pi . An additional thing to consider is the Landau damping, that often can be neglected for Langmuir waves, but generally tends to be significant for ion-acoustic waves. Thus, we must use the Landau contours when evaluating the integrals, and keep track of both the real and imaginary parts. The case of cold classical ions and a Maxwell-Boltzmann distribution for electrons was considered by Ref. (Zamanian et al. 2013), and the dispersion relation in the quasi-neutral limit was found to be where c s = m e ∕m i 1∕2 v Te is the classical ion-acoustic velocity and we have introduced the classical electron Landau damping, cl = kc s √ ∕8 √ m e ∕m i , in the cold ion limit (Boyd and Sanderson 2003). The coefficients of the real and imaginary quantum terms (0.8 and 3, respectively) are only approximate, as the final step involves a numerical integration. An important result from this calculation is that the exchange corrections now are proportional to a factor of the order unity times the quantum parameter H 2 , in contrast to the previous case of Langmuir waves. As a consequence, when studying short-scale dynamics, it is not a good approximation to keep the particle dispersive terms through the Wigner-Moyal equation, and simultaneously drop the exchange contribution.
(52) = kc s 1 + 0.8 Next, we turn our attention to the same case as above, but with degenerate electrons. In the quasi-neutral limit, the dispersion relation can then be computed as (see, e.g., Ref. (Ekman et al. 2015) for details] where = √ m e ∕3m i . Apparently, the relative magnitude of the exchange terms is even larger than for the nondegenerate case, and, accordingly, independently of the background distribution, the omission of exchange terms in the quantum regime cannot be justified for ion-acoustic waves.

Exchange effects: final discussion
As illustrated above, due to the nontrivial momentum integrals, a first principal kinetic treatment of exchange effects is complicated. Going beyond the simplest cases, like the ones studied above, would require a major numerical effort. Even so, only perturbative treatments where exchange effects are small would be possible. However, the topic is an important one, as exchange effects are not generally small compared to other quantum effects. A tentative conclusion, suggested by the above findings, would be that the exchange contribution might be limited (to the extent that it could be negligible) for high-frequency phenomena, but not for low-frequency phenomena. However, more studies need to be done to put such a conclusion on a firm ground.
An important generalization of the above treatment is to cover also electromagnetic phenomena. Some key steps toward that goal were taken in Ref. , although no concrete examples were worked out. An important issue to address in this case is the gauge invariance, as the exchange integrals tend to be dependent on the electromagnetic potentials, rather than the electromagnetic fields.
In light of the challenging nature of the kinetic exchange contribution, it would be valuable to have less complex models, allowing for computational progress also beyond the perturbative regime. A quantum hydrodynamical model fulfilling this criterion has been put forward by Ref. (Crouseilles et al. 2008). In addition to exchange effects, this model also covers the effect of electron correlations. Importantly, the hydrodynamical model is easy to use for practical calculation, as the exchange and correlation terms do not add much extra difficulty compared to pressure and particle dispersive terms. However, a drawback is that the model is based on time-independent density functional theory (DFT), and hence, the applicability to dynamical phenomena is uncertain at best. A comparison of the hydrodynamical and kinetic models, to some extent covering regimes not presented above, has been made in Ref. , and the numerical comparisons are shown in Figs. 1 and 2, borrowed with permission from the original source. The overall conclusion is that there is a reasonable agreement between the hydrodynamic and the kinetic models in the long wavelength high-frequency regime, but otherwise, this particular hydrodynamical model tends to largely underestimate the importance of exchange effects.
Page 18 of 58 Note, however, that several improved quantum hydrodynamical models have been put forward in the recent literature; see, e.g., Refs. (Moldabekov et al. 2018;Bonitz et al. 2019;. In this context, let us also point out that for solid-state plasmas in particular, a large number of functionals in DFT have been presented including exchange and correlation effects. After the early foundation of DFT was laid (Kohn and Sham 1965), much progress has been made with results going beyond the local density approximation (LDA), where it is assumed that the exchange and correlation energy at a certain position depends only on the electron number density at the same point. Thus, nowadays, DFT comes in many flavors, e.g., accounting for the different density of the two spin states (local spin density approximation, LSDA), also using functionals that depend on the current density and not just the number The exchange contribution to the susceptibility, kinetic (solid) and hydrodynamic (dashed), for ion-acoustic modes. Note that, in the kinetic case, there is also an imaginary part (dotted). The irregularities in the real part of I( ∕kv F ) are due to the numerical resolution. Reproduced from Ref. , with the permission of AIP Publishing density (current density functional theory), including effects of the local electron density gradient (referred to as the generalized-gradient approximations, GGA), and accounting for a time-dependence (applying time-dependent density functional theory TDDFT rather than just DFT). In particular, open-source software projects such as "quantum espresso" (Giannozzi et al. 2009) have been instrumental for spreading the use of DFT methods in the field of materials physics. For an overview of density functional theory in solid-state plasmas, see, e.g., Ref. (Hasnip et al. 2014).

The Pauli Hamiltonian and Spin dynamics
Next, we will generalize the treatment based on the Schrödinger equation to include the spin dynamics. In principle, the formalism is similar to that of Sect. 2. In particular, the higher order terms in the BBGKY-hierarchy (Bogolyubov 1946) are dropped, and thus, the evolution equation will be derived from the single-particle density matrix, using the von Neumann equation, followed by a Wigner transform. Nevertheless, there is an important difference compared to the previous case, as the density matrix now depends on the spin state. As a result, the single component in the Scrödinger case is replaced by a two by two matrix, where the different components represent different spin states. The changes needed to obtain the kinetic evolution equation are outlined in subsection A, and two different but equivalent systems are presented. Then, in subsection B, we use the model to derive the linear conductivity tensor in a magnetized plasma, generalizing previous results to cover the case of an an-isotropic background distribution.

Derivation of evolution equations
Including the spin dynamics, in the lowest (nonrelativisitic) approximation, we replace the Schrödinger Hamiltonian of Sect. 2 with the Pauli Hamiltonian where here denotes a vector with the 2 × 2 Pauli matrices as components. Following Ref. , we are able to construct a gauge-invariant scalar kinetic theory using a density matrix description for a spin-1/2 particle.
The basis states are � , ⟩ = � ⟩ ⊗ � ⟩ , where � ⟩ is a state with position and � ⟩ is the state with spin-up = 1 or spin-down = 2 . As a starting point for the derivation, we use the spinor state ( , , t) = ⟨ , � ⟩ which fulfill the dynamical equation i� t ( , , t) =Ĥ ( , , t) , with the Hamiltonian (54). The density matrix is now given by where, as before, p i is the probability to have a state i , but we have the additional dependence on spin state. Similarly as before, the von Neumann equation applies, i.e., the evolution equation for the density matrix is Once the density matrix has been defined, we can define the Wigner-Stratonovich transform (Stratonovic 1970) as where the phase is used to ensure gauge invariance of the resulting distribution function. This is the same transformation as in Eq. (35) with the modification that it must to be taken separately for each component of the 2-by-2 density matrix. In principle, Eq. (56) together with Eq. (57) gives us a kinetic evolution equation, and we could be content with this. However, the individual components of W have no clear physical interpretation, and it is desirable to construct functions that can be understood more intuitively. Two ways to do this have been presented in the literature, that we will describe below.
The first way is to make a spin transform [or Q-transform, as it has also been called (Scully and Wódkiewicz 1994)], which creates a single scalar function f from the Wigner matrix W . The attractive feature here is that the charge and current sources in Maxwell's equation can be calculated from a single scalar function f, which plays a role much like the distribution function of classical physics. However, there is a price to pay, since the transform extends the classical phase space to include an extra independent variable, namely the spin. Thus, our scalar function has the functional dependence f = f ( , , , t).
Following this approach, thoroughly discussed in Ref. , we here define a scalar distribution function f ( , , , t) in the extended phase space as: where W is the 2-by-2 matrix with entries W , is a vector of unit length, and tr denotes the trace over the spin indices. The spin variable can be expressed in terms of spherical coordinates = (cos s sin s , sin s sin s , cos s ) . Integrating over momentum we find that the reduced distribution function f ( , , t) gives the probability to find the particle at position with spin-up in the direction of . Similarly, defining we find that f ( , , t) gives the probability to find the particle with momentum with spin-up in the direction of . Moreover, we note that the expectation value for the spin polarization density is now given by where we stress the need for the factor 3. Here, the measure of integration for the spin is d 2 s = sin s d s d s . This follows from the form of the transformation (59) and is needed to compensate for the quantum mechanical smearing of the distribution function in spin space. Equations (56), (57), and (59) In the long-scale limit, for spatial variations much longer than the characteristic de Broglie length, the local approximations ̃ ≈ , etc, apply, and the term ∝ Δ̃ can be dropped altogether. As a result, the equation accounting for spin dynamics, but dropping short-scale physics, is given by Together with Maxwell's equations, Eq. (63) [or Eq. (68) in the long-scale limit)] provides a closed description for the spin dynamics, where the electron charge density is given by and the current density is given by Here, the two-dimensional spin integration is made over the Bloch sphere (naturally represented in spherical spin coordinates), and ŝ is the spin unit vector. As done above, the current density is naturally divided into its free contribution f and the contribution M due to the spin magnetization. The physics of the evolution equation can be understood as follows. The basic effect of short-scale particle dispersion is captured in the nonlocal variables defined in Eqs. (64)-(67), and the same physics is present already without the spin effects of the Pauli Hamiltonian, as described already in Sect. 2. The effects due to the spin that are genuinely new all survive in the long-scale limit, as seen in Eq. (68). The first terms of the equation are familiar from classical physics, and then, we have the effects of the magnetic dipole force, and the last term is the spin precession. To a large degree, Eq. (68) agrees with a semiclassical kinetic theory (Brodin et al. 2008). The main difference is that the magnetic dipole force has a slightly more complicated dependence ∝ ∇ [(̂ + ∇̂ ) ⋅ ] ⋅ ∇ f rather than the simpler semiclassical expression ∝ ∇ [(̂ ) ⋅ ] ⋅ ∇ f . The reason for this difference is discussed in Ref. .
In the next subsection, we will demonstrate how to handle the extra spin dependence in practical calculations. However, it is not necessary to use the Q-transform and introduce spin as an extra independent variable. Instead, we can proceed as Ref. Hurst et al. (2014) and split the Wigner matrix W ab into a vector part and a scalar part.
Specifically, we can split the Wigner matrix W (with entries W ) into the scalar, defined as and the vector defined as With these definitions, the Wigner matrix can be reconstructed as where the Pauli spin matrices building up the vector has been complemented with the unit matrix 0 . With the aid of (71)-(73), the equation for the Wigner matrix can be rewritten in terms of f 0 and , with the result and where the definitions for ̃ , Δ̃ , and ̃ are the same as in (64)-(66). As expected, Eqs. (74) and (75) are completely equivalent to (63). The scalar function f 0 captures the phase-space density, that is, the charge density is computed as whereas the vector gives the spin density (which in the approximation of the Pauli Hamiltonian used here coincides with the magnetization). Thus, the magnetization current is given by and the full current density to be used in Ampere's law, therefore, is Note that in the long-scale limit, the approximations ̃ → , etc., apply, in which case the long-scale version of (74) and (75) agrees with Eq. (68)). Equations (74) and (75) have been used by Ref. Hurst et al. (2014) to derive fluid equation, with the aid of moment expansions. For more applications of Eqs. (74) and (75), see, e.g., Refs. Manfredi et al. (2019); Manfredi and Hurst (2015). After presenting two different (but equivalent) systems derived from the Wigner matrix, one can ask what has been gained. After all, all the physics is already contained in the equations for W ab . However, the coupled equations for the components of W ab do not provide much help of guiding the physical intuition. The individual components do not carry physical meaning themselves, only the complete object does. By contrast, the scalar object f ( , , , t) can be applied much like a classical distribution function, but in an extended phase space. Similarly, f 0 can be viewed much like a classical distribution function, but in this case averaged over the spin state, whereas the spin properties are captured in the vector . As a consequence, the quantities of the two theories can be related by the relations Each of the two formulations has its own advantages. First, the latter formulation [Eqs. (74) and (75)] gives somewhat shorter equations. There is no explicit spin-precession term, and the magnetic dipole force is now a bit simpler [one does not need the operator (̂ + ∇̂ ) ]. As a result, these equations constitute a rather direct extension of the Vlasov equation. More importantly, not having the extra independent variable ̂ makes the latter equations more suitable for a direct numerical approach. On the other hand, Eq. (63) is attractive when aiming for analytical solutions. Since only a scalar function is involved, and the structures of the extra terms are similar to the classical ones, most analytical approaches developed for the Vlasov equation can be adopted directly, with the simple addition of an extra integration over spin space. Moreover, in the case of the long-scale version, Eq. (68), replacing the magnetic dipole force with its semiclassical correspondence, allows the equation to be written in a form consistent with classical PIC-schemes (see Ref. Crouseilles et al. 2021), in which case an efficient and attractive numerical scheme is possible.

The linear conductivity tensor
We will here limit ourselves to the case of long spatial scale length (much longer than the characteristic de Broglie length), in which case the local approximations ̃ → E , ̃ → B , etc., apply. Thus, the evolution equation for f ( , ,̂ , t) is given by Eq. (68), and the current density to be used in Ampere's law is given by (70).
In this section, we only consider the electron contribution to the current density, as the classical ion-contribution can be found in the textbook literature, see, e.g., Ref. (Swanson 2003) found. Below, we will derive the linear conductivity tensor ij in a homogeneous magnetized plasma, defined as J i = ij E j , for such a (78) system, where ij contains both the free and magnetization current densities. With the conductivity tensor known, it is straightforward to find the dispersion relations for arbitrary wave modes. We start by linearizing the kinetic Eq. (68) according to f = f 0 + f 1 and = 0 + 1 , where the subscript 0 denotes an unperturbed quantity and the subscript 1 denotes a perturbation, and we take 0 = B 0̂ . Before proceeding, let us point out a few quantum effects that may be contained already in the unperturbed distribution function: 1. Fermi-Dirac statistics This effect is well known. For ℏ 2 n 2∕3 0 ∕m e k B T much larger than unity, we will have almost complete degeneracy, whereas if the parameter is much smaller than unity, the thermodynamic background distribution can be approximated by a Maxwellian. 2. Landau quantization The quantization of perpendicular energy states becomes important in the regime of very strong magnetic fields, or very low temperatures, when ℏ | | ce | | ∕k B T → 1 , where k B is Boltzmann's constant and ce = −eB 0 ∕m e is the electron cyclotron frequency. 3. Spin splitting: The two spin states, up-and down relative to the magnetic field, have different probability distributions in spin space. As a result, the general time-independent distribution function can be written as where for a time-independent distribution function, F 0± can be arbitrary functions of (v ⟂ , v z ) , and F 0± is normalized, such that ∫ d 3 v F 0± = n 0± with n 0± being the number densities of the spin up/down states, respectively. The positive spin state here means that the spin points in the direction parallel to the magnetic field, which means that the magnetic moment points in the opposite direction. Note that with this definition, the lower energy state is the spin state with negative index, i.e., n 0− > n 0+ in case the background distribution f 0 describes a thermodynamic equilibrium.
The full quantum mechanical expression of the thermodynamic equilibrium background in an external magnetic field is derived in Ref. , and reads where L n denotes Laguerre polynomials of order n and the particle energy E n is given by We note that, in general, the degree of spin polarization in the background (n 0+ − n 0− )∕((n 0+ − n 0− )) is obtained by performing momentum integrations of the expression (79). Next, we follow Ref. Lundin and Brodin (2010) but with a slight generalization, allowing for a background distribution that is not necessarily isotropic, i.e., we also cover wave modes that are subject to Weibel type of instabilities. After linearization, Eq. (68) is written as To proceed, we make a plane wave ansatz of the perturbed parameters according to , etc. Without loss of generality, we define the wavevector as = k ⊥̂ + k ẑ , where ̂ , and ̂ are the unit vectors in the x and z-directions, respectively. We also choose to express the velocity in cylindrical coordinates The result coincides with Ref. Lundin and Brodin (2010), except that in our case is an arbitrary function, which means that the term in the integral ∝ ×̃ 1 survives the integration ∫ 2 0 ...d v . Writing out the dependence on v , s and s explicitly, the integrals can be carried out, and the result substituted into the expression for the current density, which in turn give the conductivity tensor. Apart from the extra term ∝ ×̃ 1 , and the need to avoid some other slight simplifications for an isotropic background, the calculations are similar to those given in Ref. Lundin and Brodin (2010). Therefore, we proceed directly to the final result for the conductivity tensor, which may be written as where is the classical contribution, and the spin contributions are ( )ij comes from the spin (magnetic dipole force) contribution to the free current density, whereas the terms X (sp) ( )ij and Z (sp) ( )ij come from the magnetization current. Since the conductivity tensor (87) contains all plasma currents, the general dispersion relation is obtained in the same way as in the classical case, i.e., the dispersion relation is given by detD ij = 0 , with D ij = ij (1 − k 2 c 2 ∕ 2 ) + k i k j c 2 ∕ 2 − i ij ∕ 0 . Picking the special case of an isotropic distribution, i.e., letting F 0 = F 0 (v 2 ) , it is straightforward to show that the conductivity tensor reduces to the expression derived in Ref. Lundin and Brodin (2010).
Evaluating the dispersion relation is a complicated task in its own right. Generally, the spin contribution tends to be more significant if the background magnetic fields is strong, if the plasma density is high, and if the temperature is modest. There are numerous dimensionless parameters describing this, as discussed more thoroughly in Ref. Lundin and Brodin (2010). Here, we would just like to point out a few more consequences of the conductivity tensor. Since there are certain denominators in the expression (87) proportional to − k z v z − ±( ce − cg ) , for almost perpendicular propagation (negligible k z ), certain spin terms will be magnified if the wave frequency matches the difference in the gyration frequency and the spin-precession frequency, Δ c = (g∕2 − 1) c ≈ 0.00116 c . Spin-induced wave modes with ≈ Δ c have been studied by, e.g., Refs. Brodin et al. (2008); Asenjo et al. (2012). A closely related feature due to these denominators is the spin-induced wave-particle interaction. Even if the plasma parameters are more or less classical, such that the spin terms are small, for classical wave modes with ∼ Δ c , the smallness of the overall coefficients for the spin terms can be compensated by a larger number of resonant particles. Spin-induced damping due to this mechanism have been studied, e.g., by Refs. ; Ekman et al. (2021). For an extended treatment of the electrostatic limit, see Ref. Hussain et al. (2014).
Moreover, we note that since the expression for the conductivity tensor given here allows for an-isotropic distributions, it may be used for studying instabilities of the Weibel type. In a quantum plasma context, such instabilities have been previously been studied by, e.g., Refs. Haas (2008); Rightley and Uzdensky (2018), but without accounting for the spin dynamics. In addition to the classical free energy sources, the theory presented here allows for instabilities driven by a difference in the spin temperature and the kinetic temperature. Figure 3, reprinted with permission from J. Lundin and G. Brodin, Phys. Rev. E 82, 056407 (2010), copyright 2010 by the American Physical Society, shows that a rather small difference between the spin temperature T sp and the kinetic temperature T kin is enough to drive an instability in the absence of dissipation. ij to ij plotted as a function of T sp ∕T kin for different values of ã . A positive value corresponds to damping, while a negative value gives rise to a instability. We refer to Ref. Lundin and Brodin (2010) for the detailed definitions of ã and A ij determining the instability growth rate Finally, we note that in addition to the direct application of linear wave propagation, the conductivity tensor of a plasma plays an indirect but important role also for other phenomena. A particular example is the stopping power of particles in magnetized plasmas. This problem has been thoroughly studied, e.g., in the book by Ref. Nersisyan et al. (2007), including theoretical development of the dielectric tensor for both classical and quantum plasmas.

Conservation laws, Landau quantization, and nonlinear effects
Up to now, we have aimed for a clear logical structure of the review article, building up quantum kinetic theory by starting from elementary models, gradually progressing to more advanced ones. In this section, we will leave this route to some extent, and instead offer a smorgasbord of different results, to illustrate the diversity of quantum kinetic theory. In particular, we will consider quantum kinetic conservation laws, Landau quantization in a strong magnetic field, and we will use nonlinear perturbation theory to find the spin contribution to the ponderomotive force. Doing so, to present findings of a more general nature, we will apply models that to some degree extends those presented earlier. However, all of them can be derived similarly to the schemes presented above, i.e., by finding a proper Hamiltonian, using the von Neumann equation for the density matrix, and finally making a Wigner transform (and possibly also a Q-transform). The only additional feature is the need for a Foldy-Wouthuysen transform (Foldy and Wouthuysen 1950;Silenko 2008), to isolate the electron degrees of freedom from the positron degrees of freedom.

Conservation laws
Here, we will use a relativistic quantum kinetic model that comes from separating positive and negative energy solutions of the Dirac equation by means of a Foldy-Wouthuysen (F-W) transformation (Foldy and Wouthuysen 1950;Silenko 2008). Since we are decoupling electrons and positrons, the physical condition of applicability is that pair production is negligible. Quantitatively, the electric field should be much smaller than the critical field, E ≪ E c = m 2 ∕e� and similarly for the magnetic field (we use units with c = 1 in this section). Moreover, the typical scale lengths should be long compared to the Compton wavelength ℏ∕m . Apart from the F-W transformation, the derivation of the model is similar to Sect. 4 and we refer to Ref. Ekman et al. (2017) for the details.
The evolution equation for the scalar Wigner function f is where 2 = 2 + m 2 , B = qℏ∕2m is the Bohr magneton and (88) The system is closed with Maxwell's equations, in units where c = 0 = 0 = 1 , where and are the polarization and magnetization densities, and f and f are the free charge and current densities. The source terms are slightly generalized, as compared to previous models, and given by As the theory is relativistic, the integration element dΩ of Sect. 4 is replaced by its relativistic counterpart, i.e., dΩ = d 3 pd 2 s . It follows from the evolution equation, Eq. (88), that the free charge is conserved, t f + ∇ ⋅ f = 0 , and we interpret: as the function on phase space corresponding to the velocity-it is in fact the Weyl transform of the velocity operator ̂ = i � [Ĥ,̂ ] given by the Heisenberg equation of motion. The spin-dependent term is related to the hidden momentum (Shockley and James 1967;Shockley 1968;Coleman and Van Vleck 1968;Babson et al. 2009) of systems with magnetic moments. Here, we also note an important aspect of the relativistic theory, and that the spin magnetization current is complemented by polarization currents. This is natural, of course, as a magnetic dipole moment in the rest frame of a particle correspond to both a magnetic and electric dipole moment in any other frame.
The total energy density is given by with the corresponding energy flux vector where = − . With these expression, we have a conservation law on divergence form It is straightforward to confirm the above energy conservation law by carrying out a number of partial integrations. Deriving the conservation law for momentum is somewhat more tedious, and we refer to Ref.  for the details. Here, we just present the result, which can be written in standard form in terms of energy-momentum tensors for particles and fields Here, the electromagnetic part of the energy-momentum tensor T EM ij is given by where the usual relations = + and = − applies. Similarly, the electron contribution T e ij to the energy-momentum tensor is given by where the momentum-velocity relation (96) apply. A few things can be noted. First, the corresponding conservation laws for the nonrelativistic model in Sect. 4 can be obtained as expected from the above results, i.e., by dropping terms of higher order in an expansion ∕m (letting = m(1 + p 2 ∕2) , etc.). Second, with an expression for the energy momentum tensor, in principal, we can compute the gravitational source due to quantum relativistic electrons. In practice, this is complicated by the fact that the stress tensor is not necessarily symmetric; see Ref. ) for a more thorough discussion. Finally, we note that the conserved quantities will not be modified due to our neglect of short-scale effects (of the type contained in Eqs. (64)- (67)). The reason is that all the extra terms due to short-scale effects contain higher order derivatives whose contributions to the energy-momentum tensors vanish when integrating over momentum space.

Landau quantization
When the Zeeman energy in an external magnetic field is large, i.e., comparable to the kinetic energy of particles, the phenomenon of Landau quantization becomes crucial. This means that the energy levels for motion perpendicular to a magnetic field are quantized, and also that the energy difference between the spin states is significant. In particular, this tends to occur in astrophysical scenarios, where, in extreme cases, the Zeeman energy may be comparable to the electron rest mass energy. Specifically, this happens in the vicinity of magnetars, where the magnetic field strength can be of the order 10 10 T (Harding and Lai 2006). Here, we will address the regime of fully relativistic Landau quantization, when the effect is most pronounced. A particular difficulty with relativistic Landau quantization, is that we can no longer use the inequality B B ≪ m , (or similarly for the electric field) which plays an important role when separating electron and positron states. In principle, this could be handled using the Dirac-Heisenberg-Wigner formalism (DHW) to be discussed in the next section, but this comes at the price of a considerably more complex theory. To focus solely on the problem of Landau quantization, we here takes a simpler route, and consider a strong constant magnetic field B 0 (allowing for B B 0 ∼ m or even larger), but limit the magnitude of electromagnetic field perturbations well beyond this.
A model focusing on the effect of Landau quantization may still neglect spindynamics, in case the validity condition �k 2 ∕m ≪ 1 is fulfilled. (Here, k and do not necessarily refer to plane waves; instead, they represent characteristic spatial and temporal gradients.) Also, assuming �k∇ p ≪ 1 , the previously studied short-scale effects can be dropped, and the model will be a slight generalization of the relativistic Vlasov equation, extended to account for Landau quantization in a strong (constant) magnetic field B 0 . The kinetic evolution equation derived in Ref. Al-Naseri et al. (2020) can be written as where the main difference to the (relativistic) Vlasov equation lies in the energy expression � ± , which now becomes an operator given by We note that the first two terms inside the root sign just give the classical expression. The next term with a ∓ sign gives the magnetic dipole energy for spin up and spin down, respectively. Accordingly, the upper and lower signs are described by different evolution equations for W + and W − , due to the different energies of the spinup and spin-down particles. The final nonclassical feature comes from the last term under the root sign, which makes the energy an operator instead of just an algebraic expression. We note that the full operator is defined by its Taylor expansion. While one may worry that such a Taylor expansion will not behave well due to the root sign, in fact, for any valid Wigner function, the Taylor series will be well defined and convergent. Finally, to have a closed system, we need the source terms in Maxwell's equations, which are given by summing over the spin-up and spin-down contribution, that is In principle, there are also magnetization currents that can be added to the free sources above. However, such a contribution will be negligible in comparison, provided that the condition given above for neglecting the magnetic dipole force is fulfilled. Next, we want to deduce the thermodynamic equilibrium state in a constant magnetic field. Noting that for a constant magnetic field, both the Dirac equation and the Pauli equation results in electrons obeying a quantum harmonic oscillator equation, we can make a straightforward generalization of the Pauli case . Both for the Pauli and the Dirac equations, the spatial dependence of the wave-function in Cartesian coordinates can be expressed as a Hermite polynomial times a Gaussian function (Melrose and Parle 1983) only the energy eigenvalues for the Landau quantized states are different. Specifically, applying the Dirac theory, the energy of the Landau quantized states becomes where n = 0, 1, 2, … correspond to the different Landau levels for the perpendicular contribution to the kinetic energy, the index ± represents the contribution from the two spin states, and the term proportional to p 2 z gives the continuous dependence on the parallel kinetic energy. Since the Pauli and Dirac equations for individual particle states have the same spatial dependence for the wave function, we can adopt the expression for the Wigner function from Ref.  (based on the Pauli equation) with some relatively minor adjustments.
1. Contrary to Ref. , we have made no Q-transform to introduce an independent spin variable, and thus, the spin dependence of Ref.  reduces to W ± . 2. The Wigner function of Ref.  must be expressed in terms of the momentum, i.e., m(v 2 x + v 2 y )∕2 → (p 2 x + p 2 y )∕2m. 3. The nonrelativistic energy of Ref.  is replaced by the relativistic expression (107) of the Dirac theory.
4. The normalization of the Wigner function must be adopted to fit the present case.
With these changes, the background Wigner function W TB ± for the case of electrons in thermodynamic equilibrium can be written as where n 0 = n 0+ + n 0− = ∫ (W +T + W −T )d 3 p is the electron number density of the plasma, c is the chemical potential, T is the temperature, and L n denotes the Laguerre polynomial of order n.
That the factor n (p ⟂ ) gives us the proper Wigner function for the Landau quantized eigenstates can be confirmed by an independent check. Since the expression (108) contains no dependence on the azimuthal angle in momentum space, we can write when � ± acts on n (p ⟂ ) . Computing � ± n (p ⟂ ) by Taylor-expanding the square root to infinite order, using the properties of the Laguerre polynomials, and then converting the sum back to a square root, it is straightforward to verify the relation where ce = |qB 0 | m is the electron cyclotron frequency, confirming that n (p ⟂ ) generates the proper energy eigenvalues for the perpendicular kinetic energy and the spin degrees of freedom.
While (108) gives the thermodynamic equilibrium expression W TB ± , we note that the plasma background state is not necessarily in thermodynamic equilibrium. Making use of the property (111), we note that the most general time-independent solution W 0± to (103) of physical significance can be written in the form where g n± (p z ) is a function that is normalizable, but otherwise arbitrary, and the number of particles in each Landau quantized eigenstate n n,± obeys the condition Page 36 of 58 Naturally, the expressions for W 0± and W TB ± presented here are of most significance for relativistically strong magnetic fields, when Landau quantization is pronounced. As a consequence, the above formulas will reduce to more well-known expressions when the limit � ce ∕m ≪ 1 is taken. Specifically, Eq. (108) will become a relativistically degenerate Fermi-Dirac distribution in case we let T = 0 and c = E F ≫ � ce , where E F is the Fermi energy. Alternatively, for k B T ≫ E F and k B T ≫ � ce , Eq. (108) reduces to a Synge-Juttner distribution.
To give a concrete illustration, in Fig. 4, reprinted with permission from Al-Naseri et al. (2021) copyright 2021 by the American Physical Society, we show a bar chart for the normalized number density n 0n± ∕n 0 in the different energy states, for a few values of the temperature and magnetic field, under the assumption that the density is low enough for the system to be nondegenerate, i.e., assuming T > T F .
As a result of the background dependence, Eq. (112), in the Landau quantized regime, the electrons behave as a multi-species system, where each particle species has its own rest mass, as given by Eq. (107) but with p z = 0 . This is because the separation between Landau levels is of the order of the rest mass, and all excitations by quanta with energies of that order have been neglected. If we define the effective number density of each "species" (discrete energy-state) as we see that n 0n± to a large degree will be determined by the Boltzmann factors of Eq. (108). For a study of Langmuir waves in a Landau quantized plasma, see Ref. Al-Naseri et al. (2020).
While we have here focused on the extreme case of relativistic Landau quantization, we note that the thermodynamics is much affected also in the nonrelativistic regime, provided that the Zeeman energy is of the same magnitude or larger than the characteristic kinetic energy in the background plasma. For applications of Landau quantization in the nonrelativistic regime, see e.g. (Eliezera et al. 2005;Shaukat 2017)

Nonlinear effects
Not surprisingly, quantum kinetic models can describe a great variety of nonlinear phenomena. Just like for classical plasmas, a fair share of these phenomena is induced by the ponderomotive force. For example, the ponderomotive force is the driver of plasma wake field generation [93], the key mechanism in soliton formation (Shukla et al. 1986), and the main source of nonlinear self-focusing (Kurki-Suonio et al. 1989). While the (113) Fig. 4 The normalized number density at different energy states E m for different values of the parameters = B B 0 ∕m and = k BT ∕m main features of classical and quantum kinetic models driven by the ponderomotive force are similar, nevertheless, there are important differences. For one thing, in a magnetized plasma, the classical ponderomotive force has cyclotron resonances (Karpman and Washimi 1977). In quantum kinetic theory, the classical terms are still present, but they are complemented by terms containing extra quantum resonances Stefan et al. 2011). The concept of a ponderomotive force in quantum kinetic theory is not as straightforward as in a fluid theory. Depending on definition, some low-frequency nonlinearities induced by quadratic nonlinearities may be included as a ponderomotive force term or not. To avoid ambiguities, we will focus on the regime where the phase velocity of the high-frequency wave is much higher than the thermal velocity (or characteristic velocity, in case of degeneracy effects).
As a starting point, we will base our study on the long-scale version of the model put forward in Sect. 4, that is, Eq. (68). Next, we consider circularly polarized electromagnetic waves of high frequency propagating parallel to an external magnetic field, 0 = B 0ẑ , and use the following ansatz: The amplitudes are assumed to vary much slower than the exponential phase factors, and the star denotes complex conjugates. Since the basic wave modes propagating parallel to 0 are either left-or right-circularly polarized, we have ̃ ,̃ ∝ ̂ ± î . Furthermore, all perturbations are small, such that weakly nonlinear perturbation theory is applicable, and we will focus on the ponderomotive contribution, that is the quadratically nonlinear low-frequency terms.
To calculate the weakly nonlinear low-frequency response to an incoming transverse wave packet, we make the ansatz where f 0 is the background distribution, f lf is a low-frequency part due to quadratic nonlinearities, and f 1 is a slowly modulated high-frequency wave. The background distribution will be taken to be of the form where n 0 is the equilibrium density, and the thermal velocity v T is defined as v T = √ 2k B T∕m . It should be stressed that the background distribution (117) has been picked mostly for convenience. In fact, the results below are not depending sensitively on the background, as long as the phase velocity is larger than the (116) characteristic velocity. For example, a degenerate distribution, but with the Fermi velocity smaller than the phase velocity would produce identical results. Given the ansatz (115) and (116), the aim is then to find an equation for the low-frequency part of the distribution function. From such an equation, we can then calculate the low-frequency response in the current density and magnetization, and compare with the results for a ponderomotive force, as defined from a fluid theory. First, we deduce the high-frequency linear perturbation of the distribution function, which is given by f 1 = f + ( f − ) for left-hand (right-hand) circularly polarized waves. The expressions for f ± found from Eq. (68) to linear order can be easily computed from Eq. (9) in Ref. (Lundin and Brodin 2010), and the result is Next, allowing for slow modulations and solving the equation to first order in z ∕k , t ∕ , we note that the zero-order solution applies after making the substitution → + i t and k → k − i z in Eq. (118), and then expanding to first order in the slow derivatives. Inserting the ansatz above into the evolution equation and considering the slow-time scale and keeping only up to quadratic nonlinearities, we obtain the equation where c.c. stands for complex conjugate. Here, we have also added a low-frequency electric field in the z-direction, E zlf , which has f lf as source. Equations (119) and (118) now constitute a basis for calculating the nonlinear response in the current density and magnetization.
After some algebra (described more closely in Ref. Stefan et al. 2011), based on the low-frequency part of Ampere's law, we end up with the final expression for the induced low-frequency field, which can be written in terms of the ponderomotive force f p as where f p can be divided into its classical and its spin contribution according to with the different parts given by Page 42 of 58 The last term is the classical contribution, and the second-to-last term in the square bracket is what is obtained without the weakly relativistic effects. However, if we assume that k∕ is roughly of order unity, we see that all terms in the square brackets are of the same order. This implies that when dealing with an unmagnetized plasma where spin effects are important, the spin-orbit coupling contributions must be taken into account, as well. However, the previous result for the ponderomotive force, Eq. (122), is still relevant, as spin terms tend to be more important in magnetized plasmas, and hence the simpler model (68) can still be justified.

The full Dirac theory: the DHW equations
Quantum relativistic treatments are of interest in several different contexts (Zhang et al. 2020;Elkamash et al. 2017;Shi et al. 2018). Dense astrophysical objects can have a Fermi energy approaching or exceeding the electron rest mass energy, and the strong magnetic fields of magnetars give rise to relativistic Landau quantization. Importantly, the continuous evolution of laser intensity brings a variety of quantum relativistic phenomena accessible to experimentalists. Upcoming laser facilities of interest in this context include, e.g., the extreme light infrastructure (ELI) [105], (Dunne 2009) and the European X-ray free electron laser (XFEL) [107], (Ringwald 2001). The quantum kinetic models of previous sections have all made various simplifications as compared to the full quantum relativistic theory. In particular, to avoid the mixed electron-positron states of the Dirac theory, up to now, we have imposed limitations on the electric field. Specifically, we have demanded E ≪ E cr where E cr = m 2 c 3 ∕eℏ is the critical field, to avoid the complications associated with significant pair production due to the Schwinger mechanism. In this section, however, we will take on the full complexity of the Dirac theory using the so-called Dirac-Heisenberg-Wigner (DHW) formalism (Bialynicki-Birula et al. 1991). Compared to previous sections, new features of the theory include Zitterbewegung (a rapid (speed of light) particle motion associated with interference between positive and negative energy states), vacuum polarization, and electron-positron pair creation. Also, the previous 2-by-2 Wigner matrix will be replaced by 16 components, due to the 4 components of the Dirac spinors. Nevertheless, many other aspects of the quantum kinetic theory will be familiar and we will point out how the DHW formalism can be related to the quantum kinetic approximations presented earlier.

The DHW model
The DHW model was first derived by Ref. (Bialynicki-Birula et al. 1991). Moreover, some relatively minor variations of this derivation have been published in the literature more recently, see e.g. (Hebenstreit et al. 2010;Sheng et al. 2019). As all these treatments are fully satisfactory, we will not repeat the calculations, but just point out a few of the main features.
1. The derivation is based on the Dirac equation, which gives the time evolution of Dirac four spinors, generally describing mixed electron-positron states. 2. Just like in previous theories, a gauge-invariant Wigner transformation is made, which here produce a 4-by-4 Wigner matrix, with the components depending on phase-space variables just like in standard kinetic theory. 3. The main omission is made when using the Hartree approximation where the electromagnetic field is treated as a nonquantized field. This approximation amounts to neglecting the quantum fluctuations. We will come back to the consequences of this approximation. 4. To write the equations in a physically more transparent form, the 16 components of the Wigner matrix are decomposed into 4 different four vectors, which in turn is split into temporal and spatial components. Most of these quantities have a clear physical meaning, which helps forming a physical understanding of the theory. 5. The (phase-space) current density and charge density are parts of the DHW unctions, which makes it straightforward to close the system using Maxwell's equations.
With these preliminaries, we jump directly to the DHW quations, which in units where c = ℏ = 1 can be written in the form Due to the Wigner transform, we cover short-scale quantum phenomena in much the same way as in the previous sections. This is illustrated by the appearance of the nonlocal operators in Eq. (131), which are given by Page 44 of 58 As before, the operators reduce to their local approximations (i.e., D t → ∕ t + e ⋅ ∇ p , → ∇ r + e × ∇ p , ̃ → , and ̃ → ) for scale lengths much longer than the characteristic de Broglie length.
To close the system, we need the source terms in Maxwell's equations, which are given by and Thus, v 0 is the time component of the four-vector phase-space function that gives the four current density, and is the spatial component, i.e.,the current density. Most, if not all, DHW unctions have fairly concrete interpretations which helps guiding the physical intuition. To gain a better understanding, we take a look at some of the conserved quantities of the DHW ystem (for a derivation, see Ref. Bialynicki-Birula et al. 1991). First, the total energy W is given by second, the linear momentum is and, finally, the total angular momentum is We can see in Eq. (138) that the current density can be related to the kinetic energy, as expected. However, the role of mass density is played by another function s, with no trivial relation to the charge density. As the Dirac field contains both electrons and positrons, the lack of a simple relation between the mass density and the charge density should not be surprising. Nevertheless, in the expression for momentum (139), we see that v 0 acts as a phase-space momentum density. Since electrons and positrons contribute with opposite signs to v 0 , we realize that electrons and positrons must have a different dependence on . To be concrete, for electrons and positrons moving in the same direction, we need to shift the momentum dependence for the dependent variables according to → − , as will be illustrated more explicitly below. This is not an issue when solving the DHW quations, as the DHW variables generally describe coupled electron and positron states anyway. However, this insight can be of some importance, e.g., when interpreting results, in particular when numerical calculations have been made. The shift in momentum dependence is consistent with the common interpretation of positrons as being electrons moving backwards in time.
Another observation that can be made based on the conserved quantities is that the term ∝ (1∕2) can be identified as the spin contribution to angular momentum, i.e., gives the spin density. Finally, the equation for D t 1 in Eq. (131) shows a division of the total current density into its free part, magnetization part, and polarization part. To be specific, we can deduce that 2 gives the magnetization, and 1 gives the polarization. For some further discussion of the physical interpretation of the DHW functions, see Ref. (Bialynicki-Birula et al. 1991).
Contrary to previous kinetic theories, in the DHW formalism, the kinetic variables are not zero even in vacuum; in older terminology, we would say that the vacuum is filled with the particles of the Dirac sea. However, in the absence of a spin polarizing magnetic field, the only DHW functions with nonzero vacuum values are the mass density and current density, which are given by where = √ m 2 + p 2 . The expressions above are obtained by calculating the Wigner operator for the free particle Dirac equation and taking the vacuum expectation value. Note that while the charge density is zero (due to the cancellation of the electron and positron vacuum fluctuations), the same is not true for the current density. The reason is that there are two signs that enter the picture-first, electrons and positrons have opposite signs of the charge, but second, switching electrons for positrons means letting → − , such that the vacuum contributions are additive. Nevertheless, when integrating over momentum to get the total current density, we still get zero as one would expect. The substitution → − when switching between electrons and positrons also holds for real particles as well as the vacuum contribution. In particular, for beam systems, this is important to keep in mind, as is illustrated in Fig. 5.
When adding the real particles of the Dirac field (electrons and/or positrons) into the picture, we can add distribution functions much like in the theories of the previous sections. Specifically, we could start from a background with a nonzero charge density The function f e,p ( ) can be picked as any common background distribution function from classical kinetic theory, i.e. a Maxwell-Boltzmann, Synge-Juttner, or Fermi-Dirac distribution, depending on whether the characteristic kinetic energy is relativistic and whether the particles are degenerate.
(141) Note that for a completely degenerate ( T = 0 ) Fermi-Dirac background of electrons (and no positrons f p = 0 ), the electron and vacuum contributions for the mass density and current density cancel inside the Fermi sphere. Consequently, for momenta p ≤ p F , where p F = ℏ(3 2 n 0 ) 1∕3 is the Fermi momentum, we have s = = 0 . Furthermore, note that in the presence of a strong field in the background state (e.g., a strong constant magnetic field), these fairly simple background expressions need to be modified. For example, in a case with B B 0 ∼ mc 2 , also the vacuum states will be subject to Landau quantization. Moreover, the spin density and magnetization (as described by the functions and , respectively) will no longer vanish in the background state, due to the contribution from real particles.

The electrostatic one-dimensional case
To illustrate some features of the DHW theory, we will consider the case of a onedimensional electrostatic field, i.e., = E(z, t)̂ . At the same time, the DHW functions depend on ( , z, t) , but where the momentum dependence can be reduced to two independent variables ( p ⟂ , p z ), due to the rotational symmetry.
Due to the simplified geometry, only half of the 16 scalar DHW functions will be nonzero. Moreover, only four of these variables will be independent. There are different ways of finding these nonzero variables. Here, we will just present the result, which is straightforward to verify by direct substitutions into the DHW equations, Eq. (131). For a more systematic derivation of the reduced electrostatic equations, see Ref. . As it turns out, the DHW equations can be expressed in terms of four variables 1 − 4 , related to the original DHW functions as follows: 2. The main new effect due to the quantum relativistic regime, comes from the new types of resonant denominators, associated with wave damping. In particular, even for k = 0 , we may still have a resonant denominator (corresponding to electron-positron pair-production), provided the pair-creation condition ℏ ≥ 2mc 2 is fulfilled. 3. Due to the vacuum background, the integrand is nonzero even in the absence of particles, giving raise to the effect of vacuum polarization contribution. While this term typically gives a contribution that is much smaller than that from the real particles, the given expression is subject to ultra-violet divergences that must be handled using a renormalization scheme, see, e.g., Ref. Bialynicki-Birula et al. (1991).
Naturally, the above points only give the principal features. A more thorough study of the quantum relativistic dispersion relation for Langmuir waves has to be done numerically. This will be the subject of a future paper.

Limiting cases of the DHW theory
Apart from exchange effects, the DHW equations covers all the physical phenomena presented in the previous sections. Thus, ideally, we should be able to recover all the previous models (except the parts presented in Sect. 3) as special cases of the DHW formalism. However, demonstrating the equivalence in appropriate limiting cases is somewhat nontrivial. First of all, the DHW equations do not only describe electrons, since a Foldy-Wouythuysen transformation cannot be made in the fully quantum relativistic regime. Furthermore, the equations with a spin-dependent Wigner function uses a Q-transform to get a scalar theory, which further complicates a comparison. Nevertheless, though a complete investigation is yet to be made, showing the equivalence of the DHW formalism with the models of Sect. IV (i.e., the nonrelativistic Pauli limit) is relatively straightforward. First, we note that for fields well below the critical field, there is little ambiguity whether we have electron or positron states. Considering the case of electrons only, ignoring relativistic effects, the charge and mass density are the same (due to the normalization, the constant factors involving e and m do not enter), i.e., s = v 0 . As a first exercise, let us recover the model based on the Schrödinger Hamiltonian. This implies dropping the effects of the electron spin, i.e., the spin density, magnetization, and spin polarization are zero, and thus, we let = = = 0 in the DHW equations. As a result of the above approximations, we get =̃ s =̃ v 0 , which immediately lead to a closed equation for v 0 Identifying v 0 with the Wigner function of section 4, using the definitions of the nonlocal variables to write the more explicitly, Eq. (158) gives us Next, our aim is to recover the model of Section IV, based on the Pauli Hamiltonian. As the short-scale physics associated with the nonlocal expressions (132)-(135) has already been established above, for convenience, we restrict ourselves to the case of long-scale lengths ( �∇∇ p ≪ 1 ), such that the local approximations (i.e., dropping the variables with tilde) can be used. For the nonrelativistic case, we can still put the charge and mass density equal, i.e., s = v 0 . Second, in the absence of relativistic effects (and for the given normalizations) and for the case of electrons only, the spin density and the magnetization are the same, i.e., = 2 . Third, a nonvanishing polarization due to the spin only enters in the relativistic theory, and thus, we can put 1 = 0 . Finally, the term ∝ D t in Eq. (131) is a small correction (in a quantum relativistic expansion ℏ t ∕mc 2 ), and hence, we can use the approximation a 0 = ⋅ 1 ∕m . With these simplifications as a starting point, the DHW equations can be combined to give the following evolution equation for the magnetization: Normally, we should drop the term ∝ × , since such a term is identically zero. However, recall that before introducing approximations, this term would rather be proportional to ̃ ×̃ . As it turns out, the local approximations are applicable everywhere else, as in those cases, the corrections are compared with larger surviving terms. Here, however, we need to use the full expression ∝̃ ×̃ , and evaluate the term to the first nonvanishing order in an expansion in the small parameter ℏ∇∇ p . Performing this expansion and identifying 2 with and v 0 with f 0 , after some algebra, we can confirm the exact agreement of Eq. (160) with Eq. (75).
Next, to establish agreement with the model based on the Pauli Hamiltonian, we need to re-derive Eq. (74). Using the same approximations as before, we immediately get While there is a little bit of algebra involved (since the operator contains the magnetic field), it is straightforward to show that Eq. (161) reduces to Eq. (74). Finally, we note that the current sources to be used in the DHW equations are the same as in the Pauli limit, given that polarizations currents are dropped in nonrelativistic theory. Furthermore, as the agreement between Eqs. (74)-(75) and (68) has already been established, and the short-scale physics (as described by Eq. (159)) have been recovered, for most practical purposes, the full model based on the Pauli Hamiltonian (Eq. (54)) have also been verified, although not in a strict sense.
While the above confirmation is reassuring, the main purpose of studying approximate versions of the DHW model is to find models that are easier to analyze in specific cases. The DHW equations allow for systematic expansions in numerous quantum and relativistic parameters, ℏ ∕mc 2 , ℏeF∕m 2 c 2 (where F represents electric and/or magnetic fields), ℏ∇∇ p , p/mc, and ℏk 2 ∕m to name a few. Here, and k represent general temporal and spatial scales, rather than specific frequencies and wave numbers. Depending on the ordering of the dimensionless parameters, there are several possibilities for approximate models containing different combinations of expansion parameters. A systematic search for consistent approximations of the DHW system is a project for further research.

Concluding remarks
While the aim of this review has been to cover many aspects of quantum kinetic theory, naturally, there are many interesting and important topics that we have not touched upon. While, undoubtedly, many things will be left out altogether, let us here at least partially remedy some of the omissions made, by pointing to a few relevant aspects of quantum kinetic theory that we have not covered.
Firstly, as is well known, there are close connections between hydrodynamic and kinetic theories. In particular, a common approach to derive accurate fluid theories is by making moment expansions of kinetic theories. In a quantum context, this scheme has been used to derive fluid theories from kinetic theories, e.g., by Refs. Cai et al. (2012); Haas et al. (2010aHaas et al. ( , 2010b for the model defined by the Schrödinger Hamiltonian and by Refs. Hurst et al. (2014);  for the Pauli Hamiltonian model, and also by Ref. Hurst et al. (2017) for the model where the Pauli Hamiltonian is extended by the spin-orbit term. Moreover, moment expansion for models including exchange effects has been made by, e.g., Refs. Haas (2021); Manfredi (2020).
Second, we note that for completely degenerate systems, in certain cases, kinetic models may be simplified. This happens in situations when the phase-space density is conserved, in which case the dynamics is determined by the Fermi surface. This has been explored in the so-called waterbag models of plasmas (Manfredi 2020(Manfredi , 2005 and also for the case of a semiclassical model based on the Pauli Hamiltonian . Thirdly, we stress that we here have focused on nondissipative quantum kinetic models, ignoring the effects of higher order correlation in the BBGKY-hierarchy, neglecting all the influence of collisions, as can be covered, e.g., by path integral quantum monte carlo methods (Dornheim et al. 2018;Hamann et al. 2020;Zhang et al. 2016). However, the effect of dissipation is a broad research topic in its own right, of particular importance in the strong coupling regime. For references covering this field, see e.g. (Bonitz 2016;Ichimaru 1982;Bonitz et al. 2015).
As our review has focused on quantum kinetic models, and the worked out examples has served the purpose of illustrating the applicability of the theory, many aspects of quantum kinetic theory have still been ignored; in particular, interesting theoretical aspects involving wave-particle interaction, spin-polarization dynamics, numerical computation schemes, and the interesting interplay with single-particle dynamics (Dinu et al. 2016) A BBGKY-hierarchy and the mean-field approximation The N-particle density matrix is not directly applicable in the case where we have a large number of particles, as is typically the case for plasmas. It is, therefore, necessary to make some approximations to reduce the complexity of the equations. One commonly used approximations when dealing with plasmas is the mean-field approximation, which we will consider here. Since it is more convenient for this type of considerations, we will here use the operator representation of the density matrix (as opposed to using the position representation as was done in Sect. 2. Using the bra-ket notation, we can define the density operator as where � i ⟩ are N-particle states. The evolution equation can then be written as where [⋅, ⋅] denotes the commutator. Assume that we have an Hamiltonian in the form where is the kinetic energy operator, with ̂ i being the momentum operator acting on the i:th particle. The operator in the last term V ij corresponds to particle-particle interactions, which we will later take as the Coulomb interaction. However, for now, we only need to assume that the interaction is symmetric under exchange of particles, i.e., V ij =V ji . We note that the density matrix is symmetric with respect to change of particles, that is This is true, as long as the particles are indistinguishable, and hence both for fermions and bosons. We may now define the s-particle reduced density matrix as where we, for convenience, have multiplied by the number of particles i times. Now, using the above, we can derive the so-called BBGKY-hierarchy for the reduced density matrices Here, we have assumed that N is large, so that we may use N − s ≈ N to any order s where the above equation is of any practical use. We note that the one-particle density matrix couples to the two-particle density matrix, etc. Now, we are only interested in the lowest equation which is To make progress from this, we write the two-particle density matrix as where ĝ 12 is the two-particle correlations defined by the equation above. In the lowest order approximation, we completely neglect the two-particle correlations ĝ 12 ≈ 0 . We then get the so-called mean-field, or Hartree approximation We note that the normalization of the density matrix is such that i.e., the number of particles. The diagonal elements can then be identified as the local density of particles ( , ) = n( ) . In particular, in the case of a plasma interacting via the Coulumb interaction (166)1,…i,…j,…N =̂1 ,…j,…,i,…N .
(172) Tr 1̂1 = N, and the last term is where the last step is written in terms of the mean-field potential, i.e., the potential created by all particles interacting with particle 1. Here, we have only considered particles which are indistinguishable, but we have not yet taken into account the correct symmetry. For fermions, we need to modify the procedure slightly to account for the antisymmetry of the wave function, but the basic principle is the same. We will do this modification in Sect. 3 where we consider exchange effects. The conclusion from the above considerations is that if we can assume that particle correlations are negligible, we may describe the plasma as particles moving in the mean field created by all the other particles.
Funding Open access funding provided by Umea University.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.