# Ultrafast demagnetizing fields from first principles

###### Abstract

We examine the ultrafast demagnetization process of iron-based materials, namely Fe clusters and bulk bcc Fe, with time-dependent spin-density functional theory (TDSDFT). The magnetization continuity equation is reformulated and the torque due to the spin-current divergence is written in terms of an effective time-dependent kinetic magnetic field, an object already introduced in literature before. Its time evolution, as extracted from the TDSDFT simulations, is identified as one of the main sources of the local out-of-equilibrium spin dynamics and plays a major role in the demagnetization process. Such demagnetization is particularly strong in “hot spots” where the kinetic torque is maximized. Finally, we find the rate of demagnetization in Fe to be strongly dependent on the direction of polarization of the exciting electric field and this can be linked to the out of equilibrium distribution of the kinetic field in two comparative cases.

###### pacs:

75.75.+a, 73.63.Rt, 75.60.Jk, 72.70.+m## I Introduction

The search for practical solutions for increasing the speed of manipulation of magnetic bits is essential for the progress of modern information and communication technology. It has been shown that there is an upper limit to the speed of the magnetization switching process when this is driven by a magnetic field Mag1 (); Mag2 (). An increase in power absorption beyond this limit and for higher magnetic field amplitudes push a spin system out of equilibrium into a chaotic behaviour, and the switching speed decreases. For this reason the discovery made by Beaurepaire et al. [Beau96, ] in 1996 that a ferromagnetic Ni film could be demagnetized by a femtosecond optical laser pulse attracted a lot of interest and was the seed to a new field, now called femto-magnetism.

In a standard pump-probe experiment the system is initially excited by an optical pulse (pump) and then the magnetization dynamics is monitored by analysing a second signal (probe) Rasing10 (); Kimel (). Depending on the minimal delay between the pump and the probe, one can analyse the demagnetization process at different timescales and thus observe the dissipation mechanisms active at that particular time. The interpretation of the results is, however, a complicate matter. In general for demagnetization processes observed on a timescale ranging from nanoseconds to picoseconds one considers an empirical three temperature model 3Temp (), where electrons, spins and phonons define three energy baths, all interacting with each other. In contrast, ultrafast spin dynamics, taking place within a few hundreds femtoseconds, is yet not described in terms of a single unified scheme and various models for the demagnetization process have been advanced. These include fully relativistic direct transfer of angular momentum from the light to the spins ultrarel (); ultrarel2 (), dynamical exchange splitting dynEXCH (), electron-magnon spin-flip scattering elmag (), electron-electron spin-flip scattering elel () and laser-generated superdiffusive spin currents superdiff ().

Given the complexity of the problem ab initio methods, resolved in the time domain, provide a valuable tool to probe the microscopic aspects of the ultrafast spin dynamics of real magnetic materials by means of time-dependent simulations. In this work we apply time-dependent spin density-functional-theory (TD-SDFT) TDDFT (); TDDFT2 () in its semi-relativistic, non-collinear, spin-polarized version to analyse the ultrafast laser-induced demagnetization of two ferromagnetic transition metal systems: a Fe cluster (see Fig. 1) and bulk bcc Fe. Recently, within a similar theoretical description, it has been demonstrated that the spin-orbit (SO) interaction plays a central role in the demagnetization process Krieger15 (); Tows (); Zhang (). Furthermore, it was showed by us us () that the laser-induced spin dynamics can be understood as the result of the interplay between the SO coupling potential and an effective magnetic field. The so-called kinetic magnetic field antr1 (); antr2 (), , originates from the presence of non-uniform spin currents in the system. In this work we focus on the anatomy of and we analyze in detail its role in the highly non-equilibrium process of ultrafast demagnetization.

The first formulation of the spin dynamics problem in transition metal systems was given in Refs. [antr1, ; antr2, ] by Antropov and Katsnelson, who laid down the foundation of DFT-based spin dynamics, by deriving a set of equations of motion for the local magnetization vector. In those seminal works the magnetization dynamics was analyzed at the level of the adiabatic local spin-density approximation (ALSDA), but actual applications to real out-of-equilibrium systems were not described. Our purpose is to clarify and quantify, through TDSDFT simulations at the level of the non-collinear ALSDA, the role played by in the laser-induced ultrafast spin dynamics of transition metal ferromagnets.

