Phase-space modelling of solid-state plasmas

TL;DR

This paper demonstrates the effectiveness of phase-space models in describing the complex quantum, spin, and relativistic effects in solid-state plasmas within nano-objects, offering a versatile alternative to wave function methods.

Contribution

It introduces a phase-space approach to model electron dynamics in nano-objects, incorporating quantum, spin, relativistic effects, and dissipation, with applications to linear and nonlinear responses.

Findings

Spin effects modify linear response of electron gas

Nonlinear electron dynamics in thin metal films analyzed

Phase-space models capture quantum and relativistic effects

Abstract

Conduction electrons in metallic nano-objects ($\rm 1\,nm = 10^{-9}\, m$) behave as mobile negative charges confined by a fixed positively-charged background, the atomic ions. In many respects, this electron gas displays typical plasma properties such as screening and Langmuir waves, with more or less pronounced quantum features depending on the size of the object. To study these dynamical effects, the mathematical artillery of condensed-matter theorists mainly relies on wave function $\psi(r,t)$ based methods, such as the celebrated Hartree-Fock equations. The theoretical plasma physicist, in contrast, lives and breaths in the six-dimensional phase space, where the electron gas is fully described by a probability distribution function $f(r,p,t)$ that evolves according to an appropriate kinetic equation. Here, we illustrate the power and flexibility of the phase-space approach to describe the electron dynamics in small nano-objects. Starting from classical and semiclassical scenarios, we progressively add further features that are relevant to solid-state plasmas: quantum, spin, and relativistic effects, as well as collisions and dissipation. As examples of applications, we study the spin-induced modifications to the linear response of a homogeneous electron gas and the nonlinear dynamics of the electrons confined in a thin metal films of nanometric dimensions.

Tables

Table 1. Table 1: Typical parameters for electrons in gold at room temperature. The values show that the relevant time, space, and energy scales are respectively, the femtosecond, the nanometer, and the electron-volt.

Parameter	Value	Units
$r_s$	0.16	nm
$n$	$59$	${\mathrm{nm}}^{-3}$
$T$	300	K
$2\pi /{\omega }_p$	0.46	fs
$T_F$	64 000	K
$E_F$	5.6	eV
$\hslash {\omega }_p$	9.1	eV
$v_F$	$1.4$	$\mathrm{fm}/\mathrm{ns}$
${\lambda }_{TF}$	0.1	nm
${\tau }_{ee}$	$\approx 50$	fs
${\tau }_{e-ph}$	$\approx 5$	ps
$g_q$	5.5	–
$r_s/a_0$	3	–
$v_F/c$	$0.0047$	–

Equations

$$T_F=\frac{\hslash }{2m{k}_B}{\left(3{\pi }^{2}n\right)}^{2/3}.$$

$$v_F=\sqrt{\frac{2E_F}{m}}=\frac{\hslash }{m}{\left(3{\pi }^{2}n\right)}^{1/3},$$

$$g_q=\frac{E_{int}}{E_F}\sim \frac{e^{2}m}{{\hslash }^2\epsilon n^{1/3}}\sim {\left(\frac{1}{n{\lambda }_{TF}^3}\right)}^{2/3}\sim {\left(\frac{\hslash {\omega }_p}{E_F}\right)}^2\sim \frac{r_s}{a_0},$$

$$g_c=\frac{e^2{n}^{1/3}}{{\epsilon }_0{k}_B{T}_e}.$$

$$\frac{{\nu }_{ee}}{{\omega }_p}\sim g_q^{-1/2}{\left(\frac{T_e}{T_F}\right)}^2.$$

$$C_e(T_e)\frac{\partial T_e}{\partial t}$$

$$C_l\frac{\partial T_l}{\partial t}$$

$${\Psi }^N=\Psi \left({\mathbf{r}}_1,{\mathbf{r}}_2,\cdots ,{\mathbf{r}}_n,t\right),$$

$$i\hslash \frac{\partial {\Psi }^N}{\partial t}=-\frac{{\hslash }^2}{2m}\sum _{i=1}^N\left[{\nabla }_i^2{\Psi }^N+\sum _{k=1,k\ne i}^N\frac{1}{2}\mathcal{V}\left(|{\mathbf{r}}_i-{\mathbf{r}}_k|\right){\Psi }^N\right],$$

$${\Psi }^N={\Psi }_1\left({\mathbf{r}}_1,t\right){\Psi }_2\left({\mathbf{r}}_2,t\right)\cdots {\Psi }_N\left({\mathbf{r}}_N,t\right).$$

$$i\hslash \frac{\partial {\Psi }_{\alpha }(\mathbf{r},t)}{\partial t}=-\frac{{\hslash }^2}{2m}{\nabla }^2{\Psi }_{\alpha }(\mathbf{r},t)+\underset{V_H}{\underbrace{\frac{e^2}{4\pi {ϵ}_0}\left[\sum _{{\alpha }^{\prime }=1}^N\int \frac{|{\Psi }_{{\alpha }^{\prime }}({\mathbf{r}}^{\prime },t){|}^2}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }\right]}}{\Psi }_{\alpha }(\mathbf{r},t)$$

$$-\underset{V_{sic}}{\underbrace{\frac{e^2}{4\pi {ϵ}_0}\int \frac{|{\Psi }_{\alpha }({\mathbf{r}}^{\prime },t){|}^2}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }}}{\Psi }_{\alpha }(\mathbf{r},t).$$

