# Wave-particle interactions in a spin polarized plasma

Benjamin Bakri, Nicolas Crouseilles, Paul-Antoine Hervieux, Giovanni Manfredi

► To cite this version:

Benjamin Bakri, Nicolas Crouseilles, Paul-Antoine Hervieux, Giovanni Manfredi. Wave-particle interactions in a spin polarized plasma. 2025. hal-05370955

HAL Id: hal-05370955

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

Preprint submitted on 18 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Wave-particle interactions in a spin polarized plasma

Benjamin Bakri,¹, a) Nicolas Crouseilles,², b) Paul-Antoine Hervieux,¹, c) and Giovanni Manfredi¹, d)

¹) Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504, F-67000 Strasbourg, France

²) Université de Rennes, Inria Rennes (Mingus team) and IRMAR UMR CNRS 6625, F-35042, Rennes, France

(Dated: 16 November 2025)

## ABSTRACT

We study the interplay between spin waves (magnons) and plasma waves (plasmons) in a ferromagnetic material, using an augmented Vlasov-Poisson model that includes the electron spin dynamics. The ions are fixed, but their spins can evolve in time on the unit sphere according to the Landau-Lifshitz equation, which includes nearest-neighbor magnetic interactions. The two components interact not only through the electrostatic Coulomb force, but also via magnetic-exchange interaction terms. The linear response analysis reveals the existence of a wave-particle resonance occurring at the frequency of the magnons. This resonance gives rise to significant energy exchanges between the magnons and the electrons, resulting in a rapid loss of the localized magnetism akin to the ultrafast demagnetisation observed in experiments on thin ferromagnetic films. Depending on the initial electronic spin polarization, the resonance can lead to either damping or instability of the wave. These results show that wave-particle effects, similar to those frequently encountered in plasmas physics, may play a key role in spin-polarized plasmas and electron beams.

## I. INTRODUCTION

Ultrafast magnetism refers to the study of magnetic phenomena that occur on extremely short time scales, typically in the range of femtoseconds to picoseconds. It focuses on the rapid evolution of magnetic order following external perturbations, such as laser irradiation, and provides fundamental insight into the spin dynamics at the nanoscale. A central topic within ultrafast magnetism research is the interaction between ultrashort laser pulses and magnetic materials, which enables the observation of sub-picosecond changes in the magnetic properties of the sample. The seminal work by Beaurepaire et al.¹ demonstrated that such laser excitation can induce a rapid loss of magnetization in ferromagnetic thin films. Despite significant progress, the microscopic mechanisms underlying this ultrafast demagnetization process remain the subject of ongoing investigation²,³.

Beyond such fundamental physics issues, achieving control over ultrafast demagnetization across a wide range of magnetic materials is essential for the development of future high-speed spintronic and magneto-optical devices. While early studies concentrated on ferromagnetic thin films⁴,⁵, subsequent research has extended to systems with more complex magnetic ordering, including ferrimagnets⁶ and antiferromagnets⁷. It is also known that the ultrafast demagnetization facilitates the conversion of laser pulses into terahertz (THz) electrical currents⁸,⁹. This process relies on the generation of spin

currents in response to electromagnetic pulses, which has been described in recent computational studies¹⁰. Finally, spin polarized plasmas have been considered in the contexts of nuclear fusion research¹¹,¹² and plasma wakefield acceleration, where the generation of highly spin-polarized electron gases via spin filters has been recently demonstrated experimentally¹³,¹⁴.

A widely adopted theoretical framework to describe the ultrafast demagnetization is the so-called three-temperature model¹. It is a macroscopic (non space-resolved), phenomenological model that decomposes the system into three interacting subsystems – electrons, ion lattice, and spins – each subsystem being described as a thermal bath at instantaneous thermodynamic equilibrium. Using this approach, it was shown that electron-magnon scattering (representing energy transfer between the electronic and spin subsystems) occurs on ultrashort timescales and significantly influences the initial phase of the ultrafast demagnetization¹⁵⁻¹⁷.

In the present work, we investigate the ultrafast magnetization dynamics from the viewpoint of plasma physics, using a microscopic phase-space model that goes well beyond the simple quasi-equilibrium approach adopted in the three-temperature model¹⁸⁻²⁰. In our model, we distinguish two types of magnetism: (i) itinerant magnetism, carried by the conduction electrons that are spatially delocalized (s orbitals) and (ii) fixed magnetism, carried by localized d-orbital electrons, hereafter denoted simply as “the ions”. The itinerant electrons are described by a set of Vlasov equations generalized to the case of spin-1/2 fermions, whereas the fixed ions are modeled by a Landau-Lifshitz equation²¹, which governs the evolution of the ion magnetization $S(x, t)$ in space and time. The itinerant electrons and the fixed ions interact with each other through electro-

<footnote>a) Electronic mail: benjamin.bakri@ipcms.unistra.fr</footnote>

<footnote>b) Electronic mail: nicolas.crouseilles@inria.fr</footnote>

<footnote>c) Electronic mail: paul-antoine.hervieux@ipcms.unistra.fr</footnote>

<footnote>d) Electronic mail: giovanni.manfredi@ipcms.unistra.fr</footnote>
static Coulomb forces (Poisson equation) and through magnetic-exchange interactions modeled in the RKKY (Ruderman-Kittel-Kasuya-Yosida) approximation<sup>22</sup>. A sketch of the different species intervening in the model, together with their respective interactions, is provided in figure 1. We will refer to this mathematical model as the spin-Vlasov-Poisson-Landau-Lifshitz (sVPLL) set of equations.

The fastest timescale in this system is given by the plasmon frequency $\omega_p \approx 10^4$ THz, corresponding to times around 0.1 fs, whereas oscillations of the ion spins (magnons) have much smaller frequencies $\omega_{mag} \approx 10 - 1000$ THz. Phonon modes, i.e., vibrations of the atomic lattice, occur at even slower frequencies ($\approx 1$ THz), corresponding to several picoseconds<sup>23</sup>. As we are interested in the dynamics on timescales of a few hundred femtoseconds, electron-phonon interactions will be neglected here.

An earlier version of the model, together with a description of the related computational methods, was developed in a prior work<sup>24</sup>. There, the electron gas was supposed to be weakly degenerate and therefore described using the Maxwell-Boltzmann distribution, a scenario that is appropriate for spin-polarized electron beams<sup>13,14</sup>. In contrast, in solid state systems, where the Fermi temperature is typically larger than $10^4$ K, the electron gas is highly degenerate even at room temperature and a Fermi-Dirac distribution should be used.

The present work aims at investigating the combined electron-magnon (i.e., charge-spin) collective dynamics in a realistic ferromagnetic material. We will show that resonant wave-particle effects, similar to those well-known from plasma physics, can occur in these systems. In particular, magnon waves can resonate with the electrons, as their phase velocity falls within the electron Fermi-Dirac distribution—i.e., $\omega_{mag}/k \lesssim v_F$, where k is the wave number and $v_F$ the Fermi velocity. This resonance gives rise to significant energy exchanges between the spin modes (magnons) and the electrons, resulting in a rapid loss of the localized magnetism akin to the ultrafast demagnetisation observed in the experiments<sup>1</sup>. Note that, in contrast, plasmons are not resonant, as $\omega_p/k \gg v_F$, so that Landau damping and other resonant plasma effects are excluded.

The present work is organized as follows. The details of the mathematical model are illustrated in Sec. II. The linear response and dispersion relation is analyzed in Sec. III. Numerical simulations in the nonlinear regime are shown in Sec. IV and conclusions are drawn in Sec. V.

## II. SPIN VLASOV MODEL

The system considered in this study consists of a ferromagnetic layer with a thickness of the order of a few nanometres, maintained at room temperature. For the purpose of defining the atomic and electronic parameters, we adopt the properties of pure nickel (Ni), characterized by an atomic density $n_i = 9.1 \times 10^{28} \text{ m}^{-3}$.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic view of the electron-magnon system. The electron cloud, illustrated as dark blue dots with arrows, represents the itinerant component of magnetism. In contrast, the ion spin chain, shown in red, corresponds to the localized magnetic component. Electron-electron interactions are governed by the electrostatic potential and involve the vacuum permittivity $\epsilon_0$. Ion-ion interactions arise from magnetic exchange between nearest-neighbor ions and are proportional to the exchange constant J. The coupling between electrons and ions is described by the magnetic exchange interaction in the RKKY approximation<sup>22</sup>, which links the spin orientations of the two subsystems through the coupling constant K.</figcaption></figure>