The paper is divided into four main sections. In Section II we define the various fields that couple to the spins by isolating in the continuity equation only the terms that play a major role in the dynamical process. In Section III we present the results of the calculations for Fe clusters and show that “hot spots” for demagnetization are associated with larger misalignment of the kinetic magnetic field and the local spin density. This becomes more clear through evaluation of material derivatives. A demonstration of the effect of the polarization of the electric field on the rate of demagnetization of Fe is discussed in Section IV. In Section V we show that previous observations for Fe are valid for bulk bcc Fe as well. Finally we conclude. The paper is supplemented with an Appendix where we present a detailed derivation of the spin continuity equation (A).

## Ii Theory

We consider the TDSDFT problem within the ALSDA for a spin-polarized system excited by an electric field pulse. If one neglects second-order contributions arising from the solution of the coupled Maxwell-Schrödinger system of equations, the dynamics will be governed by the usual set of time-dependent Kohn-Sham (KS) equations

(1) |

In Eq. (1) are the KS orbitals and the KS Hamiltonian, , can be expressed in the velocity gauge formulation and the minimal coupling substitution as,

(2) |

where

(3) | |||||

and

(4) |

Here represents the usual non-interacting KS potential and the full non-interacting magnetic field, , consists of the external one, , and the exchange-correlation (XC) magnetic field, . In the equations above is the electron mass, the electron charge, the speed of light, the vector potential associated to the external magnetic field, the spin operator, the Bohr magneton, the electron density and the magnetization density. Then, is decomposed into a Hartree contribution, an XC correlation one, , and into an ionic pseudo-potential . For a fully relativistic, norm-conserving pseudopotential the SO coupling enters into the KS equations in the form SOC1 ()

(5) |

In Eq. (5) the orbital momentum operator associated to the -th atomic center is , while the vectors are the associated set of spherical harmonics centered on that given atomic position. In Eq. (5) defines a generalized space-dependent SO coupling parameter providing a measure of the SO interaction strength close to the atomic site, while includes all the ionic relativistic corrections like the Darwin and the mass correction term. Within the ALSDA and are local functions in time of the electron density and magnetization, which in turn are written in terms of the time-dependent KS orbitals

(6) | |||||

(7) |

Starting from the set of time-dependent KS equations in (1) it is possible to derive an equation of motion for the magnetization, or a spin-continuity equation, in terms of the non-interacting KS observables. This reads

(8) |

where represents the non-interacting KS spin-current rank- tensor

(9) |

and the SO torque contribution reads

(10) |

The KS magnetic field is taken as in Eq. (4), which in absence of an external magnetic field reduces to . In DFT there are a set of zero-force theorems stating that the interaction between the particles cannot generate a net force levy (). In the case of the exchange-correlation magnetic field we have the exact condition , which is satisfied by the ALDA. Combining this equality with the assumption that the currents at the system boundary are negligible allows us to conclude that the only source of global spin loss is the SO coupling torque, , and that the spin lost during the temporal evolution is transferred to the orbital momentum of the system, which in turn is partially damped into the lattice (we consider frozen ions). Hence we have the relation,

(11) |

where the integration extends over the entire volume .

Within the ALDA, the exchange-correlation functional satisfies also a local variant of the zero-torque theorem SpinDFT (), which is not a property of the exact DFT functional Sharma07 (); Eich13 (); Sharma07b (). According to this condition and therefore the exchange-correlation magnetic field cannot contribute, even locally, to the magnetization dynamics. This leads us to conclude that the local magnetization dynamics is solely the result of the interplay between the spin-polarized currents and the SO torque (in reality can still contribute indirectly to the spin dynamics through a dynamical modification of the gap between up and down spin polarized bands, which in turn determines an enhancement of the spin dissipation via the spin orbit coupling channel). In order to elucidate this view further we make use of the hydrodynamical formalism applied to spin systems, which has been already introduced in References [hydro2, ; hydro3, ]. This approach needs to be slightly modified in view of the fact that we are considering an effective Kohn-Sham system and not a set of independent spin particles. In fact, as it was already pointed out in Refs. [antr1, ; antr2, ], Eq. (II) can be written in a different form (the details of the derivation are shown in the appendix A)

(12) |

