# Phase-space methods for the spin dynamics in condensed matter systems

Jérôme Hurst, Paul-Antoine Hervieux, Giovanni Manfredi

► To cite this version:

Jérôme Hurst, Paul-Antoine Hervieux, Giovanni Manfredi. Phase-space methods for the spin dynamics in condensed matter systems. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2017, 375 (2092), pp.20160199. 10.1098/rsta.2016.0199. hal-04582428

HAL Id: hal-04582428

https://hal.science/hal-04582428v1

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
## Research

<figure><img src="image_2.png" /></figure>

Check for
updates

**Cite this article:** Hurst J, Hervieux P-A, Manfredi G. 2017 Phase-space methods for the spin dynamics in condensed matter systems. *Phil. Trans. R. Soc. A* **375**: 20160199. http://dx.doi.org/10.1098/rsta.2016.0199

Accepted: 26 October 2016

One contribution of 11 to a theme issue 'Theoretical and computational studies of non-equilibrium and non-statistical dynamics in the gas phase, in the condensed phase and at interfaces'.

**Subject Areas:**

solid-state physics, quantum physics

**Keywords:**

spin and charge dynamics, ultrafast phenomena, nanostructures, Wigner quasi-probability distribution, quantum phase space

**Author for correspondence:**

Paul-Antoine Hervieux
e-mail:
paul-antoine.hervieux@ipcms.unistra.fr

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3685624.

# Phase-space methods for the spin dynamics in condensed matter systems

Jérôme Hurst, Paul-Antoine Hervieux and

Giovanni Manfredi

Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, 67000 Strasbourg, France

<figure><img src="image_3.png" /></figure>

P-AH, 0000-0002-4965-9709

Using the phase-space formulation of quantum mechanics, we derive a four-component Wigner equation for a system composed of spin-$\frac{1}{2}$ fermions (typically, electrons) including the Zeeman effect and the spin-orbit coupling. This Wigner equation is coupled to the appropriate Maxwell equations to form a self-consistent mean-field model. A set of semiclassical Vlasov equations with spin effects is obtained by expanding the full quantum model to first order in the Planck constant. The corresponding hydrodynamic equations are derived by taking velocity moments of the phase-space distribution function. A simple closure relation is proposed to obtain a closed set of hydrodynamic equations.

This article is part of the themed issue 'Theoretical and computational studies of non-equilibrium and non-statistical dynamics in the gas phase, in the condensed phase and at interfaces'.

## 1. Introduction

The formulation of quantum mechanics in the phase space was first introduced by Eugene Wigner in 1932 to study quantum corrections to classical statistical mechanics [1]. The goal was to link the wave function that appears in the Schrödinger equation to a pseudo-probability distribution function defined in the classical phase space. This pseudo-probability distribution changes in time according to an evolution equation (Wigner equation) which is somewhat similar to the classical Liouville equation. Mathematically speaking, the Wigner formulation is based on the Weyl
transformation [2,3], which is a general method to transform operators defined in the Hilbert space into phase-space functions.

As it is based on the classical phase space, the Wigner formulation is often a more intuitive approach than the standard Schrödinger equation, especially for problems where semiclassical considerations are important. For these reasons, it is used in many areas of quantum physics, including quantum optics [4], semiclassical analysis [5,6], electronic transport [7], nonlinear electron dynamics [8] and quantum plasma theory [9]. It is also the starting point to construct quantum hydrodynamic equations, which are approximate models obtained by taking velocity moments of the Wigner function. Such models were used in the past to study the electron dynamics in molecular systems [10], metal clusters and nanoparticles [11–13], thin metal films [14], quantum plasmas [15,16] and semiconductors [17].

The works cited above only considered the charge dynamics and disregarded the spin degrees of freedom. However, it is well known that spin effects (particularly the Zeeman splitting and spin-orbit coupling) can play a decisive role in nanometric systems such as semiconductor quantum dots [18,19] and diluted magnetic semiconductors [20,21]. The coupling between the spin degrees of freedom and the electron orbital motion is of the utmost importance in many experimental studies involving magnetized nano-objects. A particularly interesting example is the ultrafast demagnetization induced by a femtosecond laser pulse in ferromagnetic thin films [22], an effect that is not yet completely elucidated from the theoretical viewpoint. Recent time-dependent density functional theory (TDDFT) simulations suggest that the spin-orbit coupling plays a central role in the demagnetization process [23].

Phase-space models based on the Boltzmann equation [24]—and the corresponding fluid models [25]—were derived in the past to describe the dynamics of a gas where the constituents possess internal degrees of freedom (internal angular momentum). However, in these models, the spin is not treated *ab initio* as a fundamental quantity, but is rather incorporated into the transport equations to ensure the correct conservation properties. More recently, a few theoretical models that include the spin in the Wigner formalism were also developed. One approach [26] consists of defining a scalar probability distribution that evolves in an extended phase space, where the spin is treated as a classical two-component variable (related to the two angles on a unit-radius sphere) on the same footing as the position or the momentum. This approach was used to derive a Wigner equation that incorporates spin effects through the Zeeman interaction [26]. Semiclassical [27] and hydrodynamic [28] spin equations were also derived from those models, including other relativistic effects such as the spin-orbit coupling, the Darwin term and the relativistic mass correction.

An alternative approach keeps the 2 × 2 matrix character of the distribution function [29], so that the orbital and spin dynamics are represented by different Wigner functions. Using this approach, the corresponding Wigner equations were derived from the full Dirac theory [30]; however, their complexity makes them unsuitable for applications to condensed matter and nanophysics. A more tractable Wigner equation was derived from the Pauli (instead of Dirac) theory, but only included the Zeeman effect [31].

Both approaches (extended phase space and matrix Wigner function) are equivalent from the mathematical point of view. However, the extended phase-space approach leads to cumbersome hydrodynamic equations that are, in practice, very hard to solve, either analytically or numerically, even in the non-relativistic limit. The matrix technique, which is the one used here, separates clearly the orbital motion from the spin dynamics and leads to a simpler and more transparent hydrodynamic model.

In this paper, we go beyond our previous work [31] by including both the Zeeman effect and the spin-orbit coupling, the latter being a relativistic effect to second order in $1/c$. In terms of a semiclassical expansion, these terms are respectively first and second order in $\hbar$. Other relativistic corrections such as the Darwin term or the mass correction are neglected here, although they could be included with relative ease in our model.

First, we use a gauge-invariant formulation of the Weyl transformation [32] and the Moyal product [33] to derive a set of Wigner equations describing a system of spin-$\frac{1}{2}$ fermions.
A self-consistent mean-field model can further be obtained by coupling these Wigner equations to the set of Maxwell equations for the electromagnetic fields, where the sources (charge and current densities) are related to velocity moments of the Wigner function. A related mean-field model was obtained recently by Dixit et al. [34] in the framework of the Schrödinger–Pauli equation.

Subsequently, we derive the corresponding semiclassical limit and obtain the Vlasov equations describing the evolution of an electron gas with spin and semirelativistic effects. In this model, the orbital dynamics is treated classically, whereas the spin is represented as a fully quantum variable. Finally, the Vlasov equations will be used to derive a hierarchy of hydrodynamic equations by taking velocity moments of the probability distribution function. This is an infinite hierarchy that needs to be closed by using some additional physical hypotheses. Although this is relatively easy for spinless systems (where the closure can be obtained by assuming a suitable equation of state), things are subtler when the spin degrees of freedom are included. Here, we shall use an intuitive closure to obtain a closed set of fluid equations with spin effects.

## 2. Quantum mechanics in the phase space

The basic idea of the phase-space formulation of quantum mechanics is to associate at each operator $\hat{O}(\hat{R}, \hat{P})$, depending on the position and momentum operators $\hat{R}$ and $\hat{P}$, a function $O(r, p)$ of the classical phase-space variables $r$ and $p$. This correspondence is provided by the Weyl transformation [3,35], and is given by

$$ \hat{O}(\hat{R}, \hat{P}) \equiv \int dr dp O(r, p) \hat{F}(r, p), \qquad (2.1) $$