The corresponding interatomic spacing a is set to be equal to twice the Wigner-Seitz radius $r_s$, defined so that $4\pi n_i r_s^3/3 = 1$, yielding $a = 0.275$ nm. In nickel, two electrons per atomic site populate the conduction band ($Z = 2$). However, for simplicity, a charge state of $Z = 1$ will be assumed everywhere in this work, so that the electron and ion densities are equal, $n_i = n_e$. This assumption that does not alter the physical interpretation or conclusions of this work. The Fermi energy and Fermi velocity are respectively: $E_F = \frac{\hbar^2}{2m_e} (3\pi^2 n_e)^{2/3} = 7.41$ eV and $v_F = \sqrt{2E_F/m_e} = 1.61 \times 10^6 \text{ m/s}$, where $m_e$ and $n_e$ are the electron mass and number density.

A one-dimensional (1D) configuration with periodic boundaries is adopted, in which both position and velocity degrees of freedom are aligned along the x-axis. In contrast, the spin vector (both for the electrons and the ions) is allowed to evolve in three spatial dimensions.

### A. Evolution equations

The spin-Vlasov model<sup>23,24</sup> describes the dynamics of a spin-polarized electron gas through a set of coupled kinetic equations. In this framework, the electronic degrees of freedom are represented by two distribution functions defined over position-velocity phase space: a scalar function $f_0(x, v, t)$ and a vector function $\mathbf{f}(x, v, t)$, with $\mathbf{f} = (f_x, f_y, f_z)$. The scalar distribution $f_0$ corresponds to the charge density, while the vector distribution $\mathbf{f}$ characterizes the spin density. Both quantities evolve according to the coupled set of Vlasov equations:

$$ \begin{aligned} & \frac{\partial f_0}{\partial t} + v \frac{\partial f_0}{\partial x} - \frac{1}{m_e} \frac{\partial V_H}{\partial x} \frac{\partial f_0}{\partial v} - \frac{\mu_B}{m_e} \frac{\partial \mathbf{B}_{ion}}{\partial x} \cdot \frac{\partial \mathbf{f}}{\partial v} = 0 && (1) \\ & \frac{\partial \mathbf{f}}{\partial t} + v \frac{\partial \mathbf{f}}{\partial x} - \frac{1}{m_e} \frac{\partial V_H}{\partial x} \frac{\partial \mathbf{f}}{\partial v} - \frac{\mu_B}{m_e} \frac{\partial \mathbf{B}_{ion}}{\partial x} \frac{\partial f_0}{\partial v} - \frac{e}{m_e} \mathbf{B}_{ion} \times \mathbf{f} = 0, && (2) \end{aligned} $$
where $e > 0$ is the absolute value of the electron charge, $\hbar$ is the reduced Planck constant, $m_e$ is the electron mass, and $\mu_B = e\hbar/(2m_e)$ is the Bohr magneton. Here, $V_H$ is the Hartree potential, solution to the Poisson equation, while $\mathbf{B}_{ion}$ represents the magnetic field generated by the ions and experienced by the electrons. The above Vlasov equations are semiclassical inasmuch as the orbital motion of the electrons follows the classical trajectories, while their spin is treated as a quantum variable. The equations (1)–(2) can be derived from a fully quantum model based on the total Wigner function of the system$^{23}$.

All physical quantities can be obtained by integration of the electron distributions. For instance, the electron number density $n_e = \int f_0 dv$, magnetization $\mathbf{m} = \int \mathbf{f} dv$, and the charge and spin currents in the $x$-direction: $J_C = -e \int v f_0 dv$ and $\mathbf{J}_S = -e \int v \mathbf{f} dv$. The spin-up and spin-down distributions are related to the $f_0$ and $f_z$ components in the following way: $f_0 = f_\uparrow + f_\downarrow$ and $f_z = f_\uparrow - f_\downarrow$. The spin-resolved electron densities are given by: $n_{\uparrow,\downarrow} = n_e \pm \|\mathbf{m}\|$.

The above spin Vlasov equations describe the dynamics of the itinerant magnetism, carried by the conduction electrons that are not tied to a particular ion, but rather free to travel within the metal. In order to describe the fixed, localized magnetism, we consider an ionic spin chain where each spin $\mathbf{S}_i$ interacts with its nearest neighbors and is subject to a magnetic field $\mathbf{B}_{elec}$. This configuration can be modeled by the Heisenberg Hamiltonian

$$ \mathcal{H}_H = -\sum_i J \mathbf{S}_i \cdot \mathbf{S}_{i+1} - \mu_B \sum_i \mathbf{B}_{elec} \cdot \mathbf{S}_i, $$

where $J$ is the magnetic-exchange coupling constant, and $\mathbf{B}_{elec}$ is the magnetic field generated by the itinerant electrons. In the continuous limit, the ion spin vector becomes a continuous function of the position variable $\mathbf{S}(x,t)$, which obeys the Landau-Lifshitz equation$^{21}$

$$ \frac{\partial \mathbf{S}}{\partial t} = \frac{a^2 J}{\hbar} \mathbf{S} \times \frac{\partial^2 \mathbf{S}}{\partial x^2} - \frac{\mu_B}{\hbar} \mathbf{S} \times \mathbf{B}_{elec}, \quad (3) $$

where $a$ is the interatomic spacing. Note that $\mathbf{S}(x,t)$ is a dimensionless quantity. In the case of nickel atoms, the ion-ion magnetic-exchange coupling constant is estimated$^{10}$ to be $J = 0.022$ eV.

## B. Interactions

As illustrated in figure 1, there are three types of interactions governing the dynamics of the itinerant electrons and fixed ions: (i) the electrostatic Coulomb interactions between the electrons and the ions (coupling constant $\varepsilon_0$), (ii) the ion-ion magnetic exchange (coupling constant $J$), and (iii) the ion-electron magnetic exchange (coupling constant $K$). In the sVPLL model, all interactions are treated under the mean field approximation.

The self-consistent electrostatic potential (known as the Hartree potential in condensed matter physics) is a solution of the Poisson equation

$$ \frac{\partial^2 V_H}{\partial x^2} = \frac{e^2}{\varepsilon_0} \int f_0 dv - Z n_i, \quad (4) $$

where $n_i$ is the homogeneous ion density and $\varepsilon_0$ is the vacuum permittivity.

The electron-ion magnetic-exchange interaction is treated in the RKKY approximation$^{22}$. In this approximation, the ions experience an effective magnetic field $\mathbf{B}_{elec}$ proportional to the magnetization of the electrons

$$ \mathbf{B}_{elec}(x,t) = -\frac{K}{2\mu_B} \int \mathbf{f}(x,v,t) dv, \quad (5) $$

where $K$ is the electron-ion magnetic-exchange coupling constant. Conversely, the electrons experience an effective magnetic field $\mathbf{B}_{ions}$ generated by the magnetization of the ions

$$ \mathbf{B}_{ion}(x,t) = -\frac{K n_i}{2 \mu_B} \mathbf{S}(x,t). \quad (6) $$

For nickel atoms, the value of the coupling constant can be estimated as$^{10}$ $K = 0.014$ eV nm$^3$.

The sVPLL system of equations (1)–(2) and (3)–(6) conserves the total energy

$$ \begin{aligned} E_{\text{tot}} = & \frac{m_e}{2} \iint v^2 f_0 \, dx \, dv + \mu_B \iint \mathbf{f} \cdot \mathbf{B} \, dx \, dv \\ & + \frac{\varepsilon_0}{2} \int \left( \frac{\partial V_H}{\partial x} \right)^2 \, dx + \frac{1}{2} a^2 J Z n_i \sum_{j=\{x,y,z\}} \int \left( \frac{\partial S_j}{\partial x} \right)^2 \, dx \end{aligned} \quad (7) $$

and each component of the total magnetization

$$ \mathbf{M}_{\text{tot}} = \frac{Z}{2} \mu_B \iint \mathbf{f} \, dx \, dv + \mu_B n_i \int \mathbf{S} \, dx. \quad (8) $$

Note that, because of the 1D nature of our model, the above quantities have to be understood per unit surface.

## III. EQUILIBRIUM AND LINEAR RESPONSE

In this section, we investigate the linear response of the sVPLL system, composed of the equations (1)–(2) and (3)–(6). The procedure follows the analysis developed in our earlier work for Maxwell-Boltzmann equilibria$^{24}$.

