# Ultrafast dynamics of a spin-polarized electron plasma with magnetic ions

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

► To cite this version:

Benjamin Bakri, Nicolas Crouseilles, Paul-Antoine Hervieux, Xue Hong, Giovanni Manfredi. Ultrafast dynamics of a spin-polarized electron plasma with magnetic ions. Journal of Plasma Physics, 2025, 91 (1), pp.1-38. 10.1017/S0022377824001594. hal-04844507v2

HAL Id: hal-04844507

https://hal.science/hal-04844507v2

Submitted on 22 Jan 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
Ultrafast dynamics of a spin-polarized electron
plasma with magnetic ions

Benjamin Bakri<fnref n="1" />, Nicolas Crouseilles<fnref n="2" />, Paul-Antoine Hervieux ①,1, Xue Hong<fnref n="3" />
and Giovanni Manfredi ①,1,†

(Received 24 July 2024; revised 21 November 2024; accepted 22 November 2024)

We construct a mean-field model that describes the nonlinear dynamics of a spin-polarized electron gas interacting with fixed, positively charged ions possessing a magnetic moment that evolves in time. The mobile electrons are modelled by a four-component distribution function in the two-dimensional phase space (x, v), obeying a Vlasov-Poisson set of equations. The ions are modelled by a Landau-Lifshitz equation for their spin density, which contains ion-ion and electron-ion magnetic exchange terms. We perform a linear response study of the coupled Vlasov-Poisson-Landau-Lifshitz (VPLL) equations for the case of a Maxwell-Boltzmann equilibrium, focussing in particular on the spin dispersion relation. Conditions of stability or instability for the spin modes are identified, which depend essentially on the electron spin polarization rate η and the electron-ion magnetic coupling constant K. We also develop an Eulerian grid-based computational code for the fully nonlinear VPLL equations, based on the geometric Hamiltonian method first developed by Crouseilles et al. (J. Plasma Phys., vol. 89, no. 2, 2023, p. 905890215). This technique allows us to achieve great accuracy for the conserved quantities, such as the modulus of the ion spin vector and the total energy. Numerical tests in the linear regime are in accordance with the estimations of the linear response theory. For two-stream equilibria, we study the interplay of instabilities occurring in both the charge and the spin sectors. The set of parameters used in the simulations, with densities close to those of solids (≈10²⁹ m⁻<fnref n="3" />) and temperatures of the order of 10 eV, may be relevant to the warm dense matter regime appearing in some inertial fusion experiments.

**Keywords:** plasma waves, solid state plasma, plasma simulation

1. Introduction

The interaction of coherent electromagnetic radiation (laser light) with matter is a
well-established field within various branches of physics, particularly condensed matter

† Email address for correspondence: giovanni.manfredi@ipcms.unistra.fr

<footnote n="1">¹Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR7504,
F-67000 Strasbourg, France</footnote>

<footnote n="2">²Université de Rennes, Inria Rennes (Mingus team) and IRMAR UMR CNRS 6625, F-35042, Rennes, France</footnote>

<footnote n="3">³Department of Radiation Oncology, University of Kansas Medical Center, KA, USA</footnote>
and nanophysics, where laser pulses are often employed to study how electrons behave on extremely short time scales (femto- or attoseconds). Indeed, the most common electronic resonance found in metals – the plasmon resonance – occurs within the femtosecond time scale. This makes ultrafast laser pulses an essential tool for experimental investigations into the collective behaviour of electrons in metals.

In plasma physics, laser–plasma interactions are essential for the development of inertial fusion (triggered by powerful laser pulses) and laser–plasma accelerators (which rely on the acceleration of charged particles by plasma waves). They are also crucial in the study of warm dense matter (WDM), a state of matter that is at the frontier between solids and dense plasmas, where ultrafast non-equilibrium dynamics have been recently accessed thanks to subpicosecond laser pulses (Falk 2018).

However, in addition to their electric charge, electrons also possess an intrinsic magnetic moment, i.e. a spin. Using the electron spin as a vector to code and transfer information is at the core of the emerging field of spintronics. In nanophysics, spin effects are crucial to describe the ultrafast demagnetization observed in ferromagnetic thin films irradiated with femtosecond laser pulses (Beaurepaire et al. 1996; Bigot, Vomir & Beaurepaire 2009; Bigot & Vomir 2013). Despite intense investigations, such ultrafast demagnetization is not yet fully understood, although the spin–orbit interaction (Hinschberger & Hervieux 2012; Krieger et al. 2015, 2017), spin currents (Choi et al. 2014; Hurst et al. 2014; Schellekens et al. 2014) and superdiffusive electron transport (Battiato, Carva & Oppeneer 2010) appear to play a significant role.

The exploration of spin-dependent effects in plasma physics is a relatively new area of study. Nonetheless, it is now possible to generate and precisely control polarized electron beams with high spin polarization in laboratory settings (Wu et al. 2019, 2020; Nie et al. 2021). Theoretical studies on polarized plasmas have been revitalized in recent years (Zamanian, Marklund & Brodin 2010a; Zamanian et al. 2010b; Hurst et al. 2014; Morandi et al. 2014; Hurst, Hervieux & Manfredi 2017), although some early developments date back to the 1980s (Cowley, Kulsrud & Valeo 1986). Notably, Brodin, Holkundkar & Marklund (2013) have formulated a particle-in-cell (PIC) code that incorporates the magnetic dipole force and magnetization currents related to the electron spin. PIC methods for particles with spin have also been developed for applications in the field of laser–plasma interactions (Li et al. 2021).

Within the condensed matter and nanophysics communities, most research on ultrafast spin dynamics has relied on wavefunction-based methods, particularly time-dependent density functional theory, augmented to incorporate spin effects (spin-TDDFT) (Yin et al. 2009; Manfredi et al. 2010; Krieger et al. 2015; Sinha-Roy et al. 2020). Spin-TDDFT models have also been used to study spin effects in dense plasmas in the WDM regime (Bonitz et al. 2020).

In a recent series of papers (Crouseilles et al. 2021, 2023; Manfredi, Hervieux & Crouseilles 2023), we have proposed an alternative approach based on Wigner functions, which represent electronic quantum states through a pseudo-probability distribution in the classical phase space. The corresponding Wigner evolution equation reduces to the standard Vlasov equation of classical plasma physics. For spin-1/2 particles, such as electrons, one can construct a semi-classical model, where the orbital motion (i.e. the trajectories in the phase space) is treated classically while the spin is kept as a quantum mechanical variable. For a review of methods based on Wigner functions, see Manfredi, Hervieux & Hurst (2019).

Among these phase space models, two families can be distinguished: on the one side, Vlasov models that use a scalar distribution function on an extended phase-space (x, v, s), where x and v are the position and velocity of the electron, and s denotes the spin variable
(Brodin et al. 2008, 2011; Marklund, Zamanian & Brodin 2010; Zamanian et al. 2010a); on the other side, models using a multi-component distribution function $f_ν$, ($ν = 0, 3$) with values in the standard phase space $(x, v)$. These two approaches are almost, although not exactly, equivalent (see our detailed discussion in Crouseilles et al. 2023 for further clarifications). Hereafter, we will name these approaches respectively as ‘scalar’ and ‘vectorial’. Note that for both of them, the orbital motion is classical while the spin is a fully quantum variable. The numerical approximation of these models requires different techniques. Indeed, the scalar version involves an extended phase space of dimension 8, which naturally leads to consider PIC techniques as the method of choice (Crouseilles et al. 2021; Li 2023); in contrast, the vectorial approach is more easily amenable to grid-based methods (Crouseilles et al. 2023).

In previous works (Crouseilles et al. 2021, 2023; Manfredi et al. 2023), we had only considered the dynamics of the mobile (itinerant) electrons, whereas the ions only acted as an immobile neutralizing background. However, in ferromagnets, most of the magnetic properties are due to the fixed ions, which account for approximately 95 % of the magnetization of the material, whereas only the remaining ≈5 % can be attributed to the mobile electrons. In the present work, the ions are still fixed (because their orbital response occurs on much longer time scales), but their spin is allowed to evolve in time according to the Landau–Lifshitz (LL) equation. The latter describes the precession motion of a magnetic moment in an effective magnetic field, which can be either an external one or the field created by the spin of the itinerant electrons. In turns, the ions generate a magnetic field which acts on the spin of the electrons. The ions also interact among each other through a Heisenberg-type magnetic-exchange interaction, while the electrons feel the usual self-consistent electric field.