where a couple of new terms appear. In the equation is a material derivative, represents a single Kohn-Sham state velocity field (see appendix A), and . On the right hand-side of Eq. (II) in addition to the spin-orbit coupling torque, , we have a new term, , that describes the spin dissipation in the system due to the internal motion of the spin currents. It can be interpreted as an effective spin-current divergence object involving only transitions among different Kohn-Sham states [inter-band transitions, see Eq. (44)]. Finally the effective field is given by the sum of two terms, , with exchange-correlation field and defined as [see Eq. (53)]

(13) |

with spin vector field .

Such field has only an instrumental rôle in the equations of motion for the spin density, a very similar expression was already introduced in some previous work. In Ref. [antr1, ] it was expressed in the form , while in Ref. [antr2, ] it appears as . The interpretation of may look quite obscure at a first sight, however, in Ref. [Taka55, ; Taka57, ] it was identified as a possible source of spin wave excitations in the form of a spin-spin interaction potential.

In order to clarify this point, let us consider the Heisenberg interaction between two spins centered on atoms placed at a distance . We can assume naively, but reasonably, that the spin-spin interaction among the two spin distributions, computed at an arbitrary point in space could be expressed in the following form

(14) |

where it is more convenient for us to employ a spin field, , which describes the spin distribution in space, instead of an atom localized spin vector. Hence, defines an effective single-particle Hamiltonian. By averaging over the number of electrons in the entire space we obtain

(15) |

Then, by expanding the spin density in Taylor series up to second order in the distance and by neglecting the zeroth-order contribution (we focus our attention on the non-local term appearing in the expansion) after some straightforward rearrangement we arrive at

(16) |

which in turn becomes

(17) |

Finally, by considering a sufficiently large integration volume, the use of the divergence theorem allows to neglect all the boundary terms with consequent final expression

(18) |

which remarkable resembles the result in Eq. (13) for the kinetic magnetic field. We can therefore tentatively interpret as an effective mean-field internal magnetic field, which plays a rôle in coupling the spins at different locations in the system in the spirit of the Heisenberg spin-spin interaction.

## Iii Analyzing spin dynamics from TDSDFT simulations in Fe cluster

Here we present the results of TDSDFT calculations, performed with the Octopus code Octop1 (), where we simulate the ultrafast demagnetization process in iron-based ferromagnetic systems. In all those, at time the system is in its ground state. Then we apply an intense electric field pulse with a duration of less than fs, which initiates the dynamics. The pseudo-potentials for Fe used in the calculations are fully relativistic, norm-conserving and are generated using a Multi-Reference-Pseudo-Potential (MRPP) scheme MRPP () at the level implemented in APE APE1 (); APE2 (), which takes directly into account the semi-core states. For the XC functional we employ the ALSDA with parameterization from Perdew and Wang Wang92 (). Our simulations then consist in evolving in time the KS wave functions, i.e. in solving numerically the set of equations (1). The results are then interpreted through the magnetization continuity equation (II).

In Fig. 1 the extracted magnetization dynamics of a Fe magnetic cluster is presented. We use the LDA ground-state geometry of Fe as extracted from Ref. [Dieguez, ,gutsev, ] for which we reproduce the reported therein spin state . The nuclei are kept stationary during the dynamics. In panel (c) we observe that the total loss of the component of the total magnetization, , is exactly equal to the variation in value of its module, , since the global non-collinear contribution is negligible. This indicates that the spin is not exchanged globally between the different components of the magnetization vector, but, according to Eq. (11), it is, at least, partially transfered into the orbital momentum of the system. We note that due to the electrostatic interactions with the nuclei and due to the interaction with the laser field the rotational invariance of the electronic system is broken and the total orbital momentum is not conserved.

In Fig. 1(b) we observe that the average kinetic magnetic field (over the entire simulation box, for ) is comparable in magnitude to the exchange component. At the same time, shows a much more oscillatory behaviour compared to . In particular, While evolves smoothly in time following the action of the optical excitation, presents an abrupt variation at the on-set of the electrical pulse. This is due to the fact that the laser pulse directly excites currents, through the term , which, in turn, produces a modification of the gradients of the charge/spin density, even on a global scale since they are not conserved. Thus we observe huge variations of . can also oscillate very strongly locally, following the temporal variation of the densities, but when we measure these oscillations are averaged out given that the densities are approximately conserved over the entire simulation box. During the action of the pulse we see a tendency of the two fields to compensate each other, an effect strongly resembling the Lenz law. After the pulse, continues to oscillate dramatically with its average value that slowly increases. In contrast decreases (in absolute value) due the net dissipation of spin angular momentum.