where $\hat{F}(r, p)$ is the Wigner operator defined as

$$ \hat{F}(r, p) \equiv \frac{1}{(2\pi\hbar)^6} \int du dv \exp \left[ \frac{i}{\hbar} (u \cdot (\hat{P} - p) + v \cdot (\hat{R} - r)) \right]. \qquad (2.2) $$

The inverse of the Weyl transformation can be deduced<fnref n="1" /> from the above definition:

$$ O(r, p) = (2\pi\hbar)^3 \text{tr}[\hat{F}(r, p)\hat{O}(\hat{R}, \hat{P})], \qquad (2.3) $$

where tr denotes the trace.

Considering a system in a statistical distribution described by the density operator $\hat{\rho}$, equation (2.1) can be used to determine the mean value of an arbitrary operator:

$$ \langle \hat{O}(\hat{R}, \hat{P}) \rangle = \text{tr}[\hat{O}(\hat{R}, \hat{P})\hat{\rho}] = \int dr dp O(r, p) \text{tr}[\hat{F}(r, p)\hat{\rho}]. $$

The Wigner function is then defined as the phase-space function associated with the density operator

$$ f(r, p) = \text{tr}[\hat{F}(r, p)\hat{\rho}] = \frac{1}{(2\pi\hbar)^3} \int d\lambda \exp \left( \frac{i}{\hbar} \lambda \cdot p \right) \left\langle r - \frac{\lambda}{2} \middle| \hat{\rho} \middle| r + \frac{\lambda}{2} \right\rangle. \qquad (2.4) $$

The Wigner function obeys the following equation of motion (Wigner equation):

$$ i\hbar \frac{\partial f}{\partial t} = \{\mathcal{H}, f\}_*, \qquad (2.5) $$

where the last term is referred to as the Moyal bracket

$$ \{A(r, p), B(r, p)\}_* = 2i \sin \left[ \frac{\hbar}{2} (L \partial_i^R \partial_{p_i} - L \partial_{p_i}^R \partial_j) \right] (A(r, p), B(r, p)). \qquad (2.6) $$

