Kinetic theory for spin-1/2 particles in ultra-strong magnetic fields

When the Zeeman energy approaches the characteristic kinetic energy of electrons, Landau quantization becomes important. In the vicinity of magnetars, the Zeeman energy can even be relativistic. We start from the Dirac equation and derive a kinetic equation for electrons, focusing on the phenomenon of Landau quantization in such ultra-strong but constant magnetic fields, neglecting short-scale quantum phenomena. It turns out that the usual relativistic gamma factor of the Vlasov equation is replaced by an energy operator, depending on the spin state, and also containing momentum derivatives. Furthermore, we show that the energy eigenstates in a magnetic field can be computed as eigenfunctions of this operator. The dispersion relation for electrostatic waves in a plasma is computed, and the significance of our results is discussed.


I. INTRODUCTION
Quantum kinetic descriptions of plasmas are typically of most interest for high densities and modest temperature [1]. For this purpose, typically the the starting point is the Wigner-Moyal equation [1][2][3][4][5][6], or various generalizations thereof also accounting for physics associated with the electron spin, such as the magnetic dipole force [7,8], spin magnetization [7,9] and spin-orbit interaction [9]. Both weakly [10,11] and strongly [12,13] relativistic treatments have been presented in the recent literature.
Certain quantum phenomena depend strongly on the magnitude of the electromagnetic field, however, rather than on the density and temperature parameters. Phenomena such as radiation reaction (reviewed in, e.g., Ref. [14]) and pair-creation fall into this category. 1 Another field-dependent phenomenon is Landau quantization [16], which becomes prominent whenever the Zeeman energy due to the magnetic field is comparable to or larger than the thermal energy, or the Fermi energy, in the case of degenerate electrons.
In the atmospheres of pulsars and magnetars [17], the electron motion may become relativistic, and the magnetic field strength can be ultra-strong, i.e., the Zeeman energy may be comparable to or even larger than the electron rest mass energy [18]. Further information about pulsar properties can be gained through the emission profiles, see, e.g., Refs [19][20][21]. Theoretical studies of wave propagation relevant for strongly magnetized objects have been made both with kinetic [22] and hydrodynamic [23,24] models. However, most previous theoretical studies starting from * haidar.al-naseri@umu.se † jens.zamanian@umu.se ‡ robin.ekman@plymouth.ac.uk § gert.brodin@umu.se 1 While the classical radiation reaction, where the response is a smooth function of the orbit, is a useful approximation as long as the emitted spectrum is soft (dominated by photons with energies well below mc 2 ), the description based on QED (see e.g. [34,35]), containing discrete probabilistic contributions due to the emission of high energy quanta, is more generally applicable.
the Dirac equation (see also [11,12]) have been limited to cases where the magnetic field strength is well below the critical field B cr = m 2 c 2 / |q| .
Our objective in this work is to derive a fully relativistic kinetic model of spin-1/2 particles, applicable for ultrastrong magnetic fields, i.e. with relaxing the conditions given in previous works. Assuming that the electric field is low enough to avoid pair creation (i.e., below the critical electric field E cr = m 2 c 3 / |q| ), we can use the Foldy-Wouthuysen transformation [25,26] to separate particle states from antiparticle states in the Dirac equation. Moreover, we limit ourselves to the case where the characteristic spatial scale length of the fields is much longer than the Compton length L c = /mc.
Making a Wigner transformation of the density matrix, our approach results in an evolution equation for a 2 × 2 Wigner matrix, where the four components encode information regarding the spin states. However, the offdiagonal elements of the matrix is associated with the spin-precession dynamics which is too rapid to be resolved by the theory. Thus a further reduction is made, where only the diagonal components representing the spin-up and the spin-down states relative to the magnetic field remain.
Together with Maxwell's equations, we obtain a closed system describing the plasma dynamics. A novel feature of the model is the energy expression, where the usual gamma factor from classical relativistic theory is replaced by an operator in phase space, also depending on the magnetic field. Interestingly, the Landau-quantized states in a constant magnetic field turn out to be eigenfunctions of the energy-operator of this theory, which is helpful when computing the thermodynamic background state for the Wigner matrix. To demonstrate the usefulness of the theory, we compute the dispersion relation for Langmuir waves propagating parallel to an external magnetic field. Finally, the consequences of the theory and applications to astrophysics are discussed arXiv:2005.13916v1 [physics.plasm-ph] 28 May 2020