$$\displaystyle\left\{\begin{array}[]{lcl}\displaystyle i\hbar\frac{\partial\Psi_{\alpha}(\bm{r},t)}{\partial t}&=&\displaystyle-\frac{\hbar^{2}}{2m}\bm{\nabla}^{2}\Psi_{\alpha}(\bm{r},t)-eV_{H}(\bm{r},t)\Psi_{\alpha}(\bm{r},t),\\ \\ \displaystyle\bm{\nabla}^{2}V_{H}(\bm{r},t)&=&\displaystyle\frac{e}{\epsilon_{0}}\sum_{\alpha=1}^{N}\left|\Psi_{\alpha}(\bm{r},t)\right|^{2}.\end{array}\right.$$

$$\Psi \left({\mathbf{r}}_2,{\mathbf{r}}_1,\cdots ,{\mathbf{r}}_N,t\right)=-\Psi \left({\mathbf{r}}_1,{\mathbf{r}}_2,\cdots ,{\mathbf{r}}_N,t\right).$$

$${\Psi }^N$$

$$i\hslash \frac{\partial {\Psi }_{\alpha }(\mathbf{r},t)}{\partial t}$$

$$-\frac{e^2}{4\pi {ϵ}_0}\sum _{{\alpha }^{\prime }=1}^N{\Psi }_{{\alpha }^{\prime }}(\mathbf{r},t)\int \frac{{\Psi }_{{\alpha }^{\prime }}^{\ast }({\mathbf{r}}^{\prime },t){\Psi }_{\alpha }({\mathbf{r}}^{\prime },t)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }.$$

$$E[n]=T[n]+E_{ext}[n]+E_{int}[n]=T[n]+E_{ext}[n]+E_H[n]+E_X[n]+E_C[n],$$

$$E_H[n]=\int V_H[n]nd\mathbf{r}=\frac{e^2}{4\pi {ϵ}_0}\int \int \frac{n(\mathbf{r})n({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime }.$$

$$T[n]=\frac{3}{10}(3{\pi }^2{)}^{2/3}\frac{{\hslash }^2}{m}\int n^{5/3}d\mathbf{r}.$$

$$-\frac{{\hslash }^2}{2m}{\nabla }^2{\varphi }_j(\mathbf{r})+V_{eff}[n(\mathbf{r})]{\varphi }_j(\mathbf{r})={ϵ}_j{\varphi }_j(\mathbf{r}),$$

$$V_X=-\frac{e^2}{4\pi {ϵ}_0}{\left(\frac{3}{\pi }\right)}^{1/3}{n}^{1/3}.$$

$$-\frac{{\hslash }^2}{2m}{\nabla }^2{\varphi }_j(\mathbf{r},t)+V_{eff}[n(\mathbf{r})]{\varphi }_j(\mathbf{r},t)=i\hslash \frac{\partial {\varphi }_j(\mathbf{r},t)}{\partial t}.$$

$$\rho =\sum _{\alpha =1}^N{p}_{\alpha }\left|{\Psi }_{\alpha }⟩⟨{\Psi }_{\alpha }\right|,$$

$$f\left(\mathbf{r},\mathbf{p},t\right)$$

$$\frac{\partial f}{\partial t}+\frac{1}{m}\mathbf{p}\cdot \nabla f=\frac{ie}{\hslash }\frac{1}{{\left(2\pi \hslash \right)}^3}\int d\mathbf{\lambda }d{\mathbf{p}}^{\prime }{e}^{\frac{i\left(\mathbf{p}-{\mathbf{p}}^{\prime }\right)\cdot \mathbf{\lambda }}{\hslash }}\left[V\left({\mathbf{r}}_{+}\right)-V\left({\mathbf{r}}_{-}\right)\right]f(\mathbf{r},{\mathbf{p}}^{\prime },t),$$

$$n(\mathbf{r},t)$$

$$\mathbf{J}(\mathbf{r},t)$$

$$\frac{\partial f}{\partial t}+\frac{\mathbf{p}}{m}\cdot \nabla f+e\nabla V\cdot {\nabla }_{\mathbf{p}}f=O({\hslash }^2).$$

$$f\left(\mathbf{r},\mathbf{v},t\right)$$

$$\frac{\partial f}{\partial t}+\frac{1}{m}\left(\mathbf{\pi }+\mathbf{\Delta }\tilde{\mathbf{\pi }}\right)\cdot \nabla f-\frac{e}{m}{\left[m\tilde{\mathbf{E}}+\left(\mathbf{\pi }+\mathbf{\Delta }\tilde{\mathbf{\pi }}\right)\times \tilde{\mathbf{B}}\right]}_i{\partial }_{{\pi }_i}f=0,$$

$$\mathbf{\Delta }\tilde{\mathbf{\pi }}=-i\hslash e{\partial }_{\mathbf{\pi }}\times \left[{\int }_{-1/2}^{1/2}d\tau \tau \mathbf{B}\left(\mathbf{r}+i\hslash \tau {\partial }_{\mathbf{\pi }}\right)\right]$$

Full text

1 Introduction – Plasmas and nanophysics

The physics of metallic nano-objects has stimulated a vast amount of scientific interest in the last two decades, both for fundamental research and for potential technological applications to nanophotonics Stockman2011 ; Moreau2012 , physical chemistry Daniel2004 , and even biology and medicine Hainfeld2004 ; TatsuroEndo2006 . Metallic nano-objects are mesoscopic systems composed of a relatively small number of atoms, typically between a few tens and several millions. The typical size of these systems ranges from a few to several hundred nanometres ($\rm 1\,nm=10^{-9}\,m$), with properties that are intermediate between those of molecules and bulk solids. They can be synthesized in different geometries, ranging from spherical nanoparticles, to thin films, nanorods, cubes or pyramids (see Fig. 1).

Metallic nano-objects also present a fundamental interest as large objects that still display some quantum features Scholl2012 ; Luo2013 ; Tame2013 ; Raza2013 . Quantum effects arise because at metallic densities ($n\approx 10^{28}\,\rm m^{-3}$) the electrons are so closely packed together that their wave functions overlap even at room temperature, but also because of finite-size effects and the related presence of surfaces where the density varies significantly across very short distances (less than 1 nm).

Many of the optical properties of metallic nano-objects are mediated by their conduction electrons. The conduction electrons in a metal are delocalized over comparatively large distances, so that they can be treated, at least to first approximation, as an electron plasma confined by a fixed positively-charged background, the atomic ions. Such plasma is known to display collective effects, giving rise to typical plasma properties such as pronounced resonances near the plasma frequency $\omega_{p}=\sqrt{ne^{2}/m\epsilon_{0}}$. Due to the very large electron density, the plasma period $\tau_{p}=2\pi/\omega_{p}$ lies in the femtosecond range ($1\rm\,fs=10^{-15}\,\rm s$). Thus, it is not surprising that the study of collective electron resonances in nano-objects became possible following the development of ultrafast spectroscopy techniques (“pump-probe experiments”) in the femtosecond regime Voisin2000 ; Bigot2000 . Today, plasmonic resonances are actively investigated for applications to such diverse areas as biomedicine Salata2004 ; Loomba2013 and high-harmonic generation Butet2010 ; Singhal2010 .

Spectroscopy experiments.

Pump-probe experiments involve a sequence of two ultrafast laser pulses. The first pulse (the pump) is more intense and excites the electron dynamics, while the second (the probe) is weaker and serves as a measure the time-dependent response. A typical experimental scenario is depicted schematically in Fig. 2 for the case of a ferromagnetic nano-object, where the charge and the spin dynamics are closely intertwined. In the first few tens of femtoseconds, the electric field of the laser couples coherently to the electron charges (direct coupling between the laser magnetic field and the electron spins occurs only at very high intensities Bigot2009 ). Some of the laser energy is absorbed by the electrons, which are then driven out of equilibrium. After the laser pulse is switched off, the electrons evolve coherently under the action the self-consistent mean field. A typical response mode is the so-called surface plasmon, which is a rigid-body oscillation of the entire electron cloud around the fixed ion lattice. For an ideal spherical nanoparticle, the oscillation frequency of the surface plasmon mode is equal too $\omega_{p}/\sqrt{3}$. The surface plasmon eventually damps away by coupling to the single-particle modes, an effect that is akin to the standard Landau damping of ordinary plasmas. Such Landau damping, which occurs on a very short time scale ($10-50~{}\rm fs$), was observed experimentally in gold nanoparticles Lamprecht1999 and was studied theoretically in several works Kreibig1995 ; Molina2002 ; Shahbazyan2016 .

However, electrons possess not only a charge, but also a spin – i.e., a magnetic moment – which is a crucial property for ferromagnetic nano-objects. A remarkable result, first observed by Beaurepaire and Bigot in nickel films Beaurepaire1996 ; Bigot2009 , is the ultrafast loss of magnetization occurring during the first $100\,\rm fs$ that follows the laser pulse. Over twenty years after its discovery, there is no general agreement about the underlying causes of this ultrafast demagnetization, which was attributed to various mechanisms such as the spin-orbit interaction Krieger2015 ; Stamenova2016 (a semi-relativistic effect, as we shall see) or the superdiffusive electron transport induced by the laser field Battiato2010 .

During these initial ultrafast processes, the ionic background remains frozen and the electron energy distribution is nonthermal due to the laser excitation. For longer times ($t>50~{}\rm fs$, see Fig. 2), the laser energy is redistributed amongst the electrons through electron-electron collisions and spin-flip processes, leading to the internal thermalization of the electron gas. During this stage the electron population can still be much hotter (up to several thousand degrees) than the ion lattice, which remains close to room temperature. Finally, on the picosecond timescale, electron-phonon scattering leads to the exchange of energy between the electrons and the ion lattice and the subsequent thermalization of the full system.

Theory and modelling.

From a theoretical point of view, the description of the electronic dynamics in metallic nano-structures is a very complex challenge. Exact approaches based on the $N$-body Schrödinger equation are necessarily limited to a very small number of particles. Although such few- or even single-electron systems can nowadays be realized in the laboratory, in most practical situations a great many electrons are involved Pereira2004 ; Muller2004 . In that case, self-consistent effects arising from the Coulomb interactions (between all the electrons) play a crucial role on the dynamics. Several theoretical and computational studies, which treat the many-body dynamics in an approximate way, focused on the linear and nonlinear electron response. Earlier works were based on macroscopic phenomenological models Guillon2003 ; Aeschlimann2000 ; Rethfeld2002 that employed Boltzmann-type equations within the framework of the Fermi-liquid theory Pines1995 . Studies based on microscopic models are more recent and limited to relatively small systems, due to their considerable computational complexity. In the quantum regime, the ultrafast electron dynamics in metallic clusters was studied by Calvayrac et al. Calvayrac2000 and more recently by Teperik et al. Teperik2013 using time-dependent density functional theory. The many-particle quantum dynamics of the electron gas in a thin metal film was studied by Schwengelbeck et al. Schwengelbeck2000 within the framework of the time-dependent Hartree-Fock approximation. All the above-mentioned methods, being essentially quantum, are based on the evolution of a set of wave functions $\psi_{j}(\bm{r},t)$, each obeying a Schrödinger-like equation.

In this work, we will review a possible alternative that relies on the use of phase-space models inspired from classical plasma physics, for which the system is governed by a probability distribution function $f(\bm{r},\bm{p},t)$ that evolves according to a kinetic equation. Indeed, the semiclassical limit of the above-mentioned quantum models is the self-consistent Vlasov-Poisson system, largely employed in plasma physics.

The Vlasov-Poisson model itself was used by several authors to model semi-classically the electron dynamics in metal clusters Calvayrac2000 ; Daligault2003 ; Fomichev1999 and in thin metal films Zaretsky2004 ; Manfredi2005film . These works were later extended to the quantum regime using the Wigner phase-space description Jasiak2009 . The Wigner representation is a way to express standard quantum mechanics in a classical phase-space language and is suitable to treat both single-particle and many-particle systems. It is often more intuitive than the standard Schrödinger approach, especially for problems where semiclassical considerations are important. For these reasons, it is used in many areas of quantum physics, including quantum optics Smithey1993 , semiclassical analysis Heller1976 ; Dittrich2010 , electronic transport Bertoni1999 , nonlinear electron dynamics Jasiak2009 , and quantum plasma theory Haas2011 . It is also the starting point for the construction of 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 Brewczyk1997 , metal clusters and nanoparticles Domps1998 ; Banerjee2000 ; Manfredi2012 , thin metal films Crouseilles2008 , quantum plasmas Shukla2006 ; Shukla2010 , and semiconductors Haas2009 .

The above studies included the electron charge, but not its spin. However, it is well known that spin effects (particularly the Zeeman interaction and the spin-orbit coupling) can play a decisive role in nanometric systems such as semiconductor quantum dots Puente2000 ; Serra2001 or diluted magnetic semiconductors Morandi2009 ; Morandi2010 . The coupling between the spin degrees of freedom and the electron’s orbital motion is of the utmost importance in many experimental studies involving magnetized nano-objects, such as the above-mentioned laser-induced ultrafast demagnetization Bigot2009 .

Phase-space models based on the Boltzmann equation Thomas1970 , and the corresponding fluid models Snider1967 , 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 Zamanian2010nj consists in 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 Zamanian2010nj . Semiclassical Zamanian2010 and hydrodynamic Asenjo2012 spin equations were also derived from those models, as well as other relativistic effects.

An alternative approach is to use a matrix form for the phase-space distribution function Arnold1989 , which originates from the $2\times 2$ density matrix for spin-1/2 particles. Using this approach, the corresponding Wigner equations were derived from the full Dirac theory Bialynicki-Birula2014 . In their semiclassical limit, these equations give rise to a matrix spin-Vlasov equation, which treats the electron motion in a classical fashion while preserving the intrinsically quantum character of the spin degrees of freedoms Hurst2014 ; Hurst2017 . This approach was recently used to study the generation of spin currents in ferromagnetic thin films Hurst2018 .

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 separates clearly the orbital motion from the spin dynamics and leads to simpler and more transparent hydrodynamic models. From a computational point of view, the extended phase space method is more apt to be simulated using particle-in-cell codes, because the corresponding distribution function is transported along classical trajectories in the extended phase space (which is eight-dimensional: three positions, three velocities, and two spin angular coordinates). In contrast, the matrix Wigner function methods is more naturally amenable to grid-based Vlasov codes, because the corresponding distribution function only depends on the six variables of the ordinary phase space.

Summary.

In this review, we will first recall the basic concepts, parameters and physical mechanisms characterizing the electron plasma in a solid nano-object (Sec. 2). This will be followed by a very short review of wave-function-based methods (Hartree-Fock and density functional theory), which are the golden standard of computer simulations in condensed matter physics (Sec. 3). Then, we will introduce the Wigner phase-space representation of quantum mechanics (Sec. 4) and its extension to spin-1/2 fermions (Sec. 5). One of the attractive features of phase-space methods is that they can incorporate dissipative effects more naturally than wave-function-based methods – this is described in Sec. 6. All the above concepts are illustrated by two examples. First (Sec. 7), we study the linear response of a homogeneous electron gas including spin effects and derive the corresponding spin-dependent dispersion relation, for both Maxwell-Boltzmann and Fermi-Dirac equilibria. Second (Sec. 8), we summarize a series of computational studies on the nonlinear electron dynamics in thin metal films carried out in our research group over the last 15 years, culminating in the recent observation of spin currents in a ferromagnetic nickel film Hurst2018 .

2 Basic concepts, parameters, mechanisms

In this paper, we focus our attention on the theoretical description of the electron dynamics in metallic nano-objects. Metals are condensed-matter systems with the specificity of having a half-filled conduction band. The electrons that belong to the conduction band are not attached to a particle nucleus, but are rather delocalized in the material and behave, to first approximation, as a non-interacting electron gas. This property was exploited by Drude Drude1900 at the beginning of the twentieth century to derive approximate estimations of the electric and thermal conductivities of a metal. A more accurate understanding may be achieved by treating electrons as a one component plasma, i.e., by considering the interactions between the electrons. Because of the large density of solid-state objects, such system display prominent quantum properties, thus deserving the name of “quantum plasmas”.

In a fermion gas, quantum effects become important when the temperature of the gas is comparable to, or smaller than, the Fermi temperature, defined as:

[TABLE]

For metals, $T_{F}\approx 10^{4}\,\rm K$, and therefore conduction electrons are in the quantum regime even at room temperature. Indeed the Fermi-Dirac distribution deviates drastically from the classical Maxwell-Boltzmann distribution for temperatures much lower than the Fermi temperature. Quite often it is a sufficiently good approximation to assume that the electron temperature is equal to zero. In that case, all energy levels up to the Fermi energy $E_{F}=k_{B}T_{F}$ are occupied, whereas all levels with $E>E_{F}$ are empty, and the electron gas is said to be fully degenerate.

One also defines the Fermi velocity:

[TABLE]

and with it, the Thomas-Fermi screening length: $\lambda_{TF}=v_{F}/\omega_{p}$. This is the quantum analog of the classical Debye length $\lambda_{D}=\sqrt{\epsilon k_{B}T_{e}/(ne^{2})}$, and represents the typical distance over which the Coulomb force is screened.

In the fully degenerate regime, the coupling parameter can be written as:

[TABLE]

where $E_{int}$ is the typical Coulomb interaction energy between two electrons situated at a distance $d=n^{-1/3}$, $r_{s}$ is the Wigner-Seitz radius $r_{s}=(3/4\pi n)^{1/3}$ (i.e., the typical volume occupied on average by one electron), and $a_{0}=\hbar^{2}\varepsilon_{0}/(\pi me^{2})$ is the Bohr radius. The coupling parameter discriminates between the collisionless or weakly correlated regime ($g_{q}\ll 1$), where two-body collisions are unimportant and the system can be described in the mean-field approximation, and the collisional or strongly correlated regime ($g_{q}\approx 1$), where two-body correlations cannot be neglected. From Eq. (3), the collisionless regime can be interpreted in various ways as the regime where: (i) the number of electrons in a Thomas-Fermi volume is large, (ii) the plasmon energy $\hbar\omega_{p}$ is small compared to the Fermi energy, or (iii) the Wigner-Seitz radius is small compared to the Bohr radius. The latter ratio is the one that is usually employed in solid-state and nano-physics and is often tabulated in standard textbooks. The quantum coupling parameter $g_{q}$ should be compared to its classical counterpart:

[TABLE]

Note that a quantum plasma becomes less strongly correlated at high densities, whereas the opposite is true for a classical plasma Manfredi2005 (this is because of the density-dependence of the Fermi energy).

Finally, relativistic effects, in the quantum regime, may be quantified by the ratio $v_{F}/c$. This parameter is small ($\approx 10^{-3}$) for electrons in metals, but becomes of order unity for $n\approx 10^{35}\,\rm m^{-3}$, which is compatible with the density of compact astrophysical objects like white dwarfs and neutron stars. Nevertheless, some low-order relativistic effects – in particular, the spin-orbit coupling – may still be important for magnetic nano-objects, where the spin dynamics plays a crucial role. This will be discussed in Sec. 5.

An overview of the typical parameters for gold nano-objects is provided in Table 1.

The three dimensionless parameters $g_{c}$, $g_{q}$ and $T_{e}/T_{F}$ depend on the density and the temperature of the electron gas, and determine four different regions in the $(n,T)$ plane, represented in Fig. 3. We notice that electrons in metals are situated in the quantum and strongly coupled regime and indeed a quick estimation shows that for solid-state densities the coupling parameter is of order unity (see Table 1). Nevertheless, the mean-field approximation is still used in condensed-matter physics, in the form of the Hartree and Hartree-Fock equations, which will be reviewed briefly in Sec. 3.

This is in part due to the fact that electron-electron (e-e) collisions are mitigated in degenerate (quantum) plasmas by the so-called “Pauli blocking” effect, which is a consequence of the exclusion principle Manfredi2005 . Without going into the details, this effect stipulates that, in strongly degenerate plasmas, the electrons below the Fermi energy cannot undergo any collisions, since most quantum states are already occupied. The collision rate can be heuristically estimated as the inverse of the electron lifetime at finite temperature $k_{B}T_{e}/\hbar$, multiplied by the number of electrons above the Fermi energy $\sim T_{e}/T_{F}$, which yields: $\nu_{ee}\sim k_{B}T_{e}^{2}/(\hbar T_{F})$, or in dimensionless units:

[TABLE]

Hence, the collision rate can be small even when the coupling parameter is of order unity, provided the temperature is low enough. This contrasts with the classical case, where the collision rate scales as $\nu_{ee}/\omega_{p}\sim g_{c}^{3/2}$. This rough estimation yields an e-e collision time $\tau_{ee}=\nu_{ee}^{-1}$ of the order of the picosecond, thus much larger than the typical period of plasma oscillations. However, in current pump-probe experiments, the electrons absorb the laser energy very quickly, so that the “temperature” $T_{e}$ appearing in Eq. (5) should actually refer to the average kinetic energy of the resulting out-of-equilibrium electron gas, which can easily attain a few thousand degrees. This brings the theoretical e-e collision time down to 100 fs or less, depending on the intensity of the laser, which is consistent with the experiments Fourment2014 . At the end of this phase, the electrons have thermalized to a temperature much higher than the lattice, which has not yet had the time to interact with the electron gas.

On even longer timescales, the electrons couple to the phonons, i.e., vibrations of the ion lattice. A detailed microscopic description of this phenomenon is very complex, so it it useful to resort to empirical macroscopic models. The two-temperature model (TTM) is a simple description of electron-phonon (e-ph) interactions that considers two thermal baths for the electrons and the lattice, with temperatures $T_{e}$ and $T_{l}$, which interact through a coupling constant $G$, and obey the evolution equations:

[TABLE]

where $\kappa$ is the heat conductivity, $C_{e}$ and $C_{l}$ are respectively the electron and lattice heat capacities, $G$ is the e-ph coupling constant, and $P(t)$ is the power absorbed by the electron gas following the laser excitation. For gold, typical measured values of these parameters are Ekici2008 : $G=2\times 10^{16}\,\rm Wm^{-3}K^{-1}$, $C_{l}=2.1\times 10^{6}\,\rm Jm^{-3}K^{-1}$, and $C_{e}=70\,T_{e}\,\rm Jm^{-3}K^{-1}$. We note that $C_{e}$ is rather close to the heat capacity of an ideal Fermi gas: $C_{\rm Fermi}={\pi^{2}\over 2}nk_{B}\frac{T_{e}}{T_{F}}$ ($=62.8\,T_{e}\,\rm Jm^{-3}K^{-1}$ for the density of gold). Since $C_{l}\gg C_{e}$, the electrons respond much more quickly than the ions to the e-ph coupling. Neglecting the heat conductivity, the typical timescale is given by $\tau_{e-ph}\sim C_{e}/G\approx 1\,\rm ps$. At the end of this phase, the electron and lattice temperature are equilibrated. Finally, on even longer timescales ($\approx 1\,\rm ns$) any excess temperature is evacuated into the external environment.

The various timescales, from the ultrashort laser pulse to the lattice relaxation are represented pictorially in Fig. 4.

In summary, there exists an early stage of the laser-induced electron dynamics, lasting around 50 fs, which is essentially collisionless. This early stage can be described in the phase space by a Vlasov equation, or its quantum counterpart, the Wigner equation, as we will show in the remainder of this article. For longer times, collisions (either e-e or e-ph) should be taken into account. This is an issue where phase-space-based models have an advantage with respect to wave-function-based ones. Indeed, there is a large literature on dissipative kinetic equations for classical plasmas (Boltzmann, Fokker-Planck, Lennard-Balescu,…), which can inspire useful extensions to the quantum regime. In contrast, wave-function-based methods, being essentially Hamiltonian, face more difficulty in incorporating dissipative effects.

3 Wave-function-based methods: Hartree, Hartree-Fock, DFT

The dynamics of $N$ interacting bodies is a fascinating and challenging problem in physics. Classically, solving the exact $N$-body problem involves integrating the equations of motion of the $N$ particles using some time-stepping technique. The computational cost of such an operation grows like $N^{2}$. Sometimes, using special tricks as in hierarchical tree codes, this can be brought down to $N\log N$. For systems where the number of particles is very large, like plasmas, this is still a formidable problem, but it can be attacked using modern molecular dynamics techniques.

The quantum $N$-body problem is even worse, because the complexity of the Hilbert space of a system made of $N$ particles grows exponentially with $N$. Indeed, a quantum system is described by the $N$ body wave function

[TABLE]

where $\bm{r}_{i}$ is the position of the $i$-th particle. This evolves according to the Schrödinger equation

[TABLE]

where $\mathcal{V}(\bm{r})=\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{|\bm{r}|}$ is the Coulomb potential. The wave function lives in the $3N$-dimensional configuration space. If we use 10 points for each direction, then we need $10^{3N}$ points to code the whole wave function, which is a huge number even for small systems made of a few dozen particles.

It is clear, therefore, that some kind of approximation is mandatory if we want to model the dynamics of nano-objects that contain hundreds or thousands of electrons. In this section, we will briefly review the standard methods that are used in condensed-matter theory to describe the electron dynamics. These fall in two broad categories: (i) Hartree and Hartree-Fock methods and (ii) time-dependent density functional theory (TD-DFT). Both are based on the propagation of some reduced one-body wave functions.

3.1 The mean-field approach: Hartree equations

The approximation leading to the Hartree equations is the same that is used classically to obtain the Vlasov equation from the $N$-body Liouville equation, namely, neglecting two-body (and higher-order) correlations. Mathematically, this consists in factoring the $N$ body wave function into $N$ single-particle wave functions

[TABLE]

This approximation was first considered by Hartree in 1927 Hartree1928 , in the context of atomic physics, to describe the self-consistent effect of the atomic electrons on the Coulomb potential of the nucleus. As we saw in the preceding section, this mean-field approximation is valid only when the corresponding coupling parameter $g_{q}$ is small.

Using Eqs. (9) and (10), one obtains that the single-particle wave functions obey the following Hartree equations:

[TABLE]

where the quantum states of each particle are labelled by the wave functions

$\left\{\Psi_{\alpha}\left(\bm{r},t\right),~{}\alpha=1,\cdots,N\right\}$. The Hartree potential $V_{H}$ is the self-consistent potential created by the ensemble of all electrons. The term $V_{sic}$ is known as the self-interaction correction (SIC) and takes into account the fact that an electron should not interact with itself. The SIC is often neglected for classical plasmas, because its contribution goes like $1/N$, but in small nano-objects it has sometimes to be taken into account. For instance, the SIC correction may be important to enforce the correct asymptotic behavior of the Coulomb potential ($\sim 1/r$) at long distances Ullrich2000 . Here, we shall neglect the SIC correction, so that Eq. (11) can be rewritten as a set of Schrödinger-Poisson equations:

[TABLE]

The simplification achieved through the Hartree equations is enormous, because now we only need to solve $N$ equations for a wave function evolving in three-dimensional (3D) space. Using 10 points per dimension, we only need $10^{3}N$ points to encode the wave functions, instead of $10^{3N}$ for the full $N$-body problem.

3.2 The exchange interaction: Hartree-Fock equations

However, the Hartree equations suffer from two problems. Firstly, they violate the Pauli exclusion principle, which stipulates that two fermions cannot be in the same quantum state. Consequently, the total wave function should be antisymmetric with respect to the exchange of two particles, for instance exchanging particles 1 and 2:

[TABLE]

This property creates some special (purely quantum) correlation between the fermions, which is known as the exchange interaction. In 1930, Slater and Fock Fock1930 proposed an exact method to describe the exchange interactions. They introduced the Slater determinant Slater1929 to write the $N$-body wave function as follows:

[TABLE]

In this case, the total wave function is no longer a product of one-body wave functions, but it satisfies the antisymmetry property required by the Pauli principle. Using Eq. (17) into Eq. (8), one obtains the Hartree-Fock (HF) equations:

[TABLE]

The above HF equations differ from the Hartree equations by the last term in Eq. (18), which encodes the effect of the exchange interactions. Note that the exchange term cannot be represented as a Poisson equation like the Hartree term. It is inherently nonlocal, in the sense that it couples all wave functions together, whereas in the Hartree case the wave functions are only coupled through the partial densities $|\Psi_{\alpha}|^{2}$ appearing in Poisson’s equation (15). This fact renders the HF equations considerably more difficult to solve numerically.

The second limitation, inherent to both the Hartree and the HF methods, is that two-body correlations are not taken into account. This drawback has been overcome since the advent of density-functional theory (DFT), first for the ground state in the 1960s, and then for time-dependent problems in the 1980s.

3.3 Time-dependent density-functional theory (TD-DFT)

A big step forward in the treatment of many-body problems was achieved thanks to the results of Hohenberg and Kohn Hohenberg1964 and Kohn and Sham Kohn1965 , which can be summarize in two theorems. The first Hohenberg-Kohn theorem states that the ground-state properties of a many-electron system are exactly and uniquely determined by the electron density. This is a huge progress, as the extremely complex $N$-body problem (wave function depending on $3N$ spatial coordinates) has been reduced to finding the correct electron ground-state density (3 spatial coordinates). The second Hohenberg-Kohn theorem states that the ground-state density can be found by minimizing a certain energy functional:

[TABLE]

where $T[n]$ is the kinetic energy functional, $E_{ext}[n]$ is the external energy, and $E_{int}[n]$ is the Coulomb interaction energy. The interaction energy can be decomposed into the Hartree energy $E_{H}[n]$, the exchange energy $E_{X}[n]$, and the correlation energy $E_{C}[n]$ – the latter containing all higher-order correlations beyond the mean field (Hartree) and the Pauli principle (exchange). The above expression (19) is in principle exact, in the sense that we know that such functionals exist. However, some of the terms are not known and thus need to be approximated.

As an example, the old Thomas-Fermi theory Thomas1926 of the atomic electron gas can be viewed as an early precursor of modern DFT. In the Thomas-Fermi theory, exchange and correlations are neglected and the external potential is that of the atomic nucleus, $V_{ext}(\bm{r})=-Ze^{2}/4\pi\epsilon_{0}|\bm{r}|$. Then, we have for the external energy functional: $E_{ext}[n]=\int n(\bm{r})V_{ext}(\bm{r})d{r}$, and for the Hartree energy:

[TABLE]

Finally, Thomas and Fermi chose a semiclassical functional for the kinetic energy of an ideal fermion gas at zero temperature:

[TABLE]

Minimizing the total energy functional $E[n]$ with the constraint $\int n\,d\bm{r}=N$, yields the usual Thomas-Fermi equation for the ground-state density.

Apart from neglecting exchange and correlation effects, the main drawback of the Thomas-Fermi approach is its poor approximation of the kinetic energy functional. To improve on this approximation, Kohn and Sham Kohn1965 suggested that one uses the kinetic energy $T_{s}[n]$ of a fictitious non-interacting electron gas, with the same density $n$ as the original interacting one, which evolves in an effective potential $V_{eff}[n]=-eV_{H}[n]+V_{ext}+V_{X}[n]+V_{C}[n]$. Each term of the effective potential is obtained as a functional derivative of the corresponding energy functional: $V_{k}[n]=\delta E_{k}[n]/\delta n$. The Kohn-Sham (KS) equations are then:

[TABLE]

where the $\phi_{j}$ are the single-particle wave functions of the fictitious non-interacting system. The electron density is then obtained as: $n(\bm{r})=\sum_{j}p_{j}|\phi_{j}(\bm{r})|^{2}$, where the $p_{j}$ are occupation probabilities. The KS equations are potentially exact, provided one knows all the terms of the effective potential $V_{eff}[n]$. Unfortunately, this is not true for the exchange and correlation potentials, for which we need to find approximate expressions.

This caveat notwithstanding, the KS accomplish the truly impressive task of representing in an exact way the full $N$-body problem through a set of one-body wave functions. In addition, unlike the HF equations, the KS equations are local and only coupled to each other via the electron density $n$. Their computational complexity is thus similar to that of the Hartree equations. Thus, the KS equations can be seen as a way to make the Hartree approach exact by including the appropriate exchange and correlation functionals. As the latter functionals are not known exactly, all the art of DFT is to find the best way to approximate them. The simplest choice is the so-called local density approximation (LDA) Kohn1965 , whereby the exchange and the correlation functionals depend locally on the electron density. For instance, the LDA approximation for the exchange potential is:

[TABLE]

Many other sophisticated approximations have been developed over the years (such as the generalized gradient approximation, GGA), making DFT methods a cornerstone of computational materials science and theoretical chemistry Jones2015 .

A time-dependent version of DFT (TD-DFT) was developed by Runge and Gross Runge1984 in the 1980s. The Runge-Gross theorem stipulates that, for the same initial $N$-particle state, two external potentials differing only by a time-dependent function $c(t)$ cannot give rise to the same density $n(\bm{r},t)$. Using this theorem, one can construct the time-dependent KS equations:

[TABLE]

Compared to their static counterpart, the time-dependent KS have the additional difficulty of requiring a time-dependent approximation for the functionals. The simplest choice is to use the same functional as in the ground state, but allowing a time dependence in the density. This is known as the adiabatic local-density approximation (ALDA). Note also that, by setting $V_{X}=V_{C}=0$ in the KS equations, one recovers exactly the time-dependent Hartree equations (11).

Finally, we stress that, since the time-dependent KS equations (24) have the same mathematical form as a set of nonlinear Schrödinger equations, they can be used to construct a phase-space formalism by taking the Wigner transform of the KS wave function.

4 Quantum mechanics in the phase space

In this section, we summarize the main properties of Wigner’s phase-space formulation of quantum mechanics, which was first introduced by Eugene Wigner in 1932 to study quantum corrections to classical statistical mechanics Wigner1932 . The goal was to link the wave function that appears in the Schrödinger equation to a pseudo-probability distribution defined in the classical phase space. This pseudo-probability distribution changes in time according to an evolution equation (Wigner equation) that is somewhat similar to the classical Liouville equation.

A mathematically rigorous treatment of the Wigner formulation is based on the Weyl transformation WeylHermann1928 ; Cohen2013 , which is a general method to transform operators defined in the Hilbert space into phase-space functions. We have developed these arguments at some length in a previous publication Hurst2017 , to which we address the reader for further details. Here, we will summarize the main properties of the Wigner approach for particles in electric and magnetic fields, with and without spin.

4.1 The Wigner equation for a scalar electric potential

The density matrix of a quantum mixture of $N$ particles is given by:

[TABLE]

where $p_{\alpha}$ is the probability for one particle to be in the state $\Psi_{\alpha}$. The wave functions are supposed to obey a set of Schrödinger-like equations such as the time-dependent Hartree equations (15) or Kohn-Sham equations (24) discussed in Sec. 3.

The Wigner function of the system is defined as

[TABLE]

The Wigner function evolves in time according to the following Wigner equation

[TABLE]

where the subscripts $\bm{\pm}$ denote the shifted positions $\bm{r}_{+}=\bm{r}\pm\bm{\lambda}/2$ and $V(\bm{r},t)$ indicates a generic scalar potential, which can be the Hartree potential or the effective potential of TD-DFT.

The Wigner equation (27) is completely equivalent to the Schrödinger-like equations from which it is derived. Therefore, it provides a useful way to cast TD-DFT in a phase-space formalism without any loss of generality. In addition, the Wigner function can be used to compute all macroscopic quantities in the same way as a classical probability distribution. For instance, the particle and current densities:

[TABLE]

Finally, the Wigner approach is particularly useful to obtain the semiclassical limit. By developing the integral term in (27) up to order $O(\hbar^{2})$ we obtain

[TABLE]

The Vlasov equation is thus recovered in the formal semiclassical limit $\hbar\to 0$.

4.2 The Wigner equation in a magnetic field

So far, we did not include any magnetic effects. This approximation may be justified in some cases, for instance if we are only interested in plasmonic excitations. However, magnetic interactions cannot be escaped if we want to include the spin degrees of freedom in our treatment. The introduction of a magnetic field in the Wigner formalism is not trivial. In the presence of magnetic fields, one should use the kinetic momentum operator $\widehat{\bm{\pi}}=\widehat{\bm{p}}-q\widehat{\bm{A}}$, instead of $\widehat{\bm{p}}$ (with $q=-e$ for an electron), where $\widehat{\bm{A}}$ is the vector potential operator. However, it can be proven that simply substituting $\bm{p}$ with $\bm{\pi}\equiv m\bm{v}$ in the definition (26) will not work, as the resulting Wigner function is not gauge invariant.

A gauge-independent definition of the Wigner function was first introduced by Stratonovich Stratonovich1957 :

[TABLE]

where the momentum $\bm{p}$ was replaced by $m\bm{v}-e\int_{-1/2}^{1/2}d\tau\bm{A}\left(\bm{r}+\tau\bm{\lambda}\right)$.

After some rather convoluted algebra, one finally obtains the gauge-invariant Wigner equation for a spinless particle interacting with an electromagnetic field:

[TABLE]

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

[TABLE]

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

[TABLE]

This form of the Wigner equation was first proposed by Serimaa et al. Serimaa1986 , where the authors also discuss the case where the electromagnetic fields are quantized. In the classical limit, it is straightforward to see that: $\widetilde{\bm{E}}=\bm{E}$, $\widetilde{\bm{B}}=\bm{B}$, and $\bm{\Delta\widetilde{\pi}}=0$, so that the Wigner equations becomes

[TABLE]

thus recovering the classical Vlasov equation with the Lorentz force.

5 Spin and relativistic effects

5.1 Semi-relativistic Pauli equation

In the previous chapter, we omitted all mentions to the spin of the electrons. The spin is an intrinsic property of any elementary particle as much as the mass or the charge. It was first discovered in 1922 thanks to the experiments of Stern and Gerlach Gerlach1922 and was interpreted as an internal angular momentum of the electron. From a theoretical point of view, the spin appears naturally in the Dirac equation, which is the relativistic extension of the Schrödinger equation for spin $1/2$ particles. For an electron interacting with an external electromagnetic field, the Dirac equation reads as:

[TABLE]

where $V\left(\bm{r},t\right)$ and $\bm{A}\left(\bm{r},t\right)$ are, respectively, the scalar and vector potentials. The operators $\bm{\alpha}$ and $\beta$ are $4\times 4$ matrices

[TABLE]

where $\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector of the $2\times 2$ Pauli matrices

[TABLE]

and $\sigma_{0}$ is the $2\times 2$ identity matrix. Therefore the wave functions $\Psi^{\textrm{D}}\left(\bm{r},t\right)$ that obey to the Dirac equation (36) are four-component objects called bispinors.

The Dirac equation contains much more information than just the spin, as it deals at the same time with the dynamics of particles (electrons) and antiparticles (positrons), described respectively by the two upper (lower) components of the Dirac wave function. In our case, since we are only interested in the low-energy phenomena occurring in condensed-matter physics, we would like to discard all effects related to electron-positron pair formation.

For this purpose, it is helpful to use the unitary Foldy-Wouthuysen transformation Foldy1950 , which enables one to separate the electron and the positron dynamics in the Dirac equation. This transformation is exact in the field-free case, and leads to a semi-relativistic expansion in $1/c$ (where $c$ is the speed of light) for a particle interacting with an electromagnetic field Strange1998 .

At second order in $1/c$, the Dirac Hamiltonian transforms into the following semi-relativistic Pauli Hamiltonian Dixit2013 :

[TABLE]

where the electromagnetic fields are defined as usual as: $\bm{E}=-\nabla V-\partial_{t}\bm{A}$ and $\bm{B}=\nabla\times\bm{A}$, and $\hat{\bm{p}}=-i\hbar\nabla$. The first term on the right-hand side is the rest-mass energy of the electron; the next two terms are the standard nonrelativistic Hamiltonian in the presence of an electromagnetic field; the fourth term is the Pauli spin term (Zeeman effect); the $(\hat{\bm{p}}+e\bm{A})^{4}$ term is the first relativistic correction to the electron mass (expansion of the Lorentz factor $\gamma$ to second order); the $\nabla\cdot\bm{E}$ term is the Darwin term; and the last two terms represent the spin-orbit coupling (SOC). The wave function $\Psi={{}^{t}(}\Psi^{\uparrow},\Psi^{\downarrow})$ is a spinor, the upper and lower components describing respectively the spin-up and the spin-down electrons, and it obeys the spinorial Schrödinger equation: $i\hbar\partial_{t}\Psi=\hat{H}\,\Psi$. Higher order extensions of the Foldy-Wouthuysen expansion can be found in Hinschberger2012 .

The Zeeman and spin-orbit effects are of paramount importance to describe the magnetic properties of ferromagnetic nano-objects. For instance, there are strong indications that the SOC plays a crucial role in many experiments where the magnetisation is excited with optical pulses Krieger2015 . Thus, in the following, we will retain the Zeeman and SOC effects, but neglect the relativistic correction to the electron mass (which would generate awkward fourth-order gradients in the Schrödinger equation) and the Darwin term. The latter is a manifestation of the so-called Zitterbewegung, i.e., a quivering motion of the electron around its mean path Dixit2013 , which is due to the interference between the positive and negative energy states in the Dirac equation. This term could be reintroduced in our treatment without much difficulty.

5.2 The Wigner equation with spin

For spinless particles, the Wigner function is a scalar function related to the density matrix by Eq. (31). In the case of spin-${1/2}$ particles, both the Wigner function $\mathcal{F}$ and the density matrix $\rho$ are $2\times 2$ matrices:

[TABLE]

where $\uparrow,\downarrow$ denote respectively the spin-up and spin-down components.

It is convenient to project the matrix $\mathcal{F}$ onto the Pauli basis set Barletti2003 ; Morandi2011

[TABLE]

where

[TABLE]

Here ${\bm{\sigma}}=(\sigma_{x},\sigma_{y},\sigma_{z})$ are the Pauli matrices and $\bm{f}=(f_{x},f_{y},f_{z})$.

With this definition, the particle density $n$ and the spin polarization $\bm{S}$ of the electron gas are easily expressed as moments of the pseudo-distribution functions $f_{0}$ and $\bm{f}$:

[TABLE]

In this representation, the Wigner functions have a clear physical interpretation: $f_{0}$ is related to the electron density (in the phase space), whereas $f_{i}$ ($i=x,y,z$) is related to the spin polarization density in the direction $i$. In other words, $f_{0}$ represents the probability to find an electron at one point of the phase space at a given time, while $f_{i}$ represents the probability that the spin polarization of such electron is directed along in the $i$-th direction. One can prove the following interesting bound:

[TABLE]

Equation (45) is a direct consequence of the property of the density matrix: $Tr\left(\rho^{2}\right)\leq 1$. The equality sign holds true for a pure state or for a system where all the spins are aligned along the same direction (fully spin-polarized electron gas).

The full (quantum) evolution equation for the matrix Wigner function $\mathcal{F}$ was derived in Hurst2014 and is mathematically very complicated. A better physical insight can be gained from its semiclassical limit, i.e., keeping only terms of order $O(\hbar)$. One obtains:

[TABLE]

where the factor $\hbar$ is hidden in the definition of the Bohr magneton $\mu_{B}=e\hbar/(2m)$, which precedes all quantum corrections in Eqs. (46)-(47). Such quantum corrections couple the orbital motion with the spin terms, through the Zeeman effect or the spin-orbit interactions. There are no quantum corrections to the orbital electron dynamics, because they would only appear at second order in $\hbar$. Therefore, the orbital motion of the electrons is purely classical and determined by the Lorentz force. In contrast, the spin degree of freedom is treated as a quantum variable.

In summary, the Eqs. (46) and (47) represent the semi-relativistic (order $c^{-2}$) and semi-classical (order $\hbar$) form of the kinetic equations for particles with spin $1/2$, which we shall term the “spin-Vlasov equations” in the following.

5.3 Self-consistency

Equations (46)-(47) can be used, in a mean-field approach, to described the self-consistent spin dynamics of an ensemble of interacting charged particles. In this case, the electric and the magnetic field are solutions of the Maxwell equations:

[TABLE]

where $\bm{P}$ and $\bm{M}$ are, respectively, a polarization and a magnetization vector. These terms arise because, in the Foldy-Wouthuysen procedure, both the Dirac Hamiltonian and the wave functions are expanded in powers of $1/c$. The expansion of the wave function to order $c^{-2}$, entails that the charge and current densities are modified, as was shown using a variational method in Dixit2013 ; Manfredi2013 . The resulting source terms are the following:

[TABLE]

The spin-Vlasov equations (46)-(47), coupled to the Maxwell equations (52), constitute a self-consistent set of equations to describe the spin dynamics of an interacting electron gas in the mean-field approximation. Such mean-field approach can also be extended, in the spirit of TD-DFT (see Sec. 3.3), to include exchange and correlation effects by adding suitable potentials and fields that are functionals of the electron density Manfredi2010 .

5.4 Scalar representation in an extended phase space

There exists another way to include the spin in a phase-space Wigner formalism, which is mathematically different but physically identical to the approach described in Sec. 5.2. The basic idea Zamanian2010nj ; Marklund2010 is to define an extended phase-space distribution $g$ that depends not only on the position $\bm{r}$ and the velocity $\bm{v}$, but also on a unitary vector $\bm{s}$ (defined by two angles on the unit-radius sphere) which represents the spin direction:

[TABLE]

where $\mathcal{F}_{\beta\alpha}$ is the matrix Wigner function as in Eq. (40), $\delta_{\alpha,\beta}$ is the Kronecker delta function, and $\bm{\sigma}_{\alpha,\beta}$ represents the $(\alpha,\beta)$ element of one of the Pauli matrices $\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ defined in Eq. (38).

The distribution $g(\bm{r},\bm{v},\bm{s},t)$ is a scalar function that evolves in the extended phase space, which is therefore eight-dimensional (8D): three positions, three velocities, and two angles for the spin. This is in contrast with the approach of Sec. 5.2, where the phase space is the usual 6D one ($\bm{r},\bm{v}$), but the distribution function is a $2\times 2$ matrix. The correspondence relations between our distribution functions $f_{0}(\bm{r},\bm{v},t)$ and $f_{i}(\bm{r},\bm{v},t)$ and the scalar distribution $g(\bm{r},\bm{v},\bm{s},t)$ used by Zamanian et al. Zamanian2010nj can be written as:

[TABLE]

The semi-relativistic theory in the extended phase space, equivalent to the matrix approach of Eqs. (46) and (47), was derived in Asenjo2012 and a fully relativistic extension of this scalar theory was established later in Ekman2017 .

The semi-relativistic scalar kinetic equations reads as follows:

[TABLE]

where the summation over repeated superscripts $i=x,y,z$ is understood and we have defined the quantities

[TABLE]

The deformed magnetic field $\overline{\bm{B}}$ is the same expression that appears in the Eqs. (46)-(47) and corresponds to the first relativistic correction to the Thomas precession Thomas1926 ; Dragan2013 .

The deformed velocity $\bm{w}$ is related to the spin-orbit correction of the velocity operator. Indeed, in the Heisenberg picture, the velocity operator $\bm{\hat{V}}$ is determined by the evolution equation of the position operator $\bm{\hat{r}}$

[TABLE]

where we used the Pauli Hamiltonian (5.1). The related phase-space function is

[TABLE]

which is also the velocity appearing in the Eqs. (46)-(47).

From a computational point of view, it is interesting to note that the matrix and the scalar formalism of the spin-Vlasov equations lend themselves better to two different families of numerical methods. The scalar formalism is better adapted to particle-in-cell (PIC) methods, because the two degrees of freedom related to the spin only act as two additional labels attached to each particles. The computational cost is therefore only marginally greater than that required for spinless particle. In contrast, the scalar approach is hardly applicable to grid-based Vlasov codes, as it would require to mesh the extended 8D phase space. Grid-based codes are better suited for the matrix formalism, because the phase space is at most 6D and the only drawback compared to spinless particles is that one has to advance in time four, instead of one, distribution functions.

6 Beyond the mean field: collisions and relaxation to equilibrium

The mean-field approach described so far is accurate to treat the electron dynamics on very short timescales ($<100\rm\,fs$), see Fig. 4. On longer timescales (0.1–1ps), the laser energy is redistributed among the electrons via electron-electron (e-e) collisions, and finally delivered to the ion lattice via electron-phonon (e-ph) collisions.

One of the advantages of phase-space-based methods is that effects going beyond the mean-field approximation (e.g., two-body collisions) can be incorporated with relative ease in the governing equations. This is not the case for wave-function-based methods (Hartree-Fock, DFT), which have an essentially Hamiltonian nature reflected in the unitary propagation of the wave function.

Collisions are dynamical correlations that cannot be included in adiabatic correlation functionals in TD-DFT. There have been some attempts to include potentially dissipative effects in TD-DFT. The most accomplished of such attempts is the so-called time-dependent current-density-functional theory (TD-CDFT) Vignale97 , which uses the electron current density $\bm{j}(\bm{r},t)$ as the basic building block, instead of the density $n(\bm{r},t)$. However, the equations of TD-CDFT are mathematically very complicated and not of easy implementation in practical situations.

For phase-space methods, the construction of dissipative terms can rely on the experience acquired in plasma physics. Generally speaking, the relevant kinetic equation (Vlasov, Wigner) may be augmented by a collision term of the form:

[TABLE]

We also note that it is conceptually harder to include collisions in fully quantum models. A significant constraint is that non-unitary corrections to the Wigner equation should be written in Lindblad form Lindblad1976 , which guarantees that the evolved Wigner function always corresponds to a positive-definite density matrix. This constraint is not always satisfied for the models described below.

In the following, we will discuss three different collision integrals for e-e and e-ph collisions, based on a relaxation term, a Boltzmann-like approach, and a Fokker-Planck approach.

Relaxation methods.

To model e-e collisions, a relaxation term can be added to the right-hand side of the Vlasov or Wigner equation Manfredi2005film :

[TABLE]

where $\nu_{ee}$ is the average e-e collision rate and $f_{eq}(\bm{r},\bm{v})$ is a Fermi-Dirac distribution. The idea behind this model is that the electron distribution will eventually relax, on a time scale of the order $\nu_{ee}^{-1}$, towards a Fermi-Dirac equilibrium $f_{eq}$ with total energy equal to that of the electron distribution $f(\bm{r},\bm{v},t=0^{+})$ after the initial excitation. For electrons near the Fermi surface, the e-e collision rate can be written as Pines1995 :

[TABLE]

where $T_{e}$ is the instantaneous electron temperature. The proportionality constant can be estimated from experimental considerations or extracted from first-principles simulations Domps1998 .

Equation (64) conserves the total number of particles, but not the total momentum or energy. It can only ensure that the distribution function relaxes to the correct equilibrium distribution, but its relevance for the nonequilibrium transient is questionable. In order to conserve energy and momentum, the function $f_{eq}[n,\bm{u},T]$ should be a local equilibrium with density $n(\bm{r},t)$, mean velocity $\bm{u}(\bm{r},t)$, and temperature $T(\bm{r},t)$ computed instantaneously from the electron distribution $f(\bm{r},\bm{v},t)$. Note that the task of computing the local temperature is not trivial for a Fermi-Dirac distribution.

Boltzmann-like methods.

A Boltzmann-like electron-electron collision integral that respects Pauli’s exclusion principle was devised long ago by Ühling and Uhlenbeck UU :

[TABLE]

where $v_{12}$ is the relative velocity of the colliding particles 1 and 2, $\sigma(\Omega)$ is the differential cross section depending on the scattering angle $\Omega$, and indices 3 and 4 label the outgoing momenta, $f_{i}=f(\vec{r},\vec{p}_{i},t)$ and $\overline{f}_{i}=1-f_{i}/2$. This collision term is similar to the well-known classical Boltzmann collision term but for Pauli blocking factors $\overline{f}_{i}\overline{f}_{j}$. As seen in Sec. 2, Pauli blocking plays a major role for electronic systems. At zero electron temperature, all collisions are Pauli-blocked and the collisional mean free path of the electrons becomes infinite. However, after a strong excitation, the distribution function is no longer a Fermi-Dirac one and collisions become possible. The effect of the above e-e collision term on the semiclassical Vlasov dynamics in metal clusters was investigated numerically in Domps1998 .

Fokker-Planck methods.

By coupling to the ionic lattice vibrations (phonons), the electrons progressively relax to a thermal distribution with a temperature equal to that of the lattice $T_{l}$. This behavior can be described phenomenologically using the two-temperature model illustrated in Sec. 2. The typical relaxation time for metal nano-objects is $\tau_{e-ph}\approx 1\,\rm ps$. In addition, the lattice acts as an external environment for the electrons, leading to the loss of quantum coherence over a timescale $\tau_{dec}$ (decoherence time). The relaxation and decoherence times correspond, respectively, to the decay of diagonal and nondiagonal terms in the density matrix describing the electron population.

Such effects can be modeled, in the Wigner representation, by a classical Fokker-Planck (FP) term Jasiak2010 :

[TABLE]

where $D_{r,v}$ are diffusion coefficients in real space and velocity space respectively, $\gamma$ is the e-ph relaxation rate, and $G[f]$ is a functional that depends on the statistics and on the dimensionality of the system.

For classical particles obeying an exclusion principle, the corresponding FP equation was derived by Kaniadakis et al. Kania . For instance, $G[f]=f$ for particles obeying Maxwell-Boltzmann statistics and $G[f]=f(1-f)$ for fermions in 3D. For fermions in 1D, which is the case of the electrons in thin films that we shall consider later, we assume that the electron distribution always remains a Fermi-Dirac one in the transverse directions $(v_{y},v_{z})$. Then, integrating the collision term (66) along these two directions, one obtains the expression $G[f]=f_{0}[1-\exp(-f/f_{0})]$, where $f_{0}=\frac{3}{4}\frac{n_{0}}{v_{F}}\frac{T_{l}}{T_{F}}$, $T_{l}$ is the lattice temperature, and $f(x,v_{x},t)$ is now the 1D distribution function. The preceding expression must be changed slightly to account for the fact that the Wigner distribution can be negative: $G[f]=f_{0}[1-\exp(-|f|/f_{0})]{\rm sgn}(f)$.

It can be proven that $(\partial{f}/\partial{t})_{\rm e-ph}=0$ when the electron distribution is given by a spatially homogeneous 1D Fermi-Dirac distribution ($\mu$ is the chemical potential):

[TABLE]

provided $D_{v}$ and $\gamma$ satisfy the relation: $D_{v}=\gamma k_{B}T_{l}/m$. Thus, the FP term ensures that the electron distribution relaxes to a classical Fermi-Dirac distribution with a temperature equal to that of the lattice. The latter constitutes a perfect reservoir with infinite heat capacity, so that $T_{l}$ remains constant.

As seen in Sec. 2, the electron relaxation rate can be written as $\gamma\equiv\tau_{e-ph}^{-1}\propto G/C_{e}$, where $G$ is the electron-phonon coupling constant appearing in the two-temperature equations (6)-(7) and $C_{e}$ is electron heat capacity. For an ideal Fermi gas, the heat capacity depends on the temperature as $C_{e}(T_{e})=\pi^{2}n_{0}k_{B}(T_{e}/2T_{F})$. In a dynamical simulation, the electron temperature is a time-dependent quantity that can be computed self-consistently from the Wigner distribution function $f$.

We note that a Vlasov or Wigner equation endowed with the e-ph collision term (66) can be viewed as the microscopic counterpart of the electron equation (6) in the phenomenological two-temperature model. This can be proven easily in the case of a 1D Maxwell-Boltzmann distribution, using the diffusion operator in the phase space:

[TABLE]

The above collision term is in the Lindblad form whenever $D_{x}\,D_{v}\geq\gamma^{2}\hbar^{2}/16m^{2}$ Lindblad1976 ; Zurek2003 . Taking the first and second velocity moments of the Vlasov or Wigner equation (there are no quantum corrections to these orders), neglecting all electromagnetic fields, and defining the kinetic electron temperature in the usual way: $k_{B}T_{e}=m\langle(v-\langle v\rangle)^{2}\rangle$, where $\langle A\rangle=\int\int fAdxdv$ denotes the phase-space average, one obtains:

[TABLE]

Then, taking $D_{v}=\gamma k_{B}T_{l}/m$ and $\gamma=\tau_{e-ph}^{-1}=G/(2C_{e})$, we get

[TABLE]

which is the electron temperature equation (6) in the TTM with $\kappa=P=0$.

Finally, we would like to mention the recent work of Daligault Daligault2016 , who derived a quantum form of the Landau/Fokker-Planck collision operator.

7 Application I: Spin corrections to longitudinal plasma waves

The linear response of a homogeneous spinless electron gas has been studied theoretically for a long time, starting with the pioneering works of Vlasov vlasov1938 and Landau Landau1946 . The main linear mode is an oscillation at the plasma frequency $\omega_{p}=\sqrt{e^{2}n_{0}/m\epsilon_{0}}$ (“plasmon”) , corresponding to the collective motion of the electrons immersed in a neutralizing background of positive ions, which are assumed to be fixed because of the large mass difference between the two species.

In spherical nano-objects, a typical example of plasmonic oscillation is the localized surface plasmon (LSP), which originates from a displacement of the electron gas that creates a net charge imbalance (see Fig. 5). The resulting Coulomb force pulls the electrons back inside the system but, due to their inertia, they will travel further away, thus recreating a new Coulomb force in the opposite direction. After a few cycles, the plasma oscillations are usually damped away through Landau damping, which results from the mixing of single-electron oscillations at slightly different frequencies. LSPs have been observed in all sorts of nano-objects Scholl2012 ; Luo2013 ; Tame2013 ; Raza2013 ; Tanjia2018 and are one of the basic features of the fast-developing field of plasmonics manfredi2018 .

Here, our aim is to show a possible signature of the electron spin polarization on the plasmonic oscillations.

7.1 General linear response theory

In this section, we will derive the dielectric function of the spin-Vlasov equations coupled to the relevant Maxwell equations and then find the associated dispersion relation. In general, this is a very complicated task. Here, we simplify the problem by assuming that the equilibrium ground-state distribution is homogeneous and isotropic and that the electron spins are polarized along the $z$ direction only (collinear approximation). Therefore, the diagonal components of the Wigner matrix for the spin-up and spin-down equilibria can be written as (only one arrow is used, to simplify the notation):

[TABLE]

where $\mathcal{G}$ can be a Maxwell-Boltzmann or a Fermi-Dirac distribution and $\bm{B}=B\bm{e}_{z}$ is a uniform magnetic field parallel to the $z$ axis. From Eq. (42), one has: $f_{0}^{(0)}=f^{\uparrow(0)}+f^{\downarrow(0)}$ and $f_{z}^{(0)}=f^{\uparrow(0)}-f^{\downarrow(0)}$.

The uniform magnetic field in Eqs. (71) creates a difference between the spin-up and spin-down distributions; hence, it creates a net spin polarization at equilibrium. However, as it is homogeneous in space, it has no influence on the ensuing dynamics. The equilibrium distribution functions also verify the following properties: $\int f_{0}^{(0)}d\bm{v}=n_{0}$ and $\int f_{z}^{(0)}d\bm{v}=m_{0}$, where $n_{0}$ and $m_{0}$ are, respectively, the electron density and magnetization at equilibrium.

In order to obtain a tractable system of equations, we need to make several assumptions. First, we neglect the effect of the magnetic field on the orbital motion: thus, there is no Lorentz force $\bm{v}\times\bm{B}$ in the Vlasov equations. Second, we consider purely longitudinal plasma waves, with wave vector $\bm{k}=k\bm{e}_{x}$. Therefore, all physical quantities only depend on the spatial coordinate $x$, and the corresponding phase space $(x,v)$ is 2D, where $v\equiv v_{x}$ stands for the $x$-component of the velocity. As was mentioned before, the electron gas is polarized along the $z$ direction. With these assumptions, the spin-Vlasov equations (46)-(47), together with Maxwell’s equations (52), become:

[TABLE]

where $V_{XC}[n]$ and $B_{XC}[n]$ are the exchange and correlation functionals Maurat2009 , and $B_{z}$ is the self-consistent magnetic field.

In order to perform the linear wave analysis, we expand all quantities around the equilibrium configuration. For the distribution functions, we have

[TABLE]

and for the exchange and correlation potentials

[TABLE]

and similarly for $B_{XC}$. Inserting Eqs. (75)-(76) into Eqs. (72)-(74) and neglecting second-order terms, we obtain the linearized spin-Vlasov equations. We then follow the standard Landau procedure of Fourier transforming in space and Laplace transforming in time. This leads to the dielectric function of the system

[TABLE]

where $\omega$ is the frequency of the perturbation and $k$ the wave number, and the integrals

[TABLE]

only depend on the ground-state properties. The zeros of the dielectric function determine the dispersion relation of the system: $\omega=\omega(k)$.

Finally, one should specify the form of the exchange and correlation functionals. The minimal requirement is to neglect correlations and take an exchange functional that is spin dependent and local in space and time Gunnarsson1976 :

[TABLE]

The above functionals are the exact solutions of the Hartree-Fock equations in the case of homogeneous electronic densities. In our case, we are close to this situation, since we study perturbations around homogeneous ground states.

7.2 Fermi-Dirac ground state

In the case of electrons in metals, the spin-up and spin-down ground states obey a Fermi-Dirac distribution in one dimension, obtained after integrating the 3D Fermi-Dirac function over the transversal velocities: $f^{1d}(v_{x})=\int\int f^{3d}(\bm{v})dv_{y}dv_{z}$. Renaming $v\equiv v_{x}$, we obtain the ground state distributions:

[TABLE]

where a magnetic field $\bm{B}=B\bm{e}_{z}$ is present in order to yield a spin-polarized system. The magnetic field can be either external or internal, i.e. originating from the local magnetic exchange interaction between the localized ions and the itinerant electrons. The chemical potential $\mu$ is set to fix the electron density to $n_{0}$, i.e. $\int\left(f^{\uparrow^{(0)}}+f^{\downarrow^{(0)}}\right)dv=n_{0}$. This integration can be performed exactly only in the zero-temperature limit, leading to the following identity:

[TABLE]

where $v_{F_{\pm}}=\Re\left\{\sqrt{2\left[\mu_{T=0}\pm\mu_{B}B\right]/m}\right\}$, $\mu_{T=0}$ is the chemical potential at zero temperature, and $\Re$ denotes the real part. In the same way, we can express the magnetization at equilibrium as follows:

[TABLE]

In Fig. 6 we plot the polarization at equilibrium:

[TABLE]

as a function of the magnetic field. For a magnetic field larger than $E_{F}/\mu_{B}$, the electron gas is completely polarized. Note that, for metals, one typically has $E_{F}/\mu_{B}\approx 10^{5}\,\rm T$, which is a huge value for a magnetic field generated by external coils. However, internal magnetic fields due to the interatomic exchange interaction can reach such large values, an effect that is at the basis of ferromagnetism.

In the case of a totally unpolarized ($\eta=0$) or fully polarized ($\eta=1$) electron gas, Eq. (83) can be inverted exactly, leading to the following expressions for the chemical potential:

[TABLE]

For intermediate cases ($0<|\eta|<1$), one has to find $\mu_{T=0}$ numerically. In Fig. 6, we plot the chemical potential as a function of the applied magnetic field in the case of a Fermi-Dirac distribution with a temperature equal to $300\,\rm K$ and a density relevant to metals.

In general, the full determination of the dispersion relation has to be done numerically. However in the limit of zero temperature and long wavelengths, i.e. $k\ll|\omega|/v_{F}$, one can develop the integrals $\mathcal{I}_{0}$ and $\mathcal{I}_{z}$, defined in Eq. (78), in a power series of $kv_{F_{\pm}}/\omega$, leading to:

[TABLE]

Using these expressions together with the dielectric function defined in Eq. (77), we get the following dispersion relation:

[TABLE]

Equation (90) constitutes an extension of the celebrated Bohm-Gross dispersion relation in the case of a fully degenerated ($T=0$) spin-polarized electron gas with exchange and magnetostatic interactions. All the dispersive terms depend on the electron density $n_{0}$ and the spin polarization $\eta$. The second term on the right-hand side corresponds to the electron degeneracy pressure. The third and fourth terms represent, respectively, the exchange and the magnetostatic contributions.

In order to assess the relative importance of each term, it is useful to express Eq. (90) in non-dimensional form:

[TABLE]

where $C_{X}=2/(12\pi^{2})^{5/3}\approx 0.017$, $f(\eta)=(1+\eta)^{4/3}+(1-\eta)^{4/3}$, ${\mkern 0.75mu\mathchar 22\relax\mkern-9.75mu\lambda}_{C}=\hbar/(mc)=3.86\times 10^{-4}\rm\,nm$ is the reduced Compton wavelength, and $\lambda_{TF}=v_{F}/\omega_{p}$ is the Thomas-Fermi screening length. The first term in the curly brackets is the dominant one. The next term is small except for very large magnetic fields, of the order of $E_{F}/\mu_{B}\approx 10^{5}\,\rm T$ for metallic densities. The third term is due to the spin-dependent electron-electron exchange interaction. It is a beyond-mean-field contribution, as is evidenced from the fact that it is proportional to the quantum coupling parameter $g_{q}\sim r_{s}/a_{0}$ (see Sec. 2). The fourth term originates from the self-consistent magnetostatic contribution, see Eq. (74), and it is almost always negligible, since for metals $\lambda_{TF}\approx 0.1\,\rm nm$ (Table 1). All these corrections have an opposite sign with respect to the leading term, but they can never render the right-hand side of Eq. (91) negative. Indeed, when $r_{s}/a_{0}\gg 1$, correlation effects (here neglected) start playing a role, so that the resulting frequency is always real.

Interestingly, even in the absence of magnetic interactions, the dispersion relation still depends on the spin polarization. Indeed, in the limit of an unpolarized $(\eta=0)$ or fully polarized $(|\eta|=1)$ electron gas and using Eqs. (86)-(87), the dispersion relation becomes:

[TABLE]

This behavior is specific to the Fermi-Dirac distribution function and is absent in the case of a Maxwell-Boltzmann equilibrium, as we shall see in Sec. 7.3. A spin-polarization dependence of the longitudinal Langmuir modes was predicted in Ref. Hurst2014fluid , where it was shown, using a fluid model, that the pressure tensor is spin-dependent.

To obtain an explicit form of the dispersion relation, we have used the long wavelength approximation ($k\ll|\omega|/v_{F}$), which allowed us to avoid the singularity in the integrals of Eq. (78). As is well known, the correct treatment of the singularity leads to a complex $\omega$ and the appearance of Landau damping. However, this is a delicate mathematical issue when considering a Fermi-Dirac distribution Vladimirov2011 , due to the particular form of the Fermi-Dirac function, which has several singularities in the complex plane. We could not obtain a satisfactory solution to this problem and thus leave it open for the time being. In the next section, we will study a Maxwell-Boltzmann distribution, for which the calculations can be performed all along, at least numerically.

7.3 Maxwell-Boltzmann ground state

The plasma dispersion relation for spinless systems, obtained from the Vlasov-Poisson equations with a Maxwell-Boltzmann equilibrium, was intensively studied in the plasma physics literature since the seminal work of Landau. In that case, the dispersion relation can be fully evaluated even in the strong damping regime.

For a system of electrons endowed with spin, the Maxwell-Boltzmann distributions read as:

[TABLE]

where $v_{T}=\sqrt{2k_{B}T/m}$ is the thermal speed and $\eta=\tanh\left[\mu_{B}B/\left(k_{B}T\right)\right]$ represents the spin polarization of the electrons gas. In this case, the integrals $\mathcal{I}_{0}$ and $\mathcal{I}_{z}$ in Eq. (78) become:

[TABLE]

We use the Landau contour method Lyu2014 to compute the integral in Eq. (96). For this purpose, we introduce the plasma dispersion function

[TABLE]

where $\Gamma$ is a contour in the complex plane following the real axis at infinity and passing under the singularity $z=\varrho$. This function, originally introduced by Fried and Conte Fried1961 , is well defined for $\Im\left(\omega\right)>0$ and can be analytically continued in the lower part of the complex plane, i.e. for $\Im\left(\omega\right)<0$. Then, the integrals $\mathcal{I}_{0}$ and $\mathcal{I}_{z}$ can be expressed in terms of the plasma dispersion function:

[TABLE]

where the following property was used: $Z^{\prime}\left(\varrho\right)=-2\varrho Z\left(\varrho\right)-2.$ Using the above expressions for $\mathcal{I}_{0}$ and $\mathcal{I}_{z}$ and the local density approximation of $V_{X}$ and $B_{X}$ [see Eqs. (79)-(80)], one obtains the following dielectric function for the Maxwell-Boltzmann equilibrium:

[TABLE]

The zeros of the dielectric function provide the dispersion relation $\omega\left(k\right)$ of the system. The latter can be obtained analytically in the long wavelength limit, i.e., for $k\ll|\omega|/v_{T}$ and $\omega\in\Re$. In that case, the derivative of the plasma dispersion function can be developed in a power series:

[TABLE]

leading to the following dispersion relation:

[TABLE]

We recover almost the same expression as with the Fermi-Dirac ground state, see Eq. (90).

To obtain the full dispersion relation, we use ZEAL Kravanja2000 , a mathematical software package for computing the zeros of analytic functions. The dielectric function is parameterized by the electron density $n_{0}$ and the temperature $T$. We first consider an electron gas having a metallic density $n_{0}=5.9\times 10^{28}\,\rm m^{-3}$ (gold). Since the corresponding Fermi temperature is very high ($T_{F}=64\,000\,\rm K$), we need to use a temperature of the same order of magnitude for the Maxwell-Boltzmann approximation to be valid. Here, we take exactly $T=T_{F}$. The zeros of the dielectric function, Eq. (100), are shown in Fig. 7. They are symmetrically located with respect to the vertical axis, which is due to the symmetry properties of the plasma dispersion function. We also note that there exist an infinite number of zeros (black crosses) for which the imaginary part of the frequency is larger than the corresponding real part. These solutions are irrelevant for the dynamics, as they are rapidly damped away. Finally, there are three other zeros that have a different origin. Two of them, represented by blue circles, correspond to the plasmon response at the plasma frequency. For small value of $k$, they have a real part larger than the imaginary part, therefore they are weakly damped and govern the dynamics of the system according to the long-wavelength dispersion relation Eq. (102).

The last zero (red square) has a vanishing real part and disappears when the magnetic exchange interactions are removed. This zero corresponds to a paramagnon mode, i.e., an oscillation in antiphase of the two spin populations (up and down), with no associated charge excitation. For infinite systems, it was shown that paramagnons have a vanishing real frequency Jones1985 , which is in accordance with the result of Fig. 7. Their real frequency can become nonzero for systems with finite size Yin2009 .

Further, exchange interactions do not only lead to paramagnons, but also modify the plasmonic excitations. In Fig. 8, the plasma dispersion relation is depicted for two different physical systems. The first (Fig. 8, top panel) corresponds to the case of electrons in a metal, with the density of gold $n_{0}=5.9\times 10^{28}\,\rm m^{-3}$ and temperature $T=64\,000\,\rm K$. As already mentioned, the temperature had to be set very high for the Maxwell-Boltzmann distribution to be a valid approximation. The second system (bottom panel) corresponds to an artificial case where the electron density is about $n_{0}=10^{22}\,\rm m^{-3}$. For such a system, the Fermi temperature is close to $2\,\rm K$ so that, at room temperature (here, we take $T=300\,\rm K$), it is justified to use a Maxwell-Boltzmann equilibrium.

In the metallic case (Fig. 8, top panel), the resonant frequency is dominated by the electrostatic interactions, while the exchange interactions do not change significantly the dispersion relation. They may decrease slightly the value of the frequency for larger values of $k$, but in this case the plasma oscillations are strongly damped. Moreover, there are no appreciable changes in the dispersion relation depending on the spin polarization of the electron gas. In contrast, in the idealized low-density system (Fig. 8, bottom panel) exchange interactions play a significant role. For small wave numbers, the plasmon frequency decreases instead of increasing and we note a signature of the spin polarization, in particular a shift in $k$ for the minimum of the plasmon frequency.

8 Application II: Nonlinear electron dynamics in thin metal films

Thin metal films have long been the object of intense investigations in many subfields of nanophysics. They constitute key components in many modern technologies, ranging from integrated circuits to sensors. From a more fundamental point of view, due to the symmetry reduction and the associated electron confinement, thin metal films exhibit properties (both static and dynamic) remarkably different form their bulk counterparts. They are usually made of noble (Au, Ag, …) or simple metals (Na, K, …) and their sizes typically vary from a few nanometers to several hundred nanometers Maniyara2019 . As a notable example, the first experimental signature of ultrafast demagnetization was observed in thin nickel films in 1996 Beaurepaire1996 .

In this section, we illustrate the usefulness of phase-space methods to study the ultrafast electron dynamics in thin metal films. The particular geometry of thin films, whereby their thickness is much smaller than the transversal size, means that they can be described by a one-dimensional (1D) model, at least as far as longitudinal modes are concerned. Here, we mainly report on earlier work of ours on this topic, spanning the years from the early 2000s until today. The simulations are of increasing complexity. In Sec. 8.1, the electron motion is purely classical and follows the Vlasov equation, although the ground state is described by a Fermi-Dirac distribution. The main effect observed in these simulations is the emergence of a nonlinear ballistic mode due to bunches of electrons that travel across the film at the Fermi speed of the metal. Section 8.2 revisits this effect using a fully quantum Wigner model and it is a good example of interplay between wave-function methods (used to compute the ground state) and phase-space methods (used to compute the dynamics). The Wigner simulations allow us to identify the parameter range where the ballistic oscillations may be observed. Finally, in Sec 8.3, the description of the electron dynamics is augmented to include the spin degrees of freedom, according to the theoretical considerations presented in Sec. 5.2. The inclusion of these effects allowed us to study ballistic modes that include a magnetic component, thus exploring the possibility to generate a spin current in the film via an electromagnetic pulse.

8.1 Thin metal films: Vlasov

Several experiments have shown Brorson1962 ; Suarez1995 that electron transport in thin metal films occurs on a femtosecond timescale and involves ballistic electrons traveling at the Fermi velocity of the metal $v_{F}$. Later, a regime of low-frequency nonlinear oscillations, corresponding to ballistic electrons bouncing back and forth on the film surfaces, was measured in transient reflection experiments on thin gold films Liu2005 . These experimental findings were corroborated by the Vlasov-Poisson simulations carried out in our research group Manfredi2004 ; Manfredi2005film ; Manfredi2005optics , which are summarized in this section.

In these simulations, the electrons are initially prepared in a 1D Fermi-Dirac equilibrium at finite (typically, room) temperature:

[TABLE]

where $\epsilon(x,v_{x})=mv_{x}^{2}/2-e\phi(x)$ and $\phi$ is the electrostatic self-consistent potential, solution of the Poisson equation

[TABLE]

For the ion charge distribution, we use a jellium model with continuous density $n_{i}(x)=n_{0}/\left(1+\exp\left[\left(|x|-L/2\right)\sigma\right]\right)$, where $L$ is the thickness of the film, $n_{0}$ is the ion density of the bulk metal and $\sigma\approx 0.1\,\rm nm$ represents the typical distance over which the density falls to zero on the border of the film. The density profile is sketched in Fig. 9. The ions are supposed to be fixed during the electron motion because their typical timescale is significantly longer than that of the electrons. The effect of ionic vibrations (phonons) will be addressed in Sec. 8.2.

The electrons are subsequently excited by imposing a constant velocity shift $\delta v=0.08v_{F}$ to the initial equilibrium distribution, after which they evolve according to the 1D Vlasov-Poisson equations. This scenario is appropriate when no linear momentum is transferred parallel to the plane of the surface (i.e., $k_{\parallel}=0$), so that only longitudinal modes can be excited. As a reference case, we studied a sodium film ($r_{s}=4a_{0}$, where $a_{0}$ is the Bohr radius) with initial temperature $T_{e}=0.008\,T_{F}\approx 300$ K and thickness $L=10\lambda_{TF}\approx 12\,\rm nm$, where $\lambda_{TF}$ is the Thomas-Fermi screening length.

We consider the evolution of several energy quantities: (i) The self-consistent potential energy $E_{pot}=\frac{\varepsilon_{0}}{2}\int(\partial_{x}\phi)^{2}dx$, (ii) the center-of-mass energy $E_{cm}=\frac{1}{2}\int{\frac{j_{e}^{2}(x,t)}{n_{e}(x,t)}~{}dx}$ (where $j_{e}=\int v_{x}f_{e}dv_{x}$ is the electron current), and (iii) the thermal energy $E_{th}$. The latter is defined as: $E_{th}=E_{kin}-E_{cm}-E_{TF}$, where the total kinetic energy is $E_{kin}={m\over 2}\int v_{x}^{2}f_{e}dv_{x}$. The Thomas-Fermi energy, defined as $E_{TF}=\frac{E_{F}}{5}\int{n_{e}(x)^{5/3}dx}$, represents the kinetic energy of a state with same density but vanishing temperature and mean velocity, so that $E_{kin}(T_{e}=0)=E_{TF}$.

During an initial rapidly oscillating phase, the center-of-mass energy is almost entirely converted into thermal energy, due to Landau damping, whereas the thermal and potential energies reach a plateau. This plateau is characterized by low-frequency oscillations with a period close to $50\omega_{p}^{-1}\approx 5.3$ fs. This period is close to the time of flight of electrons traveling ballistically through the film at the Fermi velocity of the metal. The existence of such nonlinear ballistic mode is the main novel observation reported in Refs. Manfredi2004 ; Manfredi2005film . Later, we proposed to exploit this ballistic mode to enhance the absorption of infrared light in the film Manfredi2005optics .

The phase-space portraits of the electron distribution, shown in Fig. 11, clearly reveal that the perturbation starts at the film surfaces and then proceeds inward at the Fermi velocity of the metal. The structure formation at the Fermi surface, which has spread over the entire film for $\omega_{p}t>150$, is responsible for the increase of the thermal energy (and thus the electron temperature) observed in Fig. 10. As no collisional terms are present, the low-frequency oscillations of Fig. 10 are never completely damped away.

In order to include the effect of e-e collisions in a simple manner Manfredi2010 , we have augmented the Vlasov equation by a relaxation term, as described in Sec. 6. The constant in Eq. (64) is taken to be $\approx 0.4~{}{\rm fs^{-1}eV^{-2}}$, which was obtained from first-principles numerical simulations of the electron dynamics in sodium clusters Domps1998 . The electron temperature is computed instantaneously during the simulation, and plugged into the expression for the collision rate (5). When this collision model is applied to our case, the slow ballistic oscillations of Fig. 10 are still observed, but they are damped away on a timescale of the order of $500\omega_{pe}^{-1}\simeq 50\rm fs$ (see Fig. 12).

8.2 Thin metal films: Wigner

One of the most important properties of metal films, when their thickness is reduced to only a few nanometers, is the quantization of the energy levels along the direction perpendicular to the film surface (quantum size effects). The continuum that characterizes the bulk materials is replaced by a discrete spectrum of quantum states with their energies and wave functions depending on the boundary conditions Eguiluz1984 . The consequent modifications of the electronic structure with respect to a classical description lead to significant changes in their physical properties and, as shown in the following, their dynamical behaviour.

Our numerical study of the quantum electron dynamics in thin metal films is performed in two steps: first, the ground state of the electron population (possibly at finite temperature) is determined self-consistently; subsequently, the equilibrium distribution is perturbed by injecting some energy into the system. The quantum dynamics is then described by solving the time-dependent Wigner equation (27) coupled to the Poisson equation. Although the Wigner approach is a fully quantum-mechanical description, it is intrinsically ill-suited to deal with stationary states (indeed, quantization rules are overlooked in the Wigner formalism and must be imposed as additional constraints). Therefore it is more convenient to determine the ground-state from a standard density functional approach (Kohn-Sham equations) and then to construct the corresponding Wigner function from the computed Kohn-Sham wave functions.

The electron gas is assumed to be at thermal equilibrium with temperature $T_{e}$ and, as in the Vlasov approach described above in Sec. 8.1, while the ions are described by a continuous, immobile density with soft edges, as in Fig. 9. The electrons obey a Fermi-Dirac distribution at the corresponding temperature (for more details about the model, see Jasiak2009 ). The Kohn-Sham equations (22) are solved numerically using a finite-difference iterative method Eguiluz1984 . Some typical density profiles are shown in Fig. 13. In this example $N=18$ Kohn-Sham orbitals are occupied. The classical Vlasov-Poisson result (VP) yields a smooth density profile, without spatial oscillations. In the quantum regime (labeled WP), the standard Friedel oscillations are observed, particularly near the film surface (see the inset). Note that the profiles are weakly dependent on the electron temperature, because all considered cases are in the strongly degenerate regime $T_{e}\ll T_{F}$.

After computing the $N$ Kohn-Sham wave functions $\phi_{k}(x)$ and their respective occupation numbers $p_{k}$ Jasiak2009 , we can proceed to construct the ground-state Wigner function $f(x,v)$. We first compute the partial Wigner functions corresponding to each Kohn-Sham wave function $\phi_{k}$ of the ground-state

[TABLE]

The total Wigner function (Fig. 14) is given by the sum of the partial Wigner functions, weighed by the occupation numbers $p_{k}$, $f(x,v)=\sum_{k=1}^{N}p_{k}f_{k}(x,v)$.

To study the electron dynamics, we make use of the evolution equation for the Wigner function (27). In the classical limit (i.e., taking $\hbar\rightarrow 0$) this equation reduces to the Vlasov equation that was used in our previous works, summarized in Sec. 8.1. The Wigner results were obtained with a numerical code based on a uniform grid covering the relevant phase space Suh1991 . The method combines a splitting technique with fast Fourier transforms in the velocity coordinate.

Similarly to the Vlasov simulations described in Sec. 8.1, the electron dynamics is excited by perturbing the ground-state equilibrium. Again, this is done by imparting a uniform velocity shift $\delta v$ to the initial distribution. The ensuing electron dynamics also reveals some low-frequency oscillations for the thermal and potential energies, as in the Vlasov case. By plotting the observed period $T$ of the low-frequency oscillations against the value of the perturbation (see Fig. 15), we observe an interesting feature. While, for a large excitation $\delta v$, the period is close to the ballistic value $T=L/v_{F}$ (as in the Vlasov simulations, see Sec. 8.1), it diverges from this classical value below a certain excitation threshold $\delta v_{th}$

This departure from the ballistic oscillation period constitutes a clear transition between the classical and the quantum regime. The estimation of the threshold value requires the investigation of the microscopic electron dynamics in the phase space. In the classical case, we saw that the ballistic oscillations are due to bunches of nonequlibrium electrons traveling at the Fermi speed between the two surfaces of the film (Fig. 11). These bunches (i.e., vortices in the phase space) have a spatial extension roughly equal to $2\pi\lambda_{TF}$ and a width of the order of $\delta v$ in velocity space. The surface of the these vortices in the phase space (which has the dimension of an action) is thus approximately $A\approx 2\pi\lambda_{TF}m\delta v$. Quantum effects are expected to become significant when this action is of the same order as Planck’s constant, i.e. $A\approx\hbar$. For sodium, this leads to the following estimate for the threshold: $\delta v_{th}\approx 0.15v_{F}$, which is fairly close to the observed value (see Fig. 15). This is confirmed by the inspection of the classical and quantum phase-space portraits, as seen in Fig. 16. For the large excitation $\delta v=0.15v_{F}$ the phase-space vortices are clearly visible in the quantum Wigner simulations and they resemble the analog vortices of the Vlasov case. In contrast, the vortices have completely disappeared in the low-excitation case ($\delta v=0.06v_{F}$), whereas they are still present in the corresponding Vlasov portraits.

So far, we focussed exclusively on the ultrafast (up to $\approx 100\,\rm fs$) electron dynamics, while the ionic background was assumed to be frozen. In order to describe the long-time dynamics ($>1\,\rm ps$) it is mandatory to include electron-phonon interactions in the model (see Fig. 4 for a visual representation of the various timescales). In Ref. Jasiak2010 , we developed a microscopic phase-space model, based on the quantum Wigner distribution, which encompasses all relevant timescales from the femtosecond plasmon oscillations up to the phonon-mediated coupling to the ionic lattice, which occurs on a picosecond time scale. The main ingredients of the model are summarized in Sec. 6 in the paragraph devoted to Fokker-Planck methods.

For a system of quantum particles interacting with an environment, two distinct timescales are particularly important. The relaxation time measures the speed at which the energy is exchanged between the electrons and the lattice, while the decoherence time represents the typical time over which quantum correlations are lost to the external environment (here, represented by the ion lattice). In a density matrix language, the relaxation time corresponds to the decay of the diagonal terms, whereas the decoherence time is related to the nondiagonal terms of the density matrix.

Our simulations (see Ref. Jasiak2010 for details) showed that the relaxation time is of the order of $\tau_{R}\approx 2-3$ ps, for the sodium films considered here. It also increases with the amplitude of the excitation, in agreement with measurements performed on sodium clusters Maier2006 .

In our phase-space representation, the decoherence time $\tau_{D}$ is defined in a different way compared to the density matrix approach. We recall that the Wigner distribution function $f(x,v,t)$ can take negative values. The degree of “classicality” of a Wigner distribution can thus be estimated from the weight of its negative parts Kenfack2004 ; Deleglise2008 . The negative part of the Wigner function in the phase space is shown in Fig. 17 at different times, for a typical case. As time goes on, the quantum distribution progressively loses its negative values and becomes more and more classical.

In order to quantify more precisely this loss of classicality, we define the “negativity” of the distribution function as $S(t)=\int\int f_{<}(x,v,t)dxdv$, where $f_{<}=-f$ if $f<0$ and $f_{<}=0$ elsewhere ($S$ is then a positive quantity). The typical behavior of $S(t)$ is displayed in Fig. 18. The initial equilibrium distribution possesses only small negative parts, so the initial value of $S(0)$ is rather small. Soon after the excitation, when the system is driven out of equilibrium, $S$ grows very rapidly, then it decays more slowly to zero for longer times.

The decoherence time $\tau_{D}$ can be estimated by fitting $S(t)$ (after the initial sudden growth) with a simple decaying exponential function $S_{fit}(t)=S_{0}\exp(-t/\tau_{D})$. We observed that the the decoherence time $\tau_{D}$ increases with the excitation energy and decreases with the lattice temperature (Fig. 19), also in agreement with measurements on thin metal films Komori1987 . At room temperature ($T=300\,\rm K$) the decoherence time is roughly $\tau_{D}\approx 0.4\,{\rm ps}$, shorter than the relaxation time $\tau_{R}\approx 3\,\rm ps$.

8.3 Thin metal films: Spin

In this section, we extend the above investigation of the electron dynamics in thin metal films by including the spin degrees of freedom. To this end, we turn to ferromagnetic materials and study the coupled charge and spin dynamics in nickel films. We will use the self-consistent spin-Vlasov model defined in Eqs. (46)-(47). To keep the model as simple as possible, we restrict our problem to a reduced phase space ($x,v$), where the electrons are only allowed to move in the direction normal to the surfaces of the film ($x$). We only take into account electrostatic interactions through the Poisson equation and disregard effects due to exchange and correlations. The spin-orbit interaction, a higher order effect in $c^{-2}$, is also neglected. With these assumptions, the set of spin-Vlasov-Poisson equations that we solve is the following (see Sec. 5.2):

[TABLE]

where $k=(x,y,z)$ and $E_{x}=-\partial_{x}\phi$.

The typical parameters of the nickel films under consideration, expressed in both SI and atomic units (a.u.), are as follows: $r_{s}=2.6~{}\rm a.u.=0.14\,\rm nm$, $n_{0}=0.0136\rm~{}a.u.=91.8\rm~{}nm^{-3}$, $T_{p}=2\pi\omega_{p}^{-1}=15.32\rm~{}a.u.=0.37\rm~{}fs$, $v_{F}=0.74\rm~{}a.u.=1.62\rm~{}nm.fs^{-1}$, $\lambda_{TF}=1.80\rm~{}a.u.=0.095\rm~{}nm$ and $E_{\textrm{F}}=0.27\rm~{}a.u.=7.34\rm~{}eV$. We have studied films of thickness $L=5\rm~{}nm$.

In a realistic model of ferromagnetism, one should take into account the electronic structure of nickel, which comprises ten valence electrons of which eight are in the $3d$ shell and the remaining two are in the $4s$ shell. The $3d$ electrons are more localized around the ions than the $4s$ electrons. Therefore, we model the $4s$ (“itinerant”) electrons with the spin-Vlasov equations (111), whereas the localized $3d$ electrons and the corresponding ions form a motionless positively-charged background. The localized ions interact among themselves and with the itinerant electrons through magnetic exchange. The internal magnetic fields involved in these interactions can be very strong, of the order of $10^{3}-10^{4}\,\rm T$, and are responsible for the ferromagnetic properties of nickel. A study of laser-matter interaction using such realistic model of itinerant and localized magnetism was used recently to study spin-current generation in a nickel film Hurst2018 .

Here, to illustrate the main effects of the spin dynamics, we shall use a simplified “toy” model, where the internal magnetic field (generated by the magnetic ions) is replaced by a uniform external field of similar intensity directed along the $z$ axis: $\bm{B}=B\bm{e}_{z}$. This field will be responsible for the spin polarization of the electron gas. Hence, in this toy model the ions have no magnetic properties and are simply taken as a positively-charged density $n_{i}(x)$.

The ground state of the system is given by a 1D Fermi-Dirac equilibrium, see Eqs. (81)-(82). In practice, the ground state is obtained self-consistently by plugging Eqs. (81)-(82) into Poisson’s equation and solving for $\phi$ (for instance, iteratively) and then injecting the obtained potential back into the Fermi-Dirac distributions. This procedure can be performed numerically and the resulting ground state is shown in Fig. 20 for different magnetic fields. As expected, when $B=0$ the electron gas is totally unpolarized, meaning that the spin-up and spin-down densities are identical. Increasing the magnetic field induces a partial polarization of the electron gas, which becomes fully polarized when $B$ approaches $B_{F}=E_{F}/\mu_{B}$. The relation between the magnetic field and the degree of polarization is in agreement with the free electron gas model at zero temperature.

To simulate a laser pulse excitation, we shift the ground state distributions $f_{\uparrow}$ and $f_{\downarrow}$ in the velocity space by an amount $\delta v=0.05v_{F}$. This is equivalent to applying an instantaneous electric field in the direction normal to the film surfaces. We performed several simulations starting from different ground states with increasing spin polarization. To characterize the spin dynamics, we study the time evolution of the magnetic dipole

[TABLE]

as well as the electric dipole

[TABLE]

Both dipoles vanish in the ground state.

The dynamics of the magnetic dipole is shown in Fig. 21 for the different spin-polarized ground states. For an unpolarized electron gas ($B=0$), the magnetic dipole remains equal to zero, as the spin-up and spin-down distribution functions are identical for all times. For a partially polarized ground state, the magnetic dipole mainly oscillates at the ballistic frequency $\omega_{b}=2\pi v_{F}/2L=2\pi/6.17$ fs*-1*. These oscillations are similar to the ballistic oscillations observed in non-magnetic materials (see Secs. 8.1 and 8.2) and correspond to electrons traveling back and forth in the film at the Fermi speed.

The presence of a magnetic dipole implies that the spin-up and spin-down electron populations oscillate out of phase, because:

[TABLE]

Since the spin-up and spin-down components were excited in phase in the initial excitation, we deduce that their dephasing occurs during the first instants of the dynamics (roughly the first 5 fs in Fig. 21, before the magnetic dipole is clearly formed). One can show that this dephasing is a subtle effect due to the interplay of the self-consistent dynamics with the presence of strong electric fields at the boundaries of the film Hurst2018 . When the electron gas is fully polarized ($B=B_{F}$), no magnetic dipole is observed, simply because only one component (spin-up or spin-down) is present. Some faster oscillations at the plasma frequency $\omega_{p}=2\pi/0.37$ fs*-1* are also observed, but are rapidly removed by Landau damping. It is interesting to stress that, while the plasmon mode is quickly damped away, the magnetic mode at intermediate values of $B$ persists for much longer times.

Importantly, these ballistic oscillations of the magnetic dipole are still observed if one uses a more realistic model for both the magnetic material and the laser-film interaction Hurst2018 . We also stress that an oscillating magnetic dipole is equivalent to an AC spin current, since: $d\langle X\rangle_{m}/dt=\int v_{x}f_{z}dxdv_{x}=\int v_{x}(f_{\uparrow}-f_{\downarrow})dxdv_{x}=J_{\uparrow}-J_{\downarrow}\equiv J^{S}_{xz}$, where $J^{S}_{xz}$ is a spin current that propagates in the direction $x$ and is polarized along $z$ (in the general case ${\bm{J}}^{S}$ is a tensor, see Hurst2014 ). Thus the generation of a time-dependent magnetic dipole amounts to the generation of a spin current, which is an important issue for modern spintronic devices Alekhin2017 .

Finally, in Fig. 22 we show the time evolution of the electric dipole. Note that this is the same for all values of the magnetic field, because there is no backreaction of the spin dynamics on the charge dynamics in this toy model – see Eqs. (111) with a constant magnetic field. The electric dipole shows prominent plasmon oscillations that decay away through Landau damping. A weak signature of the ballistic mode can be seen between 5–15 fs.

9 Conclusions

A number of advanced methods to model the many-electron dynamics in nanometric condensed-matter systems have been developed over the years. In particular, TD-DFT has been developed to such a high level of accuracy to be capable of reproducing quantitatively many experimental results on large molecules and small nano-objects. Recent developments have extended its range of application to the spin dynamics Krieger2015 , shedding new light on old, but never properly understood, problems such as the ultrafast demagnetization occurring in ferromagnetic nano-objects illuminated with a laser pulse Beaurepaire1996 . All these methods are based on the propagation of a set of wave functions according to some Schrödinger-like equations. Notwithstanding these successes, TD-DFT and Hartree-Fock methods remain computationally very costly for large nano-objects containing more than a few hundred electrons. They are also not easily compatible with dissipative dynamics, because of their intrinsically Hamiltonian (hence, unitary) character, in spite of recent attempts to include dissipation in a mean-field formalism Dinh2018 .

Here, we have made a case for an alternative approach based on Wigner’s representation of quantum mechanics in the classical phase space. Phase-space models are certainly not a panacea and remain computationally costly for large systems. But they do have some advantages, in particular: (i) they are better adapted for the semiclassical limit and to compare to classical results, (ii) they can profit from existing numerical methods, either grid-based or particle-based, which have been developed over the years most notably in the plasma physics community, (iii) they can be adapted to include non-unitary evolutions in order to model dissipative effects, and, last but not least, (iv) the phase-space representation can be a powerful visual aid to help intuition and to display the numerical results.

The main scope of this review was to summarize the theoretical basis of phase-space methods in condensed matter physics, particularly at the nanoscale (Secs. 1-6) and to illustrate these concepts through two relevant examples, namely the propagation of linear waves in an infinite medium (Sec. 7) and the nonlinear electron dynamics in thin metal films (Se. 8). Beyond this, we would like to increase the awareness of plasma physicists for the emerging field of nanoplasmonics, i.e., the study of collective effect arising from the interaction between electromagnetic radiation and free electrons in metallic nano-objects Stockman2011 ; manfredi2018 , which is evolving very rapidly since its beginnings in the 1980s. The topics treated here show that the methods and tools of plasmas physics can be fruitfully extended and adapted to model typical plasma phenomena occurring in solid-state nanometric devices. Given these common goals and approaches, the present review hopes to encourage future exchanges between the plasmonics and plasma physics communities.