Overall, the nonlinear Vlasov–Poisson–Landau–Lifshitz (VPLL) equations describe the coupling between the itinerant magnetism generated by the mobile electrons, represented by a vector distribution function $(f_0, f')(t, x, v) ∈ R^4$, and the fixed magnetism carried by the motionless ions, represented by their local spin $S(t, x) ∈ S^2$. It can be viewed as a spin-extended version of the usual Vlasov–Poisson model with fixed ions. An earlier version of this model – employing a more rudimentary numerical technique – was used by Hurst, Hervieux & Manfredi (2018) to study spin current generation in thin nickel films. Here, we will mainly consider a parameter range relevant to WDM (Bonitz et al. 2020), with densities close to those of solids (≈10²⁹ m⁻³) and temperatures of the order of 10 eV. For these conditions, the electron plasma is weakly degenerate ($T_e ≈ T_F$), where $T_F$ is the Fermi temperature), so that its equilibrium can be characterized with relatively good accuracy by a Maxwell–Boltzmann distribution. The ions are fixed and non-degenerate.

The model is described mathematically by a set of coupled nonlinear partial differential equations (PDEs). The design of an efficient scheme for a system of PDEs is not easy and one possible strategy is to make use of a splitting algorithm. When the system under consideration enjoys a Hamiltonian structure, a systematic way to proceed relies on the Hamiltonian splitting (Crouseilles, Einkemmer & Faou 2015; Casas et al. 2017; Li, Sun & Crouseilles 2020; Crestetto et al. 2022). It turns out that the VPLL equations enjoy a Poisson structure which motivates the use of Hamiltonian time splitting. Following previous development of geometric numerical method for Vlasov-type equations (Crouseilles et al. 2015; Li et al. 2020; Crestetto et al. 2022), the Hamiltonian splitting applied to the VPLL leads to five subsystems that can be solved exactly in time, and for which efficient and high-order methods in space and velocity can be used. As a consequence, the time accuracy of the resulting scheme only depends on the splitting error (which can be made arbitrarily small using high-order composition splittings Yoshida
1990; Hairer, Lubich & Wanner 2006) and since the method is symplectic (as composition of symplectic flows), it maintains long-term accuracy on invariants such as the total energy (Hairer et al. 2006). Another interesting property that can be proven for the proposed scheme is the exact preservation of the norm of the ion spin $||S||$.

To validate the numerical results, we investigate the linearized VPLL system by deriving the pertinent dispersion relation, following Manfredi et al. (2019). When the ion-electron coupling is turned off, the dispersion relation degenerates into the standard Bohm-Gross relation for plasmons and the magnon dispersion relation for the ion spins (Eich, Pittalis & Vignale 2018). It is noteworthy that the typical plasmon time scale is approximately two orders of magnitude faster than that of magnons, which constitutes a considerable challenge for the numerical scheme. In the case of Maxwell-Boltzmann equilibria, the dispersion relations can be solved numerically using dedicated libraries, e.g. see Kravanja & Van Barel (2000). Moreover, analytical calculations are performed in the weak coupling regime. Cross-validations between the roots of the dispersion relation and the results of the nonlinear code are performed and discussed.

The rest of the paper is organized as follows. Section 2 lays the basis of the VPLL model equations and their non-dimensional form. Section 3 discusses the linear response theory and the corresponding dispersion relation. The numerical method is presented in § 4. Results of numerical simulations are presented in § 5, both for a stable Maxwell-Boltzmann equilibrium and an unstable two-stream distribution function, and compared with linear-response results obtained from the dispersion relation, particularly for damping and growth rates. Conclusions are drawn in § 6. Three appendices provide some further details on the Maxwell-Boltzmann equilibrium with spin (Appendix A), the dispersion relation (Appendix B) and the numerical splitting technique (Appendix C).

## 2. Vlasov-Poisson-Landau-Lifshitz model

We consider a generic scenario where a magnetic material (e.g. nickel) is irradiated with a strong femtosecond laser pulse, so that some or most of the electrons are extracted from the bulk and can move freely, thus constituting a mobile electron plasma. The pulse heats up the electrons to a temperature equivalent to their Fermi energy, which for nickel is $E_F \approx 10$ eV, while their density remains similar to that of the solid $n_e \approx 10^{29} \text{ m}^{-3}$. These parameters are close to those of the weakly degenerate plasmas typical of WDM (Falk 2018; Bonitz et al. 2020). During these initial instants, up to approximately 100 fs, the ions do not have time to move, and can thus be assimilated to an immobile, but magnetized, background.

Within this broad context, our purpose here is to validate our numerical code, in the linear and nonlinear regimes, for parameters that are similar to those mentioned above. Hence, we will consider a one-dimensional (1-D) model with periodic boundary conditions, and will investigate how a perturbed Maxwell-Boltzmann equilibrium evolves in time, for both the charge (plasmons) and spin (magnons) sectors. We will also analyse potentially unstable two-stream equilibria.

### 2.1. *Model equations*

The electrons are described by a four-component distribution function ($f_0, f_1(t, x, v)$) with $f = (f_1, f_2, f_3) \in \mathbb{R}^3$, which is coupled to the continuous ion spin distribution $S(t, x) = (S_1, S_2, S_3)(t, x)$. The overall system of equations, for the space variable $x \in [0, L] \subset \mathbb{R}$ and velocity variable $v \in \mathbb{R}$, is composed of a set of kinetic equations for the electron
distribution functions (Manfredi et al. 2019; Crouseilles et al. 2023),

$$ \frac{\partial f_0}{\partial t} + v \frac{\partial f_0}{\partial x} + \frac{e}{m} \frac{\partial V_H}{\partial x} \frac{\partial f_0}{\partial v} - \frac{\mu_B}{m} \frac{\partial \mathbf{B}}{\partial x} \cdot \frac{\partial \mathbf{f}}{\partial v} = 0, \qquad (2.1) $$

$$ \frac{\partial f_{\ell}}{\partial t} + v \frac{\partial f_{\ell}}{\partial x} + \frac{e}{m} \frac{\partial V_H}{\partial x} \frac{\partial f_{\ell}}{\partial v} - \frac{\mu_B}{m} \frac{\partial B_{\ell}}{\partial x} \frac{\partial f_0}{\partial v} - \frac{e}{m} (\mathbf{B} \times \mathbf{f})_{\ell} = 0, \quad \ell = 1, 2, 3, \quad (2.2) $$

and Landau–Lifshitz equation (Lakshmanan 2011) for the ion spins,

$$ \frac{\partial \mathbf{S}}{\partial t} = \frac{a^2 J}{\hbar} \left( \mathbf{S} \times \frac{\partial^2 \mathbf{S}}{\partial x^2} \right) + \frac{K}{2\hbar} \mathbf{S} \times \int \mathbf{f} \, dv, \qquad (2.3) $$

where the first term on the right-hand side is the Heisenberg ion–ion magnetic exchange, whereas the second term represents the ion–electron magnetic exchange.

The scalar distribution function $f_0(t, x, v)$ represents, as usual, the probability to find an electron in the phase space volume located around $(x, v)$, at time $t$. Its moments yield the usual macroscopic quantities, such as the density $n_e(t, x) = \int f_0(t, x, v) \, dv$. In contrast, the vector distribution function $f_\ell(t, x, v)$ represents the mean spin polarization density of the electrons in the phase space volume located around $(x, v)$ at time $t$, along the $\ell$ direction. Its first moment $M(t, x) = \int f(t, x, v) \, dv$ represents the electron spin density. For more details, see the recent review of Manfredi et al. (2019). The relationship between this $(f_0, f)$ representation and the more standard representation as a $2 \times 2$ matrix with spin-up and spin-down components is also illustrated in Appendix A.

The self-consistent electric potential (Hartree potential) $V_H(t, x)$ obeys the Poisson equation

$$ \varepsilon_0 \partial_x^2 V_H = e \int f_0 \, dv - Z e n_{\text{ion}}, \qquad (2.4) $$

and the magnetic field appearing in (2.1) and (2.2) is primarily the one created by the ions

$$ \mathbf{B}(t, x) = -\frac{Kn_{\text{ion}}\mathbf{S}(t, x)}{2\mu_B}, \qquad (2.5) $$

although external fields could also be considered. Here, $e > 0$ denotes the electron charge, $\hbar$ the Planck constant, $m$ the electron mass, $\varepsilon_0$ the permittivity of vacuum, $\mu_B = e\hbar/2m$ the Bohr magneton, $a$ the interatomic distance, $Z$ is the atomic number, $J$ and $K$ are respectively the ion–ion and electron–ion magnetic exchange constants, and $n_{\text{ion}}$ is the fixed, homogeneous ion density. The full initial condition may be denoted as $(f_0, f, V_H, S)(t=0) = (f_0^{(0)}, f^{(0)}, V_H^{(0)}, S^{(0)})$, where $\varepsilon_0 \partial_x^2 V_H^{(0)} = e \int f_0^{(0)} \, dv - Z e n_{\text{ion}}$.

Note how the $K$-terms couple the ion and electron spin dynamics: the magnetic field $\mathbf{B}$ given by (2.5) created by the ions acts on the spin part of the electron distribution functions $f$ in (2.1) and (2.2), while the electron spin density $\int f \, dv$ acts on the LL equation (2.3) for the ion spins. A schematic view of the physical system under consideration is shown in figure 1.

From a mathematical viewpoint, the model (2.1)–(2.4) enjoys a Poisson structure with the following Hamiltonian functional:

$$ \mathcal{H} = \frac{m}{2} \int v^2 f_0 \, dx \, dv + \mu_B \int f \cdot \mathbf{B} \, dx \, dv + \frac{\varepsilon_0}{2} \int (\partial_x V_H)^2 \, dx \, dv + \frac{a^2 J}{2} \int n_{\text{ion}} \sum_{\ell=1}^{3} \left( \frac{\partial S_{\ell}}{\partial x} \right)^2 \, dx. \quad (2.6) $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Schematic view of the physical system under consideration. The immobile ions (red circles) provide the main source of localized magnetism. They interact through magnetic exchange both with themselves (coupling constant J) and with the itinerant electrons, represented by green dots (coupling constant K).</figcaption></figure>

Moreover, it is possible to construct a Poisson bracket for two functionals $\mathcal{F}$ and $\mathcal{G}$ as

$$ \{\mathcal{F}, \mathcal{G}\} = \sum_{i=0}^{3} \int \frac{f_0}{m} \left[ \frac{\delta \mathcal{F}}{\delta f_i}, \frac{\delta \mathcal{G}}{\delta f_i} \right]_{xv} dx \, dv \\ + \sum_{i=1}^{3} \left( \int \frac{f_i}{m} \left[ \frac{\delta \mathcal{F}}{\delta f_0}, \frac{\delta \mathcal{G}}{\delta f_i} \right]_{xv} dx \, dv + \int \frac{f_i}{m} \left[ \frac{\delta \mathcal{F}}{\delta f_i}, \frac{\delta \mathcal{G}}{\delta f_0} \right]_{xv} dx \, dv \right) \\ + \frac{e}{\mu_B m} \int \mathcal{F} \cdot \left( \frac{\delta \mathcal{F}}{\delta \mathcal{F}} \times \frac{\delta \mathcal{G}}{\delta \mathcal{F}} \right) dx \, dv + \frac{1}{\hbar} \int \frac{\mathcal{S}}{n_{\text{ion}}} \cdot \left( \frac{\delta \mathcal{F}}{\delta \mathcal{S}} \times \frac{\delta \mathcal{G}}{\delta \mathcal{S}} \right) dx. \quad (2.7) $$

Remark 2.1. It is easy to check that the bracket (2.7) is bilinear, skew-symmetric and satisfies Leibniz’s rule, but it is not clear whether Jacobi’s identity is satisfied. Hence, this bracket is not strictly speaking a Poisson bracket; nevertheless, we will still refer to it as a Poisson bracket for the sake of simplicity.

With this notation in hand, the system (2.1)–(2.4) can be reformulated, after introducing the vector of unknowns $\mathcal{Z} = (f_0, \mathcal{F}, \mathcal{S}) \in \mathbb{R}^7$, as

$$ \frac{\partial \mathcal{Z}}{\partial t} = \{\mathcal{Z}, \mathcal{H}\}. \qquad (2.8) $$

## 2.2. Normalized dimensionless equations

We rewrite the above (2.1)–(2.4) using dimensionless units that correspond to normalizing time to the inverse of the plasmon frequency $\omega_p = \sqrt{e^2 n_e / \epsilon_0 m}$, velocities to the thermal speed $v_{\text{th}} = \sqrt{k_B T_e / m}$ and space to the Debye length $\lambda_D = v_{\text{th}} / \omega_p$, where $k_B$ is the Boltzmann constant. Hence, the electric potential is normalized to $mv_{\text{th}}^2 / e$, the electric field to $mv_{\text{th}} \omega_p / e$ and the magnetic field to $m\omega_p / e$.
Using these normalized units and defining the self-consistent electric field as $E_x = -\partial_x V_H$, the dimensionless kinetic equations read as (for simplicity of notation, we do not change the names of the dimensionless variables)

$$ \frac{\partial f_0}{\partial t} + v \frac{\partial f_0}{\partial x} - E_x \frac{\partial f_0}{\partial v} - H \frac{\partial \mathbf{B}}{\partial x} \cdot \frac{\partial \mathbf{f}}{\partial v} = 0, \qquad (2.9) $$

$$ \frac{\partial f_{\ell}}{\partial t} + v \frac{\partial f_{\ell}}{\partial x} - E_x \frac{\partial f_{\ell}}{\partial v} - H \frac{\partial B_{\ell}}{\partial x} \frac{\partial f_0}{\partial v} - (\mathbf{B} \times \mathbf{f})_{\ell} = 0, \quad \ell = 1, 2, 3, \qquad (2.10) $$

where

$$ \mathbf{B} = -\frac{\tilde{K}\mathbf{S}}{2} \qquad (2.11) $$

is the magnetic field created by the ions.

The dimensionless Planck constant,

$$ H = \frac{\hbar\omega_p}{2mv_{\text{th}}^2}, \qquad (2.12) $$

quantifies the relative importance of quantum effects with respect to thermal effects. We also note that $H$ can be written in terms of the quantum coupling parameter $\Gamma_q = \hbar\omega_p/E_F$ and the degeneracy parameter $\Theta = T_e/T_F$ as $H = \Gamma_q/(2\Theta)$. In turn, the quantum coupling parameter is related to the Wigner-Seitz radius $r_s$ through the following relationship (Bonitz et al. 2020):

$$ \Gamma_q^2 = \frac{r_s}{a_0} \frac{16}{9\pi} \left(\frac{12}{\pi}\right)^{1/3} \approx 0.88 \frac{r_s}{a_0}, \qquad (2.13) $$

where $a_0 = 4\pi\varepsilon_0\hbar^2/(me^2)$ is the Bohr radius.

The normalized LL equation becomes

$$ \frac{\partial \mathbf{S}}{\partial t} = A \left( \mathbf{S} \times \frac{\partial^2 \mathbf{S}}{\partial x^2} \right) + Z \frac{\tilde{K}}{4} \mathbf{S} \times \int \mathbf{f} \, dv, \qquad (2.14) $$

with the dimensionless magnetic exchange constants written as $A = (a^2/\lambda_D^2)(J/\hbar\omega_p)$ and $\tilde{K} = 2Kn_{\text{ion}}/\hbar\omega_p$. Finally, the dimensionless Poisson equation is

$$ -\partial_x E_x = \int f_0 \, dv - 1. \qquad (2.15) $$

The total energy in dimensionless units is given by the Hamiltonian $\mathcal{H} = \mathcal{H}_v + \mathcal{H}_E + \sum_{i=1}^3 (\mathcal{H}_{Z,i} + \mathcal{H}_{\text{spin},i})$, with

$$ \left. \begin{array}{l} \mathcal{H}_v = \frac{1}{2} \int v^2 f_0 \, dx \, dv, \qquad \mathcal{H}_E = \frac{1}{2} \int \left( \frac{\partial V_H}{\partial x} \right)^2 \, dx, \\ \\ \mathcal{H}_{Z,i} = H \int f_i B_i \, dx \, dv, \qquad \mathcal{H}_{\text{spin},i} = AH \int \left( \frac{\partial S_i}{\partial x} \right)^2 \, dx, \end{array} \right\} \qquad (2.16) $$

where the various terms correspond to the kinetic energy ($\mathcal{H}_v$), the Hartree electric energy ($\mathcal{H}_E$), the magnetic Zeeman energy ($\mathcal{H}_Z$) and the spin energy ($\mathcal{H}_{\text{spin}}$).

We consider an electron plasma in the WDM regime, with density $n_e = n_{\text{ion}} = 9.17 \times 10^{28} \, \text{m}^{-3}$ ($Z = 1$) and temperature $k_B T_e = 16.58 \, \text{eV}$. This choice yields for the
time, velocity and length scales: $\omega_p = 1.71 \times 10^{16}$ s$^{-1}$, $v_{th} = 1.71 \times 10^6$ m s$^{-1}$ and $\lambda_D = 10^{-10}$ m. As for the dimensionless parameters, we find: normalized Planck constant $H = 0.339$, quantum coupling parameter $\Gamma_q = 1.516$, Wigner-Seitz radius $r_s/a_0 = 2.60$ (corresponding to nickel) and degeneracy parameter $\Theta = 2.24$.

For the magnetic exchange coupling constants, we use values close to those of nickel (Hurst et al. 2018): $J = 0.022$ eV and $K = 0.01$ eV nm$^3$. Taking the lattice spacing $a = 2r_s = 0.275$ nm, this yields for the dimensionless parameters, $A = 0.0148$ and $\tilde{K} = 0.161$.

### 3. Linear analysis and dispersion relations

#### 3.1. Linear analysis for a generic equilibrium

To validate the model (2.9)–(2.15) in the linear response regime, we perform a linear analysis to derive the pertinent dispersion relation. First, we start with the following homogeneous stationary state:

$$ \left. \begin{array}{l} f_0^{(0)} = f_0^{(0)}(v), \quad f_3^{(0)} = f_3^{(0)}(v), \\ f_1^{(0)} = f_2^{(0)} = 0, \quad S_1^{(0)} = S_2^{(0)} = E_x^{(0)} = 0 \quad \text{and} \quad S_3^{(0)} = 1, \end{array} \right\} \qquad (3.1) $$

where the superscript ‘(0)’ stands for equilibrium. This corresponds to an ion system that is fully polarized in the $S_3$ direction, and an electron system that is partially polarized in the same direction. The degree of electron spin polarization depends on the choice of $f_3^{(0)}(v)$ and can be characterized by a single number $\eta = \int_{-\infty}^{\infty} f_3^{(0)}(v) dv$, with $\eta \in [-1, 1]$.

We then derive the linearized system and study the propagation of a perturbation around the stationary state. We thus consider solutions in the following form:

$$ f_0 = f_0^{(0)} + f_0^{(1)}, \quad f_\ell = f_\ell^{(0)} + f_\ell^{(1)}, \quad E_x = E_x^{(0)} + E_x^{(1)} \quad \text{and} \quad S_\ell = S_\ell^{(0)} + S_\ell^{(1)}. \quad (3.2a-d) $$

Inserting these solutions into the system (2.9)–(2.15) and neglecting quadratic terms leads to the following linear system:

$$ \frac{\partial f_0^{(1)}}{\partial t} + v \frac{\partial f_0^{(1)}}{\partial x} - E_x^{(1)} \frac{\partial f_0^{(0)}}{\partial v} + \frac{H \tilde{K}}{2} \frac{\partial S_3^{(1)}}{\partial x} \frac{\partial f_3^{(0)}}{\partial v} = 0, \qquad (3.3) $$

$$ \frac{\partial f_1^{(1)}}{\partial t} + v \frac{\partial f_1^{(1)}}{\partial x} + \frac{H \tilde{K}}{2} \frac{\partial S_1^{(1)}}{\partial x} \frac{\partial f_0^{(0)}}{\partial v} - \frac{\tilde{K}}{2} (f_2^{(1)} - f_3^{(0)} S_2^{(1)}) = 0, \qquad (3.4) $$