Moving from an analysis of global quantities to probing locally the spin dynamics, in Fig. 1(d) we compare the magnetization and the electron density around the atomic site at the tip of the cluster (see inset of Fig. 1(a) for the numbering labels of all the cluster atoms). We define local magnetization and charge associated to the particular atomic site as

(19) |

where the integration volume is a sphere of radius centered at site . Our results show that the loss of is not taking place just during the action of the external pulse, but it is rather distributed over the entire time evolution. This suggests that the spin-sink mechanism is not directly related to the coupling of the system to the laser field, but is rather intrinsic to the electron dynamics following the pulse. Furthermore, close to the atomic site, the temporal variation of the charge, , is much smaller in magnitude and smoother than that of . In addition for long times settles close to an average value, while continues to decrease. Hence the long-term spin dynamics is not the result of a net charge displacement from the region close to the ions to the interstitial space. These observations are valid for all the atomic sites in the cluster.

If we now consider the continuity equation for the electron density (see the Appendix A for further explanations)

(20) |

where is the material derivative of the electron density

(21) |

From Fig. 1(d) we observe that during the action of the pulse the density variation in the vicinity of the atoms appears to be very small compared to the magnetization variation. We can therefore safely assume that in this spatial region , with at the same time . From these considerations we deduce that is a reasonably good approximation for the velocity field in the vicinity of the atoms (this does not imply that the velocity field is exactly zero, but only that its effect on the spin dynamics in this particular case is negligible). The same argument is valid also for the state resolved density , given that , the contribution of the local time derivative of the Kohn-Sham state density can be neglected. By applying the latter into Eq. (II) we finally obtain a relation that could be considered approximately valid in this spatial region of the simulation box,

(22) |

where the contribution to the spin dynamics due to the velocity field term has been neglected. Note that here we have also used the condition , that is consequential to the local density approximation. In addition, the decay of , during the evolution, is not so relevant to justify a dynamical modification of the gap between up and down spin states.

In Fig. 2 we compare the behaviour of the kinetic field and of the local magnetization at two atomic sites, respectively (one of the atoms in the base plane of the bi-pyramid) and (an atom at one of the apexes). It can be seen from panel (a) that these two sites present different rates of demagnetization. In particular, at site the spin decay is considerably more prominent with respect to that observed at site . In contrast the fluctuations in are significantly more pronounced for site than for site . This can be understood from the fact that we have chosen here an electric pulse with polarization vector in the basal plane of the bi-pyramid. As such, the charge fluctuations for the atoms in the basal plane are expected to be much larger than those of the apical atoms. Finally, we note that follows similar qualitative trends as [see Fig. 2(b)]. In fact, the average change following the excitation pulse is larger for site (the one experiencing the larger demagnetization), but the fluctuations are more pronounced for site (the one experiencing the larger fluctuations in ).

The correlation between the kinetic field and the magnetization loss is also rather evident in Fig. 3. There the time-averaged variations in the -component of the two fields and are clearly comparable in magnitude and localized over the same regions of the simulation box. This demonstrates that the kinetic field can be considered as the main force driving the non-collinearity during the spin evolution. The fact that the contrast is stronger at the apex atoms (“hot spots” for demagnetization) agrees with Fig. 2(a), while the dipole-type patterns indicate how the longitudinal spin decays preserving global collinearity. The correlation between the components of and is not as evident as that for the transverse component . This is due to the fact that the and components of the field are much smaller compared to the one. Furthermore, the contribution to the spin dynamics along of the SO coupling, together with the internal dissipative term due to the spin currents, cannot be neglected.

In order to quantify the local non-collinearity we examine the evolution of the misalignment angle, , between the -axis and the direction of the magnetic fields (averaged over spheres). It can be seen in Fig. 4 that at site the averaged kinetic field and the local spin deflect very little from the quantization axis and remain rather parallel to each other. The angle that forms with the magnetization direction (the direction) is substantially negligible. Instead, at site , shows a significant deflection from the -axis after the first fs of the evolution and so does the spin, without the two being parallel to each other. It is important to notice that the angle between magnetization and starts to grow only after the action of the pulse. These results for atom and are representative for all the other sites in the base plane or outside of it, respectively. The sites located in the plane, where is mostly collinear, lose less magnetization with respect to the ones at the apices where, instead, the kinetic field shows a significant deflection from the magnetization axis and provides additional torque driving further demagnetization.