### A. Equilibrium

We decompose the electron distributions and the ion spin density into a spatially homogeneous equilibrium and a small perturbation, such as: $f_j(x, v, t) = f_j^{(0)}(v) + f_j^{(1)}(x, v, t)$, with $j \in \{0, x, y, z\}$, and $\mathbf{S}(x, t) = \mathbf{S}^{(0)} +$
$\mathbf{S}^{(1)}(x, t)$, where the superscripts (0) and (1) denote respectively the equilibrium and the perturbation. We further assume that the equilibrium corresponds to a ferromagnetic configuration with all spins aligned along the z-axis. Accordingly, the equilibrium components satisfy $f_x^{(0)} = f_y^{(0)} = S_x^{(0)} = S_y^{(0)} = 0$, and $S_z^{(0)} = 1$. The other components of the electron distribution function $f_0^{(0)}$ and $f_z^{(0)}$ are taken to be Fermi-Dirac distributions at zero temperature. The spin polarization of the electrons at equilibrium may be quantified by the parameter $\eta \in [-1, 1]$, defined as

$$\eta = \frac{\int f_z^{(0)}(v) dv}{\int f_0^{(0)}(v) dv}. \qquad (9)$$

For $\eta = 1$ (resp., $-1$) all electrons are spin polarized (resp., anti-polarized) along the z-axis; intermediate values of $\eta$ correspond to partial polarization.

In ferromagnetic materials, the Fermi temperature is typically of the order $T_F \approx 5 \times 10^4$ K. Hence, at room temperatures, the electron gas is fully degenerate and its velocity distribution is well approximated by a Fermi-Dirac distribution at zero temperature. Quantum mechanically, this means that the energy levels are all occupied up to the Fermi energy $E_F = \frac{\hbar^2}{2m_e}(3\pi^2 n_e)^{2/3}$, while they are empty for energies exceeding $E_F$. In our semiclassical approximation, this corresponds to an equilibrium velocity distribution that is constant inside the Fermi sphere, defined by $|v|^2 \le v_F^2$, where $v_F = \sqrt{2E_F/m_e}$ is the Fermi velocity.

For spin polarized electrons, the radius of the Fermi sphere differs for spin-up and spin-down electrons in the presence of an external magnetic field. In velocity space, we will denoted $v_\uparrow$ and $v_\downarrow$ the Fermi velocities of each spin component. Further, in the one-dimensional (1D) case considered here, the 1D velocity distribution along the x-axis is obtained by integrating the three-dimensional distribution over the two transverse velocity components, $v_y$ and $v_z$. The resulting 1D Fermi-Dirac distribution at zero temperature takes a parabolic form in the velocity variable$^{25}$.

Remembering that the spin-up and spin-down distributions are related to the $f_0$ and $f_z$ components in the following way: $f_0 = f_\uparrow + f_\downarrow$ and $f_z = f_\uparrow - f_\downarrow$, we can write the equilibrium 1D Fermi-Dirac distributions as

$$f_{\uparrow}^{(0)}(v) = \begin{cases} 0, & \text{for } v^2 > v_{\uparrow}^2 \\ \pi C (v_{\uparrow}^2 - v^2), & \text{for } v^2 \le v_{\uparrow}^2, \end{cases} \qquad (10)$$

$$f_{\downarrow}^{(0)}(v) = \begin{cases} 0, & \text{for } v^2 > v_{\downarrow}^2 \\ \pi C (v_{\downarrow}^2 - v^2), & \text{for } v^2 \le v_{\downarrow}^2, \end{cases} \qquad (11)$$

where $C = [m_e/(2\pi\hbar)]^3$. More details on how to obtain the above Fermi-Dirac equilibrium are provided in the Appendix A, where figure 11(a) shows the profiles of the equilibria.

The Fermi velocities $v_{\downarrow,\uparrow}$ can be written in terms of the unpolarized Fermi velocity $v_F$ and the spin polarization $\eta$ at equilibrium. We have

$$n_e = \int f_{\uparrow}^{(0)}(v) dv + \int f_{\downarrow}^{(0)}(v) dv. \qquad (12)$$

Inserting the equilibria (10)-(11) into the above expression leads to $n_e = \frac{4}{3}\pi C(v_\uparrow^3 + v_\downarrow^3)$. Using the definition of the Fermi energy and velocity in terms of the electron density, we obtain $v_\uparrow^3 + v_\downarrow^3 = 2v_F^3$. From Eq. (9), the electron polarization becomes

$$\eta = \frac{v_{\uparrow}^{3} - v_{\downarrow}^{3}}{v_{\uparrow}^{3} + v_{\downarrow}^{3}}, \qquad (13)$$

which yields, finally,

$$v_{\uparrow} = v_F(1 + \eta)^{1/3}, \quad v_{\downarrow} = v_F(1 - \eta)^{1/3}. \qquad (14)$$

## B. Linear response and dispersion relation

Expressing the perturbed quantities, such as $\mathbf{f}^{(1)}(x, v, t)$ and $\mathbf{S}^{(1)}(x, t)$ in terms of plane waves $\exp(-i\omega t + ikx)$, where $k$ is the wave vector and $\omega$ the frequency of the wave, leads to a dispersion function that can be factored as $D(\omega, k) = D_E(\omega, k) D_S(\omega, k)$. The first branch $D_E(\omega, k) = 0$ leads to the usual dispersion relation for electrostatic Vlasov-Poisson plasmas (Bohm-Gross relation), which includes Langmuir waves and Landau damping. The other branch of the dispersion relation $D_S(\omega, k) = 0$ corresponds to collective spin modes that are usually referred to as magnons in the literature$^{26}$ and will be analyzed in details in the forthcoming paragraphs.

The magnon branch of the dispersion function can be written as$^{24}$

$$D_S(\omega, k) = \omega - \frac{a^2 J k^2}{\hbar} - \frac{Z \omega_L \eta}{2} - \frac{\omega_L^2}{2Zn_i} \left( \frac{\hbar k}{2m_e} \int \frac{\partial_v f_0^{(0)}(v)}{kv - \omega + \omega_L} dv - \int \frac{f_z^{(0)}(v)}{kv - \omega + \omega_L} dv \right), \qquad (15)$$
where $\omega_L = eB_{ion}^{(0)}/m_e = Kn_i/\hbar$ is the Larmor frequency of the electrons immersed in the equilibrium magnetic field $B_{ion}^{(0)}$ generated by the ions and directed along the z-axis.

In the limit of small electron polarization $\eta$ and using a Fermi-Dirac equilibrium at zero temperature defined in Eqs. (10) and (11), the real and the imaginary parts of the frequency can be approximated as (further details are given in the Appendix B):

$$ \omega_r = \omega_S + \frac{Z}{2}\omega_L\eta - \frac{3\hbar\omega_L^2}{8E_F} + \frac{1}{4}\frac{\omega_L^2}{kv_F}\left(\frac{3\hbar\omega_L}{4E_F} - \eta\right)\ln\left|\frac{kv_F + \omega_L}{kv_F - \omega_L}\right| \quad (16) $$

and

$$ \omega_i = \frac{\pi \omega_L^2}{4 kv_F} \left( \frac{3 \hbar \omega_L}{4 E_F} - \eta \right), \quad \text{if } \frac{\omega_L}{kv_F} < 1 \\ = 0, \quad \text{otherwise.} \quad (17) $$

As mentioned earlier, we take $Z = 1$ everywhere in this work.

In the absence of electron-magnon coupling ($K = 0$), the magnon frequency is $\omega_{mag} = \omega_S = k^2 a^2 J/\hbar \approx 2.6 \times 10^{12} \text{ s}^{-1}$, which is significantly smaller than the plasmon frequency $\omega_p \approx 1.7 \times 10^{16} \text{ s}^{-1}$. When the electron-ion coupling is introduced, the ratio $\omega_p/\omega_{mag}$ becomes somewhat larger, but the following scaling remains valid:

$$ \omega_p \gg kv_F \sim \omega_L \gg \omega_{mag}. $$

Typical values of the relevant physical parameters are summarized in Table I.

The real and imaginary parts of the magnon frequency are shown in figure 2 as a function of the electron spin polarization $\eta$ (left panels) and the wave vector $k$ (right panels). The dependence of the frequency on $\eta$ is particularly interesting: below a certain positive $\eta_{th}$ the imaginary frequency is positive, corresponding to an instability, whereas above the threshold it is negative, signalling the damping of the initial perturbation.