$$ \frac{\partial f_2^{(1)}}{\partial t} + v \frac{\partial f_2^{(1)}}{\partial x} + \frac{H \tilde{K}}{2} \frac{\partial S_2^{(1)}}{\partial x} \frac{\partial f_0^{(0)}}{\partial v} + \frac{\tilde{K}}{2} (f_1^{(1)} - f_3^{(0)} S_1^{(1)}) = 0, \qquad (3.5) $$

$$ \frac{\partial f_3^{(1)}}{\partial t} + v \frac{\partial f_3^{(1)}}{\partial x} - E_x^{(1)} \frac{\partial f_3^{(0)}}{\partial v} + \frac{H \tilde{K}}{2} \frac{\partial S_3^{(1)}}{\partial x} \frac{\partial f_0^{(0)}}{\partial v} = 0, \qquad (3.6) $$

$$ \frac{\partial S_1^{(1)}}{\partial t} = -A \frac{\partial^2 S_2^{(1)}}{\partial x^2} + \frac{\tilde{K}}{4} \left( S_2^{(1)} \int f_3^{(0)} dv - \int f_2^{(1)} dv \right), \qquad (3.7) $$

$$ \frac{\partial S_2^{(1)}}{\partial t} = A \frac{\partial^2 S_1^{(1)}}{\partial x^2} - \frac{\tilde{K}}{4} \left( S_1^{(1)} \int f_3^{(0)} dv - \int f_1^{(1)} dv \right), \qquad (3.8) $$

$$ \frac{\partial S_3^{(1)}}{\partial t} = 0, \qquad (3.9) $$

$$ -\partial_x E_x^{(1)} = \int f_0^{(1)} dv. \qquad (3.10) $$
By performing Fourier (in space) and Laplace (in time) transforms of the above linear system of equations, we can derive an equation relating the frequency $\omega$ and the wavenumber $k$ (we shall further refer to $\omega_e$ for the charge branch of the dispersion relation and $\omega_s$ for the spin branch). Since $S_3$ does not depend on time, the dispersion relation for $f_0^{(1)}$ and $E_x^{(1)}$ is the same as the standard Bohm-Gross relation for unpolarized electrons, that is,

$$D_e(\omega_e, k) = 1 - \frac{1}{k} \int \frac{\partial_v f_0^{(0)}}{kv - \omega_e} dv = 0, \quad (3.11)$$

(here and in the following, velocity integrals are understood as being from $-\infty$ to $+\infty$). Hence, at the level of the linear response, the spin and charge motions are completely separated. This is an important fact, as it means that an excitation (e.g. a laser pulse) acting only on the charge density will not trigger any response in the spin dynamics. To generate spin dynamics, one needs either a strong pulse that generates nonlinear effects, or an excitation that acts directly on the spins (e.g. via the magnetic part of the laser pulse).

Next, we consider the equations for $f_1^{(1)}$, $f_2^{(1)}$, $S_1^{(1)}$ and $S_2^{(1)}$, which lead to the dispersion relation for the ion spin motion:

$$D_S(\omega_s, k) = - \left[ \omega_s - \frac{\tilde{K}^2}{8} \left( \frac{\tilde{K}Hk}{2} I_0 + I_1 \right) \right]^2 + \left[ Ak^2 + \frac{\tilde{K}\eta}{4} - \frac{\tilde{K}^2}{8} \left( HkI_3 + \frac{\tilde{K}}{2} I_2 \right) \right]^2 = 0, \quad (3.12)$$

where we have defined the integrals

$$\left. \begin{aligned} I_0 &= \int \frac{\partial_v f_0^{(0)}}{(\tilde{K}/2)^2 - (vk - \omega_s)^2} dv, & I_1 &= \int \frac{(vk - \omega_s) f_3^{(0)}}{(\tilde{K}/2)^2 - (vk - \omega_s)^2} dv, \\ I_2 &= \int \frac{f_3^{(0)}}{(\tilde{K}/2)^2 - (vk - \omega_s)^2} dv, & I_3 &= \int \frac{(vk - \omega_s) \partial_v f_0^{(0)}}{(\tilde{K}/2)^2 - (vk - \omega_s)^2} dv. \end{aligned} \right\} \quad (3.13)$$

Note that when one neglects the electron-ion coupling, i.e. $\tilde{K} = 0$, the spin branch of the dispersion relation reduces to $\omega_s = \pm Ak^2$, which is the standard magnon dispersion relation (Ashcroft & Mermin 1976). In contrast, the dispersion relation for the electrons yields, from (3.11), $\omega_e \approx \omega_p$. Taking the ratio of the magnon and plasmon frequencies yields

$$\frac{\omega_s}{\omega_e} = \frac{Ak^2}{\omega_p} \approx 8.6 \times 10^{-3}, \quad (3.14)$$

where we used the parameters given in § 2.2, i.e. $\omega_p = 1.71 \times 10^{16} \text{s}^{-1}$ and $A = 0.0148$, and considered a typical length $k^{-1} = 10 \text{ nm}$. This indicates that the time scale of magnons is approximately two orders of magnitude slower than that of plasmons. This fact has an obvious impact on the numerical simulations, as many hundreds of plasmon cycles have to be resolved before one can observe a sizeable response in the ion spins.

## 3.2. Maxwell-Boltzmann equilibrium

Now, we assume the stationary states $f_0^{(0)}$, $f_\ell^{(0)}$ to be Gaussian functions, so that $I_0$, $I_1$, $I_2$, $I_3$ can be expressed using the Fried-Comte function (Fried & Conte 1961) $Z(z) = (1/\sqrt{\pi}) \int_{\mathbb{R}} (e^{-z^2}/(t-z)) dt$, $z \in \mathbb{C}$, which can itself be expressed using the erfi function erfi(z) = $(2/\sqrt{\pi}) \int_0^z e^{t^2} dt$, $z \in \mathbb{C}$ and is tabulated in several scientific libraries.
Let consider that the following homogeneous equilibrium:

$$f_0^{(0)}(v) = \frac{1}{\sqrt{\pi}} e^{-v^2}, \quad f_3^{(0)}(v) = \eta \frac{1}{\sqrt{\pi}} e^{-v^2}, \qquad (3.15a,b)$$

where $\eta = \int f_3^{(0)} dv$ is the spin polarization rate of the electrons (see Appendix A for further details). The dispersion function $D_e$ for the charge dynamics becomes

$$D_e(\omega_e, k) = 1 + \frac{2}{k^2} \left[ 1 + \frac{\omega_e}{k} Z\left(\frac{\omega_e}{k}\right) \right], \qquad (3.16)$$

while the spin dispersion function $D_S$ is