Analogous conclusions arise from the introduction of the concept of parallel transport, commonly used in differential geometry. This requires a proper definition of the covariant derivative obtained by comparing not with , but with the value that the spin vector would have if it was translated from to while keeping the axes in the isospin space fixed,

(23) |

The connection field provides a measure of the amount of non collinearity accumulated in the translation of the spin vector from to

By using the previous expression to rewrite the first and second order spatial derivatives, the kinetic field of Eq. (13) can be divided in two components

(24) |

Here we have introduced

(25) |

which has no effects on the dynamics, having locally the same direction of the spin vector by construction, and

(26) |

According to Fig. 4, in the case of Fe, the direction in the isospin space of the connection tensor for every component can be considered in first approximation orthogonal to the direction of the spin vector , since its component along is considerably smaller than the components along and . From this we could assume and we obtain the following simplified expression for

(27) |

This represents the part of which gives rise to a non-zero torque in Eq. (22).

## Iv Directionality of the demagnetization in Fe cluster

In the previous sections we have revisited the concept of kinetic field, its derivation within DFT and its properties as a major source of torque for the spin dynamics within ALSDA. We have provided supporting evidence for the latter from TDSDFT calculations of the ultrafast demagnetizing Fe cluster under the effect of a single fs electric field pulse. Despite the conceptual clarity of as an instrumental object, very little useful physical intuition can be drawn from its definition in Eq. (13). Clearly, it is an intrinsic dynamic field that depends on the spin texture and its response to the external stimuli. It also feeds back into the dynamics of this same spin density, clearly a non-linear process. In this section we seek to extend the evidential base for the connection between the torque due to and the rate of demagnetization. Together with that we report a situation, where the direction of the polarization vector of the electric field of the laser pulse alone has a significant effect on the demagnetization of a material (the Fe cluster).

Figure 5 shows a comparison between two simulations differing only by the direction (but notably not the magnitude) of the electric field applied. In one case this is in the -direction, which is oriented along the slightly longer side of the base of the bi-pyramid Dieguez (), and in the other simulation it is along , the direction connecting the two apex atoms and . The case shows nearly 3 times faster demagnetization compared to the one. Evidently, in the former situation more energy is absorbed by the cluster (an excess of about ), we show a comparison of the total energy shift due to the pulse in Fig. 5(b).

From panel (c) we observe that the amount of non-collinearity enclosed in a relatively small radius around the apex atoms does not significatively change in the two cases, suggesting that the intra-site non-collinear component of the spin vector is already present in the ground-state configuration. What really differs is the amount of inter-sites non collinearity concentrated in the out-of-plane region. In the case, if we examine the angle between and averaged in time over the entire simulation along one particular cross-section plane (vertical through the base diagonal of the cluster as depicted in the inset), one can observe, in particular in the out-of-plane interstitial regions, a significant larger amount of non-collinearity, of the spin vector density, accumulating in the case of faster demagnetization.

Notably, there is a change in the symmetry between panels (c) and (e) in Fig. 5 – in both cases no significant spin non-collinearity arises in the plane parallel to the field connecting the atomic centers. It is, however, important to notice that the amount of intra-site spin non-collinearity for the in-plane atoms is strongly dependent on the polarization direction of the applied laser field. In fact, in the case of the temporal averaged angle between and appears much higher with respect to the one computed in the case. This quantity is mostly averaged out by the integration procedure but it is clearly visible in the panel (c) of the figure.

Although is not the only torque generator and the SO contribution is significant too, the former plays a rôle in deflecting the spins in the system in a manner that correlates with the rate of global demagnetization. Importantly, our simulations clearly show that the demagnetization process is very anisotropic and particular directions of the exciting electric field may enhance the rate of demagnetization (the study of these effects is beyond the scope of this paper and it will be explored in more details in a coming publication).

## V Demagnetization of bcc Fe

Finally, we present results of analogous simulations in bulk materials, namely in bcc Fe, with the aim of demonstrating the qualitative universal rôle played by in the ultrafast demagnetization process. We consider bcc Fe in its ferromagnetic phase with total spin in the unit cell and with 2 atoms in it.