Such threshold corresponds precisely to the electron spin polarization induced by the magnetic field generated by the ions at equilibrium, i.e., $B_{ion}^{(0)} = -Kn_i/(2\mu_B)$. A simple calculation (see Appendix A) shows that the threshold spin polarization is related to the coupling constant $K$ through the expression

$$ K = \frac{E_F}{n_i} \left[ (1 + \eta_{th})^{2/3} - (1 - \eta_{th})^{2/3} \right], \quad (18) $$

which can be written as $\eta_{th} \approx \frac{3}{4} \frac{Kn_i}{E_F}$, for $\eta_{th} \ll 1$. This value exactly cancels the imaginary part of the frequency in the approximate expression given in Eq. (17). Using the parameters of Table I, we obtain $\eta_{th} \approx 0.129$. This is the “natural” state of polarization of the electrons, in the absence of any other constraints such as an external magnetic field.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Linear response of the magnon modes. The solid orange lines represent the exact dispersion relation (15), the solid blue lines represent the approximate dispersion relations (16)-(17), and the green dots correspond to numerical simulations of the full sVPLL system close to equilibrium. The left panels show the real (a) and imaginary (c) frequencies as functions of the electron spin polarization $\eta$, varying from -1 to 1, for a fixed wave vector $k = 5.29 \text{ nm}^{-1}$. The dashed vertical line corresponds to the value $\eta_{th}$ for which the imaginary frequency vanishes. The right panels show the real (b) and imaginary (d) parts of the frequency as functions of the wave vector $k$, varying from $0.3 \text{ nm}^{-1}$ to $5.3 \text{ nm}^{-1}$, and vanishing spin polarization $\eta = 0$. Other parameters are given in Table I.</figcaption></figure>

According to figure 2, if $\eta > \eta_{th}$ (i.e., if the electrons are more strongly spin polarized than they would naturally be by the ions), then the electronic spin polarization is damped and tends to return to its initial value. In contrast, if $\eta < \eta_{th}$ an instability is observed: perturbations are amplified and bring the system away from its equilibrium. The instability rate increases with decreasing $\eta$, and is maximal when $\eta = -1$, i.e., when all the electrons are spin-polarized in the opposite direction to that of the ions (anti-polarization).

The presence of an imaginary part in the magnon frequency clearly suggests the occurrence of resonant wave-particle interactions. The resonance is visible in the expression of the full dispersion relation, Eq. (15), where the denominator vanishes for some value of the electron velocity, similarly to what happens for Vlasov-Poisson plasmas in the case of Landau damping.

Resonance occurs for a velocity $v_{res}$ such that

$$ \frac{\omega}{k} = \frac{\omega_L}{k} + v_{res}, \quad (19) $$

where $\omega/k$ denotes the phase velocity of the ion spin wave (magnon) and $\omega_L/k$ represents the phase velocity of the electron spin precessing in the magnetic field generated by the ions. This condition implies that resonance is achieved when the electron spin precesses at the same fre-
<table><caption>TABLE I. Summary of the physical parameter used in the simulations, broadly relevant to a thin nickel layer.</caption><thead><tr><th>Z</th><th>n<sub>i,e</sub> [m<sup>-3</sup>]</th><th>ω<sub>p</sub> [s<sup>-1</sup>]</th><th>E<sub>F</sub> [eV]</th><th>v<sub>F</sub> [m/s]</th><th>L<sub>F</sub> [nm]</th><th>K [eV nm<sup>3</sup>]</th><th>ω<sub>L</sub> [s<sup>-1</sup>]</th><th>J [eV]</th><th>a [nm]</th></tr></thead><tbody><tr><td>1</td><td>9.1 × 10<sup>28</sup></td><td>1.71 × 10<sup>16</sup></td><td>7.41</td><td>1.61 × 10<sup>6</sup></td><td>9.45 × 10<sup>-2</sup></td><td>0.014</td><td>1.95 × 10<sup>15</sup></td><td>0.022</td><td>0.275</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Electron velocity distribution functions at position $x = 0$, in arbitrary units. The black dashed line represents the initial Fermi-Dirac distribution $f_0(x = 0, v, t = 0)$, the solid colored lines show the velocity distributions at a later time for two different wave vectors, $k = 1.6 \text{ nm}^{-1}$ (blue line, final time $t = 36.5$ fs) and $k = 2.6 \text{ nm}^{-1}$ (orange line, final time $t = 124$ fs). The dashed vertical lines show the corresponding resonant velocities.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Temporal evolution of the magnon for three values of the electron spin polarization, $\eta = 0.9$ (green line), $\eta = -0.9$ (blue line), and $\eta = -0.1$ (orange line). The figure represents the amplitude of the lowest Fourier mode of the $S_x(x, t)$ component of the ion magnetization. The wave vector is $k = 5.29 \text{ nm}^{-1}$; other parameters are those given in Table I.</figcaption></figure>

quency as the magnon, adjusted by a Doppler shift due to the electron's velocity $v_{\text{res}}$ relative to the fixed ions. This resonance exhibits similarities to the Electron Cyclotron Resonance Heating (ECRH) mechanism observed in fusion plasmas, albeit with two key distinctions. First, the ion spin wave (magnon) assumes the role analogous to that of the externally applied electromagnetic wave in ECRH. Second, unlike ECRH where the electron magnetic moment arises from orbital motion, in this context it originates from the intrinsic spin of the electron.

The damped regime observed for $\eta > \eta_{\text{th}}$ is suggestive of an effect known as Gilbert damping in magnetic materials$^{27,28}$, whereby the magnetization vector tends to return to its initial direction when perturbed. Gilbert damping is usually introduced phenomenologically in the Landau-Lifshitz equation (3), although its microscopic origin is still unclear. Our result shows that electron-magnon interactions may constitute a possible channel contributing to the overall Gilbert damping effect.

In the unstable regime ($\eta < \eta_{\text{th}}$), the magnon amplitude increases in time until it reaches a nonlinear regime where the linear approximation is no longer valid. In this case, significant exchanges can occur between the ions and the electrons, leading to a loss of magnetization on a femtosecond timescale, which is reminiscent of the ultra-fast demagnetization long observed in the experiments$^1$. These aspects will be discussed in Section IV.

Figures 2(b) and 2(d) show the dependence of the

magnon frequency on the wave vector $k$, for unpolarized electrons with $\eta = 0$. The imaginary part of the frequency becomes nonzero above a certain wave vector corresponding to the relation $v_{\text{res}} = v_F$, where $v_{\text{res}}$ is the resonant velocity defined in Eq. (19). Indeed, for a Fermi-Dirac distribution at zero temperature [Eqs. (10)-(11)], there are no electrons with velocities larger than $v_F$, hence no wave-particle resonance is possible. Above this threshold, $\omega_i$ rises steeply because the Fermi-Dirac distribution has maximum gradient at $v = v_F$. Using $k = 1.1 \text{ nm}^{-1}$ [read from the figure 2(d)] and the parameters from Table I, one obtains $v_{\text{res}} = 1.63 \times 10^6 \text{ m/s}$, which is indeed close to $v_F = 1.61 \times 10^6 \text{ m/s}$.

The wave-particle resonance is distinctly observable in the velocity distribution at $x = 0$, as shown in figure 3 for two values of the wave vector. This result is obtained from numerical simulations of the full sVPLL system operating in the linear regime, where only a single mode is weakly excited at the initial time. The presence of resonant electrons leads to a distortion of the Fermi-Dirac equilibrium near the resonant velocity $v_{\text{res}}$, characterized by a flattening of the distribution. This behavior is reminiscent of classical plasma phenomena such as Landau damping.

The evolution of the magnon amplitude ($S_x$ component of the ion magnetization) is shown in figure 4 for three representative values of the electron spin polarization.
The cases $\eta = \pm 0.9$ are respectively damped ($\eta > 0$) and growing ($\eta < 0$), and saturate nonlinearly at about the same time, $t \approx 40$ fs. The observed damping may be linked to the Gilbert damping effect mentioned earlier$^{27}$. The case $\eta = -0.1$ is weakly unstable, as expected from the diagram on figure 2(c). The real part of the frequency is much larger for this case, in accordance with figure 2(a).