$$D_S(\omega_s, k) = - \left\{ \omega_s + \frac{c_0}{k} \left[ Z\left(\frac{\omega_s + \tilde{K}/2}{k}\right) + Z\left(\frac{\omega_s - \tilde{K}/2}{k}\right) \right] \right. \\ -c_1 \left[ Z'\left(\frac{\omega_s - \tilde{K}/2}{k}\right) - Z'\left(\frac{\omega_s + \tilde{K}/2}{k}\right) \right] \\ + \left\{ Ak^2 + d + \frac{c_0}{k} \left[ Z\left(\frac{\omega_s + \tilde{K}/2}{k}\right) - Z\left(\frac{\omega_s - \tilde{K}/2}{k}\right) \right] \right. \\ \left. +c_1 \left[ Z'\left(\frac{\omega_s - \tilde{K}/2}{k}\right) + Z'\left(\frac{\omega_s + \tilde{K}/2}{k}\right) \right] \right\}^2, \qquad (3.17)$$

with $c_0 = \tilde{K}^2\eta/16$, $c_1 = \tilde{K}^2H/16$ and $d = \tilde{K}\eta/4$. Moreover, the complex-valued function $Z$ and its derivative are given by

$$Z(z) = \sqrt{\pi} \exp(-z^2)(i - \operatorname{erfi}(z)), \quad Z'(z) = -2(zZ(z) + 1). \qquad (3.18a,b)$$

### 3.3. Analysis and computation of the spin dispersion relation

In this section, we will use another form of the dispersion function which is strictly equivalent to $D_S$ given by (3.17). Here, $D_S$ can be written as the product of two different functions $D_S = D_-D_+$ (see Appendix B.1 for further details), each of which generates the same solutions, up to a sign. In the following, we consider the function that gives rise to positive real frequencies in the limiting case $\tilde{K} = 0$, i.e.

$$D_-(\omega_s, k) = \omega_s - Ak^2 - \frac{\tilde{K}}{4} \int f_3^{(0)} dv + \frac{\tilde{K}^2}{8k} \int \frac{f_3^{(0)}}{v - \left(\frac{\omega_s - \tilde{K}/2}{k}\right)} dv \\ - \frac{H\tilde{K}^2}{8} \int \frac{\partial_v f_0^{(0)}}{v - \left(\frac{\omega_s - \tilde{K}/2}{k}\right)} dv, \qquad (3.19)$$

or, in terms of the plasma dispersion function $Z$,

$$D_-(\omega_s, k) = \omega_s - Ak^2 - \frac{\tilde{K}\eta}{4} + \frac{\tilde{K}^2\eta}{8k} Z\left(\frac{\omega_s - \tilde{K}/2}{k}\right) - \frac{H\tilde{K}^2}{8} Z'\left(\frac{\omega_s - \tilde{K}/2}{k}\right). \qquad (3.20)$$
This formulation highlights the different contributions to the magnon frequency. Let us spell out each term of the right-hand side of (3.19).

* The first two terms yield the standard dispersion relation for magnons, $\omega_s = Ak^2$.

* The next term shifts the magnon frequency due to ion precession around the magnetic field generated by electronic spins at steady state.

* The last two terms introduce corrections that are brought over by electrons that possess specific (resonant) velocities, either in their spin distribution $f_3^{(0)}$ or their charge distribution $f_0^{(0)}$ at equilibrium. This is similar to the resonant electrons that are responsible for Landau damping in spin-less plasmas.

Equation (3.20) possesses complex solutions in $\omega_s$, due to the complex-valued function Z. Physically, this means that some resonances occur in the electron population when the velocity is equal to (restoring physical dimensions for clarity) $v = \omega_s/k - \omega_L/k$, where $\omega_L = eB/m = 2\mu_B B/\hbar$ is the Larmor frequency of an electron spin in the magnetic field created by the (fully polarized) ions, $B = Kn_{ion}/(2\mu_B)$. Thus, $\omega_s/k \equiv v_s$ is the phase velocity of the ion spin wave (the magnon), whereas $\omega_L/k \equiv v_L$ is the phase velocity of the electronic spin wave propagating in the magnetized environment created by the polarized ions. The resonance occurs when the electron spin precesses at the same frequency as the magnon, shifted by Doppler effect due to the electron velocity with respect to the fixed ions. In terms of the phase velocities, this can be written as $v_s - v = v_L$.

This resonance behaves similarly to the electron cyclotron resonance heating (ECRH) effect in fusion plasmas, with two major differences. First, the ion spin wave (magnon) plays the role of the external electromagnetic wave in ECRH; second, the magnetic moment of the electrons is not orbital as in ECRH, but instead is due to the electron's intrinsic spin.

It is useful to compute the dispersion function $D_-(\omega_s, \tilde{K})$ in terms of the coupling constant $\tilde{K}$ and the frequency $\omega_s$, for a fixed value of the wavenumber $k$. Then, the solutions of the dispersion relation can be computed along a path in the $(\omega_s, \tilde{K})$ plane, by solving the equation

$$ \partial_{\tilde{K}} D_{-} d\tilde{K} + \partial_{\omega_s} D_{-} d\omega_s = 0, \qquad (3.21) $$

starting from known solutions, for instance, the one at zero coupling $\omega_s(\tilde{K} = 0) \equiv \omega_0 = Ak^2$. Solving for $\omega_s(\tilde{K})$ yields

$$ \omega_s(\tilde{K}) = \omega_0 - \int_0^{\tilde{K}} \left. \frac{\partial_{\tilde{K}} D_{-}}{\partial_{\omega_s} D_{-}} \right|_{\omega_s(\tilde{K}), \tilde{K}} d\tilde{K}. \qquad (3.22) $$

Numerically, the solution is found by starting at $\tilde{K} = 0$ and then increasing $\tilde{K}$ of small steps $d\tilde{K}$ until the desired value is reached. The derivatives of $D_-$ used in (3.22) are given in Appendix B.2.

In figures 2 and 3, we show the results obtained from (3.22) for three cases with the same wavenumber $k = 0.5$, but different electron spin polarization $\eta$. The results of the dispersion relation are compared with numerical results obtained with the fully nonlinear code with a small perturbation around the equilibrium, as detailed in § 5. For all cases, the agreement is excellent, which constitutes a cross-validation for both the numerical code and the above analytical developments.

In figure 2, we use the value of $\eta$ that is consistent with electrons at thermal equilibrium that are polarized by the magnetic field $B$ created by the magnetized ions, see (2.5) (we shall refer to this case as the 'self-consistent' case). In this case, the spin polarization is
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Magnon frequency $\omega_s$ for different normalized magnetic coupling constants $\tilde{K}$, obtained from (3.22) (continuous lines, red for the real part of the frequency, blue for the imaginary part), for $k = 0.5$ and $\eta = \tanh(H\tilde{K})$. Note that the electron spin polarization $\eta$ is different for different values of $\tilde{K}$. The dots represent numerical results obtained with the full numerical code described in the forthcoming sections. For this self-consistent case, the imaginary part remains very small with respect to the real part of the frequency.</figcaption></figure>

given by $\eta = \tanh(2\mu_B B/k_B T_e) = \tanh(H\tilde{K})$ and obviously depends on the electron-ion magnetic coupling – more details are given in Appendix A.

In contrast, in figure 3, we use two arbitrary values of the electron spin polarization, $\eta = 0.5$ and $\eta = -0.5$. The negative value means that the electrons are polarized in the opposite direction with respect to that of the self-consistent case. These values might be obtained through an external magnetic field that pre-polarizes the electrons prior to the application of a small perturbation. Nevertheless, one should keep in mind that to achieve such large spin polarizations, a very strong magnetic field would be needed, of the order of several hundred teslas.

For these values of $\eta$, the imaginary part of $\omega_s$ is significantly different from zero. In particular, for $\eta = 0.5$, there is a damping of the perturbation ($\text{Im}\omega < 0$), whereas for $\eta = -0.5$, we observe an instability ($\text{Im}\omega > 0$). This behaviour can be interpreted as follows. When $\eta > H\tilde{K} > 0$, the electron polarization has the same direction as in the self-consistent case, and hence the perturbation is damped, as the system tries to return to a state that has the ‘natural’ direction of polarization. In contrast, when $\eta < H\tilde{K}$ (and, in particular, when $\eta$ is negative), the system becomes unstable in an attempt to restore the ‘correct’ direction of polarization. When the value of $\eta$ corresponds to the self-consistent case, as in figure 2, the system is marginally stable ($\text{Im}\omega \approx 0$). Interestingly, in the self-consistent case, the first-order correction in the electron-magnon coupling $\tilde{K}$ disappears, see (3.25). Hence, figure 2 shows almost no variation of the real and imaginary parts of the magnon frequency for low values of $\tilde{K}$.

### 3.4. Weak coupling regime

From (3.20), the ion spin dispersion relation can be written as

$$ \omega_s = Ak^2 + \frac{\tilde{K}}{4}(\eta - H\tilde{K}) - Z\left(\frac{\omega_s - \omega_L}{k}\right)\left[\frac{\tilde{K}^2\eta}{8k} + \frac{\tilde{K}^2H}{4}\left(\frac{\omega_s - \omega_L}{k}\right)\right] \equiv G(\omega_s). \quad (3.23) $$

This is a transcendental equation for $\omega_s$, which cannot be solved exactly, except numerically as was done in the preceding subsection. An approximate solution to (3.23)
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Magnon frequency $\omega_s$ for different normalized magnetic coupling constants $\tilde{K}$, obtained from (3.22) (continuous lines, red for the real part of the frequency, blue for the imaginary part), for wavenumber $k = 0.5$ and electron polarizations (a) $\eta = 0.5$ and (b) $\eta = -0.5$. The dots represent numerical results obtained with the full numerical code described in the forthcoming sections. According to the value of $\eta$, the system is either (a) stable or (b) unstable.</figcaption></figure>

can be obtained iteratively, by starting with the solution for zero coupling, $\omega_0 = Ak^2$, then inserting this solution into the right-hand side of (3.23), which yields

$$ \omega_s \approx \omega_0 + \frac{\tilde{K}}{4}(\eta - H\tilde{K}) - Z \left( \frac{\omega_0 - \omega_L}{k} \right) \left[ \frac{\tilde{K}^2 \eta}{8k} + \frac{\tilde{K}^2 H}{4} \left( \frac{\omega_0 - \omega_L}{k} \right) \right], \quad (3.24) $$

which is valid for weak coupling $\tilde{K} \ll 1$. This procedure can be recast as a fixed-point problem: $\omega_s^{(\ell+1)} = G(\omega_s^{(\ell)}), \ell \in \mathbb{N}$, with $\omega_s^{(0)} = \omega_0 = Ak^2$, to obtain second- and higher-order approximations.

As the value of the dimensionless coupling constant is indeed small, $\tilde{K} \approx 0.16$, this weak-coupling approximation should hold for most cases of interest. Since $\tilde{K}/2 \equiv \omega_L/\omega_p$, physically, this approximation means that the electron Larmor frequency is much smaller than the plasmon frequency, specifically here, $\omega_L \approx 0.08\omega_p$. If we add the fact that the magnon frequency is $\omega_0 \approx 0.008\omega_p$, see (3.14), we obtain the following scaling between the three time scales that are present in this problem: $\omega_0 \ll \omega_L \ll \omega_p$.
Under such weak-coupling approximation, (3.24) simplifies to (restoring physical dimensions)

$$ \omega_s = \omega_0 + \frac{\omega_L}{2}(\eta - H\tilde{K}) \left[ 1 + \frac{2\omega_L}{v_{th}k} D\left(-\frac{\omega_L}{k}\right) - i\sqrt{\pi}\frac{\omega_L}{v_{th}k} \exp\left(-\frac{\omega_L^2}{v_{th}^2k^2}\right) \right], \quad (3.25) $$

where we used the fact that $Im Z(x) = \sqrt{\pi} e^{-x^2}$ ($x \in \mathbb{R}$) when evaluated on the real axis (i.e. $x \in \mathbb{R}$) (Fried & Conte 1961) and where $D$ is the Dawson function. By looking at the imaginary part of $\omega_s$, two regimes clearly appear. If $\eta < H\tilde{K}$, the imaginary part is positive, so that the magnetic perturbation is unstable and grows exponentially until the nonlinear regime is reached. If $\eta > H\tilde{K}$, then the perturbation is damped and disappears after a few oscillations. Interestingly, the value of $\eta$ that discriminates between these two regimes, i.e. $\eta = H\tilde{K}$, is precisely the value that corresponds to the self-consistent case, $\eta = \tanh(H\tilde{K})$, in the approximation where $\tilde{K} \ll 1$.

The form of the spin dispersion relation (3.25) reveals that all the magnetic terms in the Vlasov model (2.9) and (2.10) are important and cannot be neglected: the Zeeman terms proportional to $H$, the electron precession term proportional to $B$ (and hence to $\tilde{K}$), as well as the initial electron spin polarization $\eta$. The subtle interplay between these terms determines the stable or unstable nature of the linear response. In contrast, as we have seen, the electric charge response is completely decoupled from the spin response, at least in the linear regime. Hence, one could neglect the electric field terms in (2.9) and (2.10) (or set the initial electric perturbation to zero) and the spin response would remain unchanged. However, the plasmon oscillations would be lost.

The results for both the exact dispersion relation (3.22) and the approximate formula (3.25) are shown in figure 4 for a self-consistent case. As expected, the agreement is good for values up to $\tilde{K} \approx 1$, which cover most realistic values of the coupling constant. Finally, from (3.25), one can compute the maximum imaginary part with the parameters used in figure 4. Since $\tanh(H\tilde{K}) = H\tilde{K} - (H\tilde{K})^3/3 + O((H\tilde{K})^5)$, the imaginary part of $\omega_s$ is proportional to $\tilde{K}^5 e^{-(\tilde{K}/2k)^2}$. The maximum is then reached for $\tilde{K} = \sqrt{10k} \approx 1.58$, which is also in agreement with the exact dispersion relation.

Finally, in figure 5, we show the dependence of the magnon frequency on the wavenumber $k$, comparing the full dispersion relation with its first-order (3.24) and second-order approximations.

## 4. Numerical method

In this section, we present the numerical method used to solve the system of (2.1)–(2.4). The method is based on a Hamiltonian splitting technique, together with a phase space discretization that uses Fourier spectral approximation for the space variable $x$ and finite volumes (PSM) for the velocity direction $v$, as used by Crouseilles et al. (2023) and Li et al. (2020).

The Hamiltonian can be split into five parts:

$$ \mathcal{H} = \mathcal{H}_v + \mathcal{H}_E + \mathcal{H}_{S_1} + \mathcal{H}_{S_2} + \mathcal{H}_{S_3}, \quad (4.1) $$

where

$$ \left. \begin{aligned} \mathcal{H}_v &= \frac{1}{2} \int v^2 f_0 \, dx \, dv, & \mathcal{H}_E &= \frac{1}{2} \int \left( \frac{\partial V_H}{\partial x} \right)^2 \, dx, \\ \mathcal{H}_{S_i} &= H \int f_i B_i \, dx \, dv + A H \int \left( \frac{\partial S_i}{\partial x} \right)^2 \, dx, & i &= 1, 2, 3. \end{aligned} \right\} \quad (4.2) $$
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Magnon frequency $\omega_s$ as a function of the magnetic coupling constant $\tilde{K}$, for the self-consistent case $\eta = \tanh(\tilde{H}\tilde{K})$, and wavenumber $k = 0.5$. The solid lines represent the full dispersion relation computed numerically using (3.22), while the dashed lines are obtained with the simplified relation (3.24). Red lines refer to the real part of $\omega_s$, whereas blue lines refer to the imaginary part.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Magnon frequency $\omega_s$ as a function of the magnon wavenumber $k$, for electron polarization $\eta = 0.5$, and magnetic coupling constant $\tilde{K} = 0.16$. The solid lines represent the full dispersion relation computed numerically using (3.22) where the integral and the derivative are not with respect to $\tilde{K}$ but $k$. The other lines refer to the approximate linear theory obtained from (3.23) at first order (dashed lines, given explicitly by (3.24)), and second order (dotted lines). Red lines refer to the real part of $\omega_s$, whereas blue lines refer to its imaginary part.</figcaption></figure>

Let us remark that in this decomposition, $H_{S_i} = H_{Z,i} + H_{\text{spin},i}$, where the Zeeman energy $H_{Z,i}$ and the spin energy $H_{\text{spin},i}$ are given by (2.16). According to the Hamiltonian splitting, we get from (2.8)

$$ \frac{\partial \mathcal{Z}}{\partial t} = \{\mathcal{Z}, H_v\} + \{\mathcal{Z}, H_E\} + \{\mathcal{Z}, H_{S_1}\} + \{\mathcal{Z}, H_{S_2}\} + \{\mathcal{Z}, H_{S_3}\}, \quad (4.3) $$
which induces the five subsystems

$$ \frac{\partial \mathcal{Z}}{\partial t} = \{\mathcal{Z}, \mathcal{H}_v\}, \quad \frac{\partial \mathcal{Z}}{\partial t} = \{\mathcal{Z}, \mathcal{H}_E\}, \quad \frac{\partial \mathcal{Z}}{\partial t} = \{\mathcal{Z}, \mathcal{H}_{S_1}\}, \quad \frac{\partial \mathcal{Z}}{\partial t} = \{\mathcal{Z}, \mathcal{H}_{S_2}\}, \quad \frac{\partial \mathcal{Z}}{\partial t} = \{\mathcal{Z}, \mathcal{H}_{S_3}\}. \qquad (4.4) $$

As detailed in Appendix C, each subsystem can be solved exactly, which means that the error in time only originates from the time splitting and then can be controlled by using high-order splittings.

Denoting $\varphi_t^{\mathcal{H}_*}(\mathcal{Z}(0))$, the exact solution at time $t$ of $\partial_t \mathcal{Z} = \{\mathcal{Z}, \mathcal{H}_*\}$ (where $\star = v, E, S_1, S_2, S_3$) with the initial condition $\mathcal{Z}(t=0)$, the solution of the full model (4.3) is thus approximated by

$$ \mathcal{Z}(t) = (\Pi_{\star=v,E,S_1,S_2,S_3} \varphi_t^{\mathcal{H}_*}) \mathcal{Z}(0). \qquad (4.5) $$

This is a first-order splitting, but higher order splittings could also be derived. Since the splitting involves here five steps, we will restrict ourselves to the Strang scheme,

$$ \mathcal{Z}(t) = (\varphi_{t/2}^{\mathcal{H}_v} \circ \varphi_{t/2}^{\mathcal{H}_E} \circ \varphi_{t/2}^{\mathcal{H}_{S_1}} \circ \varphi_{t/2}^{\mathcal{H}_{S_2}} \circ \varphi_{t}^{\mathcal{H}_{S_3}} \circ \varphi_{t/2}^{\mathcal{H}_{S_2}} \circ \varphi_{t/2}^{\mathcal{H}_{S_1}} \circ \varphi_{t/2}^{\mathcal{H}_E} \circ \varphi_{t/2}^{\mathcal{H}_v}) \mathcal{Z}(0). \qquad (4.6) $$

Such Hamiltonian splittings are known to maintain long-term accuracy of the total energy. Moreover, in our case, one can also prove the scheme preserves exactly the norm of $\mathcal{S}$.

**PROPOSITION 4.1.** *The update (C13), (C20) and (C24) of the spin $\mathcal{S}$ through the Hamiltonian splitting discretization preserves the norm of the spin: $\|\mathcal{S}(\cdot, t)\| = 1$ if $\|\mathcal{S}^0(\cdot)\| = 1$.*

*Proof.* By (C13), (C20) and (C24), the vector spin $\mathcal{S}$ is updated through the multiplication of a matrix $\exp(\alpha \mathbf{J} t)$ ($\mathbf{J}$ being the symplectic matrix) which is a rotation matrix of angle $(-\alpha)$ in $\mathbb{R}^2$. Let us introduce the $3 \times 3$ matrix $\mathbf{A}$ corresponding to (C24)

$$ \mathbf{A} = \begin{pmatrix} \exp(\alpha \mathcal{H}_{S_3} \mathbf{J} t) & \mathbf{0}^T \\ \mathbf{0} & 1 \end{pmatrix}, \qquad (4.7) $$

with $\mathbf{0} = (0, 0)$ and $\alpha \mathcal{H}_{S_3} = (\tilde{K}/4) \int f_3^0 dv + A \partial_x^2 S_3^0$. We then reformulate (C24) as $\mathcal{S}(x, t) = \mathbf{A}\mathcal{S}^0(x)$ from which we easily deduce the norm is preserved. The same is true for (C13) and (C20). We finally deduce $\|\mathcal{S}(\cdot, t)\| = 1$ as long as $\|\mathcal{S}^0(\cdot)\| = 1$. $\square$

## 5. Numerical results

In this section, we present some numerical results obtained with the nonlinear code described in § 4. The results will also be compared with the analytical linear response, as detailed in § 3. In the results presented below, the numerical parameters are chosen as follows (non-dimensional units are used everywhere): number of points in space and velocity $N_x = 119$, $N_v = 1024$; time step $\Delta t = 0.1$; variable ranges in the phase space $v \in [-5, 5]$, $x \in [0, 2\pi/k]$; perturbation wavenumber $k = 0.5$.