II. THE STRONG FIELD HAMILTONIAN
To derive our theory for ultra-strong magnetic fields, we will take the Dirac Hamiltonian as our starting point. Here m is the mass,Ê = qφ(r), andÔ = α ·π, and α and β are the Dirac matrices. Furthermorep is the canonical momentum,π =p − qA(r, t), q is the charge, and φ and A are, respectively, the scalar and vector potentials. From now one we use units such that c = 1.
In their seminal paper Foldy and Wouthuysen [25] found a way to decouple the upper and lower two components of the four-spinor. For the general case, this can be done up to a given order in a suitably chosen expansion parameter. In Ref. [25] this parameter was 1/m, meaning that the expansion is valid for sufficiently small energies. In this paper we will use the results of Refs. [26,27] where a modified Foldy-Wouthuysen transformation was developed for the case where the expansion parameter is instead the scale-length of the fields. This transformation hence makes it possible to take into account arbitrarily strong fields as long as we are only concerned with variations on sufficiently long scale-lengths.
The goal of a Foldy-Wouthuysen transformation is to obtain a Hamiltonian where the upper and lower pairs of components of the four-spinor are decoupled in the regime of interest. An operator is called odd if it couples the upper and the lower pairs of components of the fourspinor, and even if it does not. The odd and even terms of the Hamiltonian (2) areÊ andÔ, respectively, satisfying [β,Ê] = 0 and {β,Ô} = 0.
In [27], Silenko found a unitary transformation of the Hamiltonian operator whereÛ whereˆ = m 2 +Ô 2 . Furthermore, using the Dirac HamiltonianĤ in (3), one obtainŝ wherê whereT = 2ˆ (m +ˆ ) andF =Ê − i ∂ t . The Hamiltonian H still has odd terms, the anti-commutators in O are linear inÔ, i.e. linear in α. However, the odd termÔ in H should be small compared toˆ , after choosing an expansion parameter. Since we are interested in the effects of an ultra-strong magnetic field, our small expansion parameters will be the electric field strength E/E cr (such that pair production is negligible) and the inverse scale-length of the fields L c /L, but we will make no assumption on the magnetic field strength. Thus a second transformation can be preformed with the following operator SinceÔ is small, only the major corrections are taken into account. Finally, the transformed Hamiltonian iŝ This Hamiltonian has no odd terms, so we are done with the transformation of the Dirac Hamiltonian. Until now, we have presented the Foldy-Wouthuysen transformation derived in Ref. [26]. Our approach now is to calculate the commutators occurring in (5) to all orders in the magnetic field. In doing this we will expand some of them in a series using Ref. [28], see the Appendix for more details. After the calculation, the Hamiltonian H in (5) becomeŝ where where µ B = q /2m is the Bohr magneton and Note that we kept all orders of the magnetic field. However we only kept up to first order of the combination of with E and ∇ r . In the expression forĤ there are still some odd operators, but sinceÔ is linear in µ B E and µ B ∂ t B, they should be smaller thanˆ . Thus it is fine to only include the minor correction of the second transformation in (7). Thus the HamiltonianĤ F W after the second transformation will have the same structure as in (8). However, the anti-commutator inĤ F W in (8) is proportional to the square ofÔ . Since we only kept up to first order of E and ∂ t B, we keep only the first two terms ofĤ F Ŵ Taking the limit of weak B-field, i.e. keeping up to first order in µ B B/m, we recover the Hamiltonian in [26]. The Hamiltonian in (12) includes the spin orbit interaction, see [12] for more details, but since E B, it is negligible compared to the magnetic interaction. We will therefore neglect the spinorbit interaction from now on, as the main idea of this paper is to study the effects of a strong magnetic field on the dynamics of a plasma. The transformed Hamiltonian is noŵ This Hamiltonian includes all orders of µ B B/m and is fully relativistic, hence is suitable to be used in deriving a kinetic equation for plasma in an environment where the magnetic field is of the order of the critical field B cr . Note that all operators in both (12) and (13) are even, thus we will let β → 1 andΣ → σ from now on.