We employ a lattice parameter , with a -points grid. In Fig. 6(c) we show the demagnetization rate of the single unit cell after it has been excited with the electric field pulse [panel (a)]. The green curve represents the demagnetization computed inside the unit cell and resembles in shape the sum of the magnetization variation, , calculated in the vicinity of the two Fe atoms, even if it is different in magnitude. This suggests that a large amount of spin is driven outside from the atomic integration region during the evolution. Similarly to the cluster, the dynamics of the onsite magnetization can be described in terms of a two-step process with an initial fast decay during the action of the external pulse, followed by a slower and noisy decrease in magnitude. The first fast decay may be attributed to the effect of the SO enhanced by the collapse of the effective field following the action of the laser pulse. In Fig. 6(b) the collapse after the first of the component of the effective field is quite clear, even if it appears to be more pronounced for the kinetic field with respect to the exchange field . Similarly to the case of Fe the role played by is dominant during the action of the pulse, but after this initial phase the dynamics is dominated by intra-band transitions and the interplay between the spin-orbit coupling and the effective field becomes dominant.

Fig. 6(d) shows the evolution of the non-collinearity of the spin vector, , and of the kinetic field, . The level of correlation among the two quantities confirms the importance of the kinetic field in the evolution of the spin non-collinearity. The long tail of spin dissipation may be explained in terms of intra-band spin-up/spin-down transitions through an Elliott-Yafet type of mechanism triggered by the scattering with the effective field ,

(28) |

represents the transition amplitude between two states with different vector and in presence of SO with different mixing of up and down spin components.

## Vi Conclusions

In conclusion, we remark the central result of our work, namely that the equation of motion for the spin dynamics within the ALSDA of TDSDFT [see Eq. (II)] can be rewritten in the form of Eq. (II), by using a formalism borrowed from magneto-hydrodynamics. Subsequently we have analyzed the properties of the so-called kinetic magnetic field and its rôle in the ultrafast demagnetization process in two different systems: a ferromagnetic Fe cluster and bulk bcc Fe. The rôle of this field is particularly significant for processes far from equilibrium, such as the ultrafast demagnetization observed in transition metals.

In both the systems studied the spin dynamics is the result of the interplay between the SO coupling and , which, in general, is strongly coupled to the external pulse and highly non-uniform in space. We have shown that the spin loss locally correlates with . Through the concept of parallel transport and the definition of a connection tensor field , we have gained further insight into the evolution of the spin texture. As describes the degree of spin rotation per infinitesimal spatial translation, it also provides a measure for the misalignment between the kinetic field and the spin texture. The regions with higher correspond to stronger local demagnetization.

Finally, the effect of the direction of the polarization vector of the electric field pulse has been studied for Fe. We have found that clusters will demagnetize about twice as fast, if the polarization vector is in the base plane and not vertical (through the apex atoms). Our analysis has shown a significant increase in the non-collinearity between and the spin density in the fast demagnetizing case. Such anisotropy, due to the electric dipole matrix elements for the valence electrons, is likely to occur in crystalline systems as well. During the application of the laser pulse, the rise of spin non-collinearity may be enhanced by the particular polarization direction of the laser pulse through the spin orbit coupling and this effect combined with the collapse of the kinetic field may explain the initial spin loss. However, in both Fe and Fe bcc the magnetization loss is more prominent after that the laser pulse has been set to zero. During this second phase of spin dissipation we need to distinguish between the spin decay observed in bcc Fe due to intra-band transitions among states with different spin up/down mixing and the spin dynamics observed in correspondence of the apex atoms in Fe that is, instead, driven by and directly related to the onsite intrinsic spin non-collinearity near the atomic sites.

###### Acknowledgements.

This work has been funded by the European Commission project CRONOS (grant no. 280879) and by Science Foundation Ireland (grant No. 14/IA/2624). We gratefully acknowledge the DJEI/DES/SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities. We also thank Prof. E.K.U. Gross for valuable discussions.## Appendix A Derivation of the hydrodynamic continuity equation

In this appendix we show in some detail how Eq. (II) can be obtained by starting from the standard TDSDFT continuity equation (II). Here we follow the hydrodynamical formalism of quantum mechanics, where a single particle with spin is considered equivalent to a non-linear vector field. In this type of hydrodynamics the quantum effects are separated as non-linear terms and are described through effective quantum potentials (see Ref. [Taka55, ]).