The initial condition is a periodic perturbation of the equilibrium $f_0^{(0)} = \mathcal{F}, f_3^{(0)} = \eta \mathcal{F}, f_1^{(0)} = f_2^{(0)} = S_1^{(0)} = S_2^{(0)} = 0, S_3^{(0)} = 1$, where $\mathcal{F}$ is a spatially homogeneous equilibrium (either a Maxwell–Boltzmann or a two-stream distribution). This equilibrium represents ions that are fully polarized in the $\ell = 3$ direction, while the electrons are partially polarized along the same direction, with a polarization rate equal to $\eta$.
<table><caption>TABLE 1. Main numerical and physical parameters of the three runs that use a Maxwell–Boltzmann (MB) equilibrium: initial electron spin polarization η, ion spin frequency ωs and initial perturbation ε. The values of ωs are those of the linear response calculation using the ZEAL code. Other values are k = 0.5, Nₓ = 119, Nᵥ = 1024 and vmax = 5.</caption><tr><td>MB1</td><td>η = tanh(CKH) ≈ 0.0547</td><td>ωs = 0.003680 - 2.739 × 10<sup>-5</sup>i</td><td>ε = 10<sup>-3</sup></td></tr><tr><td>MB2</td><td>η = 0.5</td><td>ωs = 0.02088 - 0.005253i</td><td>ε = 10<sup>-3</sup></td></tr><tr><td>MB3</td><td>η = -0.5</td><td>ωs = 0.01725 + 0.006162i</td><td>ε = 10<sup>-6</sup></td></tr></table>

After the perturbation, the initial condition is as follows:

$$ 
\left.
\begin{aligned}
f_0(t = 0^+, x, v) &= \mathcal{F}(v)(1 + \varepsilon \cos(kx)), \\
f_1(t = 0^+, x, v) &= \eta \mathcal{F}(v)\varepsilon \cos(kx), \\
f_2(t = 0^+, x, v) &= \eta \mathcal{F}(v)\varepsilon \sin(kx), \\
f_3(t = 0^+, x, v) &= \eta \mathcal{F}(v)(1 + \varepsilon \cos(kx)), \\
S_1(t = 0^+, x) &= \frac{\varepsilon}{\sqrt{1 + \varepsilon^2}} \sin(kx), \\
S_2(t = 0^+, x) &= \frac{\varepsilon}{\sqrt{1 + \varepsilon^2}} \cos(kx), \\
S_3(t = 0^+, x) &= \frac{1}{\sqrt{1 + \varepsilon^2}},
\end{aligned}
\right\} (5.1)
$$

where the amplitude of the perturbation is ε = 10⁻³. Note that the perturbation is chosen such that ||S(t = 0, x)||² = S₁²(0, x) + S₂²(0, x) + S₃²(0, x) = 1. The non-dimensional physical constants are those defined in § 2.2, i.e. A = 0.0148 (ion-ion magnetic coupling), K̃ = 0.161 (ion-electron magnetic coupling) and H = 0.339 (scaled Planck constant). The numerical results will be expressed in terms of the units defined in § 2.2. All logarithms are Napierian (base e).

### 5.1. Maxwell–Boltzmann (MB) equilibrium

Here, we consider the Maxwell–Boltzmann equilibrium (3.15a,b) that was used for the linear analysis. We will analyse three cases for different electron polarizations η. In the first case (MB1), the polarization is taken to be self-consistent with the ions, i.e. the electron polarization is due solely to the magnetic field generated by the ions, so that η = tanh(K̃H) (see Appendix A). In the remaining two cases (MB2 and MB3), the polarization will be chosen arbitrarily as η = ±0.5. This polarization may be achieved through the application of an external magnetic field. The parameters of these Maxwell–Boltzmann simulations are summarized in table 1.

MB1. The roots of the dispersion relation for charges (ωe) and spins (ωs), calculated using the ZEAL code, are the following:

$$ \omega_e = 1.225 - i0.03626, \qquad (5.2) $$

$$ \omega_s = 0.003680 - i2.739 \times 10^{-5}. \qquad (5.3) $$

We remark that: (i) the real part of ωe is close to the plasma frequency (equal to unity here), while its imaginary part is much smaller, in accordance with the Bohm–Gross dispersion
<figure><img src="image_7.png" /><figcaption>FIGURE 6. MB1 simulation. Time history of the square root of the electric energy $H_E^{1/2}$ (given by (4.2)), (a) in semi-log scale and (b) corresponding frequency spectrum. The red straight line represents the linear damping rate given in (5.2). Time history of the absolute value of the real part of the first Fourier mode of the ion spin $\hat{S}_1(k, t)$ in (c) semi-log scale and (d) corresponding frequency spectrum. The red straight line corresponds to zero damping, see (5.3). The arrows in the spectral plots correspond to the results of linear response theory.</figcaption></figure>

relation; (ii) the real part of $\omega_s$ is much smaller than the plasma frequency, in accordance with (3.14), while its imaginary part is even smaller, signifying the almost absence of spin damping.

In figure 6, we plot the time evolution of some physical quantities associated to the electron charge in panels (a,b) and ion spin in panels (c,d). The Coulomb electric energy decays exponentially with a rate $Im\omega_e$ very close to the one predicted by the linear response analysis (Landau damping). The real part of the frequency is also very close to the analytical prediction of (5.2), with an additional factor of 2 due to the modulus.

In figure 6(c,d), we show the evolution of the absolute value of the real part of the first Fourier mode of the ion spin $S_1(x, t)$, i.e. $\hat{S}_1(k, t)$, with $k = 0.5$ in this case. In agreement with (5.3), this mode is virtually undamped (the red line is horizontal and corresponds to zero damping). The corresponding frequency spectrum peaks in the vicinity of the theoretical magnon frequency $Re\omega_s$. Note that, due to the great disparity between the magnon and the plasmon frequencies, only a few ($\approx 6$) magnon frequencies could be observed, resulting in a limited accuracy for the magnon spectrum.

In addition to the good agreement with the linear theory for $\omega_e$ and $\omega_s$, we also emphasize that the modulus of the ion spin vector $||S(t, x)||$ is preserved up to machine accuracy and that the (relative) total energy is preserved up to $10^{-7}$.
<figure><img src="image_8.png" /></figure>

FIGURE 7. MB2 simulation ($\eta = 0.5$). Time history of the absolute value of the real part of the first Fourier mode of the ion spin $\hat{S}_1(k, t)$ (a) in semi-log scale and (b) corresponding frequency spectrum. The slope of the red straight line is -0.005186, very close to the linear response result given in table 1. The peak of the frequency spectrum also matches the linear result $\text{Re}\,\omega_s = 0.02088$ (indicated by an arrow on the plot) with good accuracy. Panels (c,d) show the same quantities for the electronic spin mode $\hat{M}_1(k, t)$. The real and imaginary parts of the frequency are the same as for the ion spins.

MB2. For this second test, we consider an initial condition with an electron spin polarization rate $\eta = 0.5$. This can be achieved through an external magnetic field $B_3^{\text{ext}}$ directed along the same direction as the ion polarization. The positive value of $\eta$ corresponds to the ‘natural’ polarization direction for the electrons, parallel to that of the ions and oriented in the same way, as in the self-consistent case. Hence, we expect this equilibrium to be magnetically stable.

As was mentioned earlier, the charge dynamics is decoupled from the spin dynamics in the linear regime, and hence the electric response (not shown here) is the same as that of figure 6, displaying plasmonic oscillations and Landau damping.

The spin response is depicted in figure 7, where we show the first Fourier moments of the ion and electron spins and their frequency spectra. In this case, a clear damping of the magnon mode is observed, which is in good agreement with the roots of the dispersion relation: $\omega_s = 0.02088 - i0.005253$, which is to be compared with the damping rate obtained from the simulation, $\gamma = -0.005186$. The real part of the frequency, see figure 7(b), shows a peak near $\text{Re}\,\omega_s \approx 0.02$, also in good accordance with the linear response result.

The electron spin density $M$, shown in figure 7(c,d), follows the same evolution as the ions, with very similar frequency and damping rate.
<figure><img src="image_9.png" /><figcaption>FIGURE 8. MB3 simulation ($\eta = -0.5$). Time history of the absolute value of the real part of the first Fourier mode of the ion spin $\hat{S}_1(k, t)$ (a) in semi-log scale and (b) corresponding frequency spectrum. The slope of the red straight line is 0.00607, very close to the linear response result given in table 1. The peak of the frequency spectrum also matches the linear result $Re\omega_s = 0.01725$ (indicated by an arrow on the plot) with good accuracy. Panels (c,d) show the same quantities for the electronic spin mode $\hat{M}_1(k, t)$. The real and imaginary parts of the frequency are the same as for the ion spin.</figcaption></figure>

MB3. Here, we consider an electron gas which is initially polarized in the opposite direction to the one corresponding to the self-consistent case. In this case, the polarization rate is negative and we take $\eta = -0.5$. Since the electron polarization is opposite to the self-consistent scenario, we expect the system to be unstable, as it attempts to restore the ‘natural’ direction of polarization.

In figure 8(a,b), we plot the evolution of the first Fourier mode of the ion spin and its frequency spectrum. The real part of the frequency and the instability rate are very close to the linear response result $\omega_s = 0.01725 + i0.006162$. After approximately $2000\omega_p^{-1}$, the instability saturates nonlinearly. The electric field evolution is the same as in figure 6(a).

The electron spin density $M$, shown in figure 8(c,d), follows the same evolution as the ions, with very similar frequency and instability rate.

## 5.2. Two-stream (TS) equilibrium

In this subsection, we consider a two-stream equilibrium for the initial electron distribution

$$ \mathcal{F}(v) = \frac{1}{2\sqrt{\pi}}(e^{-(v-u)^2} + e^{-(v+u)^2}). \qquad (5.4) $$
<table><caption>TABLE 2. Main numerical and physical parameters of the runs that use a two-stream equilibrium. Other values are $k = 0.2$ and $N_x = 129$.</caption><tr><td>TS1</td><td>η = tanh(CKH) ≈ 0.0547</td><td>u = 1.4</td><td>N<sub>v</sub> = 1024, v<sub>max</sub> = 8, ε = 10<sup>-3</sup></td></tr><tr><td>TS2</td><td>η = tanh(CKH) ≈ 0.0547</td><td>u = 3</td><td>N<sub>v</sub> = 1536, v<sub>max</sub> = 12, ε = 10<sup>-6</sup></td></tr><tr><td>TS3</td><td>η = -0.5</td><td>u = 3</td><td>N<sub>v</sub> = 512, v<sub>max</sub> = 14, ε = 10<sup>-6</sup></td></tr><tr><td>TS4</td><td>η = -0.5</td><td>u = 1.4</td><td>N<sub>v</sub> = 512, v<sub>max</sub> = 14, ε = 10<sup>-6</sup></td></tr></table>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. TS1 simulation. (a) Time evolution of the square root of the electric energy $\mathcal{H}_E^{1/2}$ (given by (4.2)) for short times $t \in [0, 50]$. (b) Time evolution of the absolute value of the real part of the fundamental mode of the electron spin $\hat{M}_1$. (c) Time evolution of the absolute value of the real part of the fundamental mode of the ion spin $\hat{S}_1$. The red straight lines have slopes equal to zero for the electric energy and $-6.3 \times 10^{-4}$ for $\hat{M}_1$ and $\hat{S}_1$.</figcaption></figure>

This equilibrium can be either stable or unstable for the charge dynamics, depending on the value of the stream velocity $u$. In the numerical runs reported below, we have chosen $u = 1.4$, which corresponds to a stable case (run TS1), and $u = 3$ which corresponds to an unstable case (runs TS2 and TS3). In TS1 and TS2, we use the self-consistent value for the electron spin polarization, $\eta = \tanh(CKH) \approx 0.0547$, while in TS3 and TS4, we force a spin instability by setting $\eta = -0.5$.

The parameters of these runs are summarized in table 2.

TS1. In this case, the stream velocity is weak ($u = 1.4$) so that the charge sector of the dynamics is basically undamped, as seen in figure 9(a) for the electric field. The spin sector is more interesting, both for the ions and the electrons, which are rather strongly damped
<figure><img src="image_11.png" /><figcaption>FIGURE 10. TS2 simulation. (a) Time evolution of the square root of the electric energy $H_E^{1/2}$ (given by (4.2)) for short times $t \in [0, 100]$. (b) Time evolution of the absolute value of the real part of the fundamental mode of the electron spin $\hat{M}_1$. (c) Time evolution of the absolute value of the real part of the fundamental mode of the ion spin $\hat{S}_1$. The red straight lines have slopes equal to 0.2845 for the electric energy and $4 \times 10^{-5}$ for $\hat{S}_1$.</figcaption></figure>

at a rate $\approx 6.3 \times 10^{-4}$. This is in contrast with the corresponding Maxwell-Boltzmann simulation (MB1, figure 6) where the spin mode was very weakly damped. Although the wavenumber is not the same ($k = 0.5$ for MB1 and $k = 0.2$ for TS1), it appears that the equilibrium profile has a strong impact on the stability properties of the ion magnon mode.

TS2. This run uses the same parameters as TS1, except that the stream velocity is larger, $u = 3$. We also changed the magnitude of the initial perturbation, now set to $\epsilon = 10^{-6}$, to get a longer-lasting linear phase. Linear theory predicts an instability in the charge sector, with growth rate equal to 0.2845, which is confirmed by the numerical data shown in figure 10(a). The ion spin sector displays a very weak instability, with an observed growth rate $\approx 4 \times 10^{-5}$. The electron spin remains at very low amplitude all along the simulation time.

TS3. Here, we wish to consider a case where an instability is expected both in the charge and in the spin sectors. Therefore, we take the same value $u = 3$ for the stream velocity as in TS2, and an electron spin polarization $\eta = -0.5$, which leads to the instability of the magnon mode in MB3. The results are plotted in figure 11 and show that the evolution of the electric field is almost the same as in the case TS2. This is natural, as the linear response of the charge sector is independent of $\eta$ (nonetheless, one may have expected some differences after the nonlinear regime is attained, around $\omega_p t = 50$, but in practice, the two curves are very similar, although not identical). Interestingly, the electron $\hat{M}_1$ and
<figure><img src="image_12.png" /><figcaption>FIGURE 11. TS3 simulation. (a) Time evolution of the square root of the electric energy $H_E^{1/2}$ (given by (4.2)) for short times $t \in [0, 100]$; the red straight line has slope equal to 0.2845. (b) Time evolution of the absolute value of the real part of the fundamental mode of the electron spin $\hat{M}_1$. (c) Time evolution of the absolute value of the real part of the fundamental mode of the ion spin $\hat{S}_1$.</figcaption></figure>