III. GAUGE-INVARIANT WIGNER FUNCTION
Now we want to derive a kinetic equation using the Hamiltonian (13). To do that, we start with the evolution equation for the density matrixρ αβ which is given by the Von Neumann equation Our goal now is to transform this equation into a kinetic equation for the Wigner quasi-distribution function. Considering the gauge-invariant Wigner function derived by Stratonovich [29] W αβ (r, p, t) = 1 2π We express (14) in terms of ρ(r, r ) and use the identities [29] F These identities can be utilized for a function F that only depends onπ. However, our Hamiltonian in (13) depends on both B(r) andπ, thus we divide the magnetic field as where B 0 is a constant strong magnetic field and δB(r) is a varying magnetic field. Furthermore, a Taylor series around B = B 0 can be done, see Ref. [30] for more details. However, the Taylor series gets more complicated for the higher order terms, thus a restriction on δB(r) needs to be done. Considering the case where µ B δB/m 1, note that we did consider µ B E/m 1 in the previous section, the Hamiltonian becomeŝ While this Hamiltonian does not contain δB explicitly, as will be seen below, the perturbed magnetic field will still be contained in the Lorentz force. Using this Hamiltonian and keeping up to first order in ∇ r as we did in the derivation of the Hamiltonian, the kinetic equation is Note that is a function of σ and that the momentum derivatives act on everything to the right of the operator. Thus, in the second and fourth terms of (20) 1/ acts also on p. In order to get a scalar theory, we can Taylorexpand around σ where I is the identity matrix and Next, we note that if initially W αβ has no off-diagonal elements, let us say that W 11 = W + , W 22 = W − , the evolution (20) for W + and W − will decouple into separate equations for the spin-up and spin-down populations, as defined relative to B 0 . While limiting ourselves to such initial conditions may seem unwarranted, the only thing left out by this restriction is the spin precession dynamics. However, since the time-scale for spin-precession is the inverse Compton frequency, we note that the present theory, based on the assumption ∂/∂t c/L c , is not designed to resolve the spin precession dynamics anyway. Hence, from now on, we will be using the above representation for W αβ , in which case (20) decouples into the scalar equations for W + and W − as follows: The kinetic equation in (24) is our main result in this work. The new effects of this kinetic equation are hiding in ± . Firstly, we have all orders of the spin magnetic moment, compared to previous models [11,12] where only the first order correction is included. Moreover, we have momentum derivatives in , which turn to be energy operators, see Section V for more details. (24) describes the dynamics of an ensemble of spin-1/2 particles in an ultra-strong magnetic field in the meanfield approximation. In this approximation, the electric and magnetic fields are generated by the sources via where ρ f and j f are the free charge and current density respectively Under the assumptions we have made, the bound sources arising from the spin are negligible, but these can be and have been included in other models, e.g. Refs. [7,11]

IV. CONSERVATION LAWS
To check the validity of the derived model, we derive the conservation law of energy and the mass continuity. Starting with the mass continuity, the number density of spin-up (spin-down) particles n ± can be given by n ± = d 3 p W ± . To show that this quantity is conserved we take the time derivative of it and use (24). Since is independent of r, it is trivial to show that The number densities are separately conserved because transitions between spin-up and spin-down states require absorption or emission of quanta with energies on the order of m.
Moving to the conservation of energy, the total energy density is We want now to show that the energy is conserved, taking the time derivative of E tot , and using Maxwell's equations together with the kinetic equation, we get where K is the energy flux This is precisely what one would expect: the Poynting vector for the fields, and the kinetic energy flux for the particles as required in a relativistic theory. Since we have neglected polarization and magnetization, there is no Abraham-Minkowski dilemma in this model; see Ref. [13] and references therein for a related discussion.