Finally, we note that the charge branch of the dispersion relation, $D_E(\omega, k) = 0$, does not exhibit any resonant behavior. This is due to the fact that, in solid-state plasmas, the plasma frequency vastly exceeds the magnon frequency (see Table I), so that the phase velocity $\omega_p/k$ is much larger than the Fermi velocity of the metal and no wave-particle interactions can occur. Therefore, the resonant wave-particle interactions observed in this system are purely magnetic in nature and do not arise in spin-less plasmas. Nonetheless, the underlying resonant mechanisms bear resemblance to those encountered in classical plasma phenomena such as Landau damping.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Time evolution of the Fourier modes $k_n = nk_{min} = n2\pi/L$ of the ion magnetization $S_x(x, t)$ along the x-axis. The amplitude of the modes is given by the color code shown at the right of the figure. At time $t = 0$, the ion spins are weakly excited at the mode $k_1 = 0.264 \text{ nm}^{-1}$, which is linearly stable. At later times ($t \approx 50$ fs) several unstable higher harmonics with $n > 1$ appear, and for $t > 200$ fs the Fourier spectrum becomes very broad. The electron spin polarization is $\eta = 0$.</figcaption></figure>

## IV. NONLINEAR REGIME

When the initial configuration is unstable, the perturbation grows exponentially and quickly attains the nonlinear regime, where the linear response analysis of the preceding section no longer applies. To explore the nonlinear regime, we have solved numerically the full sVPLLS system using an Eulerian code based on a Hamiltonian splitting method, which was described in detail in our earlier work$^{24}$.

The boundary conditions are periodic in $x$, with period $L$ and corresponding wave vector $k_{min} = 2\pi/L$. The initial condition is a homogeneous Fermi-Dirac distribution as in Eqs. (10)-(11), perturbed sinusoidally at a certain wave vector $k$, usually the lowest mode $k_{min}$:

$$ 
\begin{aligned}
f_0(x, v, t = 0) &= f_0^{(0)}(v) [1 + \varepsilon \cos(kx)], \\
f_z(x, v, t = 0) &= f_z^{(0)}(v), \\
f_x(x, v, t = 0) &= f_0^{(0)}(v) \varepsilon \cos(kx), \\
f_y(x, v, t = 0) &= f_0^{(0)}(v) \varepsilon \sin(kx),
\end{aligned}
$$

where $f_0^{(0)} = f_{\uparrow}^{(0)} + f_{\downarrow}^{(0)}$, $f_z^{(0)} = f_{\uparrow}^{(0)} - f_{\downarrow}^{(0)}$, and $\varepsilon \ll 1$ is a small perturbation. The above initial condition describes an electron gas that is spin polarized along the z-direction, with small perturbations in the transverse plane $(x, y)$. The charge distribution $f_0$ is also perturbed with the same amplitude and wave vector.

We consider a similar initial condition for the ion spin distribution $S(x, t)$

$$ 
\begin{aligned}
S_z(x, t = 0) &= \sqrt{1 - \varepsilon^2}, \\
S_x(x, t = 0) &= \varepsilon \cos(kx), \\
S_y(x, t = 0) &= \varepsilon \sin(kx),
\end{aligned}
$$

whereby the ions are fully magnetized along the z-axis, while their spin direction is slightly perturbed in the transverse plane. We emphasize that the perturbation preserves the norm of the spin vector, with $\|S(x, t)\| = 1$ at each spatial point.

Although one single mode is excited in the initial condition, many more Fourier modes ($k_n = nk_{min}$, with $n = 1, 2, \dots$) arise at later times in an unstable configuration. The couplings between these modes are responsible for the nonlinear saturation observed, for instance, in figure 4. Figure 5 shows the evolution of the various modes $k_n$ as the instability develops. At $t = 0$, only the mode $k_1 = 0.264 \text{ nm}^{-1}$ is present, but this mode is stable ($\omega_i = 0$) because its phase velocity exceeds $v_F$. However, higher order modes $k_n$ may be unstable, as is evident from figure 2(d). These unstable modes, although not present in the initial condition, are excited via mode couplings, and start growing exponentially around $t = 50$ fs. Note that the growth rate decreases with $k$ [see figure 2(d)], so that higher-order modes become significant at later times, as is apparent from figure 5. After $t \approx 200$ fs, a wide spectrum of modes is excited and the system is clearly in a strongly nonlinear regime.

Figures 6 and 7 show respectively the charge distribution $f_0(x, v, t)$ and the spin distribution $f_x(x, v, t)$ at different times, for an unstable case with $\eta = -0.2$ and $Kn_i/(\hbar\omega_p) = 0.08$. The fundamental mode $k = 2\pi/L = 1.59 \text{ nm}^{-1}$ is excited at $t = 0$ and grows exponentially due to the wave-particle interactions near the resonant velocity $v_{res} = -0.58v_F$. After $t = 66$ fs – see figures 6(c) and 7(c) – nonlinear effects become dominant and the distribution function develops phase-space vortices, as is usually observed in plasma instabilities.

Figure 8 shows the transfer of energy from the electron-
<figure><img src="image_7.png" /><figcaption>FIG. 6. Electron phase-space distribution function $f_0(x, v, t)$ at different times, $t = 0$ fs (a), $t = 37$ fs (b), $t = 66$ fs (c), $t = 220$ fs (d), $t = 242$ fs (e), and $t = 271$ fs (f). Other parameters are: $L = 3.96$ nm, $k = 2\pi/L = 1.59 \text{ nm}^{-1}$, $Kn_i/(\hbar\omega_p) = 0.08$, $\eta = -0.2$, and perturbation $\epsilon = 10^{-4}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Same as figure 6 for the electron spin distribution $f_x(x, v, t)$. The values of $f_x$ in the first panel (a) have been multiplied by a factor 200 for better readability.</figcaption></figure>

ion magnetic exchange energy $E_{mag} = -\frac{Kn_i}{2} \iint \mathbf{f}(x, v, t) \cdot \mathbf{S}(x, t)dx dv$ to the electrons' kinetic energy $E_{kin} = \frac{m_e}{2} \iint v^2 f_0(x, v, t) dx dv$. This transfer of energy, due to the wave-particle interactions detailed in the preceding sections, is particularly evident during the nonlinear phase of the evolution, highlighted in green on the figure.

Finally, figure 9 shows the loss of magnetization that is induced by the wave-particle interactions. Initially, the ions are fully magnetized along the z-directions. At later times, their average magnetization over the entire domain is given by $\langle S_z \rangle(t) = \int S_z(x, t) dx/L$. Following the instability, part of their magnetization is transferred
<figure><img src="image_9.png" /><figcaption>FIG. 8. Time evolution of the electron kinetic energy $E_{kin} = (m_e/2) \iint v^2 f_0 dxdv$ (orange curve) and electron-ion magnetic energy $E_{mag} = -(Kn_i/2) \iint \mathbf{f} \cdot \mathbf{S} dxdv$ (blue curve), for an initial electron spin polarization $\eta = -0.6$ and excitation wave number $k = 2.12 \text{ nm}^{-1}$. The energies are expressed in arbitrary units. The curves are shifted to improve readability: the kinetic energy is set to zero at $t = 0$, while the magnetic energy is set to 0.5 at the final time. The inset shows the amplitude of the first Fourier mode of the ion magnetization $S_x(x, t)$ (i.e., the magnon amplitude), on a semi-logarithmic scale. The nonlinear regime is highlighted by the shaded green background.</figcaption></figure>

to the electrons, whose average spin density is $\langle f_z \rangle(t) = \frac{Z}{2} \iint f_z(x, v, t) dxdv / \iint f_0(x, v, t) dxdv$. Such transfer of magnetization is evident from the blue curves on figure 9, which represents a case with $\eta = 0$, i.e., electrons that are initially unpolarized. When the electrons are anti-polarized—i.e., polarized in the direction opposite to that of the ions (orange and green curves)—the loss of ion magnetization is even larger, and approaches 30% for $\eta = -0.4$.

We note that such loss of magnetization occurs on a few tens of femtoseconds, in accordance with the experimental findings$^{1-3}$. Therefore, the wave-particle mechanism proposed in this work may be a key element in the still elusive explanation of ultrafast demagnetization occurring in ferromagnetic materials when irradiated with laser light.

According to figure 9, the loss of magnetization is greater the larger the difference between the initial electron polarization $\eta$ and the threshold value $\eta_{th} = \frac{3}{4} Kn_i E_F$ for which the imaginary part of the frequency (growth rate) vanishes. In order to explain this behavior, we consider a simple “quasi-linear” scenario. Due to the instability, the ion magnetization decreases, while the electrons’ increases (and possibly changes sign if the initial $\eta$ is negative). We postulate that the instability stops when the instantaneous growth rate vanishes, i.e. when, from Eq. (17),