ion $\hat{S}_1$ spins are initially stable until $\omega_p t \approx 2500$, i.e. well into the nonlinear regime, and only become unstable later. Their growth rate is much smaller than the one associated with the charges.

The phase space portraits at the end of the simulation are displayed in figure 12 for the four distributions $f_0(t, x, v)$ and $f_\ell(t, x, v)$, $\ell = 1, 2, 3$. Typically, for this type of instability, the two-stream structure has been destroyed in the nonlinear regime and a single vortex centred at $v = 0$ can be observed. The vortex is present not only in the charge distribution $f_0$, but also in the spin distributions $f$.

TS4. Finally, we repeat the same simulation as TS3, but for a smaller stream velocity $u = 1.4$, so that there is no instability in the charge sector (see figure 13). In this case, the usual magnon instability ($\eta < 0$) develops immediately, in contrast to the preceding TS3 case. Although it is difficult to draw definite conclusions, it is clear that the onset, or otherwise, of a charge instability interacts strongly with the development of a magnon instability. This is further evidence that the charge and spin sectors are closely intertwined and need to be both included in the model for an accurate description of the magnonic dynamics.

## 6. Conclusion

In this work, we have built on previous developments (Crouseilles et al. 2021, 2023; Manfredi et al. 2023) to construct a fully kinetic 1-D model of the interaction between
<figure><img src="image_13.png" /></figure>

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

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

<figure><img src="image_16.png" /><figcaption>FIGURE 12. TS3 simulation. Contour plots of the distribution functions in the (x, v) phase space at the final time $\omega_{pt} = 10^4$: (a) $f_0$; (b) $f_1$; (c) $f_2$ and (d) $f_3$.</figcaption></figure>

the charge and the spin dynamics in a material with intrinsic magnetization (ferromagnet). The electron dynamics is described by a four-component phase space distribution function $f_0(t, x, v)$, $f_\ell(t, x, v)$, $\ell = 1, 2, 3$, where $f_0$ is related to the electron charge and $f_\ell$ to the electron spin polarization in the $\ell$ direction. The fixed ions are modelled by the Landau–Lifshitz equation for the magnetization $S(t, x)$. The electron charges interact through the self-consistent electric field, solution of the Poisson equation. The electron and ion spins interact through the magnetic exchange, whose magnitude is controlled by the coupling constant $K$. Finally, the ion spins interact among themselves via the ion–ion magnetic exchange, with coupling constant $J$.

This model can be seen as an extension of the standard Vlasov–Poisson equations for mobile electrons and fixed ions, taking into account the electron spin and allowing for a spin dynamics for the ions.

We first focused on the linear response of this system when the equilibrium is a Maxwell–Boltzmann function. The full dispersion relation is rather complex, but can be split into a charge sector and a spin sector. The former is independent of the spin and leads to the standard Bohm–Gross relation. The spin sector was analysed more in detail, particularly the occurrence of damping and instability when the ion–electron magnetic coupling constant and the electron spin polarization at equilibrium are varied. Interestingly, we observe damping when the electron spin polarization is directed along the ‘natural’ direction of magnetization (the one dictated by the magnetic field generated by the ions) and instability when it is directed opposite to it.
<figure><img src="image_17.png" /><figcaption>FIGURE 13. TS4 simulation. (a) Time evolution of the square root of the electric energy $H_E^{1/2}$ (given by (4.2)) for short times $t \in [0, 100]$. (b) Time evolution of the absolute value of the real part of the fundamental mode of the electron spin $\hat{M}_1$. (c) Time evolution of the absolute value of the real part of the fundamental mode of the ion spin $\hat{S}_1$. The red straight lines have slope equal to 0.004.</figcaption></figure>

Next, we built a computational code based on the Hamiltonian splitting method first developed by Crouseilles et al. (2023) and Li et al. (2020). This is an Eulerian grid-based method that solves simultaneously the coupled Vlasov-Poisson-Landau-Lifshitz equations. This technique allowed us to achieve great accuracy for the conserved quantities: the modulus of the ion spin vector $||S(t, x)||$ is preserved up to machine accuracy and the (relative) total energy is preserved up to $\approx 10^{-7}$.

We have used the code to validate the estimations of the linear response theory, with very good agreement between the two approaches for Maxwell-Boltzmann equilibria. We also tested it on two-stream equilibria, which may lead to instability in the charge sector, depending on the streams' relative velocities. Particularly interesting was the case where an instability in the charge sector leads to a much delayed instability in the spin sector, which develops well after the charge dynamics has saturated nonlinearly. This is further evidence of the close interaction between the charge and spin sectors in the coupled plasmon-magnon dynamics.

The Maxwell-Boltzmann equilibria and parameter range used in this work, with densities close to those of solids ($\approx 10^{29} \text{ m}^{-3}$) and temperatures of the order of 10 eV, are relevant to the warm dense matter (WDM) regime (Bonitz et al. 2020) that appears, among others, in inertial fusion experiments. For these conditions, the electron plasma is weakly degenerate ($T_e \approx T_F$), so that it can be characterized with relatively good accuracy by a MB distribution. The ions are fixed and non-degenerate. In this WDM regime,
ultrafast non-equilibrium dynamics has been recently observed thanks to subpicosecond laser pulses (Falk 2018). At these very short time scales, and for magnetic materials, the electron and ion spin polarization may not yet be lost, and impact the early instants of the dynamics.

However, MB distributions are not relevant to condensed-matter systems – for which the Fermi temperature is well above the room temperature – and the latter should therefore be described by a Fermi–Dirac (FD) equilibrium. Calculations of the dispersion relation for FD distributions are notoriously more involved than for MB distributions, particularly in the finite-temperature case. These developments are left for future work.

### Acknowledgements

The authors thank Aurélien Manchon for useful discussions and advice. Editor Dmitri Uzdensky thanks the referees for their advice in evaluating this article.

### Funding

This work is 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 was partially funded by the French National Research Agency (ANR) through the Programme d'Investissement d'Avenir under contract ANR-11-LABX-0058-NIE and ANR-17-EURE-0024 within the Investissement d'Avenir program ANR-10-IDEX-0002-02. N.C. and X.H. would like to thank the Centre Henri Lebesgue, program ANR-11-LABX-0020-0 and the support from Brittany's regional council.

### Declaration of interest

The authors report no conflict of interest.

### Author contributions

N.C. and X.H. developed the computational method and wrote the code; B.B., G.M. and P.A.H. derived the theory; B.B. and N.C. performed the analytical and numerical calculations. All authors contributed equally to analysing data and reaching conclusions, and in writing the paper.

## Appendix A. Spin-polarized equilibrium

To compute stationary states, it is more convenient to go back to the standard representation of the Wigner function (Manfredi et al. 2019):

$$ \mathcal{F} = \begin{pmatrix} \mathcal{F}_{++} & \mathcal{F}_{+-} \\ \mathcal{F}_{-+} & \mathcal{F}_{--} \end{pmatrix}, \qquad (A1) $$

where + (−) stands for spin-up (spin-down) with respect to the direction $\ell = 3$. The relationship between this representation and the Pauli representation used in the main text is the following: $f_\ell = \text{tr}(\sigma_\ell \mathcal{F})$, $f_0 = \text{tr}(\mathcal{F})$, where $\sigma_\ell (\ell = 1, 2, 3)$ are the Pauli matrices:

$$ \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}, \quad \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \qquad (A2a-c) $$

For a spatially homogeneous equilibrium, the terms corresponding to the self-consistent electric energy $\mathcal{H}_E$ and the spin energy $\mathcal{H}_{\text{spin}}$ vanish from the expression of the Hamiltonian (2.16). In the above basis, the Hamiltonian is a diagonal $2 \times 2$ matrix
diag($\mathcal{H}_+$, $\mathcal{H}_-$), where $\mathcal{H}_\pm = (m/2)v^2 \pm \mu_B B_3$ is the signed sum of the kinetic and Zeeman energies and $B_3$ is the magnetic field generated by the (fully polarized ions), see (2.11). In our dimensionless units, $B_3 = -\tilde{K}/2$, and we get for the Hamiltonian, $\mathcal{H}_\pm = v^2 \mp H\tilde{K}$.

For a stationary state, the distribution function must be a function of the Hamiltonian, i.e. in the Maxwell-Boltzmann case, $\mathcal{F} = C \exp(-\beta\mathcal{H})$, where $C$ is a normalization constant. Hence, the distribution function is also diagonal, with $\mathcal{F}_{++} = C \exp(-\beta\mathcal{H}_+)$, and similarly for $\mathcal{F}_{--}$, where $\beta = 1/(k_B T_e) = 1$ in our units.

Going back to the Pauli basis used in the main text, we obtain

$$ \left. \begin{aligned} f_0(v) &= \mathcal{F}_{++} + \mathcal{F}_{--} = 2C e^{-v^2} \cosh(H\tilde{K}), \\ f_3(v) &= \mathcal{F}_{++} - \mathcal{F}_{--} = 2C e^{-v^2} \sinh(H\tilde{K}). \end{aligned} \right\} \quad (A3) $$

With the normalization $\int f_0(v) dv = 1$, we get $C = 1/(2\sqrt{\pi} \cosh(H\tilde{K}))$. As a consequence, the equilibrium distribution function becomes

$$ f_0(v) = \frac{e^{-v^2}}{\sqrt{\pi}}, \quad f_3(v) = \eta \frac{e^{-v^2}}{\sqrt{\pi}}, \qquad (A4a,b) $$

with $\eta = \tanh(H\tilde{K})$, which is identical to (3.15a,b) in the main text.

Finally, if the magnetic field in the Hamiltonians $\mathcal{H}_\pm$ is not the one generated self-consistently by the ions, but instead an external one $B_3^{\text{ext}}$, then the electron spin polarization is $\eta = \tanh(2\mu_B B_3^{\text{ext}}/k_B T_e)$ and can take any values in $[-1, 1]$. Note that $\eta > 0$ corresponds to a case where the ion spin $S(t, x)$ and electrons spin $M(t, x) = (\hbar/2) \int f(t, x, v) dv$ are aligned along the same direction, which is a stable ferromagnetic equilibrium. In contrast, when $\eta < 0$, the ion and electron spins point in opposite directions, leading to an unstable equilibrium. This is confirmed by the simulations reported in § 5.1.

## Appendix B. Dispersion relation details

In this appendix, some details are given about the analytical dispersion relation. In particular, a new form of the dispersion function $D_S(\omega, k)$ is presented and its derivatives are computed explicitly.

### B.1. Alternative form of the dispersion relation

The dispersion relation (3.17) writes as (in this appendix, we denote the magnon frequency $\omega$ instead of $\omega_s$, for simplicity of notation)

$$ \begin{aligned} D_S(\omega, k) = & - \left[ \omega + \frac{c_0}{k} \left[ Z\left(\frac{\omega + \tilde{K}/2}{k}\right) + Z\left(\frac{\omega - \tilde{K}/2}{k}\right) \right] \right. \\ & \left. - c_1 \left[ Z'\left(\frac{\omega - \tilde{K}/2}{k}\right) - Z'\left(\frac{\omega + \tilde{K}/2}{k}\right) \right] \right]^2 \\ & + \left[ Ak^2 + d + \frac{c_0}{k} \left[ Z\left(\frac{\omega + \tilde{K}/2}{k}\right) - Z\left(\frac{\omega - \tilde{K}/2}{k}\right) \right] \right. \\ & \left. + c_1 \left[ Z'\left(\frac{\omega - \tilde{K}/2}{k}\right) + Z'\left(\frac{\omega + \tilde{K}/2}{k}\right) \right] \right]^2, \end{aligned} \quad (B1) $$
with $c_0 = \tilde{K}^2\eta/16$, $c_1 = \tilde{K}^2H/16$, $d = \tilde{K}\eta/4$ (recall that $\eta = \int f_3^{(0)} dv$). Factorizing leads to

$$ D_S(\omega, k) = - \left[ \omega + Ak^2 + d + \frac{2c_0}{k} Z\left(\frac{\omega + \tilde{K}/2}{k}\right) + 2c_1 Z'\left(\frac{\omega + \tilde{K}/2}{k}\right) \right] \\ \times \left[ \omega - Ak^2 - d + \frac{2c_0}{k} Z\left(\frac{\omega - \tilde{K}/2}{k}\right) - 2c_1 Z'\left(\frac{\omega - \tilde{K}/2}{k}\right) \right]. \quad (B2) $$

Naming $D_+$ the first term on the right-hand side and $D_-$ the second term, $D_-(-\omega^*, k)$ can be computed, where the asterisk denotes the complex conjugate:

$$ D_{-}(-\omega^{*}, k) = -\omega^{*} - Ak^{2} - d + \frac{2c_{0}}{k}Z\left(\frac{-\omega^{*} - \tilde{K}/2}{k}\right) - 2c_{1}Z'\left(\frac{-\omega^{*} - \tilde{K}/2}{k}\right) \\ = -\omega^{*} - Ak^{2} - d + \frac{2c_{0}}{k}Z\left(-\left(\frac{\omega + \tilde{K}/2}{k}\right)^{*}\right) - 2c_{1}Z'\left(-\left(\frac{\omega + \tilde{K}/2}{k}\right)^{*}\right). \quad (B3) $$

Now, some symmetries in $Z(-z^*)$ and $Z'(-z^*)$ can be used (Fried & Conte 1961): $Z(-z^*) = -Z^*(z)$ so $Z'(-z^*) = -2(1 - z^*Z(-z^*)) = -2(1 + z^*Z^*(z)) = Z^*(z)$. Finally, $D_-(-\omega^*, k)$ is expressed as