The superscripts L and R mean that the derivative acts only on the left or on the right term in the parenthesis ($i=x, y, z$ and here and in the following we use Einstein's summation convention). The Wigner equation is the analogue of the density matrix evolution equation in the operator representation of quantum mechanics: $\hbar\partial_t\hat{\rho} = [\hat{H}, \hat{\rho}]$, sometimes called the von Neumann

<footnote n="1">¹For the demonstration, we use the following properties of the Wigner operator: tr[$\hat{F}(r, p)\hat{F}(r', p')$] = (1/(2πħ)³)δ(p - p')δ(r - r').</footnote>
equation. The Moyal brackets can be easily developed as a power series of $\hbar$, which makes the Wigner formulation particularly interesting to study the semiclassical limit. The lowest-order term leads to the standard Poisson brackets and to the equations of classical mechanics.

In equation (2.5), $\mathcal{H}$ is the phase-space function associated with the Hamiltonian operator $\hat{\mathcal{H}}$ of the system, and they are related by equation (2.1). In order to determine the phase-space function of any arbitrary operator $\hat{\mathcal{O}}(\hat{\mathbf{R}}, \hat{\mathbf{P}})$, one should apply the Weyl correspondence rule [3,32], defined as follows: (i) first symmetrize the operator $\hat{\mathcal{O}}(\hat{\mathbf{R}}, \hat{\mathbf{P}})$ with respect to the position and the momentum operators $\hat{\mathbf{R}}$ and $\hat{\mathbf{P}}$, (ii) then replace $\hat{\mathcal{O}}(\hat{\mathbf{R}}, \hat{\mathbf{P}})$ by their associated classical variables. For instance, for the operator $\hat{P}_x \hat{X}$, one finds

$$ \hat{P}_x \hat{X} = \frac{1}{2} (\hat{P}_x \hat{X} + \hat{X} \hat{P}_x) - \frac{i\hbar}{2} \rightarrow xp_x - \frac{i\hbar}{2}, \qquad (2.7) $$

where use has been made of the commutator $[\hat{P}_x, \hat{X}] = \hbar/\mathrm{i}$. We note that the Weyl correspondence defined above is not unique, and one could have defined other rules leading to a different phase-space function, such as the Husimi representation [36].

In the case of a spinless particle moving in a scalar potential $V(r)$, the Wigner evolution (equation (2.5)) reads as

$$ \frac{\partial f(r, p, t)}{\partial t} + \frac{1}{m} p \cdot \nabla f(r, p, t) = -\frac{i}{\hbar} \frac{1}{(2\pi\hbar)^3} \int d\lambda \, d\lambda' \exp \left[ \frac{i(p-p') \cdot \lambda}{\hbar} \right] \\ \times \left[ V\left(r+\frac{\lambda}{2}\right) - V\left(r-\frac{\lambda}{2}\right) \right] f(r, p', t). \qquad (2.8) $$

However, complications arise when we want to include magnetic interactions. It is well known that in the presence of magnetic fields one should use the kinetic momentum operator $\hat{\mathbf{H}} = \hat{\mathbf{P}} - q\hat{\mathbf{A}}$ instead of $\hat{\mathbf{P}}$ ($q = -e, e > 0$ for an electron). This situation cannot be addressed by simply replacing $\hat{\mathbf{P}}$ with $\hat{\mathbf{H}}$ in the Weyl transformation. Indeed, it can be easily proven that with such substitution the Wigner function, equation (2.4), is not gauge invariant. As spin effects, such as the Zemann interaction or the spin-orbit coupling, strongly depend on the magnetic field, it is of paramount importance to work with a gauge-invariant formulation of the Weyl transformation. A gauge-independent definition of the Wigner function was first introduced by Stratonovich [37]:

$$ f(r, v, t) = \left(\frac{m}{2\pi\hbar}\right)^3 \int d\lambda \exp\left[\frac{i\lambda}{\hbar} \cdot (mv + q) \int_{-1/2}^{1/2} d\tau A(r + \tau\lambda)\right] \left\langle r - \frac{\lambda}{2} \middle| \hat{\rho} \middle| r + \frac{\lambda}{2} \right\rangle, \qquad (2.9) $$

where the momentum $p$ was replaced by $mv + q \int_{-1/2}^{1/2} d\tau A(r + \tau\lambda)$.

To be consistent with this new definition of the Wigner function, one should also modify the Weyl correspondence rule [32]. The procedure is identical, except that one should use $\hat{\mathbf{H}}$ instead of $\hat{\mathbf{P}}$. The main difference is that one must also symmetrize operators with respect to the different components of $\hat{\mathbf{H}}$, because they do not commute, i.e. $[\hat{\mathbf{H}}_i, \hat{\mathbf{H}}_j] = i\hbar q \epsilon_{ijk} \hat{\mathbf{B}}_k(\hat{\mathbf{R}})$, where $\epsilon_{ijk}$ is the Levi-Civita symbol and $\hat{\mathbf{B}}$ is the magnetic field operator ($\hat{\mathbf{B}} = \nabla \times \hat{\mathbf{A}}$). The classical phase-space variable associated with the kinetic momentum operator is the linear momentum $\hat{\mathbf{H}} \rightarrow \pi = mv$.

The Moyal product defined in equation (2.6) is also modified in the presence of magnetic fields. A gauge-invariant Moyal product was derived by Müller [38], and reads

$$ A(r, \pi) \star C(r, \pi) = \exp \left[ i\hbar\mathcal{L} + ie \sum_{n=1}^{\infty} \hbar^n \mathcal{L}_n \right] (A(r, \pi), C(r, \pi)), \qquad (2.10) $$

where $\mathcal{L}$ is the operator corresponding to the free magnetic field case:

$$ \mathcal{L}(A(r, \pi), C(r, \pi)) = \frac{1}{2} (\partial_i^R \partial_{\pi_i} - \partial_j^L \partial_{\pi_j}) (A(r, \pi), C(r, \pi)) \qquad (2.11) $$
and $\mathcal{L}_n$ is a new operator that depends on the magnetic field:

$$ \mathcal{L}_n (A(\mathbf{r}, \boldsymbol{\pi}), C(\mathbf{r}, \boldsymbol{\pi})) = \left(\frac{i}{2}\right)^{n+1} \frac{\epsilon_{jlr}}{(n+1)^2 n!} \sum_{i_1 \dots i_{n-1}=1}^{3} \left(\frac{\partial^{n-1}}{\partial r_{i_1} \dots \partial r_{i_{n-1}}} B_r\right) L_{\partial \pi_j} R_{\partial \pi_l} \sum_{p=1}^{n} \binom{n+1}{p} g(n,p) \frac{\partial^{p-1}}{\partial \pi_{i_1} \dots \partial \pi_{i_{p-1}}} \frac{\partial^{n-p}}{\partial \pi_{i_p} \dots \partial \pi_{i_{n-1}}} (A(\mathbf{r}, \boldsymbol{\pi}), C(\mathbf{r}, \boldsymbol{\pi})), \quad (2.12) $$

with $g(n,p) = [(1 - (-1)^p)(n+1) - (1 - (-1)^{n+1})p]$, $(r_1, r_2, r_3) = (x, y, z)$ and $(\pi_1, \pi_2, \pi_3) = (\pi_x, \pi_y, \pi_z)$. This new definition of the Moyal product makes the calculation of the evolution equation much more cumbersome than in the unmagnetized case. Its great advantage is that it ensures that the final equations of motion are gauge invariant.

### 3. Derivation of the spin Wigner model

We consider an ensemble of fermions in the presence of an electromagnetic field $E, B$. We denote the Schrödinger wave function of the $\mu$th particle state by

$$ \Psi_{\mu}(\mathbf{r}, t) = \Psi_{\mu}^{\uparrow}(\mathbf{r}, t)|\uparrow\rangle + \Psi_{\mu}^{\downarrow}(\mathbf{r}, t)|\downarrow\rangle, \quad (3.1) $$

where $\Psi_{\mu}^{\uparrow}(\mathbf{r}, t)$ and $\Psi_{\mu}^{\downarrow}(\mathbf{r}, t)$ are respectively the spin-up and spin-down components of the wave function. The evolution of the system is governed by the Pauli–Schrödinger equation

$$ i\hbar \frac{\partial \Psi_{\mu}(\mathbf{r}, t)}{\partial t} = \hat{\mathcal{H}} \Psi_{\mu}(\mathbf{r}, t) \quad (3.2) $$

and

$$ \hat{\mathcal{H}} = \left( \frac{\hat{\mathbf{\Pi}}^2}{2m} - e\hat{V} \right) \sigma_0 + \left[ \mu_B \hat{B} + \frac{\mu_B}{4mc^2} (\hat{E} \times \hat{\mathbf{\Pi}} - \hat{\mathbf{\Pi}} \times \hat{E}) \right] \cdot \boldsymbol{\sigma}. \quad (3.3) $$

Here, $\mu_B = e\hbar/2m$ is the Bohr magneton, $\sigma$ is the vector made of the $2 \times 2$ Pauli matrices, and $\sigma_0$ is the $2 \times 2$ identity matrix. $V$ and $A$ are, respectively, the scalar and vector potential. Equation (3.3) can be derived from the Dirac equation by means of a Foldy–Wouthuysen transformation [39,40]. This semirelativistic development leads to plenty of terms that couple the spin to the charge dynamics. In this work, we keep terms only up to second order in $1/c$, where $c$ is the speed of light in vacuum, namely the Zeeman interaction (order $1/c^0$) and the spin–orbit coupling (order $1/c^2$). We neglect, however, the Darwin term and the relativistic mass correction, which are also second-order effects.

Without spin, the Wigner function is a scalar function related to the density matrix $\rho$ through equation (2.9). This definition can be generalized as follows to take into account the spin degrees of freedom:

$$ F(\mathbf{r}, v, t) = \left(\frac{1}{2\pi\hbar}\right)^3 \int d\lambda \exp \left[ \frac{i\lambda}{\hbar} \cdot \left( mv + q \int_{-1/2}^{1/2} d\tau A(\mathbf{r} + \tau\mathbf{\lambda}) \right) \right] \rho \left( \mathbf{r} - \frac{\lambda}{2}, \mathbf{r} + \frac{\lambda}{2}, t \right), \quad (3.4) $$

where, for particles with spin $\frac{1}{2}$, $F$ is a $2 \times 2$ matrix. The elements of the density matrix $\rho^{\eta\eta'}(\mathbf{r}, \mathbf{r}', t)$, where $\eta = \uparrow, \downarrow$, are given by

$$ \rho^{\eta\eta'}(\mathbf{r}, \mathbf{r}') = \sum_{\mu} \Psi_{\mu}^{\eta}(\mathbf{r}, t) \Psi_{\mu}^{\eta*}(\mathbf{r}', t). \quad (3.5) $$

In order to study the macroscopic properties of the system, it is convenient to project $F$ onto the Pauli basis set [41,42]

$$ F = \frac{1}{2}\sigma_0 f_0 + \frac{1}{\hbar}f \cdot \boldsymbol{\sigma}, \quad (3.6) $$

where

$$ f_0 = \text{tr}\{F\} = f^{\uparrow\uparrow} + f^{\downarrow\downarrow}, \quad f = \frac{\hbar}{2} \text{tr}(F\sigma) \quad (3.7) $$
and tr denotes the trace. With this definition, the particle density $n$ and the spin polarization $S$ of the electron gas are easily expressed by the moments of the pseudo-distribution functions $f_0$ and $f$:

$$n(r,t) = \sum_{\mu} |\Psi_{\mu}^{\dagger}(r,t)|^2 = \int f_0(r,v,t) \, dv \qquad (3.8)$$

and

$$S(r,t) = \frac{\hbar}{2} \sum_{\mu} \psi_{\mu}^{\dagger}(r,t) \sigma \psi_{\mu}(r,t) = \int f(r,v,t) \, dv. \qquad (3.9)$$

In this representation, the Wigner functions have a clear physical interpretation: $f_0$ is related to the total electron density (in phase space), whereas $f_k$ ($k = x, y, z$) is related to the spin polarization in the direction $k$. In other words, $f_0$ represents the probability of finding an electron at one point of the phase space at a given time, while $f_k$ represents the probability of that electron having a spin-polarization probability in the direction $k$.

There exist different ways to include the spin in the Wigner formalism other than the one described above. For instance, Brodin and co-workers [26] introduced an extended phase space $(r, v, s)$, where $s$ is a unitary vector that defines the spin direction. The corresponding probability distribution is a scalar function of the extended phase-space variables. This is in contrast with our approach, where the spin is treated as a fully quantum variable (evolving in a two-dimensional Hilbert space). Nevertheless, the two approaches are equivalent, as shown in [31]. The correspondence relations between our distribution functions $f_0(r, v, t)$ and $f(r, v, t)$ and the scalar distribution used by Zamanian et al. [26] $f_Z(r, v, s, t)$ read as

$$f_0 = \int f_Z \, d^2s, \quad f = 3 \int s f_Z \, d^2s. \qquad (3.10)$$

Let us now turn to the evolution equation for the Wigner functions $f_0(r, v, t)$ and $f_k(r, v, t)$. After some tedious calculations, developed in the electronic supplementary material, equation (2.5) leads to the following Wigner equations:

$$\begin{aligned} \frac{\partial f_0}{\partial t} + \frac{1}{m}(\boldsymbol{\pi} + \Delta\tilde{\boldsymbol{\pi}}) \cdot \nabla f_0 - \frac{e}{m}[m\tilde{E} + (\boldsymbol{\pi} + \Delta\tilde{\boldsymbol{\pi}}) \times \tilde{B}]_i \partial_{\pi_i} f_0 \\ - \mu_B \nabla \left( \tilde{B} - \frac{1}{2mc^2} \boldsymbol{\pi} \times \tilde{E} \right)_i \cdot \nabla_{\boldsymbol{\pi}} f_i + \frac{\mu_B}{4mc^2} [(E_+ + E_-) \times \nabla] \cdot f \\ - \frac{\mu_B e}{2mc^2} [\tilde{E} \times [\tilde{B} \times \nabla_{\boldsymbol{\pi}}]] \cdot f - \frac{\mu_B}{2mc^2} \frac{i}{\hbar} [\Delta\tilde{\boldsymbol{\pi}} \times (E_+ - E_-)] \cdot f = 0, \end{aligned} \qquad (3.11)$$

$$\begin{aligned} \frac{\partial f_k}{\partial t} + \frac{1}{m}(\boldsymbol{\pi} + \Delta\tilde{\boldsymbol{\pi}}) \cdot \nabla f_k - \frac{e}{m}[m\tilde{E} + (\boldsymbol{\pi} + \Delta\tilde{\boldsymbol{\pi}}) \times \tilde{B}]_i \partial_{\pi_i} f_k \\ - \mu_B \nabla(\tilde{B} - \frac{1}{2mc^2}\boldsymbol{\pi} \times \tilde{E})_k \cdot \nabla_{\boldsymbol{\pi}} f_0 + \frac{\mu_B}{4mc^2}[(E_+ + E_-) \times \nabla]_{f0} \\ - \frac{\mu_B e}{2mc^2}[\tilde{E} \times [\tilde{B} \times \nabla_{\boldsymbol{\pi}}]]_{kf0} - \frac{\mu_B}{2mc^2} \frac{i}{\hbar} [\Delta\tilde{\boldsymbol{\pi}} \times (E_+ - E_-)]_{kf0} \\ - \frac{e}{2m} \left[ \left(B_+ + B_- - \frac{1}{2mc^2}(\boldsymbol{\pi} + \Delta\tilde{\boldsymbol{\pi}}) \times (E_+ + E_-)\right) \times f \right]_k \\ + \frac{\mu_B}{2mc^2} \frac{i}{2} [((E_+ - E_-) \times (\nabla - e\tilde{B} \times \nabla_{\boldsymbol{\pi}})) \times f]_k = 0, \end{aligned} \qquad (3.12)$$

where $\Delta\tilde{\boldsymbol{\pi}}$ depends on the magnetic field and corresponds to a quantum shift of the velocity

$$\Delta\tilde{\boldsymbol{\pi}} = -i\hbar e \nabla_{\boldsymbol{\pi}} \times \left[ \int_{-1/2}^{1/2} d\tau \tau B(r + i\hbar\tau \nabla_{\boldsymbol{\pi}}) \right] \qquad (3.13)$$

and $\tilde{E}$ and $\tilde{B}$ are written in terms of the electric and magnetic fields

$$\tilde{E} = \int_{-1/2}^{1/2} d\tau E(r + i\hbar\tau \nabla_{\pi}), \quad \tilde{B} = \int_{-1/2}^{1/2} d\tau B(r + i\hbar\tau \nabla_{\pi}). \qquad (3.14)$$
The subscripts $\pm$ mean that the corresponding quantity is evaluated at a shifted position $r \pm i\hbar\nabla_{\pi}/2$. This particularly illuminating form of the Wigner equations was proposed by Serimaa et al. [32] for the case of a charged particle without spin evolving in an external electromagnetic field.

Equations (3.11) and (3.12) can be viewed as a generalization of those obtained in our previous work [31], where only the Zeeman interaction was included. The latter has two effects: the first is to couple the spin to the orbital dynamics through the gradient of the magnetic field (terms $\mu_B \nabla \tilde{B}_k \cdot \nabla_{\pi}$ in the equations); the second effect is the spin precession around an effective magnetic field [terms $(\mathbf{B}_+ + \mathbf{B}_-) \times \mathbf{f}$]. In addition, many new terms appear owing to the spin-orbit interaction, which can be easily identified because they are proportional to $1/c^2$. Some of these terms couple the spin to the orbital dynamics, whereas others provide corrections to the spin precession or the Lorentz force. The physical origin of all these terms will appear clearly in the next session, when we discuss the semiclassical limit of the Wigner equations.

Equations (3.11) and (3.12) can be used, in the context of a mean-field approach, to describe the self-consistent spin dynamics of an ensemble of interacting electrons. In this case, the electric and the magnetic fields are solutions of the Maxwell equations:

$$ 
\left.
\begin{aligned}
\nabla \cdot \mathbf{E} &= \frac{\rho}{\epsilon_0} - \frac{\nabla \cdot \mathbf{P}}{\epsilon_0}, \\
\nabla \cdot \mathbf{B} &= 0, \\
\nabla \times \mathbf{E} &= -\frac{\partial \mathbf{B}}{\partial t}, \\
\nabla \times \mathbf{B} &= \mu_0 j + \mu_0 \epsilon_0 \frac{\partial \mathbf{E}}{\partial t} + \mu_0 \frac{\partial \mathbf{P}}{\partial t} + \mu_0 \nabla \times \mathbf{M},
\end{aligned}
\right\} \quad (3.15)
$$

and

where we introduced some relativistic corrections to the source terms, namely a spin magnetization $M$, a spin polarization $P$ and a new contribution to the current density; see equation (3.17). These corrections appear when one considers an expansion in $1/c$ of the Dirac-Maxwell equations and are consistent with the Hamiltonian of equation (3.3), as was shown recently using a Lagrangian method [34,43]. Using equation (2.9), we can transpose their results to our formulation, which yields

$$ \rho = -e \int f_0 dv, \quad (3.16) $$

$$ j = -e \left[ \int v f_0 \, dv + \frac{E \times M}{2mc^2} \right], \quad (3.17) $$

$$ M = -\mu_B \int f \, dv \quad (3.18) $$

$$ P = -\frac{\mu_B}{2c^2} \int v \times f \, dv. \quad (3.19) $$

and

This mean-field approach could, in principle, be extended, in the spirit of density functional theory, to include exchange and correlation effects by adding suitable potentials and fields that are functionals of the electron density [44].

## 4. Semiclassical limit and spin Vlasov model

The form of equations (3.11) and (3.12) is particularly useful to study the semiclassical limit of the model. Indeed, we can easily expand $\tilde{E}$, $\tilde{B}$ and $\Delta\tilde{\pi}$ as a power series of $\hbar$

$$ 
\begin{aligned}
\tilde{E} &= \sum_{n=0}^{\infty} \left(\frac{\hbar}{2m}\right)^{2n} \frac{(-1)^n}{(2n+1)!} \sum_{i_1 \dots i_{2n}=1}^{3} \left(\frac{\partial^{2n}}{\partial r_{i_1} \dots \partial r_{i_{2n}}} E\right) \frac{\partial^{2n}}{\partial v_{i_1} \dots \partial v_{i_{2n}}} \\
&= E - \frac{\hbar^2}{12m^2} \sum_{i_1, i_2=1}^{3} \frac{\partial^2 E}{\partial r_{i_1} \partial r_{i_2}} \frac{\partial^2}{\partial v_{i_1} \partial v_{i_2}} + \mathcal{O}(\hbar^4),
\end{aligned}
\quad (4.1)
$$
$$
\begin{align}
\tilde{B} &= \sum_{n=0}^{\infty} \left(\frac{\hbar}{2m}\right)^{2n} \frac{(-1)^n}{(2n+1)!} \sum_{i_1 \dots i_{2n}=1}^{3} \left(\frac{\partial^{2n}}{\partial r_{i_1} \dots \partial r_{i_{2n}}} \mathbf{B}\right) \frac{\partial^{2n}}{\partial v_{i_1} \dots \partial v_{i_{2n}}} \nonumber \\
&= \mathbf{B} - \frac{\hbar^2}{12m^2} \sum_{i_1, i_2=1}^{3} \frac{\partial^2 \mathbf{B}}{\partial r_{i_1} \partial r_{i_2}} \frac{\partial^2}{\partial v_{i_1} \partial v_{i_2}} + \mathcal{O}(\hbar^4), \tag{4.2}
\end{align}
$$

$$
\begin{align}
\Delta \tilde{\pi} &= m \mu_B \sum_{n=0}^{\infty} \left(\frac{\hbar}{2m}\right)^{2n+1} \frac{(-1)^n (2n+2)}{(2n+3)!} \sum_{i_1 \dots i_{2n}=1}^{3} \left(\frac{\partial^{2n+1}}{\partial r_{i_1} \dots \partial r_{i_{2n+1}}} \mathbf{B}\right) \frac{\partial^{2n+1}}{\partial v_{i_1} \dots \partial v_{i_{2n+1}}} \nonumber \\
&= \sum_{i=1}^{3} \frac{\mu_B \hbar}{6} \frac{\partial \mathbf{B}}{\partial r_i} \frac{\partial}{\partial v_i} + \mu_B \mathcal{O}(\hbar^3). \tag{4.3}
\end{align}
$$

From these semiclassical expansions, we note that the velocity shift $\Delta\tilde{\pi}$ has a purely quantum origin, because the leading term in the expansion is of first order in $\hbar$. Therefore, it has no classical counterpart. In the case of $\tilde{E}$ and $\tilde{B}$, the leading term in the expansion simply corresponds to the classical electric or magnetic field.

To zeroth order, the equations for $f_0$ and $f_i$ decouple, so that one can study the particle motion irrespective of the spin degrees of freedom. To first order in $\hbar$, equations (3.11) and (3.12) become

$$
\begin{equation}
\begin{split}
& \frac{\partial f_0}{\partial t} + \mathbf{v} \cdot \nabla f_0 - \frac{e}{m} (\mathbf{E} + \mathbf{v} \times \mathbf{B}) \cdot \nabla \mathbf{v} f_0 + \frac{\mu_B}{2mc^2} (\mathbf{E} \times \nabla) f_i \\
& \qquad - \frac{\mu_B}{m} \nabla \left[ B_i - \frac{1}{2c^2} (\mathbf{v} \times \mathbf{E})_i \right] \cdot \nabla \mathbf{v} f_i - \frac{\mu_B e}{2m^2 c^2} [\mathbf{E} \times (\mathbf{B} \times \nabla \mathbf{v})]_{if_i} = 0,
\end{split}
\tag{4.4}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
& \frac{\partial f_i}{\partial t} + \mathbf{v} \cdot \nabla f_i - \frac{e}{m} (\mathbf{E} + \mathbf{v} \times \mathbf{B}) \cdot \nabla \mathbf{v} f_i + \frac{\mu_B}{2mc^2} (\mathbf{E} \times \nabla) f_0 \\
& \qquad - \frac{\mu_B}{m} \nabla \left[ B_i - \frac{1}{2c^2} (\mathbf{v} \times \mathbf{E})_i \right] \cdot \nabla \mathbf{v} f_0 - \frac{\mu_B e}{2m^2 c^2} [\mathbf{E} \times (\mathbf{B} \times \nabla \mathbf{v})]_{if_0} \\
& \qquad - \frac{2\mu_B}{\hbar} \left\{ \left[ \mathbf{B} - \frac{1}{2c^2} (\mathbf{v} \times \mathbf{E}) \right] \times \mathbf{f} \right\}_i = 0,
\end{split}
\tag{4.5}
\end{equation}
$$

where the factor $\hbar$ is hidden in the definition of the Bohr magneton $\mu_B = e\hbar/(2m)$. The quantum corrections in equations (4.4) and (4.5) couple the orbital and the spin dynamics through the Zeeman and spin-orbit interactions. There are no quantum corrections to the orbital electronic dynamics because they would appear only at the second order in $\hbar$. For instance, the Darwin term would not appear in the above equations (even if we had included it in the full Wigner equations) because it is a correction of order $\hbar^2$ to the orbital motion of the electron. In summary, equations (4.4) and (4.5) represent a semiclassical model where the orbital dynamics is classical (hence the familiar Lorentz force terms), whereas the spin is treated as a fully quantum variable (two-dimensional Hilbert space).

In equations (4.4) and (4.5), the Zeeman interaction gives two contributions: (i) the term $\mu_B \nabla B_i \cdot \nabla v$, which represents the force exerted on a magnetic dipole by an inhomogeneous magnetic field and is at the basis of Stern–Gerlach-type experiments, and (ii) the term $f \times B$, which describes the precession of the spin around the magnetic field lines.

The spin-orbit interaction yields a correction to the magnetic field $B \rightarrow B - (v \times E)/2c^2$, which is the first-order correction in the non-relativistic limit of the Thomas precession [45,46]. The other terms are related to the spin-orbit correction of the velocity operator. Indeed, in the Heisenberg picture, the velocity operator $\hat{V}$ is determined by the evolution equation of the position operator

$$
\hat{V} = \frac{1}{i\hbar} [\hat{R}, \hat{\mathcal{H}}] = \frac{\hat{\Pi}}{m} - \frac{\mu_B}{2mc^2} \hat{E} \times \hat{\sigma}, \qquad (4.6)
$$

where we used the Hamiltonian defined in equation (3.3). The associated phase-space function is determined by the Weyl correspondence rule and reads as

$$
\mathbf{V} = \mathbf{v} - \frac{\mu_B}{2mc^2} \mathbf{E} \times \boldsymbol{\sigma}. \qquad (4.7)
$$
This is the phase-space function that can be used to calculate the average velocity or the charge current. Therefore, the particles are transported with a modified velocity. The term $(E \times \nabla)_i$ in equation (4.4) and (4.5) is a direct consequence of this effect, whereas the term $[E \times (B \times \nabla_v)]_i$ corresponds to the same velocity correction in the Lorentz force $v \times B$.

The spin Vlasov equations (4.4) and (4.5) are correct to second order in $1/c$ and can be viewed as a generalization of those obtained in [31], where only the Zeeman interaction was included. An alternative form of these equations was obtained by Asenjo et al. [28] in the extended phase-space formalism.

The Maxwell equations (3.15), combined with the spin Vlasov equations (4.4) and (4.5), form a self-consistent model for the charge and the spin dynamics of a system of interacting particles. One can show that the following quantities are conserved during the time evolution:

$$M_{\text{tot}} = m \int f_0 \, dv \, dr, \qquad (4.8)$$

$$P_{\text{tot}} = m \int v f_0 \, dv \, dr + \int D \times B \, dr, \qquad (4.9)$$

$$E_{\text{tot}} = \frac{m}{2} \int v^2 f_0 \, dv \, dr + \mu_B \int f \cdot B \, dv \, dr + \frac{\epsilon_0}{2} \int E^2 \, dr + \frac{1}{2\mu_0} \int B^2 \, dr \qquad (4.10)$$

$$\text{and} \qquad J_{\text{tot}} = m \int (r \times v) f_0 \, dr \, dv + \frac{\hbar}{2} \int f \, dr \, dv + \int r \times (D \times B) \, dr, \qquad (4.11)$$

where $D = \epsilon_0 E + P$ and $H = B - \mu_0 M$. The conserved quantities are the total mass $M_{\text{tot}}$, the total linear momentum $P_{\text{tot}}$ (sum of the particles and fields momenta), the total energy $E_{\text{tot}}$ (sum of the kinetic, Zeeman and the electromagnetic field energies), and the total angular momentum $J_{\text{tot}}$ (sum of the orbital, spin and electromagnetic field angular momenta).

For simulation purposes, the spin Vlasov equations (4.4) and (4.5) are much easier to solve numerically than the corresponding Wigner equations (3.11) and (3.12), mainly because the former are local in space while the latter are not. The Vlasov approximation is valid when quantum effects in the orbital dynamics are small. From equation (2.6), it appears that the semiclassical expansion is valid when $\hbar/(mL_0v_0) \ll 1$, where $L_0$ and $v_0$ are typical length and velocity scales. For a degenerate electron gas with density $n$, the typical velocity is the Fermi speed $v_F = \hbar(3\pi^2n)^{1/3}/m$. Inserting into the previous inequality, we obtain the validity condition $L_0n^{1/3} \gg 1$, which means that the typical length scale must be larger than the interparticle distance $d = n^{-1/3}$. For this reason, the semiclassical limit is also referred to as the long-wavelength approximation. All in all, the above spin Vlasov equations constitute a valuable tool to simulate the charge and spin dynamics in condensed matter systems, particularly semiconductor and metallic nano-objects.

## 5. Hydrodynamic model with spin-orbit coupling

In this section, starting from equations (4.4) and (4.5), we derive the hydrodynamic evolution equations by taking velocity moments of the phase-space distribution functions. In addition to the particle density and spin polarization (equations (3.8) and (3.9)), we define the following macroscopic quantities:

$$u = \frac{1}{n} \int v f_0 \, dv, \qquad (5.1)$$

$$J_{i\alpha}^{S} = \int v f_{i\alpha} \, dv, \qquad (5.2)$$

$$P_{ij} = m \int w_i w_j f_0 \, dv, \qquad (5.3)$$

$$\Pi_{ij\alpha} = m \int v_i v_j f_{\alpha} \, dv \qquad (5.4)$$

$$Q_{ijk} = m \int w_i w_j w_k f_0 \, dv, \qquad (5.5)$$

and
where we separated the mean fluid velocity $\mathbf{u}$ from the velocity fluctuations $\mathbf{w} \equiv \mathbf{v} - \mathbf{u}$. Here, $P_{ij}$ and $Q_{ijk}$ are respectively the pressure and the generalized energy flux tensors. They coincide with the analogous definitions for spinless fluids with probability distribution function $f_0$. The spin-velocity tensor $J_{i\alpha}^S$ represents the fluid current along the $i$th direction of the $\alpha$th spin polarization vector, whereas $\Pi_{ij\alpha}$ represents the corresponding spin-pressure tensor.<fnref n="2" /> The evolution equations for the above fluid quantities are obtained by the straightforward integration of equations (4.4) and (4.5) with respect to the velocity variable. One obtains

$$ \frac{\partial n}{\partial t} + \nabla_r \cdot (n\bar{\mathbf{u}}) = 0, \qquad (5.6) $$

$$ \frac{\partial S_{\alpha}}{\partial t} + \partial_i \bar{J}_{i\alpha}^S + \frac{e}{m}(S \times B)_\alpha + \frac{e}{2mc^2} \epsilon_{jk\alpha} \epsilon_{rlj} E_{jl} J_{rk}^S = 0, \qquad (5.7) $$

$$ \frac{\partial u_i}{\partial t} + u_j(\nabla_j u_i) + \frac{1}{nm}\nabla_j P_{ij} + \frac{e}{m}[E_i + (\bar{\mathbf{u}} \times B)_i] + \frac{e}{nm^2}S_{\alpha}(\partial_i B_{\alpha}) \\ + \frac{\mu_B}{2mc^2n}\epsilon_{jkl}[u_i \partial_j(E_k S_l) + E_j(\partial_k J_{il}^S) - (\partial_i E_k)J_{jl}^S - (\partial_j E_k)J_{il}^S] = 0, \qquad (5.8) $$

$$ \frac{\partial J_{i\alpha}^S}{\partial t} + \partial_j \Pi_{ij\alpha} + \frac{eE_i}{m}S_{\alpha} + \frac{e}{m}\epsilon_{jki}B_k \bar{J}_{j\alpha}^S + \frac{e}{m}\epsilon_{jk\alpha}B_k J_{ij}^S + \frac{\mu_B\hbar}{2m}(\partial_i B_{\alpha})n \\ + \frac{\mu_B}{2mc^2}\epsilon_{kl\alpha}\partial_l(E_k n u_i) - \frac{\mu_B}{2mc^2}\epsilon_{kl\alpha}(\partial_i E_l)n u_k + \frac{\mu_B}{\hbar c^2}\epsilon_{kl\alpha}\epsilon_{rsk}E_s \Pi_{irl}^S = 0 \qquad (5.9) $$

$$ \begin{aligned} \frac{\partial P_{ij}}{\partial t} + u_k \partial_k P_{ij} + P_{jk} \partial_k u_i + P_{ik} \partial_k u_j + P_{ij} \partial_k u_k + \partial_k Q_{ijk} + \frac{e}{m} [\epsilon_{klj} P_{jk} + \epsilon_{klj} P_{ik}] B_l \\ + \frac{\mu_B}{m} [\partial_i B_k (J_{jk}^S - u_j S_k) + \partial_j B_k (J_{ik}^S - u_i S_k)] + \frac{\mu_B}{2mc^2} \epsilon_{rsl} \partial_s [E_r (\Pi_{ijl}^S - u_i u_j S_l)] \\ + \frac{\mu_B}{2mc^2} \epsilon_{rkp} E_r [\epsilon_{kli} (J_{jp}^S - u_j S_p) + \epsilon_{klj} (J_{ip}^S - u_i S_p)] B_l \\ - \frac{\mu_B}{2mc^2} \epsilon_{rsl} [\partial_i E_s (\Pi_{jrl}^S - u_j J_{rl}^S) + \partial_j E_s (\Pi_{irl}^S - u_i J_{rl}^S)] \\ - \frac{\mu_B}{2mc^2} u_i \epsilon_{rsl} \partial_s [E_r (J_{il}^S - u_j S_l)] - \frac{\mu_B}{2mc^2} u_j \epsilon_{rsl} \partial_s [E_r (J_{il}^S - u_i S_l)] = 0, \end{aligned} \qquad (5.10) $$

where we introduced a new average velocity and a new spin current

$$ \bar{\mathbf{u}} = \mathbf{u} - \frac{\mu_B}{2mc^2n} E \times S, \quad \bar{J}_{ij}^S = J_{ij}^S + \frac{\mu_B}{2mc^2} \epsilon_{ijk} E_k n. \qquad (5.11) $$

The above corrections reflect the modification of the velocity owing to the spin-orbit coupling. Indeed, the average velocity can be immediately obtained from the velocity phase-space function, equation (4.7), yielding

$$ \bar{\mathbf{u}} = \frac{1}{nm} \mathrm{tr} \left[ \int \mathbf{V}(\mathbf{r}, \boldsymbol{\pi}) F \, d\boldsymbol{\pi} \right] = \mathbf{u} - \frac{\mu_B}{2mc^2n} E \times S, \qquad (5.12) $$

where $F$ is the $2 \times 2$ distribution function defined in equation (3.6). The same holds for the spin current operator, which is defined as follows:

$$ \hat{J}_{ij}^S = \hat{v}_i \sigma_j = \frac{\hat{\Pi}_i}{m} \sigma_j - \frac{\mu_B}{4mc^2} [(\hat{E} \times \boldsymbol{\sigma})_i \sigma_j + \sigma_j (\hat{E} \times \boldsymbol{\sigma})_i], \qquad (5.13) $$

where we symmetrized the operator so that it is Hermitian. Then the associated phase-space function,

$$ \hat{J}_{ij}^S(\mathbf{r}, \boldsymbol{\pi}) = \frac{\hat{\Pi}_i}{m} \sigma_j - \frac{\mu_B}{4mc^2} [(\hat{E} \times \boldsymbol{\sigma})_i \sigma_j + \sigma_j (\hat{E} \times \boldsymbol{\sigma})_i], \qquad (5.14) $$

can be used to determine the spin current

$$ \bar{J}_{ij}^S = \frac{1}{m} \mathrm{tr} \left[ \int \hat{J}_{ij}^S(\mathbf{r}, \boldsymbol{\pi}) F \, d\boldsymbol{\pi} \right] = \hat{J}_{ij}^S + \frac{\mu_B}{2mc^2} \epsilon_{ijk} E_k n. \qquad (5.15) $$

<footnote n="2">²Strictly speaking, a pressure tensor should be defined in terms of the velocity fluctuations $w_i w_j$, but this would unduly complicate the notation. Thus, we stick to the above definition of $\Pi_{ij\alpha}$ while still using the term 'pressure' for this quantity.</footnote>
As is always the case for hydrodynamic models, some further hypotheses are needed to close the above set of equations (5.6)–(5.10). A particularly interesting strategy, based on the maximum entropy method (MEP), was used in a previous work [31] to close the set of hydrodynamic equations in the presence of the sole Zeeman interaction. Unfortunately, when one adds the spin-orbit interaction, the MEP does not provide any conclusive analytical results (the difficulty arises from the fact that the spin-orbit interaction couples all the components of the velocity). However, an intuitive closure can be found by inspecting the evolution equation (5.10) for the pressure tensor. There, most spin-dependent terms cancel if we set

$$J_{i\alpha}^{S} = u_{i}S_{\alpha} \quad \text{and} \quad \Pi_{ij\alpha}^{S} = u_{i}J_{j\alpha}^{S}. \qquad (5.16)$$

The physical interpretation of the above equations is that the spin of a particle is simply transported along the mean fluid velocity. This is of course an approximation that amounts to neglecting some spin-velocity correlations [27]. With this assumption, equation (5.9) and the definition of the spin-pressure $\Pi_{ij\alpha}$ are no longer necessary. The system of fluid equations simplifies to

$$\frac{\partial n}{\partial t} + \nabla_r \cdot (n\bar{u}) = 0, \qquad (5.17)$$

$$\frac{\partial S_{\alpha}}{\partial t} + \partial_i(u_i S_{\alpha}) - \frac{\mu_B}{2mc^2}(\nabla \times nE)_{\alpha} + \frac{e}{m}\left[S \times \left(B - \frac{1}{2c^2}u \times E\right)\right]_{\alpha} = 0, \qquad (5.18)$$

$$\frac{\partial u_i}{\partial t} + u_j(\nabla_j u_i) + \frac{1}{nm}\nabla_j P_{ij} + \frac{e}{m}[E_i + (\bar{u} \times B)_i] + \frac{e}{nm^2}S_{\alpha}(\partial_i B_{\alpha}) + \frac{\mu_B}{2mc^2n}\epsilon_{jkl}[E_j(\partial_k u_i) - u_k(\partial_i E_j)]S_l = 0 \qquad (5.19)$$

$$\text{and} \quad \frac{\partial P_{ij}}{\partial t} + u_k \partial_k P_{ij} + P_{jk} \partial_k u_i + P_{ik} \partial_k u_j + P_{ij} \partial_k u_k + \partial_k Q_{ijk} + \frac{e}{m} [\epsilon_{kli} P_{jk} + \epsilon_{klj} P_{ik}] B_l = 0. \qquad (5.20)$$

In order to complete the closure procedure, one can proceed in the same way as is usually done for spinless fluids, for instance by supposing that the system is isotropic and adiabatic. The isotropy condition imposes that $P_{ij} = (P/3)\delta_{ij}$, where $\delta_{ij}$ is the Kronecker delta, whereas the adiabaticity condition requires that the heat flux $Q_{ijk}$ vanish. In that case, the pressure takes the usual form of the equation of state of an adiabatic system, i.e. $P = \text{const.} \times n^{(D+2)/D}$ (where $D$ is the dimensionality of the system), which replaces equation (5.20). In summary, equations (5.17)–(5.19), together with the preceding expression for the pressure, constitute a closed system of hydrodynamic equations with spin-orbit effects.

## 6. Conclusion

Phase-space methods can be applied to condensed matter and nanophysics to model the electron dynamics in either the quantum or the semiclassical regime. Several studies were performed in the past but neglected spin effects [8,44,47]. In this paper, we showed that phase-space methods can be conveniently generalized to include the spin dynamics at different orders. In an earlier work, we had developed a phase-space model that includes the lowest order spin term (the Zeeman effect), but neglects all relativistic corrections (spin-orbit coupling, Darwin term, mass correction, ...). Here, we considered the case where both the Zeeman and the spin-orbit interaction are present (other relativistic corrections could be added with relative ease). The spin-orbit interaction plays an important role, for instance, in ultrafast spectroscopy experiments on magnetic nano-objects, where the electron spin is known to interact with the incident laser field and with the self-consistent field generated by the electron gas.

We first derived a four-component Wigner equation to describe the quantum dynamics of a system of spin-$\frac{1}{2}$ particles. These equations, together with the appropriate Maxwell equations,
form a fully quantum self-consistent model to study the spin and charge dynamics in the mean-field approximation. This model is not limited to the linear response, but can deal with nonlinear effects, which are often important, particularly for large incident laser powers.

Next, using a semiclassical expansion to first order in $ħ$, we obtained a four-component Vlasov equation. The orbital part of the motion is classical, i.e. the particles follow the classical phase-space trajectories, while the spin degrees of freedom are treated in a fully quantum fashion (two-dimensional Hilbert space). These spin Vlasov equations constitute a good approximation of the quantum electron dynamics for wavelengths larger than the typical inter-electron distance.

The corresponding hydrodynamic equations were derived by taking velocity moments of the phase-space distribution functions. The spin-orbit interaction modifies considerably our earlier hydrodynamic equations [31], where the only spin effect was the Zeeman interaction. We proposed a simple, intuitive closure for the hydrodynamic equations whereby the spin is simply transported along by the fluid velocity of the electrons.

The present models (Vlasov and hydrodynamic) constitute two valuable tools to study the intertwined spin and charge dynamics in condensed matter systems and nano-objects. The challenge now is to implement these models in performing numerical codes and to study the electron dynamics in realistic nanoscale systems that are relevant to current experiments.

**Authors' contributions.** J.H. performed the analytical calculations and drafted the manuscript. P.-A.H. and G.M. conceived the study, checked the calculations, and contributed to the writing of the manuscript.

**Competing interests.** The authors declare no competing interests.

**Funding.** We thank the Agence Nationale de la Recherche, project Labex 'Nanostructures in Interaction with their Environment', for financial support.

## References

1. Wigner E. 1932 On the quantum correction for thermodynamic equilibrium. *Phys. Rev.* **40**, 749–759. (doi:10.1103/PhysRev.40.749)

2. Case WB. 2008 Wigner functions and Weyl transforms for pedestrians. *Am. J. Phys.* **76**, 937–946. (doi:10.1119/1.2957889)

3. Cohen L. 2013 *The Weyl operator and its generalization*. Basle, Switzerland: Springer.

4. Smithey DT, Beck M, Raymer MG, Faridani A. 1993 Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomography: application to squeezed states and the vacuum. *Phys. Rev. Lett.* **70**, 1244–1247. (doi:10.1103/PhysRevLett.70.1244)

5. Dittrich T, Gòmez EA, Pachòn LA. 2010 Semiclassical propagation of Wigner functions. *J. Chem. Phys.* **132**, 214102. (doi:10.1063/1.3425881)

6. Heller EJ. 1976 Wigner phase space method: analysis for semiclassical applications. *J. Chem. Phys.* **65**, 1289–1298. (doi:10.1063/1.433238)

7. Bertoni A, Bordone P, Brunetti R, Jacoboni C. 1999 The Wigner function for electron transport in mesoscopic systems. *J. Phys. Condens. Matter* **31**, 5999–6012. (doi:10.1088/0953-8984/11/31/308)

8. Jasiak R, Manfredi G, Hervieux P-A, Haefele M. 2009 Quantum-classical transition in the electron dynamics of thin metal films. *New J. Phys.* **11**, 063042. (doi:10.1088/1367-2630/11/6/063042)

9. Haas F. 2011 *Quantum plasmas—an hydrodynamic approach*. Heidelberg, Germany: Springer.

10. Brewczyk M, Rzażewski K, Clark CW. 1997 Multielectron dissociative ionization of molecules by intense laser radiation. *Phys. Rev. Lett.* **78**, 191–194. (doi:10.1103/PhysRevLett.78.191)

11. Banerjee A, Harbola MK. 2000 Hydrodynamic approach to time-dependent density functional theory; response properties of metal clusters. *J. Chem. Phys.* **113**, 5614–5623. (doi:10.1063/1.1290610)

12. Domps A, Reinhard P-G, Suraud E. 1998 Theoretical estimation of the importance of two-electron collisions for relaxation in metal clusters. *Phys. Rev. Lett.* **81**, 5524–5527. (doi:10.1103/PhysRevLett.81.5524)

13. Manfredi G, Hervieux P-A, Haas F. 2012 Nonlinear dynamics of electron-positron clusters. *New J. Phys.* **14**, 75012. (doi:10.1088/1367-2630/14/7/075012)
14. Crouseilles N, Hervieux P-A, Manfredi G. 2008 Quantum hydrodynamic model for the nonlinear electron dynamics in thin metal films. *Phys. Rev. B* **78**, 155412. (doi:10.1103/PhysRevB.78.155412)

15. Shukla PK, Eliasson B. 2006 Formation and dynamics of dark solitons and vortices in quantum electron plasmas. *Phys. Rev. Lett.* **96**, 245001. (doi:10.1103/PhysRevLett.96.245001)

16. Shukla PK, Eliasson B. 2010 Nonlinear aspects of quantum plasma physics. *Physics-Uspekhi* **53**, 51–76. (doi:10.3367/UFNe.0180.201001b.0055)

17. Haas F, Manfredi G, Shukla PK, Hervieux P-A. 2009 Breather mode in the many-electron dynamics of semiconductor quantum wells. *Phys. Rev. B* **80**, 073301. (doi:10.1103/PhysRevB.80.073301)

18. Puente A, Casas M, Serra L. 2000 A semiclassical approach to the ground state and density oscillations of quantum dots. *Physica E* **8**, 387–397. (doi:10.1016/S1386-9477(99)00042-9)

19. Serra L, Puente A. 2001 Magnetic Thomas–Fermi–Weizsäcker model for quantum dots: a comparison with Kohn–Sham ground states. *Eur. Phys. J. D* **14**, 77–81. (doi:10.1007/s100530170237)

20. Morandi O, Hervieux P-A, Manfredi G. 2009 Ultrafast magnetization dynamics in diluted magnetic semiconductors. *New J. Phys.* **11**, 73010. (doi:10.1088/1367-2630/11/7/073010)

21. Morandi O, Hervieux P-A, Manfredi G. 2010 Time-dependent model for diluted magnetic semiconductors including band structure and confinement effects. *Phys. Rev. B* **81**, 155309. (doi:10.1103/PhysRevB.81.155309)

22. Bigot J-Y, Vomir M, Beaurepaire E. 2009 Coherent ultrafast magnetism induced by femtosecond laser pulses. *Nat. Phys.* **5**, 515–520. (doi:10.1038/nphys1285)

23. Krieger K, Dewhurst JK, Elliott P, Sharma S, Gross EKU. 2015 Laser-induced demagnetization at ultrashort time scales: predictions of TDDFT. *J. Chem. Theory Comput.* **11**, 4870–4874. (doi:10.1021/acs.jctc.5b00621)

24. Thomas MW, Snider RF. 1970 Boltzmann equation and angular momentum conservation. *J. Stat. Phys.* **2**, 61–81. (doi:10.1007/BF01009711)

25. Snider RF, Lewchuk KS. 1967 Irreversible thermodynamics of a fluid system with spin. *J. Chem. Phys.* **46**, 3163–3172. (doi:10.1063/1.1841187)

26. Zamanian J, Marklund M, Brodin G. 2010 Scalar quantum kinetic theory for spin-$\frac{1}{2}$ particles: mean field theory. *New J. Phys.* **12**, 043019. (doi:10.1088/1367-2630/12/4/043019)

27. Zamanian J, Stefan M, Marklund M, Brodin G. 2010 From extended phase space dynamics to fluid theory. *Phys. Plasmas* **17**, 102109. (doi:10.1063/1.3496053)

28. Asenjo FA, Zamanian J, Marklund M, Brodin G, Johansson P. 2012 Semi-relativistic effects in spin-$\frac{1}{2}$ quantum plasmas. *New J. Phys.* **14**, 073042. (doi:10.1088/1367-2630/14/7/073042)

29. Arnold A, Steinriick H. 1989 The electromagnetic Wigner equation for an electron with spin. *J. Appl. Math. Phys.* **40**, 793–615. (doi:10.1007/BF00945803)

30. Bialynicki-Birula I. 2014 Relativistic Wigner functions. *EPJ Web Conf.* **78**, 1001. (doi:10.1051/epjconf/20147801001)

31. Hurst J, Morandi O, Manfredi G, Hervieux P-A. 2014 Semiclassical Vlasov and fluid models for an electron gas with spin effects. *Eur. Phys. J. D* **68**, 176. (doi:10.1140/epjd/e2014-50205-5)

32. Serimaa OT, Javanainen J, Varrò S. 1986 Gauge-independent Wigner functions: general formulation. *Phys. Rev. A* **33**, 2913–2927. (doi:10.1103/PhysRevA.33.2913)

33. Royal JE, Bartlett MS, Wigner E, McCoy NH, Dirac PAM, Groenewold HJ, Feynman RP. 1949 Quantum mechanics as a statistical theory. *Math. Proc. Camb. Philos. Soc.* **45**, 99–124. (doi:10.1017/S0305004100000487)

34. Dixit A, Hinsberger Y, Zamanian J, Manfredi G, Hervieux P-A. 2013 Lagrangian approach to the semirelativistic electron dynamics in the mean-field approximation. *Phys. Rev. A* **88**, 032117. (doi:10.1103/PhysRevA.88.032117)

35. Weyl H. 1927 Quantenmechanik und Gruppentheorie. *Z. Phys.* **46**, 1–46. (doi:10.1007/BF02055756)

36. Husimi K. 1940 Some formal properties of the density matrix. *Proc. Phys. Math. Soc. Jpn* **22**, 264–314.

37. Stratonovich RL. 1956 A gauge invariant analog of the Wigner distribution. *Sov. Phys. D* **1**, 414–418.

38. Müller M. 1999 Product rule for gauge invariant Weyl symbols and its application to the semiclassical description of guiding centre motion. *J. Phys. A, Math. Gen.* **32**, 1035–1052. (doi:10.1088/0305-4470/32/6/014)
39. Hinschberger Y, Hervieux P-A. 2012 Foldy-Wouthuysen transformation applied to the interaction of an electron with ultrafast electromagnetic fields. *Phys. Lett. A* **376**, 813–819. (doi:10.1016/j.physleta.2012.01.023)

40. Strange P. 1998 *Relativistic quantum mechanics: with applications in condensed matter and atomic physics*. Cambridge, UK: Cambridge University Press.

41. Barletti L. 2003 Wigner envelope functions for electron transport in semiconductor devices. *Transp. Theory Stat. Phys.* **32**, 253–277. (doi:10.1081/TT-120024764)

42. Morandi O, Schürrer F. 2011 Wigner model for quantum transport in graphene. *J. Phys. A, Math. Theor.* **44**, 265301. (doi:10.1088/1751-8113/44/26/265301)

43. Manfredi G. 2013 Non-relativistic limits of Maxwell's equations. *Eur. J. Phys.* **34**, 859–871. (doi:10.1088/0143-0807/34/4/859)

44. Manfredi G, Hervieux P-A, Yin Y, Crouseilles N. 2010 *Collective electron dynamics in metallic and semiconductor nanostructures*, pp. 1–44. Berlin, Germany: Springer. (doi:10.1007/978-3-642-04650-6_1).

45. Thomas LH. 1926 Motion of the spinning electron. *Nature* **117**, 514. (doi:10.1038/117514a0)

46. Dragan A, Odrzygóźdź T. 2013 A half-page derivation of the Thomas precession. *Am. J. Phys.* **81**, 631–632. (doi:10.1119/1.4807564)

47. Calvayrac F, Reinhard P-G, Suraud E, Ullrich CA. 2000 Nonlinear electron dynamics in metal clusters. *Phys. Rep.* **337**, 493–578. (doi:10.1016/S0370-1573(00)00043-0)