$$ \frac{3}{4} Kn_i S_{z,f} - \eta_f = 0, \qquad (20) $$

where $S_{z,f}$ is the final value of the ion magnetization

<figure><img src="image_10.png" /><figcaption>FIG. 9. Time evolution of the ion magnetization in the z-direction $\langle S_z \rangle(t) = \int S_z(x, t) dx/L$, expressed as a fraction of the initial magnetization, for three values of the initial electron polarization: $\eta = 0$ (solid blue line), $\eta = -0.2$ (solid orange line), and $\eta = -0.4$ (solid green line). The dashed blue line represents the electrons’ average polarization $\langle f_z \rangle(t) = \frac{Z}{2} \iint f_z(x, v, t) dxdv / \iint f_0(x, v, t) dxdv$, for $\eta = 0$, upshifted by 0.9 to improve the readability of the diagram. The excited wave number is the lowest allowed by the periodic conditions: $k = k_{min} = 1.1 \text{ nm}^{-1}$. In this case, the electron-ion coupling constant $K$ is defined by $Kn_i = 0.3 \hbar \omega_p$; it differs from the value given in Table I for nickel and is closer to the one that can be observed in cobalt samples$^{29}$.</figcaption></figure>

along the z-axis, and $\eta_f$ is the final value of the electron polarization, defined by generalizing Eq. (9) to nonequilibrium situations.

The conservation of angular momentum along $z$-see Eq. (8) -implies that

$$ S_{z,f} + \frac{1}{2}\eta_f = 1 + \frac{1}{2}\eta. \qquad (21) $$

By combining Eqs. (20) and (21), we obtain the final magnetization as a function of the initial electron polarization $\eta$

$$ S_{z,f} = \frac{\eta + 2}{\eta_{\text{th}} + 2}, \quad \text{for } -1 \le \eta \le \eta_{\text{th}}, \qquad (22) $$

and $S_{z,f} = 1$, for $1 \ge \eta > \eta_{th}$.

Figure 10 shows the final ion magnetization as a function of $\eta$, for various wave numbers $kL_F$, where $L_F = v_F/\omega_p$. All final magnetization values lie above the curve defined by Eq. (22), which serves as a lower bound, representing the maximum possible demagnetization achievable in this scenario. Such lower bound is attained for small wave numbers $kL_F \ll 1$ (large systems), whereas for $kL_F \approx 1$ the data points are shifted upwards, although the slope of the curve remains the same as in Eq. (22). Since $L_F \approx 0.1 \text{ nm}$, which is of the same order as the interatomic spacing $a$ (see Table I), wave numbers for which $kL_F \approx 1$ are not very relevant in practice, and one can consider that the expected demagnetization is the one given by Eq. (22).
<figure><img src="image_11.png" /><figcaption>FIG. 10. Final value of the ion magnetization $S_{z,f}$ as a function of the initial electron polarization $\eta$. The colored dots represent simulation results for different values of the wave number $kL_F$, where $L_F = v_F/\omega_p$. The solid black line is the theoretical result of Eq. (22).</figcaption></figure>

As concluding remarks, one should not forget that our model is Hamiltonian and does not include any dissipative effects. Therefore, the magnetization can only be exchanged between the two species – the itinerant electrons and the localized ions – but cannot be evacuated out of the system. By incorporating an appropriate dissipative channel, such as electron-phonon couplings (which have a typical timescale of a few picoseconds), it may be possible to describe the full demagnetization process, although this would require longer and more time consuming simulations.

It is also important to emphasize that the electron-magnon interactions described here occur even in the absence of spin-orbit coupling (SOC), which is not included in our model. SOC is a relativistic effect that couples the electron spin to an internal or external electric field, which possesses a magnetic component in the reference frame of the itinerant electron. Many existing models of the ultrafast magnetization dynamics rely on SOC as a primary source of demagnetization$^{15,30}$. Our approach shows that resonant electron-magnon interactions can lead to demagnetization on a femtosecond time scale, even in the absence of SOC.

V. CONCLUSION

In this work, we have explored the complex dynamics of wave-particle interactions in a spin-polarized plasma, focusing on the coupling between magnons and plasmons within a ferromagnetic material. By extending the Vlasov-Poisson framework to include spin degrees of freedom and coupling it with the Landau-Lifshitz equation for ion magnetization, we developed a comprehensive phase-space approach that captures both itinerant and localized magnetism.

Our linear response analysis revealed the existence

of a magnon-electron resonance, a phenomenon analogous to classical plasma effects such as Landau damping and electron-cyclotron resonance heating. This resonance occurs when the phase velocity of the magnon wave matches the electron spin precession frequency, adjusted by Doppler shift. The resulting energy exchange leads to either damping or instability of the spin wave, depending on the initial electron spin polarization. Notably, the threshold polarization at which the system transitions from damping to instability corresponds to the natural self-consistent equilibrium induced by the ion-generated magnetic field.

Numerical simulations confirmed these theoretical predictions, showing that in the unstable regime, the system rapidly evolves into a nonlinear state characterized by significant magnetization loss on femtosecond timescales. This behavior closely mirrors the ultrafast demagnetization observed experimentally in ferromagnetic thin films$^{1,3}$. Furthermore, the damping regime offers a microscopic interpretation of Gilbert damping, suggesting that electron-magnon resonant interactions may contribute to this effect.

Importantly, our findings demonstrate that wave-particle resonances — long studied in plasma physics — play an important role in spin-polarized systems and may offer new insights into ultrafast magnetization dynamics. This opens avenues for future research in spintronics, magneto-optical devices, and plasma-based acceleration schemes, where controlling spin-wave interactions could lead to novel effects and enhanced performance.

ACKNOWLEDGMENTS

We thank Aurélien Manchon and Nicolas Bergeard for several helpful discussions. This work was supported by France 2030 government investment plan managed by the French National Research Agency under grant reference PEPR SPIN – [SPINTHEORY] ANR-22-EXSP-0009. This work of the Interdisciplinary Thematic Institute QMat, as part of the ITI 2021 2028 program of the University of Strasbourg, CNRS and Inserm, was supported by IdEx Unistra (ANR 10 IDEX 0002), and by SFRI STRAT'US project (ANR 20 SFRI 0012) and EUR QMAT ANR-17-EURE-0024 under the framework of the French Investments for the Future Program. The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsac calcul/Big Data.
# Appendix A: Equilibrium

## 1. Electron gas without spin

In the case of a homogeneous equilibrium, the Hartree potential vanishes, so that the Hamiltonian of a spinless electron gas is simply

$$H = \frac{m_e(v_x^2 + v_y^2 + v_z^2)}{2}. \qquad (A1)$$

The equilibrium distribution $\mathcal{F}$ follows is Fermi-Dirac distribution. At zero temperature, this implies that all energy states below the Fermi energy $E_F$ are occupied, while all states above $E_F$ are empty. The distribution is then

$$\mathcal{F}(H) = \frac{2\pi C}{m_e} \Theta(H - E_F), \qquad (A2)$$

where $\Theta$ is the Heaviside step function and $C = [m_e/(2\pi\hbar)]^3$ is a normalization constant that enforces the total number of electrons.

The 3D velocity distribution $f^{3D}(v_x, v_y, v_z)$ associated to $\mathcal{F}$ is constant inside the Fermi sphere – defined by $|v|^2 \le v_F^2$, where $v_F = \sqrt{2E_F/m_e}$ is the Fermi velocity – and zero outside:

$$f^{3D}(v_x, v_y, v_z) = \begin{cases} \pi C, & \text{if } \|v\|^2 < v_F^2 \\ 0, & \text{otherwise.} \end{cases} \qquad (A3)$$

The 1D distribution is obtained by integration $f^{1D}(v_x) = \iint f^{3D} dy dz$, which yields

$$f^{1D}(v_x) = \begin{cases} \pi C(v_F^2 - v_x^2), & \text{if } v_x^2 < v_F^2 \\ 0, & \text{otherwise.} \end{cases} \qquad (A4)$$

## 2. Polarized electron gas

When the electron gas is polarized by a magnetic field $\mathbf{B}$ (either external or internal), the spin degrees of freedom cannot be neglected. The magnetic field splits the electron distribution into two different populations, with spins either parallel or antiparallel to the field. The Hamiltonian is a diagonal $2 \times 2$ matrix $\mathcal{H} = \text{diag}(H_\uparrow, H_\downarrow)$, with