$$ D_{-}(-\omega^{*}, k) = - \left[ \omega + Ak^2 + d + \frac{2c_0}{k} Z\left(\frac{\omega + \tilde{K}/2}{k}\right) + 2c_1 Z'\left(\frac{\omega + \tilde{K}/2}{k}\right) \right]^* = -D_{+}^*(\omega, k). \quad (B4) $$

Then, we get $D_S(\omega, k) = D_-(-\omega^*, k)D_-(\omega, k)$. Hence, if $\omega$ satisfies $D_S(\omega, k) = 0$, then $D_S(-\omega^*, k)$ also vanishes. Therefore, we will consider

$$ D_{-}(\omega, k) = \omega - Ak^2 - d + \frac{2c_0}{k} Z\left(\frac{\omega - \tilde{K}/2}{k}\right) - 2c_1 Z'\left(\frac{\omega - \tilde{K}/2}{k}\right) \quad (B5) $$

as the dispersion relation instead of $D_S$. Since $Z'(z) = -2(1 + zZ(z))$,

$$ D_{-}(\omega, k) = \omega - Ak^2 - d + 4c_1 + Z(z) \left[ \frac{2c_0}{k} + 4c_1z \right], \quad (B6) $$

with $z = (\omega - \tilde{K}/2)/k$. Using the expressions of $d$, $c_0$, $c_1$ in terms of $\tilde{K}$, we obtain

$$ D_{-}(\omega, k) = \omega - Ak^2 - \frac{\tilde{K}\eta}{4} + \frac{\tilde{K}^2H}{4} + Z(z) \left[ \frac{\tilde{K}^2\eta}{8k} + \frac{\tilde{K}^2H}{4}z \right], \quad (B7) $$

which can also be interpreted as a function of $(\omega, \tilde{K})$ for a constant value of $k$.
## B.2. Computation of the derivatives of $D_-$

The partial derivatives of $D_-(\omega, \tilde{K})$ given by (B7) with respect to $\omega$ and $\tilde{K}$ can be computed as follows (with $\eta = \tanh(H\tilde{K})$ and $z = (\omega - \tilde{K}/2)/k$):

$$ \left. \begin{aligned} \frac{\partial D_{-}}{\partial \tilde{K}} &= -\frac{\eta}{4} - \frac{\tilde{K}H(1-\eta^2)}{4} + \frac{\tilde{K}H}{2} + \frac{\tilde{K}^2\eta}{8k^2} + \frac{\tilde{K}^2Hz}{4k} \\ &\quad + Z(z) \left[ \frac{\tilde{K}^2\eta z}{8k^2} + \frac{\tilde{K}^2Hz^2}{4k} + \frac{\tilde{K}\eta}{4k} + \frac{\tilde{K}^2H(1-\eta^2)}{8k} + \frac{\tilde{K}Hz}{2} - \frac{\tilde{K}^2H}{8k} \right] \\ \frac{\partial D_{-}}{\partial \omega} &= 1 - \frac{\tilde{K}^2\eta}{4k^2} - \frac{\tilde{K}^2Hz}{2k} + Z(z) \left[ -\frac{\tilde{K}^2\eta z}{4k^2} - \frac{\tilde{K}^2Hz^2}{2k} + \frac{\tilde{K}^2H}{4k} \right] \end{aligned} \right\} . \quad (B8) $$

# Appendix C. Time splitting

In this appendix, we give the details of the time solution of the different subsystems induced by the Hamiltonian splitting, as detailed in § 4. Regarding the space approximation, Fourier spectral methods are used, so that the linear transport operators (for the Vlasov part) and the elliptic operators (for the Poisson equation) reduce to a simple multiplication in the Fourier space. In the velocity direction, the linear transport operators in the Vlasov equations are approximated by using a semi-Lagrangian method based on finite volumes (see Crouseilles, Mehrenberger & Sonnendrücker (2010) for more details). Finally, all the integrals in velocity space are approximated by standard rectangle quadratures.

## C.1. Subsystem for $\mathcal{H}_v$

The subsystem $\partial \mathcal{Z}/\partial t = \{\mathcal{Z}, \mathcal{H}_v\}$ associated to $\mathcal{H}_v = \frac{1}{2} \int v^2 f_0 \, dx \, dv$ is

$$ \left. \begin{aligned} \frac{\partial f_0}{\partial t} &= \{f_0, \mathcal{H}_v\} = -v \frac{\partial f_0}{\partial x} \\ \frac{\partial f_\ell}{\partial t} &= \{f_\ell, \mathcal{H}_v\} = -v \frac{\partial f_\ell}{\partial x}, \quad \ell = 1, 2, 3 \\ \frac{\partial \mathcal{S}}{\partial t} &= \{\mathcal{S}, \mathcal{H}_v\} = \mathbf{0} \\ \partial_x^2 V_H &= \int f_0 \, dv - 1 \end{aligned} \right\} . \quad (C1) $$

We denote the initial value as $(f_0^0(x, v), f^0(x, v), \mathcal{S}^0(x))$ at time $t=0$. The solution at time $t$ of this subsystem can be written explicitly:

$$ f_0(t, x, v) = f_0^0(x - vt, v), \quad f(t, x, v) = f^0(x - vt, v), \quad \mathcal{S}(t, x) = \mathcal{S}^0(x). \quad (C2a-c) $$

## C.2. Subsystem for $\mathcal{H}_E$

The subsystem $\partial \mathcal{Z}/\partial t = \{\mathcal{Z}, \mathcal{H}_E\}$ associated to $\mathcal{H}_E = \frac{1}{2} \int (\partial V_H / \partial x)^2 \, dx$ is

$$ \left. \begin{aligned} \frac{\partial f_0}{\partial t} &= \{f_0, \mathcal{H}_E\} = -\frac{\partial V_H}{\partial x} \frac{\partial f_0}{\partial v} \\ \frac{\partial f_\ell}{\partial t} &= \{f_\ell, \mathcal{H}_E\} = -\frac{\partial V_H}{\partial x} \frac{\partial f_\ell}{\partial v}, \quad \ell = 1, 2, 3 \\ \frac{\partial \mathcal{S}}{\partial t} &= \{\mathcal{S}, \mathcal{H}_E\} = \mathbf{0} \end{aligned} \right\} . \quad (C3) $$
With the initial value $(f_0^0(x, v), f^0(x, v), S^0(x))$ at time $t=0$, the solution at time $t$ is as follows:

$$f_0(t, x, v) = f_0^0 \left( x, v - t \frac{\partial V_H}{\partial x}(x) \right), \quad f(t, x, v) = f^0 \left( x, v - t \frac{\partial V_H}{\partial x}(x) \right), \quad S(t, x) = S^0(x). \tag{C4a-c}$$

### C.3. Subsystem for $\mathcal{H}_{S_1}$

The subsystem $\partial \mathcal{Z}/\partial t = \{\mathcal{Z}, \mathcal{H}_{S_1}\}$ associated to $\mathcal{H}_{S_1} = H \int f_1 B_1 dx dv + AH \int (\partial S_1 / \partial x)^2 dx$ is

$$\left.
\begin{aligned}
\frac{\partial f_0}{\partial t} &= \{f_0, \mathcal{H}_{S_1}\} = H \frac{\partial B_1}{\partial x} \frac{\partial f_1}{\partial v} \\
\frac{\partial f_1}{\partial t} &= \{f_1, \mathcal{H}_{S_1}\} = H \frac{\partial B_1}{\partial x} \frac{\partial f_0}{\partial v} \\
\frac{\partial f_2}{\partial t} &= \{f_2, \mathcal{H}_{S_1}\} = -B_1 f_3 \\
\frac{\partial f_3}{\partial t} &= \{f_3, \mathcal{H}_{S_1}\} = B_1 f_2 \\
\frac{\partial S_1}{\partial t} &= \{S_1, \mathcal{H}_{S_1}\} = 0 \\
\frac{\partial S_2}{\partial t} &= \{S_2, \mathcal{H}_{S_1}\} = \frac{\tilde{K}}{4} S_3 \int f_1 \, dv + A S_3 \partial_x^2 S_1 \\
\frac{\partial S_3}{\partial t} &= \{S_3, \mathcal{H}_{S_1}\} = -\frac{\tilde{K}}{4} S_2 \int f_1 \, dv - A S_2 \partial_x^2 S_1
\end{aligned}
\right\}, \tag{C5}$$

with the initial value $(f_0^0(x, v), f^0(x, v), S^0(x))$ at time $t=0$ and $B^0 = -\tilde{K}S^0/2$. By using $S_1 = S_1^0$, $B_1 = B_1^0$ and $\int f_1 dv = \int f_1^0 dv$, we reformulate (C5) as

$$\partial_t \begin{pmatrix} f_0 \\ f_1 \end{pmatrix} - H \frac{\partial B_1^0}{\partial x} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \partial_v \begin{pmatrix} f_0 \\ f_1 \end{pmatrix} = 0, \tag{C6}$$

$$\partial_t \begin{pmatrix} f_2 \\ f_3 \end{pmatrix} + B_1^0 J \begin{pmatrix} f_2 \\ f_3 \end{pmatrix} = 0, \tag{C7}$$

$$\partial_t \begin{pmatrix} S_2 \\ S_3 \end{pmatrix} - \left( \frac{\tilde{K}}{4} \int f_1^0 dv + A \partial_x^2 S_1^0 \right) J \begin{pmatrix} S_2 \\ S_3 \end{pmatrix} = 0, \tag{C8}$$

where $J$ denotes the symplectic matrix

$$J = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}. \tag{C9}$$

By the eigen-decomposition

$$\begin{pmatrix} \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \tag{C10}$$
(C6) can diagonalized to get two transport equations that can be solved exactly in time

$$ \partial_t \begin{pmatrix} \frac{1}{2}f_0 + \frac{1}{2}f_1 \\ \frac{1}{2}f_0 - \frac{1}{2}f_1 \end{pmatrix} - H \frac{\partial B_1^0}{\partial x} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \partial_v \begin{pmatrix} \frac{1}{2}f_0 + \frac{1}{2}f_1 \\ \frac{1}{2}f_0 - \frac{1}{2}f_1 \end{pmatrix} = 0. \quad (C11) $$

The exact solution for (C7) is

$$ \begin{pmatrix} f_2 \\ f_3 \end{pmatrix} (t, x, v) = \exp(-B_1^0 \mathbf{J} t) \begin{pmatrix} f_2^0(x, v) \\ f_3^0(x, v) \end{pmatrix}, \quad \text{with } \exp(\mathbf{J} s) = \begin{pmatrix} \cos(s) & \sin(s) \\ -\sin(s) & \cos(s) \end{pmatrix}. \quad (C12) $$

Similarly, we can get the exact solution for last system (C8)

$$ \begin{pmatrix} S_2 \\ S_3 \end{pmatrix} (t, x) = \exp \left( \left( \frac{\tilde{K}}{4} \int f_1^0 dv + A \partial_x^2 S_1^0 \right) \mathbf{J} t \right) \begin{pmatrix} S_2^0(x) \\ S_3^0(x) \end{pmatrix}. \quad (C13) $$

### C.4. Subsystem for $\mathcal{H}_{S_2}$

The subsystem $\partial \mathcal{Z}/\partial t = \{\mathcal{Z}, \mathcal{H}_{S_2}\}$ associated to $\mathcal{H}_{S_2} = H \int f_2 B_2 dx dv + A H \int (\partial S_2 / \partial x)^2 dx$ is

$$ \left. \begin{aligned} \frac{\partial f_0}{\partial t} &= \{f_0, \mathcal{H}_{S_2}\} = H \frac{\partial B_2}{\partial x} \frac{\partial f_2}{\partial v} \\ \frac{\partial f_1}{\partial t} &= \{f_1, \mathcal{H}_{S_2}\} = B_2 f_3 \\ \frac{\partial f_2}{\partial t} &= \{f_2, \mathcal{H}_{S_2}\} = H \frac{\partial B_2}{\partial x} \frac{\partial f_0}{\partial v} \\ \frac{\partial f_3}{\partial t} &= \{f_3, \mathcal{H}_{S_2}\} = -B_2 f_1 \end{aligned} \right\} , \quad (C14) \\ \frac{\partial S_1}{\partial t} = \{S_1, \mathcal{H}_{S_2}\} = -\frac{\tilde{K}}{4} S_3 \int f_2 dv - A S_3 \partial_x^2 S_2 \\ \frac{\partial S_2}{\partial t} = \{S_2, \mathcal{H}_{S_2}\} = 0 \\ \frac{\partial S_3}{\partial t} = \{S_3, \mathcal{H}_{S_2}\} = \frac{\tilde{K}}{4} S_1 \int f_2 dv + A S_1 \partial_x^2 S_2 \right\} $$

with the initial value $(f_0^0(x, v), f_1^0(x, v), S^0(x))$ at time $t=0$ and $B^0 = -\tilde{K}S^0/2$. This subsystem is very similar to the $\mathcal{H}_{S_1}$ one, and hence, as was done previously, we reformulate the equations by using $S_2 = S_2^0$, $B_2 = B_2^0$ and $\int f_2 dv = \int f_2^0 dv$,

$$ \partial_t \begin{pmatrix} f_0 \\ f_2 \end{pmatrix} - H \frac{\partial B_2^0}{\partial x} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \partial_v \begin{pmatrix} f_0 \\ f_2 \end{pmatrix} = 0, \quad (C15) $$

$$ \partial_t \begin{pmatrix} f_1 \\ f_3 \end{pmatrix} - B_2^0 \mathbf{J} \begin{pmatrix} f_1 \\ f_3 \end{pmatrix} = 0, \quad (C16) $$

$$ \partial_t \begin{pmatrix} S_1 \\ S_3 \end{pmatrix} + \left( \frac{\tilde{K}}{4} \int f_2^0 dv + A \partial_x^2 S_2^0 \right) \mathbf{J} \begin{pmatrix} S_1 \\ S_3 \end{pmatrix} = 0. \quad (C17) $$
As in the step $\mathcal{H}_{s_1}$, we have two transport equations from (C15) that can be solved exactly:

$$ \partial_t \begin{pmatrix} \frac{1}{2}f_0 + \frac{1}{2}f_2 \\ \frac{1}{2}f_0 - \frac{1}{2}f_2 \end{pmatrix} - H \frac{\partial B_2^0}{\partial x} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \partial_v \begin{pmatrix} \frac{1}{2}f_0 + \frac{1}{2}f_2 \\ \frac{1}{2}f_0 - \frac{1}{2}f_2 \end{pmatrix} = 0. \quad (C18) $$

Moreover, the exact solutions for the systems (C16) and (C17) are respectively

$$ \begin{pmatrix} f_1 \\ f_3 \end{pmatrix} (t, x, v) = \exp(B_2^0 \mathbf{J} t) \begin{pmatrix} f_1^0(x, v) \\ f_3^0(x, v) \end{pmatrix}, \quad \text{with } \exp(\mathbf{J}s) = \begin{pmatrix} \cos(s) & \sin(s) \\ -\sin(s) & \cos(s) \end{pmatrix}, \quad (C19) $$