The formalism is based on the assumption that it is possible to describe the dynamical evolution of a single particle immersed in an external vector potential, , through the so-called Madelung decomposition (see Ref. [Holl93, ]) of the system wave function, in which the amplitude is translated into the probability density and the gradient of the phase determines the velocity field. A hydrodynamical description of the wave function was also obtained in Ref. [Schon54, ] starting from the ordinary interpretation of quantum mechanics and by introducing an operator for the charge density and the current density.

The formalism was also later extended to the semi-relativistic description (Pauli approximation) of a single particle in an external electro-magnetic field in Ref. [hydro2, ]. However, while in all the previous studies the main objective was to derive the single-particle dynamics of the spin plasma, only recently the study of the collective dynamical properties of the quantum plasma started to attract some interests (see Ref. [hydro3, ] and [Brod11, ]).

In deriving the Eq. (II) for the spin density in the Kohn-Sham system we introduce also of the electron density, , and the velocity field, . The equation of motion for the velocity field are not explicitly written here for the reasons explained in section III.

The electron density is written in terms of the Kohn-Sham wave functions as

(29) |

while the spin density is

(30) |

and the covariant velocity field appears as

(31) |

By making use of the Kohn-Sham equations (1) the charge continuity equation can be written straightforwardly as

(32) |

while for spin it is written in terms of and as

(33) |

By evaluating explicitly the spatial partial derivative of the spin vector and by multiplying it with the component we obtain the following equality

(34) |

We now need to focus our attention on the first term on the right-hand-side of Eq. (34), by using the following notation

(35) |

that leads to

(36) |

The anti-symmetric part, , of the tensor , defined as may be written as

(37) |

By making use of the following relation between Pauli matrices [see Ref. Taka55, ]

(38) |

we obtain the following final expression for that can be splitted in two parts as follows

(39) |

where we have introduced the tensor

(40) |

The procedure that we have followed up to now is formally exact. Then, in order to simplify the previous expression we substitute the Kohn-Sham ratio with its average over the various occupied states . From the fact that, to a good degree of approximation, with total number of particles in the system, we will consider from now on to be spatially homogeneous and constant in time. Then Eq. (39) becomes

(41) |

In order to simplify the formalism we employ the notation . Immediately from Eq. (A) follows that

(42) |

where and define, respectively, the single Kohn-Sham state magnetization and velocity field. By employing the properties of the Levi-Civita tensor we have

(43) |

where we have introduced the new tensor quantity

(44) |

Here and the other many-particle objects are defined as

(45) | |||||

(46) | |||||

(47) |

Finally, from Eq. (43) the divergence of the spin current tensor may be rewritten as

(48) |

Then, by substituting Eq. (48) into Eq. (33) we obtain

(49) |

Finally, by decomposing the magnetization into its single-particle components, , we can define the magnetization material derivative as follows

(50) |

with the spin continuity equation that becomes

(51) |

or

(52) |

where we have introduced an effective magnetic field

(53) |

The continuity equation for the electron density instead follows immediately from Eq. (32) through the definition of velocity field

(54) |

It should be noted that the kinetic field, , written in Eq. (53) is expressed in term of the density and spin density, that are observables of the many-body system. This means that the kinetic field obtained within the Kohn-Sham formalism is identical to its many-body counterpart.

## References

