On the conservation of the angular momentum in ultrafast spin dynamics
OOn the conservation of the angular momentum in ultrafast spin dynamics
Jacopo Simoni ∗ and Stefano Sanvito Molecular Foundry, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA School of Physics, Trinity College, Dublin 2, Ireland
The total angular momentum of a close system is a conserved quantity, which should remainconstant in time for any excitation experiment once the pumping signal has extinguished. Suchconservation, however, is never satisfied in practice in any real-time first principles description of thedemagnetization process. Furthermore, there is a growing experimental evidence that the same takesplace in experiments. The missing angular momentum is usually associated to lattice vibrations,which are not measured experimentally and are never considered in real-time simulations. Here wecritically analyse the issue and conclude that current state-of-the-art simulations violate angularmomentum conservation already at the electronic level of description. This shortcoming originatesfrom an oversimplified description of the spin-orbit coupling, which includes atomic contributionsbut neglects completely that of itinerant electrons. We corroborate our findings with time-dependentsimulations using model tight-binding Hamiltonians, and show that indeed such conservation canbe re-introduced by an appropriate choice of spin-orbit coupling. The consequences of our findingson recent experiments are also discussed.
PACS numbers: 75.75.+a, 73.63.Rt, 75.60.Jk, 72.70.+m
I. INTRODUCTION
Thanks to recent advances in femtosecond laser tech-nologies the possibility of achieving control over the mag-netization dynamics at time scales of the order of 100 fs isnow within our reach. After the discovery of the ultrafastoptical demagnetization of a Nickel film irradiated witha sub-picosecond laser pulse in 1996 , several additionalexperiments have shown a rich variety of laser-inducedphenomena in magnetic compounds, including ultrafastdemagnetization , spin-reorientation , and the mod-ifications of the magnetic structure .Such race towards the control of the magnetization dy-namics at the femtosecond time scale is driven by boththe fundamental physical investigation of the underlyingmechanism leading to ultrafast spin dynamics and tech-nological innovation in the fields of high-speed magneticrecording and spin electronics .The now-standard pump-probe experimental protocolrepresents a valid technique for probing the magnetiza-tion dynamics of a sample at short time scales. Here,the system is first excited by the application of an opti-cal pulse (pump), and then the magnetization dynamicsis reconstructed by probing it with a second small per-turbing signal (probe) . Depending on the time de-lay between the pump and the probe one can accuratelytrace the magnetization trajectory with a few femtosec-onds resolution. Then, the results are typically ratio-nalised by choosing among the most relevant dissipationmechanisms that are at play at the given time scale. Inparticular, the different types of magnetization dynamicshave been usually classified within two main groups.Fast magnetization processes are those observed ona time scale ranging between a few nanoseconds to a hun-dred picoseconds. These are usually described by meansof the so-called three-temperatures model, where elec-trons, phonons and spins form three different thermalbaths that are brought out of equilibrium by the appli-cation of the laser pulse. The subsystems are able toexchange energy and they thermalise to achieve a final equilibrium state.An ultrafast process, instead, is active at a muchshorter time scale, approximately of the order of 100 fs.In this case there is much less general agreement onthe ultimate cause behind the observed spin dynam-ics. In recent years different magnetization dynamicsmodels have been proposed in literature, including fullyrelativistic direct transfer of angular momentum fromthe laser field to the spins , electron-magnon spin-flip scattering , electron-electron spin-flip scattering ,Elliott-Yafet mechanism and laser induced superdiffu-sive spin currents . Among all these different schemesonly the last one does not require the spin orbit cou-pling to play a dominant role in the demagnetization pro-cess. Crucially, it was experimentally observed that dur-ing the ultrafast demagnetization the behaviour of theelectronic orbital momentum alone cannot account forthe spin decay rate in magnets . Thus, several theoret-ical works suggested that at least part of the spinmomentum should be directly transferred to the atomicdegrees of freedom. However, the question has not beensettled yet, and there is still some debate on the impor-tance of such angular-momentum transfer channel andits relation to the ultrafast demagnetization process.In our contribution we analyse these issues in greatdetail, in particular focussing on the conservation of thetotal electronic angular momentum. We find that thestandard treatment of the spin-orbit interaction, basedon the sum of atomic contributions, is at the origin ofthe violation of the electronic angular momentum con-servation. Relaxing some of the approximations, and inparticular the assumption that in a multi-atom environ-ment the spin-orbit interaction could still be written asif the atoms were isolated, leads to a more complete ex-pression for the spin-orbit interaction. This allows one tosatisfy the conservation law. Our results are put to thetest with a simple tight-binding model and confirmed bytime-dependent simulations.In section II we will discuss the problem of the totalorbital momentum conservation, in particular, with ref- a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b erence to the ab initio methods that are more commonlyused to simulate these processes. In section III we explainwhy these methods cannot conserve the total system’s an-gular momentum at the present level of development andshow how to modify the spin orbit coupling operator inorder to enforce this conservation law. In section IV welook at some results obtained in the case of a very simpletight binding model simulation for atomic clusters andwe compare the results with the ones obtained by usingthe standard models. In section V we conclude. II. NON CONSERVATION OF THE TOTALANGULAR MOMENTUM IN AB-INITIO SPINDYNAMICS
The general assumption underpinning any currentwork in the field of ultrafast magnetism is that the totalangular momentum of the system, ˆ J , is conserved duringthe dynamics. This means that the following equationholds ∆ (cid:104) ˆ J (cid:105) = ∆ (cid:104) ˆ S (cid:105) + ∆ (cid:104) ˆ L e (cid:105) + ∆ (cid:104) ˆ L atom (cid:105) + ∆ (cid:104) ˆ L ph (cid:105) = 0 , (1)where ˆ S is the spin operator, ˆ L e the electronic orbital mo-mentum operator, ˆ L atom the atomic orbital momentumand ˆ L ph the orbital momentum carried by the electro-magnetic field interacting with the material. In Eq. (1)the symbol (cid:104) ˆ O (cid:105) = (cid:104) ˆ O (cid:105) t represents the expectation valueof the vector operator ˆ O at time t . Under this assumptionthe nuclear spin degrees of freedom are not considered,since they do not contribute appreciably to the entire spindynamics due to the fact that the hyperfine interactionis small compared to the other interactions at play.When simulating the magnetization dynamics in realtime, none of the methods routinely used to interpretthe experiments is able to satisfy Eq. (1). These include ab initio approaches, such as real-time time-dependentdensity functional theory (rtTDDFT) , semi- ab - initio schemes, such as time-dependent tight-binding (TDTB)models or quantum chemistry methods . In gen-eral, this deficiency is not due to a lack of numericalaccuracy, instead it represents a well-known limitationof the underpinning theoretical formalism. The mostpromising among these methods, both in terms of lackof free parameters and ease of computation, is certainlyrtTDDFT . In rtTDDFT only the electronic subsystemis evolved quantum mechanically by solving a set of singleparticle Schr¨odinger-like equations, known as the Kohn-Sham equations . By solving exactly the rtTDDFTproblem it is possible to reproduce the temporal depen-dence of both the electron and the magnetization den-sity. However, the interaction between the electrons, theatoms and the laser field are treated by mean of effectiveexternal scalar potentials. This means that, in practice,the photons and the atomic degrees of freedom are notself-consistently evolved during the dynamics. As a con-sequence, only ∆ (cid:104) ˆ S (cid:105) + ∆ (cid:104) ˆ L e − KS (cid:105) is accessible from theknowledge of the spin and charge densities, where ˆ L e − KS is the electron KS angular momentum . Unfortunately,as one can deduce from Eq. (1), ∆ (cid:104) ˆ S (cid:105) +∆ (cid:104) ˆ L e − KS (cid:105) is not aconstant of motion. Hence, if we insist on this level of de-scription, we will not be able to answer to the question on what is the relevant channel for angular momentum dis-sipation during the demagnetization process. The sameconsiderations are of course valid also for TDTB models,which are based on a similar approach.In the remaining part of this section we will look athow ∆ (cid:104) ˆ S (cid:105) and ∆ (cid:104) ˆ L e (cid:105) evolve in time by solving the set ofrtTDDFT equations at the level of implementation pro-vided by the OCTOPUS code . OCTOPUS expandsthe wave function and the operators in real space over anumerical grid and approximates the electron-ion inter-action with norm-conserving pseudo potentials incor-porating relativistic effects through a spin-orbit couplingterm . In particular, the pseudo potential, ˆ V ps ( r ), hasthe following general form ˆ V ps ( r ) = (cid:88) l l (cid:88) m = − l V ps l ( r ) | l, m (cid:105) (cid:104) l, m | ,V ps l ( r ) = ¯ V ion l ( r ) + V SO l ( r )4 + λV SO l ( r )ˆ L · ˆ S , (2)where ˆ L is the orbital momentum operator of the atomand the expansion is performed over its eigenstates | l, m (cid:105) (spherical harmonics). In Eq. (2) the scalar componentof the potential, ¯ V ion l ( r ), describes the effect of the massshift and the Darwin term, while V SO l ( r ) sets the rangeof the spin-orbit coupling that, for convenience, we canrescale by a factor λ . We start by considering the sim- ( e V f s ) < S z >, 0×V SO < S z >, 1×V SO < S z >, 2×V SO < J z >, 0×V SO < J z >, 1×V SO < J z >, 2×V SO t (au) ( e V f s )
1e 1 < L >, 0×V SO < L >, 1×V SO < L >, 2×V SO < S >, 0×V SO < S >, 1×V SO < S >, 2×V SO ( e V / Å ) E x (t) FIG. 1: (Color online) Time evolution of different observablesfor an isolated Fe atom excited by the application of an ex-ternal laser field E x ( t ) along the x axis (dashed green line,upper panel). The comparison between (cid:104) ˆ S z (cid:105) (solid line) and (cid:104) ˆ J z (cid:105) (dotted line) is shown in the upper panel for differentspin orbit coupling strengths λV SO l ( λ = 0 , , (cid:104) ˆ L (cid:105) (solid line) and ∆ (cid:104) ˆ S (cid:105) (dottedline) obtained with the same values of λ shown in the upperpanel. The units for (cid:104) ˆ S z (cid:105) and (cid:104) ˆ J z (cid:105) are eV · fs (left-hand sidescale), while that for the field is eV/˚A (right hand-side scale). plest possible case of an iron atom isolated in vacuum andexcited by the application of an external laser pulse. Thisexample, of course, has little practical relevance, but ithas important implications for our analysis. In the caseof a single isolated atom, the (cid:104) ˆ L atom (cid:105) term is identicallyzero, so that, after the electric field has extinguished, thequantity (cid:104) ˆ J (cid:105) = (cid:104) ˆ S (cid:105) + (cid:104) ˆ L e (cid:105) is conserved. This, however,is not the case, as one can see from the time dependenttraces of Fig. (1). We find, in fact, that the total orbitalmomentum (cid:104) ˆ J z (cid:105) still oscillates even after the laser pulsehas vanished. Importantly, we notice a difference be-tween the calculations performed with and without spinorbit coupling. If the spin orbit strength, λ , is set to zero(red lines) (cid:104) ˆ S z (cid:105) does not change in time, whereas (cid:104) ˆ J z (cid:105) is affected only during the application of the pulse. For λ = 1 ,
2, instead, (cid:104) ˆ J z (cid:105) is not conserved even after the ap-plication of the pulse with the fluctuations getting largeras the spin-orbit strength increases. The lower panel ofFig. (1) confirms these findings by showing that (cid:104) ˆ L (cid:105) re-mains constant over long times only in the absence ofspin-orbit coupling. In contrast, for λ (cid:54) = 0 the modulesquared of the orbital momentum is not conserved dur-ing the evolution, even at times where the laser pulse isnot present. The next level of complexity is achieved ( e V f s ) < S z >, fixed atoms< S z >, Ehrenfest< L z >, fixed atoms < L z >, Ehrenfest L atz t (au) ( e V f s )
1e 1 < L >, fixed atoms< L >, Ehrenfest < S >, fixed atoms< S >, Ehrenfest ( e V / Å ) E x (t) FIG. 2: (Color online) Time evolution of different observablesfor the iron dimer excited by the application of an externallaser field E x ( t ) along the x axis (bond axis, dashed green line,lower panel) in presence of spin-orbit coupling. The compari-son between (cid:104) ˆ S z (cid:105) (solid line) and (cid:104) ˆ L e z (cid:105) (dotted line) is shownin the upper panel. It is obtained by fixing the atomic po-sitions during the dynamics and by performing a Ehrenfestmolecular dynamics calculation. In this last case also the to-tal atomic orbital momentum is shown (dashed black line). Inthe lower panel we show instead (cid:104) ˆ L (cid:105) (solid line) and ∆ (cid:104) ˆ S (cid:105) (dotted line) obtained for fixed atoms and for Ehrenfest dy-namics calculations. The units for (cid:104) ˆ S z (cid:105) and (cid:104) ˆ J z (cid:105) are eV · fs(left-hand side scale), while that for the field is eV/˚A (right-hand side scale). by looking at a Fe ferromagnetic dimer. In this case we first obtain the electronic ground state at the optimizedgeometry and then we let the system evolve in time un-der the effect of the same laser field of Fig. (1), polarizedalong the bond axis of the dimer. In the simulation theatomic degrees of freedom are explicitly included by per-forming the rtTDDFT evolution together with Ehrenfestmolecular dynamics for the ions In Fig. (2) we com-pare the temporal evolution obtained without allowingthe atomic motion to that obtained from Ehrenfest dy-namics. In the upper panel we look at (cid:104) ˆ L e z (cid:105) and (cid:104) ˆ S z (cid:105) and we note that the atomic motion plays no role in theirdynamics (the curves obtained with Ehrenfest dynamicsoverlap perfectly with the ones obtained by keeping theatoms fixed). The time scale is, in fact, too short andthe atoms do not have enough time to move. In the caseof fixed atomic coordinates the quantity (cid:104) ˆ L e z (cid:105) + (cid:104) ˆ S z (cid:105) isnot conserved, as expected. However, also the quantity (cid:104) ˆ L e z (cid:105) + (cid:104) ˆ S z (cid:105) + L atom z is not conserved during the Ehrenfestdynamics. In particular, the atomic orbital momentumis characterized by huge oscillations persisting also afterthat the pulse amplitude sets to zero. In the lower panelof Fig. (2) it is confirmed that the spin lost during the dy-namics is not transfered to the orbital degrees of freedom.The same conclusions are valid for both atomic clustersand bulk systems, and increasing the size of the clusterdoes not help to reinstate angular momentum conserva-tion (see also Ref. 30).In the next section we will critically review the ap-proximations taken in our time-dependent simulationsand show that the explicit form of spin-orbit potentialused is the source of the non-conservation of the angu-lar momentum. Such oversimplified description of thespin-orbit coupling can be corrected and the conserva-tion re-instated. III. BREAKDOWN OF THE CONSERVATIONLAW AND SOLUTION
For the remaining of the paper we will focus on finitesize systems, however, these arguments can be general-ized also to the case of extended, periodic, systems.
A. case 1) ˆ H SO = 0 The first case we analyse corresponds to the situationwhere there is no spin-orbit interaction. As a conse-quence, our attention is focussed on the conservation ofthe total angular momentum only (we do not consider thespin, given that it commutes with the system’s Hamilto-nian, it is a constant of motion). We start by writing theHamiltonian ˆ H of a system of N e electrons in the fieldgenerated by N i identical ions of mass M (we make thischoice to keep the treatment easier). Here the ions aretreated also as quantum particles, there is no additionalexternal field. The Hamiltonian thus writes,ˆ H = N i (cid:88) a =1 ˆ P a M + ˆ V ii + ˆ H ie (ˆ r , ˆ p , ˆ R ) , (3)where (ˆ r , ˆ p ) defines the set of electronic positions andlinear momentum operators, ˆ R = ( ˆ R ,x , ˆ R ,y , . . . ˆ R N i ,z )is the set of 3 N i atomic coordinate operators, ˆ V ii is theion-ion interaction and ˆ H ie represents the full electronicHamiltonian. This latter includes the electron kineticenergy term, the electron-electron interaction and theelectron-ion interaction. By using d ˆ R a /dt = ˆ P a /M , itis easy to show that the time derivative of the orbitalmomentum associated to ion a is ddt ˆ L a = − ˆ R a × ∇ R a (cid:104) ˆ H ie + ˆ V ii (cid:105) , (4)by writing ˆ H ie = ˆ H + ˆ U ee , where ˆ U ee is the electron-electron interaction potential and ˆ H = (cid:80) N e i=1 ˆ p i / m e + (cid:80) N e i=1 (cid:80) N i a =1 v ie ( | ˆ r i − ˆ R a | ), with v ie ( | ˆ r i − ˆ R a | ) being theelectrostatic potential exerted by the a -th ion on the i-thelectron, it is easy to show that ddt ˆ L a = − ˆ R a × (cid:20) N e (cid:88) i =1 ∇ R a v ie ( | ˆ r i − ˆ R a | ) + ∇ R a ˆ V ii (cid:21) , = − N e (cid:88) i =1 ˆ r i × ∇ R a v ie ( | ˆ r i − ˆ R a | ) − ˆ R a × ∇ R a ˆ V ii , (5)by combining Eq. (5) with the time derivative of the ioniclinear momentum operator, d ˆ P a /dt = −∇ R a (cid:2) ˆ H ie + ˆ V ii (cid:3) ,we obtain the following important relation ddt ˆ L a = ddt (cid:26) ˆ L a − N e (cid:88) i =1 ˆ r i × ˆ P a (cid:27) == − m e N e (cid:88) i =1 ˆ p i × ˆ P a + N e (cid:88) i =1 ˆ r i × ∇ R a ˆ V ii − ˆ R a × ∇ R a ˆ V ii , (6)from which we can conclude that (cid:80) N i a =1 ˆ L a is a constantof motion for the system’s Hamiltonian ˆ H , ddt N i (cid:88) a =1 ˆ L a = 0 , (7)(further details are given in A). In addition, we canrewrite (cid:80) N i a =1 ˆ L a = ˆ L atom + ˆ L e with Eq. (1) being man-ifestly satisfied in the case of a spin unpolarized sys-tem with no externally applied electromagnetic field. Itmay be also useful to separate the electronic contribu-tion to the atomic angular momentum from the ionicone, by defining ˆ L e a = − (cid:80) N e i =1 ˆ r i × ˆ P a we can finallywrite ˆ L a = ˆ L a + ˆ L e a . a. Proof of the inequality ˆ L e a (cid:54) = ˆ L a :Let us define the angular momentum ˆ L a as that asso-ciated to the atom a in the absence of other nuclei, thisoperator commutes with the Hamiltonian of the isolatedatom: ˆ H a ie = (cid:80) N e i =1 ˆ p i / m e + (cid:80) N e i =1 v ie ( | ˆ r i − ˆ R a | ) + ˆ U ee .In general it is easy to show that for every multi-atomsystem the inequality [ (cid:80) N i a =1 ˆ L a , ˆ H ie ] (cid:54) = 0 is valid due tothe broken spherical symmetry of the system (namely tothe term (cid:80) N e i =1 (cid:80) N i b (cid:54) = a v ie ( | ˆ r i − ˆ R b | )).From the previous considerations we have also[ (cid:80) N i a =1 ˆ L a , ˆ H ie + ˆ V ii ] = 0, while for the purely electronic part we write − i [ ˆ L e a , ˆ H ie ] / (cid:126) = ˆ R a × ∇ R a ˆ H ie . In gen-eral it can be shown that − i [ˆ L a , ˆ H ie ] / (cid:126) (cid:54) = ˆ R a × ∇ R a ˆ H ie (all the details are given in B) and, as a consequence,ˆ L e a (cid:54) = ˆ L a . In conclusion the electronic orbital momen-tum operator corresponding to ˆ L e = (cid:80) N i a =1 ˆ L e a cannot beidentified with the sum of the operators ˆ L a for isolatedatoms. b. Isolated atom case In the case of a single isolatedatom we have [ˆ L a , ˆ H a ie ] = 0 due to spherical symmetry,in addition ∇ R a ˆ H a ie = 0, and, as a consequence ˆ L e a =ˆ L a + c ˆ I where we can freely set c = 0. c. Non conservation of ˆ L at + ˆ L elec in TDTB In TDTB models the electronic orbital momentum op-erator around each atom is usually set by definition toˆ L e a = ˆ L a . The time derivative of the expectation value ofthe electronic orbital momentum operator around atom a is (see B) ddt L a = − i (cid:126) (cid:68)(cid:104) ˆ L a , ˆ H ie + ˆ V ii (cid:105)(cid:69) = − i (cid:126) (cid:68)(cid:104) ˆ L a , ˆ H ie (cid:105)(cid:69) (cid:54) = R a × < ∇ R a ˆ H ie > + R a × ∇ R a V ii = − ddt L atom a , (8)that leads to ddt N i (cid:88) a =1 L a + ddt L atom (cid:54) = 0 . (9)As a consequence these models cannot conserve the to-tal orbital momentum even, if the atoms are allowed tomove. Note that This result is independent on the choiceof ˆ H ie while ˆ U ee is usually parametrized by means of amean field approximation. B. case 2) ˆ H SO (cid:54) = 0 Here we consider the case of finite spin orbit cou-pling, so that the spin degrees of freedom need to bere-introduced in the discussion. We will not generalizethe treatment to the full Dirac formalism, since it is un-necessary for our conclusions, and so we will simply addthe spin-orbit coupling term to the general Hamiltonianof Eq. (3). The spin-orbit coupling operator associatedto atom a may be written in the following general formˆ H SO a = − e (cid:126) m e c N e (cid:88) i =1 ˆ σ i · (cid:2) ˜ E a (ˆ r i , t ) × ˆ p i (cid:3) , (10)where ˜ E a (ˆ r , t ) is the effective screened electric field dueto atom a and experienced by the electrons. From now onwe will treat the electron-electron interaction in a mean-field way, so that the expression for the effective field canbe written as, ˜ E a (ˆ r , t ) = E ext ( t ) /N i + ˜ E mf a (ˆ r − ˆ R a ), where˜ E mf a is short range due to the screening of the conductionelectrons, while E ext is the applied external field assumedto be spatially homogeneous. By summing over all theatoms the spin-orbit Hamiltonian becomesˆ H SO = ˆ H SO − ext + N i (cid:88) a =1 ˆ˜ H SO a , where, (11)ˆ H SO − ext = − e (cid:126) m e c (cid:20) E ext ( t ) · N e (cid:88) i =1 ˆ p i × ˆ σ i (cid:21) , ˆ˜ H SO a = − e (cid:126) m e c N e (cid:88) i =1 ˆ σ i · (cid:104) ˜ E mf a (ˆ r i − ˆ R a ) × ˆ p i (cid:105) . (12)The first contribution due to the externally applied elec-tric field is not of great interest here since it will act onlyduring the application of the laser pulse, while we aremainly interested in the long term dynamics of the sys-tem, therefore, we will focus on the second term only.Without any loss of generality we can rewrite the mean-field screened electric field as e ˜ E mf a = −∇ R a ˜ v ie ( | r − R a | ),where ˜ v ie is the screened electron-ion potential. The spin-orbit coupling operator then becomesˆ H a SO = i m e c (cid:126) N e (cid:88) i =1 ˆ S i · (cid:26)(cid:2) ˆ P a , ˜ v ie ( | ˆ r i − ˆ R a | ) (cid:3) × ˆ p i (cid:27) = 12 m e c N e (cid:88) i =1 (cid:88) ljk (cid:15) ljk ˆ S li ˆ P ja (cid:26) − i (cid:126) (cid:2) ˆ p ki , ˜ v ie ( | r − R a | ) (cid:3)(cid:27) ++ i m e c (cid:126) N e (cid:88) i =1 (cid:88) ljk (cid:15) ljk ˆ S li (cid:20) ˆ P ja ˆ p ki , ˜ v ie ( | ˆ r i − ˆ R a | ) (cid:21) . (13)The second term is exactly zero since the electronic andatomic momentum commute. Thus, the spin-orbit cou-pling experienced by a single electron becomesˆ H a SO = − ∂ r ˜ v ie ( r )2 m e c | ˆ r − ˆ R a | ˆ S · (cid:8) ( ˆ R a − ˆ r ) × ˆ P a (cid:9) = − ∂ r ˜ v ie ( r )2 m e c | ˆ r − ˆ R a | ˆ S · ˆ L a (14)where ˆ L a is the orbital momentum operator introducedin Eq. (6). It is evident that (cid:80) N i a =1 L a is not a constant ofmotion of the new Hamiltonian ˆ H (cid:48) = ˆ H + ˆ H SO , with ˆ H given by Eq. (3). In particular from Eq. (14), by writingˆ H SO = (cid:80) a f a ( r ) ˆ L a · ˆ S , it is easy to show that[ ˆ L a , ˆ H SO ] = ( ˆ R a − ˆ r ) × N i (cid:88) b =1 (cid:2) ˆ P a , f b (ˆ r ) (cid:3) ˆ L b · ˆ S ++ N i (cid:88) b =1 f b (ˆ r ) (cid:2) ˆ L a , ˆ L b (cid:3) · ˆ S , (15)where the first term on the right hand side is exactly zero.Furthermore, given the fact that ˆ L a is a generator of thespatial rotation group and satisfies the commutation rela-tions [ ˆ L a,x , ˆ L b,y ] = i (cid:126) δ a,b (cid:80) z ε xyz ˆ L a,z , we easily obtain d ˆ L a /dt = f a ( r )ˆ S × ˆ L a . Since d ˆ S /dt = (cid:80) N i a =1 f a ( r ) ˆ L a × ˆ S we conclude that d ˆ J /dt = 0 with ˆ J = (cid:80) N i a =1 ˆ L a + ˆ S . a. Non conservation of the orbital momentum in noncollinear TDDFT The non conservation of the total an-gular momentum in the rtTDDFT simulations presentedin Figs. (1 and 2) and reported previously in reference 30is due to the fact that the system’s spin orbit coupling operator (introduced via the pseudopotential approxima-tion) is given by ˆ H SO ( r ) = (cid:80) N i a =1 ˆ V SOps ( r − R a ) and it is afunction of the isolated atom orbital momentum operatorˆ L . As a consequence d ˆ L a /dt (cid:54) = f a ( r )ˆ S × ˆ L a and ˆ J is nota constant of motion. Note that also (cid:80) N i a =1 ˆ L a + ˆ S is nota constant of motion, since (cid:80) N i a =1 ˆ L a does not commutewith the crystal field Hamiltonian ˆ H . IV. RESULTS
In this section we consider a set of different atomicclusters and we analyze their temporal evolution by usingthe correct orbital-momentum-conserving spin-orbit cou-pling [see Eq. 14]. The results are then compared with theapproximated, non-orbital-momentum-conserving evolu-tion. For the analysis we use a TB approximation, whoseHamiltonian writes as followsˆ H ( t ) = N i (cid:88) a =1 P a M + ˆ˜ H [ R ]ie ( t ) + V ii ( R ) + N i (cid:88) a =1 Z ∗ a e R a · E ( t ) , (16)ˆ˜ H [ R ]ie ( t ) = N i (cid:88) a =1 (cid:88) i (cid:15) ai ˆ c † ai ˆ c ai − (cid:88) a
Here we compare the temporal evolution of an atomicdimer under the effect of the two different spin orbit cou-pling operators introduced in the previous section. Thebasis set is built of the s and p orbitals of the two atoms.The model here for explanatory reasons does not haveto be physically realistic and we will focus on more re-alistic systems in future works. In addition the modeldepends on a set of free parameters that are set arbitrar-ily and varied in order to look at their influence on thedynamics. These are hopping terms t ss , t sp and t pp , thespin orbit strengths λ a , the on-site energies (cid:15) ai and theelectric dipole parameters. The atoms are assumed fromnow on to be identical. The two spin-orbit operatorsconsidered here areˆ H = N i (cid:88) a =1 λ a ˆ L a · ˆ S , (18)ˆ H SO = N i (cid:88) a =1 λ a ˆ L a · ˆ S . (19) < S z > ( e V f s )
1e 1 < S z > ( e V f s )
1e 1
1e 1
1e 1 t (fs) t ss =0.01 eV t pp =0.05 eV =0.5 eVt ss =0.03 eV t pp =0.07 eV =0.1 eVt ss =0.03 eV t pp =0.07 eV =0.5 eV FIG. 3: (Color online) Time evolution of different observablesfor the atomic dimer excited by the application of an exter-nal laser field E x ( t ) slightly tilted with respect to the bondaxis, with the effect of the atomic motion being neglected. Inthe upper plots we use a spin orbit strength λ = 0 . t ss = 0 .
01 eV, t pp = 0 .
05 eV, we con-sider the expectation values < . . . > obtained by evolving thesystem’s Hamiltonian with spin-orbit coupling from Eq. (19)and < . . . > obtained by using the Hamiltonian with approx-imated spin-orbit from Eq. (18). In the plot on the left-handside we compare the expectation value of the spin, ˆ S z , forthe two calculations, while in the plots on the right-hand sidewe compare instead the expectation values of ˆ L z , ˆ J z for thetwo calculations with < ˆ L z > . For the middle plots we per-form the same calculations with parameters λ = 0 . t ss = 0 .
03 eV, t pp = 0 .
07 eV, while in the lower plots we use λ = 0 . t ss = 0 .
03 eV, t pp = 0 .
07 eV.
Note that the first is of the type commonly used intight-binding calculations and similar to that employedin rtTDDFT, while the second represents its generaliza-tion, as discussed in the previous section. In Fig. (3)we evolve the electronic system under the action of theHamiltonian (17), by approximating the spin-orbit cou-pling operator with Eq. (18) or with the complete spin-orbit coupling of Eq. (19). The different sets of parame-ters used in the calculations are provided in the figure’scaption, and we assume the atoms to be fixed in theirinitial positions. We distinguish the evolution of the ob-servables obtained with spin-orbit coupling (19) from theones obtained with spin-orbit (18) by using respectively < . . . > and < . . . > for the expectation values.As expected < ˆ J z > = < ˆ L z > + < ˆ S z > is a con-stant of motion after that the external electric-field pulsehas vanished (see the solid blue line in the three plotson the right hand side of Fig. (3) after t = 10 f s ), while (cid:104) ˆ J z (cid:105) = (cid:104) ˆ L z (cid:105) + (cid:104) ˆ S z (cid:105) is not (see dotted blue lines).been set to zero (see solid blue line in the three plotson the right after t = 10 f s ), while < ˆ J z > = < ˆ L z > + < ˆ S z > is not (see dotted blue lines). Thus, the total orbital momentum is, in general, not conserved un-der the Hamiltonian including the approximated spin-orbit coupling interaction. The expectation value of thetotal orbital momentum (cid:104) ˆ L z (cid:105) is also not conserved forboth types of evolution (see red dotted and solid lines).The dynamics of (cid:104) ˆ L z (cid:105) is oscillatory with the frequencythat depends on the strength of the spin-orbit interac-tion, while the dynamics of < ˆ L z > appears to be not asstrongly dependent on the spin-orbit interaction.The plots on the left hand side of Fig. (3) show the tem-poral evolution of the two observables (cid:104) ˆ S z (cid:105) and (cid:104) ˆ S z (cid:105) .Their evolution is, not surprisingly, quite different, due tothe fact that the operator ˆ L is dependent on the strengthof the hopping terms, since [ ˆ H , ˆ L ] = 0 must always besatisfied. The spin decay appears more pronounced forstronger spin-orbit interaction, although, for the samevalue of λ it seems that the effect of the hopping termsshould not be neglected either. In Fig. (4) we study the t (fs) < S z > ( e V f s ) ×10 < S z > < S z > FIG. 4: (Color online) Time evolution of the observables (cid:104) ˆ S z (cid:105) (solid lines) and (cid:104) ˆ S z (cid:105) (dashed lines) for the different modelsof spin-orbit coupling in the case of a three-atoms system.The spin-orbit strength is fixed to a value λ = 0 . t ss =0 .
03 eV while t sp (cid:39) t pp have values (cid:39) .
12 eV (black curves), (cid:39) .
22 eV (red curves) and (cid:39) .
32 eV (blue curves). Theexternal electric field is applied along one of the bond axes. time evolution of a three-atom system initiated by anexternally applied laser pulse and driven by two differ-ent models of spin-orbit coupling. The figure shows theevolution of the spin expectation values (cid:104) ˆ S z (cid:105) and (cid:104) ˆ S z (cid:105) when the spin-orbit coupling strength is fixed and thehopping terms are set to the values given in the caption.This choice of the parameters has no particular mean-ing from a physical point of view, and our purpose is toshow how the choice may affect the orbital momentumoperator ˆ L and in turn the spin-orbit coupling, whileits strength λ is kept fixed. The evolution under thetwo spin-orbit operators is clearly different, in particu-lar we observe that the non approximated one [Eq. (19)]leads to a greater loss of spin magnetization comparedto the results obtained with the approximated evolution[Eq. (18)]. The effect in the case of real magnetic systemis still to be investigated, but already at this level it isclear that the two spin-orbit operators may lead to differ-ent rates of spin decay thanks to the different nature ofthe orbital operators ˆ L and ˆ L . The fact that ˆ L describesin addition to the localized orbital properties also the or-bital properties of the delocalized electronic states leadsto an additional contribution to the spin-orbit operatorwith a non trivial influence on the dynamical evolutionof the magnetic system. B. Coupling between electrons and molecularvibrations
In this section we briefly consider the effect of molecu-lar vibrations on the spin dynamics. In order to accountfor this effect our tight-binding Hamiltonian is modifiedto writeˆ˜ H [ R ]ie ( t ) = N i (cid:88) a =1 (cid:88) i (cid:15) ai ˆ c † ai ˆ c ai − (cid:88) a
03 ˚ A (greendotted line) and compared with the evolution in absenceof vibrations for the exact (red line) and the approxi-mated (black line) spin-orbit interaction. The effect doesnot seem particularly relevant at these time scales withonly a noticeable change in the amplitude of the oscilla-tions, while the evolution under the approximated spin-orbit interaction leads to huge modifications in the spindynamics, as we have already observed. t (fs) < S z > ( e V f s ) ×10 < S z > < S z > < S z > (1)< S z > (2) FIG. 5: (Color online) Time evolution of the observables (cid:104) ˆ S z (cid:105) (solid red lines) and (cid:104) ˆ S z (cid:105) (solid black lines) for the differentmodels of spin-orbit coupling and without considering molec-ular vibrations in the case of a five-atoms system. The othertwo lines represent (cid:104) ˆ S z (cid:105) obtained by using two different am-plitudes of atomic vibrations, namely 0.03 ˚A [green dottedline (2)] and 0.1 ˚A [blue dashed line (1)]. V. CONCLUSIONS AND FUTURE WORKS
In this work we have discussed how to properly accountfor the conservation of the total orbital momentum instate-of-the-art ab initio simulations of ultrafast demag-netization processes. We have identified the main prob-lem behind the observed breaking of the full rotationalinvariance to be associated with the standard implemen-tation of the spin-orbit coupling operator. We have thenexplained why an approximated spin-orbit coupling op-erator, built as a function of the isolated atom orbitalmomentum, ˆ L , may lead to difficulties in the estimationof the spin and orbital momentum temporal evolution.A proper analysis of the ultrafast demagnetization pro-cess should not be limited to the correct description ofthe electron-atom (phonon) channel of energy and or-bital momentum dissipation, but should be based alsoon a more accurate description of the spin-orbit couplinginteraction. In section IV we have shown how, evenfor simple models, an oversimplified description of thespin-orbit coupling may affect the temporal evolution ofthe spin magnetization. These results are preliminaryand based on simple model tight-binding Hamiltonians.At this point they only aim to show the importance ofthese new effects. In the future the model will be ex-tended to describe more realistic systems (bulk mate-rials and transition metal atomic clusters). This newdescription of the ultrafast spin dynamics has some in-teresting consequences and it may explain the apparentnon-conservation of the orbital momentum observed inexperiments . The measure of the electronic orbital mo-mentum is, in fact, performed by means of X-ray mag-netic circular dichroism (XMCD), which provides sepa-rate information on the dynamics of the spin and theelectronic orbital momentum. However, it is the expec-tation value of ˆ L , which is not a globally conserved quan-tity, to be directly accessible through these methods. Intransition metals, ˆ L a , is expected to significantly differfrom ˆ L a around each atom, suggesting that the missingorbital momentum observed in experiments is not neces-sarily transfered to the nuclei but could remain hiddenin the form of delocalized electronic orbital momentum. VI. ACKNOWLEDGMENTS
JS was supported by the Molecular Foundry, a DOEOffice of Science User Facility supported by the Office ofScience of the U.S. Department of Energy under ContractNo. DE-AC02-05CH11231. SS thanks the Irish ResearchCouncil (IRCLA/2019/127) for financial support.
Appendix A: Proof of Eq. (7)
We take Eq. (6) and we sum over all the ions a on theleft and the right hand side. ddt N i (cid:88) a =1 ˆ L a = − m e N e (cid:88) i =1 ˆ p i × N i (cid:88) a =1 ˆ P a + N e (cid:88) i =1 ˆ r i × N i (cid:88) a =1 ∇ R a ˆ V ii −− N i (cid:88) a =1 ˆ R a × ∇ R a ˆ V ii , (A1)The first term on the right hand side is exactly zero dueto the conservation of the linear momentum, (cid:80) N e i =1 ˆ p i = − (cid:80) N i a =1 ˆ P a . The second term is also zero. In Gaussianunits we write N e (cid:88) i =1 ˆ r i × N i (cid:88) a =1 ∇ R a ˆ V ii = N e (cid:88) i =1 ˆ r i × N i (cid:88) a =1 (cid:88) a (cid:48) (cid:54) = a ∇ R a ( Ze ) | R a − R a (cid:48) | = N e (cid:88) i =1 ˆ r i × N i (cid:88) a =1 (cid:88) a (cid:54) = a (cid:48) (cid:20) − ( Ze ) ˆ R a | R a − R a (cid:48) | + ( Ze ) ˆ R a (cid:48) | R a − R a (cid:48) | (cid:21) = 0 . (A2)Analogously for the third term we have N i (cid:88) a =1 ˆ R a × ∇ R a ˆ V ii = N i (cid:88) a =1 (cid:88) a (cid:48) (cid:54) = a ( Ze ) ˆ R a × ˆ R a (cid:48) | R a − R a (cid:48) | = ( Ze ) N i (cid:88) a =1 (cid:88) a (cid:48) (cid:54) = a (cid:20) ˆ R a × ˆ R a (cid:48) | R a − R a (cid:48) | − ˆ R a × ˆ R a (cid:48) | R a − R a (cid:48) | (cid:21) = 0 , (A3)that finally proves Eq. (7), namely, that the operator (cid:80) N i a =1 ˆ L a is a constant of motion for the system’s Hamil-tonian ˆ H . x (0 , y ( − d/ ,
0) ( d/ , R R δθ FIG. 6: (Color online) picture of the homonuclear dimer thatwe are using as a model system for the analysis of the or-bital momentum conservation. In its initial configuration thebonding axis is along directed along x , in its successive con-figuration the entire dimer is rotated by an infinitesimal angle δθ around the out-of-plane z axis. Appendix B: Proof of inequality (8)
Here we explicitly verify the inequality (8) for the elec-tronic orbital momentum. The simplest possible casewe may consider here is that of a homonuclear diatomicmolecule (see Fig. (6), the nature of the atoms is notrelevant), with wave functions expanded over a minimalbasis of s , p x , p y and p z orbitals with no spin polariza-tion. Here we do not want to describe a realistic situ-ation, just prove conceptually that the orbital momen-tum of the atomic plus electronic subsystem cannot beconserved in general if we employ the operators ˆ L a todescribe the electronic orbital momentum. According toFig. (6) we assume that the dimer is uniformly rotatingaround the z axis (out-of plane) by an infinitesimal angle δθ , the length of the bond is d , the Hamiltonian and or-bital momentum operators are expanded over the atomicbasis set {| s (cid:105) , | p − (cid:105) , | p (cid:105) , | p (cid:105) , | s (cid:105) , | p − (cid:105) , | p (cid:105) , | p (cid:105)} , in the new rotated configuration we canwriteˆ L z = (cid:32) L ,z L ,z (cid:33) = − − (B1)that corresponds to the total electronic orbital momen-tum operator, obtained as a direct sum of the orbitalmomentum operators corresponding to the two isolatedatomic sites, ˆ L = ˆ L ⊕ ˆ L . The hamiltonian ˆ H ie of thesystem, in turn, depends on 4 free parameters, the on-siteorbital energies (cid:15) s , (cid:15) p (we assume degeneracy of the p or-bitals for simplicity) and the hopping coefficients dependon a set of primitive integrals t ss , t sp , t σpp and t πpp thathere are left arbitrary but in general are a function of thedistance between the two atoms of the dimer molecule.These integrals t nlm a ; n (cid:48) l (cid:48) m (cid:48) a ij are diagonal in the angularmomentum projection m a along the bond axis and do notdepend on its sign , t nlm a ; n (cid:48) l (cid:48) m (cid:48) a ij = t nl | m a | ; n (cid:48) l (cid:48) | m a | ij δ m a ,m (cid:48) a .These can be written through the following general ex-pression t nlm a ,n (cid:48) l (cid:48) m a ij = (cid:90) V d r ψ ∗ nlm a ( r − R i ) h ij ( r ) ψ n (cid:48) l (cid:48) m a ( r − R j ) , (B2)where ψ nlm a is an atomic wave function and h ij ( r ) = − (cid:126) ∇ / m e + v ie ( r − R i ) + v ie ( r − R j ), the free hop-ping parameters are given by the following two-centersintegrals t ss = (cid:104) n = 2 , l = 0 , m a = 0 | ˆ h | n = 2 , l (cid:48) = 0 , m a = 0 (cid:105) ,t sp = (cid:104) n = 2 , l = 0 , m a = 0 | ˆ h | n = 2 , l (cid:48) = 1 , m a = 0 (cid:105) ,t σpp = (cid:104) n = 2 , l = 1 , m a = 0 | ˆ h | n = 2 , l (cid:48) = 1 , m a = 0 (cid:105) ,t πpp = (cid:104) n = 2 , l = 1 , m a = 1 | ˆ h | n = 2 , l (cid:48) = 1 , m a = 1 (cid:105) . In addition we also have the exact relation t nlm a ,n (cid:48) l (cid:48) m a ij =( − l + l (cid:48) t nlm a ,n (cid:48) l (cid:48) m a ji that leads to the final expression forthe dimer’s Hamiltonianˆ H ie = (cid:15) s t ss λ − tsp √ t sp λ z − λ + tsp √ (cid:15) p − λ − tsp √ tσpp + tπpp − λ z ∆ tpp λzλ − ∆ tpp √ − λ − ∆ tpp (cid:15) p − t sp λ z λzλ − ∆ tpp √ λ z ∆ t pp + t πpp − λzλ +∆ tpp √ (cid:15) p λ + tsp √ − λ − ∆ tpp − λzλ +∆ tpp √ tσpp + tπpp − λ z ∆ tpp t ss − λ + √ t sp − t sp λ z λ −√ t sp (cid:15) s λ + √ t sp tσpp + tπpp − λ z ∆ tpp λzλ + √ ∆ t pp − λ − ∆ t pp (cid:15) p t sp λ z λzλ + √ ∆ t pp λ z ∆ t pp + t πpp − λzλ −√ ∆ t pp (cid:15) p − λ −√ t sp − λ ∆ t pp − λzλ −√ ∆ t pp tσpp + tπpp − λ z ∆ tpp (cid:15) p where we have used ∆ t pp = t σpp − t πpp , the dimer bondlength is d= | R − R | with λ x,y,z = ( R − R ) x,y,z / | R − R | and λ ± = λ x ± iλ y .In the rotated configuration we have ( R − R ) x = d cos δθ , ( R − R ) y = d sin δθ and we can substitute inthe Hamiltonian λ ± = 1 − δθ ± iδθ , (B3) λ z = 0 , (B4)the new atomic coordinatets according to Fig. (6) is givenby R = d (cos δθ, sin δθ, R = − d (cos δθ, sin δθ, L = − R × ∇ R (cid:104) δ ˆ H ie (cid:105) = − R × ∇ R δθ · ∂ δθ (cid:104) δ ˆ H ie (cid:105) ˙ L = − R × ∇ R (cid:104) δ ˆ H ie (cid:105) = − R × ∇ R δθ · ∂ δθ (cid:104) δ ˆ H ie (cid:105) and the total orbital momentum time derivative becomes˙ L atom = ˙ L + ˙ L = − ∂ δθ (cid:104) ˆ H ie (cid:105) · ( R − R ) × ∇ R δθ , = − e z · ∂ δθ (cid:104) ˆ H ie (cid:105) . (B5)In order to simplify the calculation we proceed by assum-ing that we have a single electron in a state given by thesuperposition of | s (cid:105) around atom 1 and | p − (cid:105) aroundatom 2; | Ψ (cid:105) = [1 , , , , , , , / √ δθ∂ δθ (cid:104) ˆ H ie (cid:105) = − t sp √ δθ , ˙ L atom = 2 √ · ˆ e z t sp δθ . (B6)The time derivative of the expectation value of the elec-tronic orbital momentum is given by ( (cid:126) is set to 1) ddt (cid:104) ˆ L z (cid:105) = i (cid:104) [ ˆ H ie , ˆ L z ] (cid:105) , (B7)that explicitly computed by using the wave function | Ψ (cid:105) previously defined gives ddt (cid:104) ˆ L z (cid:105) = − t sp √ δθ . (B8)The time derivative of the total orbital momentum˙ L atom + ddt (cid:104) ˆ L z (cid:105) along the z axis is˙ L atom z + ddt (cid:104) ˆ L z (cid:105) = 3 √ t sp δθ , (B9)that proves Eq. (7) in the main text.0 ∗ Contact email address: [email protected] E. Beaurepaire, J.-C. Merle, A. Daunois and J.-Y. Bigot,Phys. Rev. Lett. , 4250 (1996). J.-Y. Bigot, L. Guidoni, E. Beaurepaire and P.N. Saeta,Phys. Rev. Lett. , 077401 (2004). J. Hohlfeld, E. Matthias, R. Knorren and K.H. Benne-mann, Phys. Rev. Lett. , 4861 (1997). J. G¨udde, U. Conrad, V. J¨ahnke, J. Hohlfeld andE. Matthias, Phys. Rev. B , R6608 (1999). A. Scholl, L. Baumgarten, R. Jacquemin and W. Eber-hardt, Phys. Rev. Lett. , 5146 (1997). G. Ju, A.V. Nurmikko, R.F.C. Farrow, R.F. Marks,M.J. Carey and B.A. Gurney, Phys. Rev. Lett. , 3705(1999). B. Koopmans, M. van Kampen, J.T. Kohlhepp andW.J.M. de Jonge, Phys. Rev. Lett. , 844 (2000). T. Ogasawara, K. Ohgushi, Y. Tomioka, K.S. Takahashi,H. Okamoto, M. Kawasaki and Y. Tokura, Phys. Rev. Lett. , 087202 (2005). A.V. Kimel, A. Kirilyuk, A. Tsvetkov, R.V. Pisarev andT. Rasing, Nature , 850 (2004). M. Vomir, L.H.F. Andrade, L. Guidoni, E. Beaurepaireand J.-Y. Bigot, Phys. Rev. Lett. , 237601 (2005). G. Ju, J. Hohlfeld, B. Bergman, R.J.M. van de Veerdonk,O.N. Mryasov, J.Y. Kim, X. Wu, D. Weller and B. Koop-mans, Phys. Rev. Lett. , 197403 (2004). J.U. Thiele, M. Buess and C.H. Back, Appl. Phys. Lett. , 2857 (2004). A.V. Kimel, A. Kirilyuk, P.A. Usachev, R.V. Pisarev,A.M. Balbashov and T. Rasing, Nature , 655 (2005). C.D. Stanciu, F. Hansteen, A.V. Kimel, A. Kirilyuk,A. Tsukamoto, A. Itoh and T. Rasing, Phys. Rev. Lett. , 047601 (2007). J. Berezovsky, M.H. Mikkelsen, N.G. Stoltz, L.A. Coldrenand D.D. Awschalom, Science , 349 (2008). A.V. Kimel, A. Kirilyuk and T. Rasing, Laser Photon. Rev. , 275 (2007). A. Kirilyuk, A.V. Kimel and T. Rasing, Rev. Mod. Phys. , 2731 (2010). M. F¨ahnle and C. Illg, J. Phys.: Condens. Matter ,493201 (2011). , Y. Hinschberger and P.-A. Hervieux, Phys. Lett. A ,813 (2012). H. Vonesch and J.Y. Bigot, Phys. Rev. B , 180407(R)(2012). E. Carpene, E. Mancini, C. Dallera, M. Brenna, E. Puppinand S. De Silvestri, Phys. Rev. B , 174422 (2008). M. Krauss, T. Roth, S. Alebrand, D. Steil, M. Cinchetti,M. Aeschlimann and H.C. Schneider, Phys. Rev. B ,180407(R) (2009). K. Carva, M. Battiato and P.M. Oppeneer, Phys. Rev.Lett. , 207201 (2011). M. Battiato, K. Carva and P.M. Oppeneer, Phys. Rev. B , 024404 (2012). M. Hennecke, I. Radu, R. Abrudan, T. Kachel, K. Holl-dack, R. Mitzner, A. Tsukamoto and S. Eisebitt, Phys.Rev. Lett. , 157202 (2019). D. Steiauf, C. Illg and M. F¨ahnle, J. of Phys.: Conf. Series , 042024 (2010). M. F¨ahnle, T. Tsatsoulis, C. Illg, M. Haag, B.Y. M¨ullerand L. Zhang, J. Supercond. Nov. Magn. , 1381 (2017). M. F¨ahnle, M. Haag and C. Illg, J. Magn. Magn. Mater. , 45 (2013). K. Krieger, J.K. Dewhurst, P. Elliott, S. Sharma andE.K.U. Gross, J. Chem. Theory Comput. , 4870 (2015). M. Stamenova, J. Simoni and S. Sanvito, Phys. Rev. B ,014423 (2016). W. T¨ows and G.M. Pastor, Phys. Rev. Lett. , 217204(2015). D. Chaudhuri, G. Lefkidis and W. H¨ubner, Phys. Rev. B , 184413 (2017). E.K.U. Gross and W. Kohn, Adv. Quant. Chem. , 255(1990). E. Runge and E.K.U. Gross, Phys. Rev. Lett. , 997(1984). ˆ L KSelectrons is the orbital momentum of the Kohn-Sham sys-tem and it should always be distinguished from the truemany body electronic orbital momentum, L electrons . How-ever, in order to simplify our analysis, we will assumethroughout the paper that two quantities are equal to agood degree of approximation. This assumption, in fact,will not play any role in the coming considerations. M.A.L. Marques, A. Castro, G.F. Bertsch and A. Rubio,Comp. Phys. Comm. , 60 (2003). D.R. Hamann, M. Schl¨uter and C. Chiang, Phys. Rev.Lett. , 1494 (1979). L. Kleinman, Phys. Rev. B , 2630 (1980). M.J.T. Oliveira and F. Nogueira, Comp. Phys. Comm. , 524 (2008). X. Andrade, A. Castro, D. Zueco, J.L. Alonso,P. Echenique, F. Falceto and A. Rubio, J. Chem. TheoryComput. , 728 (2009). P. Koskinen, V. M¨akinen, J. Com. Mat. Sci. , 237 (2009)., 237 (2009).