and

$$ \begin{pmatrix} S_1 \\ S_3 \end{pmatrix} (t, x) = \exp \left( - \left( \frac{\tilde{K}}{4} \int f_2^0 \, dv + A \partial_x^2 S_2^0 \right) \mathbf{J} t \right) \begin{pmatrix} S_1^0(x) \\ S_3^0(x) \end{pmatrix}. \quad (C20) $$

### C.5. Subsystem for $\mathcal{H}_{s_3}$

The subsystem $\partial \mathcal{Z}/\partial t = \{\mathcal{Z}, \mathcal{H}_{s_3}\}$ associated to $\mathcal{H}_{s_3} = H \int f_3 B_3 \, dx \, dv + AH \int (\partial S_3 / \partial x)^2 \, dx$ is

$$ \left. \begin{aligned} \frac{\partial f_0}{\partial t} &= \{f_0, \mathcal{H}_{s_3}\} = H \frac{\partial B_3}{\partial x} \frac{\partial f_3}{\partial v} \\ \frac{\partial f_1}{\partial t} &= \{f_1, \mathcal{H}_{s_3}\} = -B_3 f_2 \\ \frac{\partial f_2}{\partial t} &= \{f_2, \mathcal{H}_{s_3}\} = B_3 f_1 \\ \frac{\partial f_3}{\partial t} &= \{f_3, \mathcal{H}_{s_3}\} = H \frac{\partial B_3}{\partial x} \frac{\partial f_0}{\partial v} \\ \frac{\partial S_1}{\partial t} &= \{S_1, \mathcal{H}_{s_3}\} = \frac{\tilde{K}}{4} S_2 \int f_3 \, dv + A S_2 \partial_x^2 S_3 \\ \frac{\partial S_2}{\partial t} &= \{S_2, \mathcal{H}_{s_3}\} = -\frac{\tilde{K}}{4} S_1 \int f_3 \, dv - A S_1 \partial_x^2 S_3 \\ \frac{\partial S_3}{\partial t} &= \{S_3, \mathcal{H}_{s_3}\} = 0 \end{aligned} \right\} , \quad (C21) $$

with the initial value $(f_0^0(x, v), f_1^0(x, v), S_1^0(x))$ at time $t=0$. This subsystem is also very similar to the $\mathcal{H}_{s_1}$ one, and hence, as was done previously, we reformulate the equations by using $S_3 = S_3^0$, $B_3 = B_3^0$ and $\int f_3 \, dv = \int f_3^0 \, dv$. The update of $(f_0, f_3)$ is performed by solving the following transport equation:

$$ \partial_t \begin{pmatrix} \frac{1}{2}f_0 + \frac{1}{2}f_3 \\ \frac{1}{2}f_0 - \frac{1}{2}f_3 \end{pmatrix} - H \frac{\partial B_3^0}{\partial x} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \partial_v \begin{pmatrix} \frac{1}{2}f_0 + \frac{1}{2}f_3 \\ \frac{1}{2}f_0 - \frac{1}{2}f_3 \end{pmatrix} = 0. \quad (C22) $$

The exact solution for $(f_1, f_2)$ is

$$ \begin{pmatrix} f_1 \\ f_2 \end{pmatrix} (t, x, v) = \exp(-B_3^0 \mathbf{J} t) \begin{pmatrix} f_1^0(x, v) \\ f_2^0(x, v) \end{pmatrix}, \quad \text{with } \exp(\mathbf{J}s) = \begin{pmatrix} \cos(s) & \sin(s) \\ -\sin(s) & \cos(s) \end{pmatrix}, \quad (C23) $$
and for $(S_1, S_2)$, we have

$$ \left(\frac{S_1}{S_2}\right)(t, x) = \exp\left(\left(\frac{\tilde{K}}{4} \int f_3^0 dv + A \partial_x^2 S_3^0\right) \mathbf{J} t\right) \left(\frac{S_1^0(x)}{S_2^0(x)}\right). \quad (C24) $$

REFERENCES

ASHCROFT, N.W. & MERMIN, N.D. 1976 *Solid State Physics*. Saunders College Publishing.

BATTIATO, M., CARVA, K. & OPPENEER, P.M. 2010 Superdiffusive spin transport as a mechanism of ultrafast demagnetization. *Phys. Rev. Lett.* **105** (2), 027203.

BEAUREPAIRE, E., MERLE, J.-C., DAUNOIS, A. & BIGOT, J.-Y. 1996 Ultrafast spin dynamics in ferromagnetic Nickel. *Phys. Rev. Lett.* **76** (22), 4250.

BIGOT, J.-Y. & VOMIR, M. 2013 Ultrafast magnetization dynamics of nanostructures. *Ann. Phys.* **525** (1–2), 2–30.

BIGOT, J.-Y., VOMIR, M. & BEAUREPAIRE, E. 2009 Coherent ultrafast magnetism induced by femtosecond laser pulses. *Nat. Phys.* **5** (7), 515–520.

BONITZ, M., DORNHEIM, T., MOLDBEKOV, Z.A., ZHANG, S., HAMANN, P., KÄHLERT, H., FILINOV, A., RAMAKRISHNA, K. & VORBERGER, J. 2020 Ab initio simulation of warm dense matter. *Phys. Plasmas* **27** (4), 042710.

BRODIN, G., HOLKUNDKAR, A. & MARKLUND, M. 2013 Particle-in-cell simulations of electron spin effects in plasmas. *J. Plasma Phys.* **79** (4), 377–382.

BRODIN, G., MARKLUND, M., ZAMANIAN, J., ERICSSON, Å. & MANA, P.L. 2008 Effects of the g factor in semiclassical kinetic plasma theory. *Phys. Rev. Lett.* **101** (24), 245002.

BRODIN, G., MARKLUND, M., ZAMANIAN, J. & STEFAN, M. 2011 Spin and magnetization effects in plasmas. *Plasma Phys. Control. Fusion* **53** (7), 074013.

CASAS, F., CROUSEILLES, N., FAOU, E. & MEHRENBERGER, M. 2017 Hamiltonian splitting for the Vlasov–Poisson equations. *Numer. Math.* **135**, 769–801.

CHOI, G.-M., MIN, B.-C., LEE, K.-J. & CAHILL, D.G. 2014 Spin current generated by thermally driven ultrafast demagnetization. *Nat. Commun.* **5** (1), 1–8.

COWLEY, S.C., KULSRUD, R.M. & VALEO, E. 1986 A kinetic equation for spin-polarized plasmas. *Phys. Fluids* **29** (2), 430–441.

CRESTETTO, A., CROUSEILLES, N., LI, Y. & MASSOT, J. 2022 Comparison of high-order Eulerian methods for electron hybrid model. *J. Comput. Phys.* **451**, 110857.

CROUSEILLES, N., EINKEMMER, L. & FAOU, E. 2015 Hamiltonian splitting for the Vlasov–Maxwell equations. *J. Comput. Phys.* **283**, 224–240.

CROUSEILLES, N., HERVIEUX, P.-A., HONG, X. & MANFREDI, G. 2023 Vlasov-Maxwell equations with spin effects. *J. Plasma Phys.* **89** (2), 905890215.

CROUSEILLES, N., HERVIEUX, P.-A., LI, Y., MANFREDI, G. & SUN, Y. 2021 Geometric particle-in-cell methods for the Vlasov-Maxwell equations with spin effects. *J. Plasma Phys.* **87** (3).

CROUSEILLES, N., MEHRENBERGER, M. & SONNENDRÜCKER, E. 2010 Conservative semi-Lagrangian schemes for Vlasov equations. *J. Comput. Phys.* **229** (6), 1927–1953.

EICH, F.G., PITTALIS, S. & VIGNALE, G. 2018 A shortcut to gradient-corrected magnon dispersion: exchange-only case. *Il Nuovo Cimento* **91** (8), 173.

FALK, K. 2018 Experimental methods for warm dense matter research. *High Power Laser Sci. Engng* **6**, e59.

FRIED, B.D. & CONTE, S.D. 1961 *The Plasma Dispersion Function: The Hilbert Transform of the Gaussian*. Academic Press.

HAIRER, E., LUBICH, C. & WANNER, G. 2006 *Geometric Numerical Integration*, 2nd edn. Springer Series in Computational Mathematics, vol. 31. Springer.

HINSCHBERGER, Y. & HERVIEUX, P.-A. 2012 Foldy-Wouthuysen transformation applied to the interaction of an electron with ultrafast electromagnetic fields. *Phys. Lett. A* **376** (6), 813–819.

HURST, J., HERVIEUX, P.-A. & MANFREDI, G. 2017 Phase-space methods for the spin dynamics in condensed matter systems. *Phil. Trans. R. Soc. A Math. Phys. Engng Sci.* **375**, 20160199.
HURST, J., HERVIEUX, P.-A. & MANFREDI, G. 2018 Spin current generation by ultrafast laser pulses in ferromagnetic Nickel films. *Phys. Rev. B* **97**, 014424.

HURST, J., MORANDI, O., MANFREDI, G. & HERVIEUX, P.-A. 2014 Semiclassical Vlasov and fluid models for an electron gas with spin effects. *Eur. Phys. J. D* **68** (6), 176.

KRAVANJA, P. & VAN BAREL, M. 2000 *Zeros of Analytic Functions*, pp. 1–59. Springer.

KRIEGER, K., DEWHURST, J.K., ELLIOTT, P., SHARMA, S. & GROSS, E.K.U. 2015 Laser-induced demagnetization at ultrashort time scales: predictions of TDDFT. *J. Chem. Theory Comput.* **11** (10), 4870–4874.

KRIEGER, K., ELLIOTT, P., MÜLLER, T., SINGH, N., DEWHURST, J.K., GROSS, E.K.U. & SHARMA, S. 2017 Ultrafast demagnetization in bulk versus thin films: an ab initio study. *J. Phys.: Condens. Matter* **29** (22), 224001.

LAKSHMANAN, M. 2011 The fascinating world of the Landau-Lifshitz-Gilbert equation: an overview. *Phil. Trans. R. Soc. A: Math. Phys. Engng Sci.* **369** (1939), 1280–1300.

LI, F., DECYK, V.K., MILLER, K.G., TABLEMAN, A., TSUNG, F.S., VRANIC, M., FONSECA, R.A. & MORI, W.B. 2021 Accurately simulating nine-dimensional phase space of relativistic particles in strong fields. *J. Comput. Phys.* **438**, 110367.

LI, Y. 2023 Energy conserving particle-in-cell methods for relativistic Vlasov-Maxwell equations of laser-plasma interaction. *J. Comput. Phys.* **473**, 111733.

LI, Y., SUN, Y. & CROUSEILLES, N. 2020 Numerical simulations of one laser-plasma model based on Poisson structure. *J. Comput. Phys.* **405**, 109172.

MANFREDI, G., HERVIEUX, P.-A. & CROUSEILLES, N. 2023 Spin effects in ultrafast laser-plasma interactions. *Eur. Phys. J. Spec. Top.* **232** (13), 2277–2283.

MANFREDI, G., HERVIEUX, P.-A. & HURST, J. 2019 Phase-space modeling of solid-state plasmas. *Rev. Mod. Plasma Phys.* **3** (1), 1–55.

MANFREDI, G., HERVIEUX, P.-A., YIN, Y. & CROUSEILLES, N. 2010 *Collective Electron Dynamics in Metallic and Semiconductor Nanostructures*, pp. 1–44. Springer.

MARKLUND, M., ZAMANIAN, J. & BRODIN, G. 2010 Spin kinetic theory-quantum kinetic theory in extended phase space. *Transp. Theory Stat. Phys.* **39** (5-7), 502–523.

MORANDI, O., ZAMANIAN, J., MANFREDI, G. & HERVIEUX, P.-A. 2014 Quantum-relativistic hydrodynamic model for a spin-polarized electron gas interacting with light. *Phys. Rev. E* **90** (1), 013103.

NIE, Z., LI, F., MORALES, F., PATCHKOVSKII, S., SMIRNOVA, O., AN, W., NAMBU, N., MATTEO, D., MARSH, K.A., TSUNG, F., MORI, W.B. & JOSHI, C. 2021 In Situ generation of high-energy spin-polarized electrons in a beam-driven plasma wakefield accelerator. *Phys. Rev. Lett.* **126**, 054801.

SCHELLEKENS, A.J., KUIPER, K.C., DE WIT, R.R.J.C. & KOOPMANS, B. 2014 Ultrafast spin-transfer torque driven by femtosecond pulsed-laser excitation. *Nat. Commun.* **5** (1), 1–7.

SINHA-ROY, R., HURST, J., MANFREDI, G. & HERVIEUX, P.-A. 2020 Driving orbital magnetism in metallic nanoparticles through circularly polarized light: a real-time TDDFT study. *ACS Photonics* **7** (9), 2429–2439.

WU, Y., JI, L., GENG, X., THOMAS, J., BÜSCHER, M., PUKHOV, A., HÜTZEN, A., ZHANG, L., SHEN, B. & LI, R. 2020 Spin filter for polarized electron acceleration in plasma wakefields. *Phys. Rev. Appl.* **13**, 044064.

WU, Y., JI, L., GENG, X., YU, Q., WANG, N., FENG, B., GUO, Z., WANG, W., QIN, C., YAN, X., et al. 2019 Polarized electron-beam acceleration driven by vortex laser pulse. *New J. Phys.* **11**, 073052.

YIN, Y., HERVIEUX, P.-A., JALABERT, R.A., MANFREDI, G., MAURAT, E. & WEINMANN, D. 2009 Spin-dependent dipole excitation in alkali-metal nanoparticles. *Phys. Rev. B* **80** (11), 115416.

YOSHIDA, H. 1990 Construction of higher order symplectic integrators. *Phys. Lett. A* **150**, 262–268.

ZAMANIAN, J., MARKLUND, M. & BRODIN, G. 2010a Scalar quantum kinetic theory for spin-1/2 particles: mean field theory. *New J. Phys.* **12** (4), 043019.

ZAMANIAN, J., STEFAN, M., MARKLUND, M. & BRODIN, G. 2010b From extended phase space dynamics to fluid theory. *Phys. Plasmas* **17** (10), 102109.