V. BACKGROUND WIGNER FUNCTION IN A CONSTANT MAGNETIC FIELD
In principle we can compute the time-indpendent solutions for W in a constant magnetic field B = B 0ẑ by solving the Dirac-equation for this geometry, making a sum over different particle states, and then perform the Foldy-Wouthuysen and Wigner transformations of Section III. Except for the Foldy-Wouthuysen transformation, this was done in a covariant approach in Ref. [31]. However, here we will take a shorter route to arrive at the same results. 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 trivial generalization of the Pauli case [7]. Both for the Pauli and the Dirac equations, the spatial dependence of the wavefunction in Cartesian coordinates can be expressed as a Hermite polynomial times a Gaussian function [32] only the energy eigenvalues for the Landau quantized states are different. Specifically, applying the Dirac theory, the energy of the Landau quantized states become where n = 0, 1, 2, . . . corresponds 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 parallell 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. [7] (based on the Pauli equation) with some relatively minor adjustments.
1. Contrary to Ref. [7], we have made no Q-transform to introduce an independent spin variable, and thus the spin-dependence of Ref. [7] reduces to W ± .
2. The Wigner function of Ref. [7] must be expressed in terms of the momentum, i.e., m(v 2 3. The non-relativistic energy of Ref. [7] is replaced by the relativistic expression (32) 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 T B ± for the case of electrons in thermodynamic equilibrium can be written 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 (33) contains no dependence on the azimuthal angle in momentum space, we can write when ± acts on φ n (p ⊥ ). Computing ± φ n (p ⊥ ) by Taylorexpanding 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 ± φ n (p ⊥ ) = m 1 + (2n + 1 ± 1) where ω ce = |qB0| 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 (33) gives the thermodynamic equilibrium expression W T B ± , we note that the plasma background state is not necessarily in thermodynamic equilibrium. Making use of the property (36), we note that the most general time-independent solution W 0± to (24) 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 Naturally, the expressions for W 0± and W T B ± 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, (33) 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 , (33) reduces to a Synge-Juttner distribution.
To give a concrete illustration, in Fig. 1 we have made 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 will be demonstrated in the next section, to a large degree the electrons behave as a multi-species system, where each particle species has its own rest mass, as given by (32) but with p z = 0. This is because the separation between Landau levels is on the order of the rest mass, and all excitation 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± essentially will be determined by the Boltzmann factors of (33). However, for the cases where our model (24) is of most interest, the magnetic field B 0 is strong enough to make relativistic Landau quantization prominent. Thus, in the next section and for the remainder of this paper, we will consider background distributions W 0± where simplifications based on ω ce /m 1 do not apply.

VI. LINEAR WAVES
The operator in the square root in (24) gives the impression that it is very technical and complex to apply the model in studying, e.g., waves in plasma. In this section we consider electrostatic waves in a homogeneous plasma by using (24). We consider the wave vector k = ke z and express the momentum p in cylindrical coordinates p ⊥ , ϕ p , and p z . To linearize (24), we separate variables according to W ± = W 0± (p ⊥ , p z ) + W 1± (z, p ⊥ , p z , t), E = E 1 e z and B = B 0 e z , where the subscripts 0 and 1 denote unperturbed and perturbed quantities respectively. Moreover, the perturbed quantities follow the wave plane ansatz The unperturbed Wigner function W 0± is given by the thermal background Wigner function in (33). Since the operators in ± have been shown to be energy eigenvalues when acting on W 0± ; we act on both sides by ω − kp z / ± −1 , such that the perturbed distribution function becomes where in the second equality we used that W 0± is a sum of eigenfunctions of ± , and the summation is over the Landau levels indexed by n. Using Poisson's equation, the dielectric tensor for the electrostatic case is Note that if we set to zero in E n± in the dielectric tensor, then the denominator in the the second term of the dielectric tensor is the same as for the relativistic Vlasov equation [33]. The background distribution W 0n± will be divided into its eigenstates depending on the temperature and the magnetic field, see Fig. 1. As a result, the dispersion relation (42) is that of a relativistic multi-species plasma where each species has its own rest mass, E n± .

VII. SUMMARY AND DISCUSSION
In the present paper we have derived a kinetic model for plasmas immersed in a relativistically strong magnetic field, i.e., with a field strength of the order of the critical field. Based on a Foldy-Woythausen transformation and a Wigner transformation, an evolution equation for the spin-up and spin-down components W ± has been found.
Besides having two components, the main difference to a classical relativistic model is that the gamma-factor in that theory is replaced by an operator containing momentum derivatives. Since our theory is formulated in phase space, we stress that such operators are fundamentally different from the operators of Hilbert space, and to the best of our knowledge, this effect has not been seen in any previous quantum kinetic models. An immediate effect of the energy being an operator, comes when studying the background Wigner function. Classically, or in less advanced quantum mechanical models, the evolution equation does not predict the detailed momentum dependence for a given Landau level. Thus the background expression has been put in by hand, and one has to return to the starting point of the theory (e.g., the Pauli or Dirac equations for single particles) to find proper expressions. Here, however, the eigenvalue equation ε W 0 = E n,± W 0 determines the background state for a given Landau level (bar a constant for the number density), and there is no need to go outside the kinetic theory itself to find proper initial conditions.
To avoid some technical difficulties related to the operator orderings, we have here divided the magnetic field into an ultra-strong but constant part (B 0ẑ ) and a fluctuating part δB. This approach allows for the treatment of large classes of problems in magnetar atmospheres, for example linear and nonlinear wave propagating in homogeneous backgrounds, even up to relativistic wave amplitudes, as long as the flutuating part fulfills µ B |δB| m. Of course, a full modeling of the magnetar surroundings, covering the dipole nature of the background source field, is beyond the scope of such a theory. Moreover, in order to focus on the physics due to ultra-strong magnetic fields, the present theory excludes effects such as the magnetic dipole force, the spin magnetization, and the spin-orbit interaction included in some previous models [10,11], which can be justified for the long scale lengths and moderate frequencies that we focus on here. In this context, however, it should be noted that omission of the spin-orbit interaction is closely related to a correction term of the free current density (see e.g. Eq. (21) of Ref. [10]). While the additional term kept by Ref. [10] is a small correction for the conditions studied in this paper, it can contribute with currents perpendicular to B 0 that may be of importance for certain problems, specifically for geometries where the (otherwise larger) parallel currents vanish. This and other possible extensions of the current theory is a project for future research.
To illustrate the usefulness of the present theory, we have computed the dispersion relation for Langmuir waves in a strong magnetic field for a relativistic temperature. To a large extent, we find that the electrons behave as if they are divided into different species. More concretely, each Landau level of the background plasma contributes to the susceptibility with a term similar to the classical relativistic expression, but with its own effective mass m n± = m (1 + (2n + 1 ± 1) ω ce /m) 1/2 . We expect this result to generalize to some other problems, but not be completely general, as for certain problems, the difference between the standard gamma factor and the energy expression of the current theory will be apparent. A more complete study of the effects due to ultra-strong magnetic fields is a project for future research.

Appendix A: Commutators
The Hamiltonian in (5) contains some commutators of functions of operators. To calculate these commutators, we need to expand them in a series. We present here some of the calculations that were done in order to obtain the result in (9).
We can now rewrite (5) aŝ Looking firstly at the commutator ofˆ andF (the same result can be used to the commutator ofT andF), we have To calculate the commutator [ˆ , φ(r)], we expand the functions in the commutators in a series [28] [ˆ , φ(r)] = − ∞ k=1 (i ) k k!ˆ (k) φ k (r) where in the last equality, higher derivative terms were neglected in accordance with the long scale approximation. Thus, we have ˆ ,F = i qπ · E.
Next, we calculate the commutator where in the last equality, we neglected the commutator of and E j in accordance with the long-scale approximation. For the the commutator ofˆ andπ j , we need to expand it in a series. However this time, we have a commutator of functions that depend on bothp andr. Using the result from [28], we expand the commutator in a series where in the last equality we only kept up to k = 1 in the summation since higher order terms vanish in the long-scale approximation. Finally, we have ˆ , ˆ ,F = (i q) 2π 2 · (E × B).
We did not need to do any approximation in calculating this commutator since it is linear inÔ. Finally, we calculate where in the last equality we have used Using (A6) and (A9) to (A12) in (A3), we get the Hamil-tonianĤ in (9).