$$H_\uparrow = \frac{m_e(v_x^2 + v_y^2 + v_z^2)}{2} - \mu_B B \qquad (A5)$$

$$H_\downarrow = \frac{m_e(v_x^2 + v_y^2 + v_z^2)}{2} + \mu_B B. \qquad (A6)$$

Each population follows its own Fermi-Dirac distribution

$$\mathcal{F}_{\uparrow,\downarrow}(H_{\uparrow,\downarrow}) = \frac{2\pi C}{m_e} \Theta(H_{\uparrow,\downarrow} - \mu), \qquad (A7)$$

where the chemical potential $\mu$ is the energy needed to add one electron to the system. The chemical potential

is adjusted to correspond to the number of particles in the system, and it plays the same role as $E_F$ in the non-polarized case.

In terms of the velocity, there are now two Fermi spheres, one for each population, whose radii are $v_\uparrow$ and $v_\downarrow$. From Eq. (A7), the energies associated with these velocities are given by $\frac{m_e}{2}v_{\uparrow,\downarrow}^2 = \mu \pm \mu_B B$. By eliminating $\mu$, one obtains

$$2\mu_B B = \frac{m_e(v_\uparrow^2 - v_\downarrow^2)}{2}. \qquad (A8)$$

Integrating over $v_y$ and $v_z$, we finally obtain the 1D Fermi-Dirac distributions at zero temperature

$$f_{\uparrow,\downarrow}^{(0)}(v_x) = \begin{cases} \pi C (v_{\uparrow,\downarrow}^2 - v_x^2), & \text{if } v_x^2 < v_{\uparrow,\downarrow}^2 \\ 0, & \text{otherwise,} \end{cases} \qquad (A9)$$

which are identical to the Eqs. (11) in the main text. An example of these equilibrium distributions is shown in figure 11(a).

Finally if, in Eq. (A8), we take $B$ to be the magnetic field generated by the fully magnetized ions $\mu_B B = -Kn_i/2$, and use Eq. (14), we obtain

$$K = \frac{E_F}{n_i} \left[ (1 + \eta)^{2/3} - (1 - \eta^{2/3}) \right],$$

which is the same as Eq. (18) for $\eta = \eta_{th}$.

## 3. Low-polarization equilibrium ($\eta \ll 1$)

For low values of $\eta$, the electron charge density $f_0^{(0)} = f_\uparrow^{(0)} + f_\downarrow^{(0)}$ is almost identical to that of an unpolarized gas. The $z$-component of the spin density $f_z^{(0)} = f_\uparrow^{(0)} - f_\downarrow^{(0)}$ is given by the difference of two truncated parabolas with same curvature and boundaries $v_\uparrow \approx v_\downarrow \approx v_F$, see Eq. (A9). Therefore, it can be approximated by a constant $f_z^{(0)} = A$, within the interval $v \in [-v_F, v_F]$. The constant $A$ is determined by the normalization condition $\int_{-v_F}^{v_F} f_z^{(0)} dv = n_e \eta$, so that $A = \eta n_e / (2v_F)$. This approximation is illustrated in figure 11.

# Appendix B: Approximate magnon dispersion relation

In this Appendix, we provide some details on the derivation of the approximate dispersion relation given in Eqs. (16)-(17).

We first define the following functions, which are useful for computing the dispersion relation of a Fermi gas at zero temperature

$$W_n(\gamma) = \int_{-1}^{1} \frac{v^n}{v - \gamma} dv.$$
<figure><img src="image_12.png" /><figcaption>FIG. 11. Equilibrium distribution of an electron gas with spin, for a polarization $\eta = 0.7$. Top panel (a): Spin-up and spin-down Fermi-Dirac distributions at zero temperature. $v_↑$ and $v_↓$ are the Fermi velocities for each spin population. Bottom panel (b): The solid blue line represents the exact spin distribution in the z-direction $f_z = f_↑ - f_↓$. The dashed red line is an approximation valid for low values of $\eta$. The error (shaded green area) is localized near the Fermi velocity and comes from the difference $v_↑ - v_↓ \approx 2\eta/3$, see Eq. (14). Here, we used a large polarization $\eta = 0.7$ to better illustrate the approximation.</figcaption></figure>

These functions are the Fermi-Dirac counterpart of the Fried-Conte function $Z(\gamma)$ for a Maxwell-Boltzmann distribution$^{31}$.

In particular, the first two functions $W_0$ and $W_1$ are relevant here, as they intervene in the integrals of Eq. (15), after some normalization that can be recovered easily using the scaling property

$$ \int_{-a}^{a} \frac{v^n}{v - \gamma} dv = a^n W_n(\gamma/a). $$

Using the definition of the complex logarithmic function, one can write: $W_0(\gamma) = \ln(1 - \gamma) - \ln(-1 - \gamma)$, and $W_1(\gamma) = 2 + \gamma W_0(\gamma)$.

With the above mathematical tools and making use of the low-$\eta$ approximation illustrated in figure 11, we can

compute the integrals appearing in Eq. (15):

$$ \begin{aligned} \int \frac{\partial_v f_0^{(0)}(v)}{v - \frac{\omega - \omega_L}{k}} dv &= -4\pi C v_F \left[ 2 + \frac{\omega - \omega_L}{kv_F} W_0 \left( \frac{\omega - \omega_L}{kv_F} \right) \right], \\ \int \frac{f_z^{(0)}(v)}{v - \frac{\omega - \omega_L}{k}} dv &= \pi C \frac{4\eta}{3} v_F^2 W_0 \left( \frac{\omega - \omega_L}{kv_F} \right). \end{aligned} $$

When these expressions are inserted into Eq. (15), one gets

$$ \begin{aligned} D(\omega, k) = \omega &- \frac{a^2 J}{\hbar} k^2 - \frac{Z \omega_L \eta}{2} + \frac{4 \pi \omega_L^2}{Z n_i} C v_F \frac{\hbar}{2 m_e} \\ &+ \frac{2 \pi \omega_L^2}{Z n_i} C v_F \left( \frac{\omega - \omega_L}{k v_F} \frac{\hbar}{2 m_e} + \frac{\eta v_F}{3 k} \right) W_0 \left( \frac{\omega - \omega_L}{k v_F} \right). \end{aligned} $$

Injecting the value of the constant $C = \frac{3Zn_i}{8\pi v_F^3} = \left(\frac{m_e}{2\pi\hbar}\right)^3$ leads to

$$ \begin{aligned} D(\omega, k) = \omega &- \frac{a^2 J}{\hbar} k^2 - \frac{Z \omega_L \eta}{2} + \frac{3 \hbar \omega_L^2}{8 E_F} \\ &+ \frac{\omega_L^2}{4} \left( \frac{\omega - \omega_L}{k v_F} \frac{3 \hbar}{4 E_F} + \frac{\eta}{k v_F} \right) W_0 \left( \frac{\omega - \omega_L}{k v_F} \right) \end{aligned} $$

The difficulty now lies in simplifying the expression $W_0(\frac{\omega - \omega_L}{kv_F})$. Using the scaling

$$ \omega_p \gg kv_F \sim \omega_L \gg \omega_S, $$

which is satisfied for typical sizes $k^{-1} > 1$ nm (see also Table I), one can neglect $\omega$ compared to $\omega_L$. In that case, the imaginary part of the argument of $W_0$ vanishes and one can write

$$ W_0 \left( \frac{\omega - \omega_L}{kv_F} \right) \approx W_0 \left( \frac{-\omega_L}{kv_F} \right) = i\pi + \ln \left( \frac{1 + \frac{\omega_L}{kv_F}}{1 - \frac{\omega_L}{kv_F}} \right). $$

The dispersion function then becomes

$$ \begin{aligned} D(\omega, k) = \omega &- \frac{a^2 J}{\hbar} k^2 - \frac{Z \omega_L \eta}{2} + \frac{3 \hbar \omega_L^2}{8 E_F} \\ &+ \frac{\omega_L}{4} \frac{\omega_L}{kv_F} \left( -\frac{3}{4} \frac{\hbar \omega_L}{E_F} + \eta \right) \left( i\pi + \ln \left| \frac{1 + \frac{\omega_L}{kv_F}}{1 - \frac{\omega_L}{kv_F}} \right| \right). \end{aligned} $$