- (1) I. Tudosa, C. Stamm, A.B. Kashuba, F. King, H.C. Siegmann, J. Stohr, G. Ju, B. Lu and D. Weller, Nature 428, 831 (2004).
- (2) Th. Gerrits, H.A.M. Van den Berg, J. Hohlfeld, L. Bär and Th. Rasing, Nature 418, 509 (2002).
- (3) E. Beaurepaire, J.-C. Merle, A. Daunois and J.-Y. Bigot, Phys. Rev. Lett. 76, 4250 (1996).
- (4) A.V. Kimel, A. Kirilyuk and Th. Rasing, Laser & Photon. Rev. 1, 275 (2007).
- (5) A. Kirilyuk, A.V. Kimel and Th. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
- (6) B. Koopmans, G. Malinowski, F. Dalla Longa, D. Steiauf, M. Fähnle, T. Roth, M. Cinchetti and M. Aeschlimann, Nature Materials 9, 259 (2010).
- (7) Y. Hinschberger and P.-A. Hervieux, Phys. Lett. A 376, 813 (2012).
- (8) H. Vonesch and J.Y. Bigot, Phys. Rev. B 85, 180407(R) (2012).
- (9) B.Y. Mueller, A. Baral, S. Vollmar, M. Cinchetti, M. Aeschlimann, H.C. Schneider and B. Rethfeld, Phys. Rev. Lett. 111, 167204 (2013).
- (10) E. Carpene, E. Mancini, C. Dallera, M. Brenna, E. Puppin and S. De Silvestri, Phys. Rev. B 78, 174422 (2008).
- (11) M. Krauss, T. Roth, S. Alebrand, D. Steil, M. Cinchetti, M. Aeschlimann and H. C.Schneider, Phys, Rev. B 80, 180407(R) (2009).
- (12) M. Battiato, K. Carva and P.M. Oppeneer, Phys. Rev. B 86 024404 (2012).
- (13) E. Runge and E.K.U. Gross, Phys. Rev. Lett. 52, 997 (1984).
- (14) E.K.U. Gross and W. Kohn, Adv. Quantum Chem. 21, 255 (1990).
- (15) K. Krieger. J.K. Dewhurst, P. Elliott, S. Sharma and E.K.U. Gross, J. Chem. Theory Comput. 11, 4870 (2015).
- (16) W. Töws and G.M. Pastor, Phys. Rev. Lett. 115, 217204 (2015).
- (17) G.P. Zhang, Phys. Rev. Lett. 101, 187203 (2008).
- (18) M. Stamenova, J. Simoni and S. Sanvito, Phys. Rev. B 94, 014423 (2016).
- (19) V. Antropov, J. Appl. Phys. 97, 10A704 (2005).
- (20) M.I. Katsnelson and V.P. Antropov, Phys. Rev. B 67, 140406(R) (2003).
- (21) L. Fernandez-Seivane, M.A. Oliveira, S. Sanvito and J. Ferrer, J. Phys: Condens. Matter 18, 7999 (2006).
- (22) M. Levy and J.P. Perdew, Phys. Rev. A 32, 2010 (1985).
- (23) K. Capelle, G. Vignale and B.L. Györffy, Phys. Rev. Lett. 87, 206403 (2001).
- (24) T. Koide, Phys. Rev. C 87, 034902 (2013).
- (25) T. Takabayasi, Prog. Theor. Phys. 14, 283 (1955).
- (26) T. Takabayasi and J.-P. Vigier, Prog. Theor. Phys. 18, 573 (1957).
- (27) P.R. Holland, The quantum theory of motion, (Cambridge: Cambridge University Press), (1993).
- (28) M. Marklund and G. Brodin, Phys. Rev. Lett. 98, 025001 (2007).
- (29) G. Brodin, M. Marklund, J. Zamanian and M. Stefan, Plasma Phys. Control. Fusion 53, 074013 (2011).
- (30) M. Schönberg, Il Nuovo Cimento Vol.XII, N. (1954).
- (31) A. Castro, H. Appel, M. Oliveira, C.A. Rozzi, et al., Phys. Stat. Sol. B 243, 2465 (2006).
- (32) C.L. Reis, J.M. Pacheco and J.L. Martins, Phys. Rev. B 68, 155111 (2003).
- (33) M.J.T. Oliveira and F. Nogueira, Comp. Phys. Comm. 178, 524 (2008).
- (34) M.A.L. Marques, M.J.T. Oliveira and T. Burnus, Comp. Phys. Comm. 183, 2272 (2012).
- (35) J. Zhu, X.W. Wang and S.G. Louie, Phys. Rev. B 45, 8887 (1992).
- (36) O. Dieguez, M.M.G. Alemany, C. Ray, P. Ordejon and C.W. Beuschlicher, Phys. Rev. B 63, 205407 (2001).
- (37) G. Gutsev et al., J. Chem. Phys. 107, 7013 (2003).
- (38) S. Sharma, J.K. Dewhurst, C. Ambrosch-Draxl et al., Phys. Rev. Lett. 98, 196405 (2007).
- (39) F.G. Eich and E. K.U.Gross, Phys. Rev. Lett. 111, 156401 (2013).
- (40) S. Sharma, S. Pittalis, S. Kurth et Al., Phys. Rev. B 76, 100401 (2007).