Setting $D(\omega, k) = 0$, leads to the approximate expressions for the real and imaginary parts of the magnon frequency, as given in Eqs. (16) and (17) of the main text.
<footnote n="1">¹E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, “Ultrafast spin dynamics in ferromagnetic nickel.” *Physical Review Letters* (1996), 10.1103/physrevlett.76.4250.</footnote>

<footnote n="2">²K. Krieger, J. K. Dewhurst, P. Elliott, S. Sharma, and E. K. U. Gross, “Laser-induced demagnetization at ultrashort time scales: Predictions of toddf,” *Journal of Chemical Theory and Computation* (2015), 10.1021/acs.jctc.5b00621.</footnote>

<footnote n="3">³J.-Y. Bigot, M. Vomir, and E. Beaurepaire, “Coherent ultrafast magnetism induced by femtosecond laser pulses,” *Nature Physics* (2009), 10.1038/nphys1285.</footnote>

<footnote n="4">⁴B. Koopmans, G. Malinowski, F. D. Longa, D. Steiauf, M. Fähnle, T. Roth, L. Floreano, M. Cinchetti, and M. Aeschlimann, “Explaining the paradoxical diversity of ultrafast laser-induced demagnetization,” *Nature Materials* (2010), 10.1038/nmat2593.</footnote>

<footnote n="5">⁵G. Malinowski, F. D. Longa, J. Rietjens, P. V. Paluskar, R. Huijink, R. Huijink, H. Swagten, and B. B. Koopmans, “Control of speed and efficiency of ultrafast demagnetization by direct transfer of spin angular momentum,” *Nature Physics* (2008), 10.1038/nphys1092.</footnote>

<footnote n="6">⁶D. Gupta, M. Pankratova, M. Riepp, M. Pereiro, B. Sanyal, S. Ershadrad, M. Hehn, N. Pontius, C. Schüßler-Langeheine, R. Abrudan, N. Bergeard, A. Bergman, O. Eriksson, and C. Boeglin, “Tuning ultrafast demagnetization with ultrashort spin polarized currents in multi-sublattice ferrimagnets,” *Nature Communications* (2025), 10.1038/s41467-025-58411-3.</footnote>

<footnote n="7">⁷I. Radu, K. Vahaplar, C. Stamm, T. Kachel, N. Pontius, H. A. Dürr, T. Ostler, J. Barker, R. F. L. Evans, R. W. Chantrell, A. T. A. Itoh, A. Kirilyuk, T. Rasing, and A. Kimel, “Transient ferromagnetic-like state mediating ultrafast reversal of antiferromagnetically coupled spins,” *Nature* (2011), 10.1038/nature09901.</footnote>

<footnote n="8">⁸J. Walowski and M. Muenzenberg, “Perspective: Ultrafast magnetism and thz spintronics,” *ChemInform* (2016), 10.1002/chin.201648262.</footnote>

<footnote n="9">⁹W.-T. Lu, Y. Zhao, M. Battiato, Y. Wu, Y. Wu, and Z. Yuan, “Interface reflectivity of a superdiffusive spin current in ultrafast demagnetization and terahertz emission,” *arXiv: Mesoscale and Nanoscale Physics* (2019), 10.1103/physrevb.101.014435.</footnote>

<footnote n="10">¹⁰J. Hurst, P.-A. Hervieux, and G. Manfredi, “Spin current generation by ultrafast laser pulses in ferromagnetic nickel films,” *Phys. Rev. B* 97, 014424 (2018).</footnote>

<footnote n="11">¹¹S. C. Cowley, R. M. Kulsrud, and E. Valeo, “A kinetic equation for spin-polarized plasmas,” *Physics of Fluids* (1986), 10.1063/1.865726.</footnote>

<footnote n="12">¹²R. M. Kulsrud, E. Valeo, and S. C. Cowley, “Physics of spin-polarized plasmas,” *Nuclear Fusion* (1986), 10.1088/0029-5515/26/11/001.</footnote>

<footnote n="13">¹³Y. Wu, L. Ji, X. Geng, Q. Yu, N. Wang, B. Feng, Z. Guo, W. Wang, C. Qin, X. Yan, L. Zhang, J. Thomas, A. Hützen, M. Büscher, T. P. Rakitzis, A. Pukhov, B. Shen, and R. Li, “Polarized electron-beam acceleration driven by vortex laser pulses,” *New Journal of Physics* 21, 073052 (2019).</footnote>

<footnote n="14">¹⁴Y. Wu, L. Ji, X. Geng, J. Thomas, M. Büscher, A. Pukhov, A. Hützen, L. Zhang, B. Shen, and R. Li, “Spin filter for polarized electron acceleration in plasma wakefields,” *Phys. Rev. Appl.* 13, 044064 (2020).</footnote>

<footnote n="15">¹⁵A. Manchon, Q. Li, L. Xu, and S. Zhang, “Theory of laser-induced demagnetization at high temperatures,” *Phys. Rev. B* 85, 064408 (2012).</footnote>

<footnote n="16">¹⁶E. G. Tveten, A. Brataas, and Y. Tserkovnyak, “Electron-magnon scattering in magnetic heterostructures far out of equilibrium,” *Phys. Rev. B* 92, 180412 (2015).</footnote>

<footnote n="17">¹⁷M. Weißenhofer and P. Oppeneer, “Ultrafast demagnetization through femtosecond generation of non-thermal magnons,” *Advanced Physics Research* 4 (2024), 10.1002/apxr.202300103.</footnote>

<footnote n="18">¹⁸J. Hurst, O. Morandi, G. Manfredi, and P.-A. Hervieux, “Semiclassical Vlasov and fluid models for an electron gas with spin effects,” *European Physical Journal D* (2014), 10.1140/epjd/e2014-50205-5.</footnote>

<footnote n="19">¹⁹G. Manfredi, P.-A. Hervieux, and N. Crouseilles, “Spin effects in ultrafast laser-plasma interactions,” *European Physical Journal-Special Topics* (2022), 10.1140/epjs/s11734-022-00669-5.</footnote>

<footnote n="20">²⁰N. Crouseilles, P. Hervieux, X. Hong, and G. Manfredi, “Vlasov-Maxwell equations with spin effects,” *Journal of Plasma Physics* (2023), 10.1017/s0022377823000314.</footnote>

<footnote n="21">²¹M. Lakshmanan, “The fascinating world of the Landau-Lifshitz-Gilbert equation: an overview,” *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 369, 1280-1300 (2011).</footnote>

<footnote n="22">²²M. A. Ruderman and C. Kittel, “Indirect exchange coupling of nuclear magnetic moments by conduction electrons,” *Phys. Rev.* 96, 99-102 (1954).</footnote>

<footnote n="23">²³G. Manfredi, P.-A. Hervieux, and J. Hurst, “Phase-space modeling of solid-state plasmas,” *Reviews of Modern Plasma Physics* (2019), 10.1007/s41614-019-0034-0.</footnote>

<footnote n="24">²⁴B. Bakri, N. Crouseilles, P.-A. Hervieux, X. Hong, and G. Manfredi, “Ultrafast dynamics of a spin-polarized electron plasma with magnetic ions,” *Journal of Plasma Physics* 91 (2025), 10.1017/s0022377824001594.</footnote>

<footnote n="25">²⁵G. Manfredi and P.-A. Hervieux, “Finite-size and nonlinear effects on the ultrafast electron transport in thin metal films,” *Phys. Rev. B* 72, 155421 (2005).</footnote>

<footnote n="26">²⁶A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, “Magnon spintronics,” *Nature physics* 11, 453-461 (2015).</footnote>

<footnote n="27">²⁷T. Gilbert, “A phenomenological theory of damping in ferromagnetic materials,” *IEEE Transactions on Magnetics* 40, 3443-3449 (2004).</footnote>

<footnote n="28">²⁸M. C. Hickey and J. S. Moodera, “Origin of intrinsic gilbert damping,” *Phys. Rev. Lett.* 102, 137601 (2009).</footnote>

<footnote n="29">²⁹C. Kittel and P. McEuen, *Introduction to solid state physics* (John Wiley & Sons, 2018).</footnote>

<footnote n="30">³⁰M. Stamenova, J. Simoni, and S. Sanvito, “Role of spin-orbit interaction in the ultrafast demagnetization of small iron clusters,” *Phys. Rev. B* 94, 014423 (2016).</footnote>

<footnote n="31">³¹B. D. Fried and S. D. Conte, *The plasma dispersion function: the Hilbert transform of the Gaussian* (Academic press, 2015).</footnote>