## Abstract

Despite more than 20 years of development, the underlying physics of the laser-induced demagnetization process is still debated. We present a fast, real-time time-dependent density functional theory (rt-TDDFT) algorithm together with the phenomenological atomic Landau-Lifshitz-Gilbert model to investigate this problem. Our Hamiltonian considers noncollinear magnetic moment, spin-orbit coupling (SOC), electron-electron, electron-phonon, and electron-light interactions. The algorithm for time evolution achieves hundreds of times of speedup enabling calculation of large systems. Our simulations yield a demagnetization rate similar to experiments. We found that (i) the angular momentum flow from light to the system is not essential and the spin Zeeman effect is negligible. (ii) The phonon can play a role but is not essential. (iii) The initial spin disorder and the self-consistent update of the electron-electron interaction play dominant roles and enhance the demagnetization to the experimentally observed rate. The spin disorder connects the electronic structure theory with the phenomenological three-temperature model.

## INTRODUCTION

The interaction of femtosecond laser pulses with magnetically ordered materials has attracted a large number of studies (*1*–*7*) ever since the discovery of subpicosecond demagnetization about 20 years ago (*8*). Not only does it have potential applications in information technology, but also the basic physics in such an ultrafast demagnetization process is intriguing and has attracted many fundamental studies in this field (*9*–*16*). Besides simple demagnetization, many new related phenomena have been discovered along the way, including optically generated spin current (*17*, *18*), laser-induced spin reorientation (*19*–*23*), and spin flip in complex structures (*24*, *25*). Despite this progress, the basic mechanism and the fundamental process underneath the demagnetization phenomenon especially from the electronic structure point of view are still under debate (*26*–*31*). There are a few long-standing problems regarding the demagnetization mechanism (*13*), including the following: What is the channel of angular momentum flow that leads to the change of magnetic moment? Is spin-orbit coupling (SOC) strong enough to cause such ultrafast demagnetization? What is the original reservoir of the angular momentum for demagnetization: photon, phonon, or electron orbital? How do we properly describe the demagnetization process? Should it be described as a thermodynamic process represented by three temperatures for spin, orbital, and phonon degree of freedom, or a coherent precession process due to the direct light-spin interaction? Is the electron-electron interaction important, and can it be described by independent single-particle picture or perturbation theory? To understand these problems, a variety of models have been proposed, including electron-magnon coupling (*32*, *33*), Elliott-Yafet phonon scattering (*34*, *35*), spin-flip Coulomb scattering (*36*, *37*), direct photon-induced spin flips (*12*), relativistic spin-light interaction (*11*), and superdiffusive spin transport (*16*, *38*–*40*). Because of the complexity of the problem, many of these analytical models are based on single-particle pictures without electron-electron interaction and correlation. In the electron-magnon coupling picture (*32*, *33*), the electron excitation will induce magnon, which then reduces the magnetic moment. In the Elliott-Yafet phonon scattering (*34*, *35*) and spin-flip Coulomb scattering explanation (*36*), the noninteracting single-particle scattering by phonon and charge causes the change of spin via SOC. In the direct photon-induced spin flip picture or the relativistic spin-light interaction picture (*11*, *12*), a Stoner-like excitation directly changes the electron spin because the spin is not a good quantum number in the presence of SOC. Last, in the superdiffusive spin transport picture (*16*, *38*, *39*), the total spin of the system is not changed. Instead, the electron is excited into the itinerary orbital with the majority spin and is then quickly diffused into surrounding regions, hence depleting the magnetic moment at the original excitation spot. Despite all these pictures, a consensus is yet to be reached (*13*, *41*). By studying the laser fluence dependence and scattering mechanisms in the demagnetization process, Roth *et al.* (*41*) argued that demagnetization is a thermal process, instead of a coherent process between light and spin, or superdiffusive transport. Carva and Oppeneer *et al*. (*42*–*44*) have shown that the perturbation treatment of the well-known Elliott-Yafet scattering picture (without considering the self-consistent electrostatic potential here) (*34*, *35*) is not strong enough to explain the ultrafast demagnetization. Zhang *et al*. (*12*, *45*) and Oppeneer *et al*. (*46*) debated whether spin flip during optical transitions via SOC is sufficient to explain the observed large demagnetization. There are also different conclusions for the central question of spin angular momentum reservoirs (e.g., whether it is photon, electron orbital, or phonon) and how the total angular momentum is conserved (*32*, *34*, *41*, *45*, *47*–*52*). While it is possible that some of these pictures describe some aspects of the demagnetization process, a consensus from an ab initio electronic structure point of view is yet to emerge.

In recent years, there is some consensus in using phenomenological micromagnetic simulations based on the Landau-Lifshitz-Gilbert (LLG) model to describe the demagnetization process (*53*, *54*). Such simulations are based on some prior assumptions. For example, usually three energy reservoirs—phonon, electron, and spin—are assumed each with their own temperature. The heat exchange between these reservoirs, the critical damping constant, and the thermal noise amplitude that controls the demagnetization in the LLG model are parameterized to fit with the experiment. The basic picture of such simulation is to assume that the temperature is above the Curie temperature due to the heat injection from the laser field, which leads to demagnetization. However, the question regarding the flow of angular momentum is not addressed because the LLG model uses thermal noise as the source of angular momentum. Thus, although there is some phenomenological understanding of the problem, there is still a lack of understanding from the electronic structure ab initio point of view, and a connection between the quantum mechanical description and the phenomenological description is still lacking. Given this situation, it is thus desirable to simulate the process directly using ab initio numerical methods (*39*, *43*, *45*, *49*–*52*, *55*), which describe the system without any assumptions. Zhang *et al.* (*45*, *49*, *50*) performed first-principles and time-dependent simulations for Ni and NiO under laser illumination. Gross and co-workers (*51*, *55*) developed a real-time time-dependent density functional theory (rt-TDDFT) code (“real-time” here means using explicit time propagation of electron wave functions compared to the perturbation TDDFT treatment) and simulated a small bulk system (one or four atoms) for a short time (20 fs) without taking into account the phonon degree of freedom. Such direct rt-TDDFT simulation can be advantageous compared to many analytical/theoretical studies because the electron-electron interaction, which is very important for magnetic systems but difficult to be described with analytical models (*13*), is included throughout the simulation. Unfortunately, the current direct ab initio simulations suffer a major drawback: The simulated demagnetization rate is one order of magnitude smaller than the experimentally observed one. This casts some doubt on whether the real physics is simulated. For example, only about 1% demagnetization is observed with a laser fluence of 11.5 mJ/cm^{2} in a time-dependent Liouville equation calculation (*56*). This demagnetization rate can be increased to about 10% when exchange-correlation functional allows the inclusion of spin polarization change (*49*), while experimental 50% (or higher) demagnetization is observed with similar laser fluence (*8*, *10*). In an rt-TDDFT simulation, 43% demagnetization is reached only when the laser fluence reaches 934.8 mJ/cm^{2} (*55*), which can totally change the electronic structure of the system. At this point, it is not clear what has caused such a major discrepancy. Because of the extreme time-consuming nature of the ab initio real-time simulations, there are many limitations in the current studies. For example, very often, a rather small system (e.g., one-atom primary cell) is used, and the phonon degree of freedom is also ignored.

In this work, we use a new algorithm to carry out rt-TDDFT simulations. This algorithm markedly increases the time step of the time evolution integration, e.g., from a typical 10^{−5} fs to 0.05 fs, hence achieving hundreds of times of speedup and enables us to calculate bigger systems. Our rt-TDDFT simulations have considered the effect of noncollinear magnetic moment, light-spin, light-orbital, spin-orbit, electron-electron, and electron-phonon (by nuclear movements) interactions as well as the effect of finite-temperature spin disorder in larger systems (up to 64 atoms). These advantages, together with the LLG model, enable us to yield the experimental demagnetization rate and reveal roles of various effects. An isolated system calculation also reveals the flow and conservation of different types of angular momentum. Through our study on ferromagnetic Ni systems, we conclude that (i) the angular momentum flow from photon to the system is not essential. If any, such flow goes to the electron orbital, not directly to the spin. The direct photon-spin interaction through the spin Zeeman effect is negligibly small. (ii) The phonon can play a role, but is not notable for demagnetization, as including the phonon effect alone cannot yield the experimental rate of demagnetization. (iii) There are major roles for the initial spin disorder at room temperature and the self-consistent update of the electron-electron interaction. Initial disorder significantly enhances the demagnetization and helps the system to reach the experimental demagnetization rate. The spin disorder also connects the electronic structure theory with the phenomenological three-temperature model. The importance of the self-consistency of the electron-electron interaction casts some doubts about many current analytical explanations with non–self-consistent treatments. (iv) All the angular momentum needed for the spin demagnetization comes from the electron orbital via SOC. Overall, we have the following picture of the demagnetization process: The absorption of light induces excitation of the electron orbital, which further excites the spin degree of freedom via SOC. Such excitation is much like the thermal noise in the LLG, hence assigning an effective high temperature to the spin degree of freedom. In such a picture, the demagnetization is a random and collective process. As a result, the initial spin disorder/randomness (in a lesser degree, the phonon can also provide some randomness) is essential to describe this process. Without it, e.g., in a single-atom unit-cell description, the simulation can only describe a deterministic precession driven by an external force. As a result, a single-atom unit-cell simulation might miss an essential point of the phenomena.

## RESULTS

In our rt-TDDFT method, the time-dependent Kohn-Sham equation has the following formalism* _{ik}* contains spin-up and spin-down components;

*i*is the band index and

*k*is the

*k*-point index; and

*V*

_{ion},

*V*

_{H},

*V*

_{XC}, and

*V*

_{SOC}are the ionic potential, Hartree potential, exchange-correlation potential, and SOC term, respectively.

*V*

_{SOC}is based on a nonlocal pseudopotential term representing the relativistic effects of the core level. More specifically,

*V*

_{SOC}= ∑

_{R,J}

*D*∣

_{J}*J*> <

_{R}*J*∣, where

_{R}*J*is a reference state constructed from the Ni

_{R}*d*atomic orbital and spin with total (orbital + spin) angular momentum

*J*and

*D*is a coefficient for this

_{J}*J*projector. σ is the Pauli matrix.

*A*(

*r*,

*t*) and

*B*(

*r*,

*t*) represent the magnetic vector potential and magnetic field in the laser, respectively. We apply the dipole approximation and ignore the spatial dependence of

*A*(

*r*,

*t*) and

*B*(

*r*,

*t*), because the laser wavelength is much longer than the size of our supercell (

*55*). The 1/2σ ⋅

*B*(

*r*,

*t*) term represents the spin Zeeman interaction between magnetic field and spin magnetic moment. This is the only direct interaction term between light and the spin degree of freedom. The time dependence of the laser is represented by a sin oscillation multiplied by a Gaussian envelop (see the Supplementary Materials). Because of the spin-up and spin-down components of the wave function, the charge density ρ(

*r*) is represented by a 2 × 2 matrix:

*o*(

*i*,

*k*) is the occupation number of the state

*i*at the

*k*-point

*k*[i.e., ψ

*(*

_{ik}*r*)]. Note that the Hamiltonian in Eq. 1 is also a matrix, and the exchange-correlation potential is given by

_{1}(

*r*) and ρ

_{2}(

*r*) are the “spin-up” and “spin-down” magnitudes of charge matrix ρ

_{αβ}(

*r*,

*t*) at its principle direction

*m*(

*r*,

*t*) at point

*r*. Local spin density approximation is used to calculate the derivative of

*E*

_{XC}with respect to ρ

_{1}(

*r*) and ρ

_{2}(

*r*).

To solve the time-evolving Eq. 1, we expand the wave function during the time interval of *t*_{1} to *t*_{2} in terms of the adiabatic eigenstates ϕ* _{lk}*(

*r*,

*t*

_{1}) at

*t*

_{1}

*(*

_{lk}*r*,

*t*

_{1}) is the eigenstate of the Hamiltonian of Eq. 1 at

*t*

_{1}and is given by our planewave pseudopotential DFT calculation code (PEtot) (

*57*). Then, Eq. 1 can be transformed to

*t*

_{1}to

*t*

_{2}. Specifically, we calculate the charge matrix at

*t*

_{2}from ψ

*(*

_{ik}*r*,

*t*

_{2}) and calculate

*H*(

^{k}*t*

_{2}). Then, we approximate

*H*(

^{k}*t*) by linear interpolation

*H*(

^{k}*t*) to integrate

*t*

_{1}to

*t*

_{2}by Eq. 4 to obtain ψ

*(*

_{ik}*r*,

*t*

_{2}), and recalculate the charge density matrix at

*t*

_{2}. This process is repeated until

*H*(

^{k}*t*

_{2}) is converged. Equation 3 is then updated using ϕ

*(*

_{lk}*r*,

*t*

_{2}) as the basis set to integrate the wave function from

*t*

_{2}to

*t*

_{3}. The Hellmann-Feynman forces are calculated in Ehrenfest dynamics to move the ions with classical Verlet algorithm. The above algorithm allows us to use a relatively large Δ

*t*=

*t*

_{2}−

*t*

_{1}(∼0.05 fs), while the integration of Eq. 4 from

*t*

_{1}to

*t*

_{2}is carried out using a much smaller time step (10

^{−5}fs). This enables hundreds of times of speedup compared to the conventional rt-TDDFT procedure and enables us to calculate larger systems, which have not been studied before. For example, the test on an eight-atom Ni bulk for 100-fs simulation lasts less than 1 day and the same calculation on a unit-cell bulk only lasts 1 hour.

Using the above rt-TDDFT method, we first studied the laser-induced demagnetization of a bulk Ni crystal. A Gaussian-enveloped linearly polarized laser light with a fluence of 22.5 mJ/cm^{2} is used. The wavelength of light is 550 nm, and the pulse duration is 60 fs (time between the start and the end of the pulse; see the Supplementary Materials). We used a 2 × 2 × 2 supercell of eight Ni atoms with initial thermal atomic movements. The 4 × 4 × 4 *k*-points are used without symmetry. Norm-conserving pseudopotential with an energy cutoff of 653 eV is used. A room temperature is used to set up the initial ionic velocity and atomic positions through a short ground-state molecular dynamics simulation. Figure 1A shows the laser pulse, Fig. 1B shows the total energy of the Ni system, while Fig. 1C shows the average magnetic moment per atom. One can see that, after the laser pulse, the *z*-direction magnetic moment has dropped from 0.7 μB per atom to about 0.64 μB per atom, representing an 8.3% demagnetization. Compared with previous real-time simulations (*55*, *56*), this 8.3% demagnetization is larger for the same light fluence, probably due to the phonon effect, as will be discussed later. Nevertheless, the 8.3% demagnetization is still much smaller than the experimental demagnetization of >50% for the similar laser fluence. We next use rt-TDDFT to analyze what contributes to the demagnetization, what is the channel of angular momentum transfer, and why the demagnetization is too small.

We first plot the density of state (DOS) for unoccupied holes in the valence bands and occupied electrons in the conduction bands during the demagnetization process in Fig. 1D (the total DOS is shown in the Supplementary Materials). The occupied electron DOS is defined as

We next study the effects of the phonon. To do this, we have turned off the ionic movement (fix the atoms in their ideal crystal position). The results are shown in Fig. 2A. One can see that removing the ionic movement will reduce the demagnetization by about 40%, which might explain why our demagnetization is larger than previous calculations. This demonstrates that the phonon does play a role. However, the phonon effect is not large enough to explain the experiment, thus is not dominant. Our conclusion is also in accordance with the analytical spin-flip Eliashberg function calculations of Carva *et al.* (*42*, *44*), where phonon-mediated spin-flip scattering is not strong enough to be accountable for ultrafast demagnetization.

We further study the role of light polarization. Thus far, a linear polarization laser has been used, with the polarization direction perpendicular to the Ni magnetization direction. We then replaced it with circular polarization and kept the other settings the same (e.g., ions are not frozen, the light propagation direction is along the magnetization, and light polarization is perpendicular to the magnetization). As shown in Fig. 2A, the left circular light reduces the demagnetization (the result of the right circular light is similar; see the Supplementary Materials). This difference might be due to (i) the light-spin interaction that directly transfers the photon angular momentum to spin or (ii) the light-orbital interaction that changes the electron excitations. To show the possible effect of direct angular momentum flow from photon to spin, we turned off SOC while keeping the σ ⋅ *B* term in Eq. 1. We found that, after SOC is turned off, there is no demagnetization for both linear and circular lights. Similarly, we found that, if we keep SOC but remove the σ ⋅ *B*(*r*, *t*) term in Eq. 1, the amount of demagnetization is almost unchanged. These show that the direct light-spin interaction (as described by the σ ⋅ *B* Zeeman term) plays a negligible role, which is in accordance with previous analytical works (*45*). It is the light-orbital interaction that transfers energy and angular momentum to the electron orbital, which then affects the spin via SOC. The negligible role of the spin Zeeman term can also be understood in terms that the spin motion cannot follow the rapidly oscillating magnetic field of light.

Because many of the analytical explanations center around the angular momentum transfer channels (*32*, *34*, *41*, *45*, *48*–*52*), it is interesting to calculate the spin, electron orbital, and ionic angular momentums separately and study the transfer between them. Unfortunately, in an infinite periodic crystal, both the electron orbital and ionic angular momentum are not well defined. To overcome this, we have studied a two-atom Ni dimer. The center of mass *R*_{0} is used to evaluate the angular momentum, with the electron orbital angular momentum defined as *L*_{e} = ∑* _{i}*〈ψ

*∣(*

_{i}*r*−

*R*

_{0}) ×

*i*∇ ∣ψ

*〉. The ion angular momentum is defined as*

_{i}*L*

_{ion}= ∑

*(*

_{j}*R*−

_{j}*R*

_{0}) ×

*M*, where

_{j}V_{j}*R*,

_{j}*M*, and

_{j}*V*are the position, mass, and velocity of the

_{j}*j*

^{th}ion, respectively. We have performed the rt-TDDFT simulations with an ultrashort laser of 0.4-fs duration. Small time step (10

^{−2}fs) and very high convergence criteria (10

^{−15}eV) in each time step have been used to evaluate the ionic angular momentum precisely, because the Ni ion mass is 10

^{5}times of the electron mass. The results for the spin, electron orbital, and ion angular momentums are shown in Fig. 2B. We see that, because of the interaction with light (electron excitation), the orbital angular momentum varies markedly within the duration of light (before 0.4 fs). Meanwhile, the spin angular momentum only changes slightly (not due to the direct spin-light Zeeman interaction, but due to SOC), and the ionic momentum also changes slightly. After the light is turned off (after 0.4 fs), the total angular momentum is conserved, with a big exchange of the angular momentum between electron orbital and ionic degree of freedom; meanwhile, some of the angular momentum has been transferred to the spin. This clearly shows that electron-ion can exchange angular momentum in a very fast rate (so is the exchange between electron orbital and spin). We also note that the final total angular momentum is not zero, indicating that the whole system absorbs a small angular momentum from the laser. Because our Hamiltonian does not trace the light angular momentum, such absorbed angular momentum is inferred from the conservation law.

From the above results, we can conclude that (i) the direct light-spin Zeeman interaction is negligible. (ii) The strong interaction comes from the light-orbital. Different light polarization induces different light absorption. (iii) Spin demagnetization needs the channel of SOC, and SOC is strong enough to cause an angular momentum exchange rate of about 3 × 10^{−3}μB/fs (judged from Fig. 2B), which is in the same order as in the experimental demagnetization (*1*, *31*, *41*). (iv) The angular momentum exchange between the orbital and the ion can be very fast due to the strong electron-nuclear Coulomb interaction. As a matter of fact, in a periodic crystal, one can consider the orbital-ion system as an infinite angular momentum reservoir in terms of angular momentum conservation. Thus, the relevant question should not be where the angular momentum comes from. Instead, it should be how the orbital-ion system transfers its angular momentum to the spin via SOC. The overall angular momentum flow picture is depicted in the schematic diagram of Fig. 2C.

Note that the effect of superdiffusive spin transport (*16*, *38*, *39*) is not taken into account in the above calculations. This is because the laser field that we used has no spatial dependence (because the laser wavelength is much longer than the supercell size). There are some arguments (*16*, *38*, *39*) and counterarguments (*27*) using superdiffusive spin transport from the laser-excited region to the unexcited region to explain the spin loss of the system. This is not the view adopted by our current study.

Having clarified the issue of angular momentum flow, we now turn our attention to the central question: Why is the theoretically simulated demagnetization rate much smaller than the experimental results? While we get about 8.3% demagnetization (Fig. 1C), the corresponding experimental value is about 50% (or higher) for the same laser fluence. There are two types of possible demagnetization pictures. In the first type, the magnetic moments of all atoms change (or rotate) in unison, so a unit cell is able to describe the whole process. In the second type, the magnetic moments of different atoms change differently in a thermally disordered fashion, much like what happens in a paramagnetic system in the LLG model. In this case, it is essential to describe the system with a supercell and with some initial randomness. A supercell without initial randomness and without taking random phonon movement into account will only behave like a single-atom cell.

In the above calculations, we have explored the effect of randomness caused by the phonon. However, as we can see from Fig. 2A, although such phonon randomness helps to enhance the demagnetization, its magnitude is still not large enough. Another more direct randomness is the spin random orientation at room temperature before light illumination. Although at room temperature the bulk Ni is in its ferromagnetic phase, its individual atomic spin orientations are not exactly pointing to the *z* direction. Significant orientational disorder exists. To estimate the amplitude of this disorder, we have carried out atomic LLG simulations with the nearest neighboring ferromagnetic interaction coefficient *J _{ij}* of 17.23 meV (

*53*) for the bulk Ni. Monte Carlo (MC) simulation with Metropolis sampling is used to compute the statistical average, and Langevin dynamics with the Heun time integration algorithm (

*53*) is used to perform the real-time simulation (see the Supplementary Materials). As shown in Fig. 3A, we see that the average angle (Δθ

_{LLG}) between two nearest neighboring spins increases with temperature (

*T*). At zero temperature, spins are completely aligned and Δθ

_{LLG}= 0. At a higher temperature, Δθ

_{LLG}notably increases, indicating remarkable spin disorder. This disorder results in the decline of the magnetization

*m*

_{LLG}, as the temperature increases until the magnetization goes to zero above the Curie temperature (

*T*

_{c}, at about 627 K). Note that the classical MC simulation well estimates

*T*

_{c}but overestimates the decrease of magnetization at a low temperature. This is a well-known issue of classical spin models (

*58*). As shown in Fig. 3A, the experimental magnetization curve (

*m*

_{exp}) is flatter at a low temperature (

*59*). Using this

*m*

_{exp}and the mapping relation between

*m*

_{LLG}and Δθ

_{LLG}, we can compute the experimental Δθ

_{exp}(see Materials and Methods). It can be seen that Δθ

_{exp}also increases sharply with temperature. Its value is about 22.3° at 300 K, indicating a strong spin randomness. In addition, from Fig. 3B, one can see that the spin-spin correlation function decays until about four lattice constants. Thus, to adequately describe the random behavior of such systems, one should use a large supercell, with a size of perhaps four to eight (or larger) lattice constants. Unfortunately, our current rt-TDDFT scheme still cannot handle such large systems. To simulate even larger systems, we further improve our algorithm with a fixed basis set expansion.

In this scheme, in Eq. 4, after integration from *t*_{1} to *t*_{2}, instead of replacing the basis set from ϕ* _{lk}*(

*t*

_{1}) to ϕ

*(*

_{lk}*t*

_{2}) to carry out the integration from

*t*

_{2}to

*t*

_{3}, we will continue to use the original basis ϕ

*(*

_{lk}*t*

_{1}) throughout the simulation. If the ϕ

*(*

_{lk}*t*

_{1}) basis set is large enough and the simulation time is short, then this should be exact (see the Supplementary Materials). The use of a fixed basis set allows us to update the wave function, charge density matrix, and potential in real space based on the coefficient

*t*=

*t*

_{2}−

*t*

_{1}to be sufficiently small (10

^{−4}fs), so charge density self-consistent iteration is no longer necessary. To further speed up the calculation with many large systems, we have used a strong but short light pulse (4-fs duration) with total fluence similar to the experiment (see the Supplementary Materials). Although this pulse is shorter than the typical experimental pulse, we like to use this to investigate the effects of initial magnetic disorder in a comparative study.

Using the above fixed basis rt-TDDFT (to be called FB-rt-TDDFT thereafter), we have studied the demagnetization under different initial spin disorder. To prepare the initial disordered electronic structure of the system, we start with the ferromagnetic ground state and then introduce some random disorder [Δ*m*(*r*)] on its atomic spin vector while keeping the scalar charge density unchanged (see the Supplementary Materials). The resulting Hamiltonian will yield a set of eigenstate {ϕ* _{lk}*(

*r*,0)} and its corresponding charge density matrix

*m*(

*r*) is tuned, and the disorder in ρ

_{αβ}(

*r*,0) results in a desired average Δθ between nearest atomic spins, as described in Fig. 3A. This {ϕ

*(*

_{lk}*r*, 0), ρ

_{αβ}(

*r*, 0)} constitutes our initial disordered electronic state. This state is not the ground state, and its energy at Δθ = 40° is about 30 meV/atom higher than the ground state at Δθ = 0°. This is a relatively low-energy excitation in the same order as the room temperature

*k*(hence can be considered as the spin thermal excitation). To focus on the laser-induced demagnetization of this disordered system, not to be obscured by its own spin fluctuation, we have added an effective correction term in the Hamiltonian (see the Supplementary Materials), so the time-dependent evolution of {ϕ

_{B}T*(*

_{lk}*r*,

*t*), ρ

_{αβ}(

*r*,

*t*)} without the external light incident is a steady state (no change of charge density matrix and spin). As a result, the demagnetization in the following discussion is only induced by light illumination. Note that the change of the total spin with time is evaluated on the basis of the electronic wave functions of the whole system, which follows the time-dependent TDDFT equation (we do not have a specific equation, like the torque in the LLG model, for the evolution of individual atoms’ spin).

We first test the effect of initial disorder on the 4 × 4 × 4 system with 64 atoms. We have introduced different initial random orientation of spins for a series of independent simulations. For each case, we applied different random disorder [Δ*m*(*r*)] to prepare the initial electronic states (see the Supplementary Materials), and the average canting angle (Δθ) between two nearest neighboring spins is taken as the fingerprint. The demagnetization results for 16 different cases are shown in Fig. 3C. It can be seen that there is a marked effect of the initial spin disorder. When there is no initial disorder (Δθ = 0), the amplitude of the demagnetization is about 10%, similar to the above rt-TDDFT simulation. However, when the initial disorder reaches 22.3° (corresponding to the room temperature shown in Fig. 3A), the system has about 20 to 25% demagnetization. This rate can be even bigger (>40%) in a larger disorder. To illustrate the effect of the electron-electron interaction during the demagnetization process (a many-electron correlation effect), we have also artificially fixed the *V*_{H} and *V*_{XC} terms in Eq. 1 during the simulation (fix el-el). In this case, the demagnetization will only be contributed by the single-electron light-orbital, light-spin, and SOC interactions, while the many-electron correlation effect is turned off. This is the situation in many conventional analytical perturbation treatments (e.g., Stoner excitation, Elliott-Yafet phonon scattering, and light-spin direct interaction). As shown in Fig. 3C, we found that the demagnetization is significantly reduced if the self-consistent update of the electron-electron interaction is turned off. Without spin disorder and such self-consistence, the demagnetization is less than 2%. Even with spin disorder, the demagnetization is still significantly lower than the case with electron-electron self-consistence and is less than 10%. This reveals the important role of the update of the electron-electron interaction and the collective behavior in the demagnetization process.

The above 4 × 4 × 4 system might not be large enough to fully describe the long-range spin-spin correlation as shown in Fig. 3B. To further study this, we have simulated an 8 × 1 × 1 system following the same procedure. The results with 120 cases of different degree of initial spin disorder are shown in Fig. 3D. As one can see, the effect of spin disorder is indeed slightly larger than the 4 × 4 × 4 case, and at Δθ = 22.3°, the system reaches about 40% demagnetization, similar to that of the experimental observations. We have also confirmed that the demagnetization is significantly reduced when the contribution of the electron-electron interaction is turned off.

The effect of the initial spin disorder can also be observed in the magnitude of spin-orbit torque (SOT; the expectation value of the operator σ × *L*). Its *z* component value along the magnetization direction is zero at *t* = 0 for the 4 × 4 × 4 system without the initial spin disorder. This can be understood in the point that σ and *L* are parallel in the ferromagnetic ground state. This value becomes 1.66 × 10^{−5} μB/fs (per atom) for the case with Δθ = 26.1° and 3.01 × 10^{−5} μB/fs for the case with Δθ = 37.7°. It indicates that the spin randomness leads to a canting angle between σ and *L*. Such canting will be enlarged by the laser excitation. After the laser excitation, the SOT value increases to 1.68 × 10^{−3}, 1.06 × 10^{−2}, and 1.91 × 10^{−2} μB/fs for the non-disordered ground-state case and the two disordered cases, respectively. Their difference shows the role of spin disorder.

Because the degree of spin disorder corresponds to different equilibrium temperature (Fig. 3A), our results predict that there is strong relation between the demagnetization process and the initial material temperature. This is observed in some recent experiments where the demagnetization rate was shown to increase with the initial equilibrium temperature (*41*, *60*). The phenomenological microscopic three-temperature model (M3TM) was used to explain their experiments. From the entropy point of view, the nature of spin temperature is disorder. Our simulations thus provide the ab initio insights into this problem.

To understand why the initial magnetic disorder enhances the demagnetization from the electronic structure point of view, we can write down the rate of the total magnetic moment change as* _{ik}*(0) basis set, and

*A*(

_{ij}*k*) = 〈ϕ

*(0) ∣ [σ,*

_{ik}*H*] ∣ ϕ

*(0)〉 is the spin change matrix. Note that if we ignore the small σ ⋅*

_{jk}*B*(

*r*,

*t*) term in Eq. 1 as justified from our previous tests, then the [σ,

*H*] term will be independent of time. This leads to a time-independent

*A*(

_{ij}*k*) that is only related to the initial wave functions.

We have shown the *A _{ij}* matrix and

*D*(

_{ij}*t*) matrix for the 4 × 4 × 4 system at

*t*= 2.5 fs (during the laser excitation) in Fig. 4. Two cases are compared, one without the initial spin disorder and the other with the spin disorder (Δθ = 37.7°). Only the

*z*component (ferromagnetic direction) of

*A*is shown. We can see that only a limited number of off-diagonal

_{ij}*A*elements are nonzero in the case without disorder, while the

_{ij}*A*elements of the disordered case are widely spread and are overall much bigger. The

_{ij}*D*elements of the disordered case are relatively more uniform than the case without disorder, but they are overall similar, representing similar light adsorption. The main difference comes from

_{ij}*A*and the resulting multiplication ∑

_{ij}*⋅*

_{ij}A_{ij}*D*. One can view

_{ij}*A*as a torque exerted on the spin of a given atom due to the electron excitation between states

_{ij}*i*and

*j*. The density matrix

*D*(

_{ij}*t*) represents such excitation due to the laser light exposure. The initial spin disorder significantly increases the torque in a disordered fashion, resulting in the acceleration of demagnetization.

*A _{ij}* can also be analyzed by considering the symmetry of [σ,

*H*]. Here, we will focus on the spin operator on one atom. Then, the SOC term in Eq. 1 can be approximated as ασ ⋅

*L*, where

*L*is the angular momentum operator of this atom and α is a coefficient. As a result, [σ

*,*

_{z}*H*] = α(σ ×

*L*)

*= α(σ*

_{z}*− σ*

_{x}L_{y}*). If we focus on the*

_{y}L_{x}*p*,

_{x}*p*, and

_{y}*p*symmetry components (not necessarily the atomic orbital character) within the wave function ϕ

_{z}*(0), then the orbital coupling due to*

_{ik}*L*and

_{y}*L*are <

_{x}*p*∣

_{z}*L*∣

_{y}*p*> and <

_{x}*p*∣

_{z}*L*∣

_{x}*p*>, respectively. If the occupied ϕ

_{y}*(0) has a high symmetry (for high-symmetry*

_{ik}*k*-point) with a spin pointing strictly to the

*z*direction, then it will be likely to have some

*p*component, but might not have

_{z}*p*and

_{x}*p*components due to symmetry. Then,

_{y}*A*of such two high-symmetry states will be zero. Furthermore, in such an ordered case, the supercell wave functions that belong to different primary cell

_{ij}*k*-points will also not couple. All these indicate that there could be many zeros of

*A*for an ordered system. For a disordered system, such forbidden coupling will become allowed, resulting in a more spreading

_{ij}*A*. The disorder also makes

_{ij}*D*more uniform (Fig. 4D), instead of more pattern as in Fig. 4C, which is probably a result of the

_{ij}*k*-point selection rule in the excitation of the ordered system. All these lead to a larger overlap of

*A*and

_{ij}*D*in the disordered system. The randomness in both

_{ij}*A*and

_{ij}*D*supports a picture of a random torque to the spin degree of freedom by the orbital excitation, which mimics the finite-temperature thermal noise in the LLG model. The thermal noise serves as an effective spin temperature as in the M3TM to drive the demagnetization. Thus, this connects our electronic structure description to the phenomenological micromagnetic descriptions (M3TM and LLG) for the demagnetization process.

_{ij}## CONCLUSIONS

In conclusion, we used fast rt-TDDFT methods, together with the phenomenological LLG model, to study the ultrafast spin dynamics induced by the laser. Our rt-TDDFT simulations yield a demagnetization rate similar to the experiment, a goal not achieved in previous TDDFT works. We found that (i) the angular momentum flow from light to the system is not important. Such flow goes to the electron orbital, while the direct light-spin interaction is negligible. (ii) The phonon can play a role, but it is not the most critical one. (iii) The initial spin disorder and the electron-electron interaction play dominant roles. Such initial disorder significantly enhances the demagnetization and helps the system to reach the degree of demagnetization observed experimentally. The spin disorder also connects the electronic structure theory with the phenomenological micromagnetic descriptions (M3TM and LLG). The importance of the self-consistent update of the electron-electron interaction casts some doubts about previous analytical models with non–self-consistent treatments. (iv) All the angular momentum needed for the spin demagnetization comes from the orbital-ion system via SOC. Overall, we have the following picture: The laser induces the excitation of the electron orbital. Such an electron orbital excited state exerts a random torque on the spin via SOC. This random torque represented by the *A* and *D* matrix in Fig. 4 serves as an effective temperature for the spin degree of freedom and results in the collective demagnetization process.

## MATERIALS AND METHODS

All the rt-TDDFT and FB-rt-TDDFT calculations were performed by our noncollinear magnetic version of PEtot code (*57*). The electron wave functions were expanded by a plane-wave basis set with an energy cutoff of 653 eV. Each wave function has two components (spin-up and spin-down), and the Hamiltonian operators were extended to 2 × 2 matrices (see the Supplementary Materials). We tested four types of Ni systems in these TDDFT calculations. For the 2 × 2 × 2 bulk in Figs. 1 and 2A , we used 4 × 4 × 4 *k*-points without symmetry for the summation over the Brillouin zone. For the two-atom dimer in Fig. 2B, Γ point was used. For the 4 × 4 × 4 bulk in Figs. 3C and 4, 2 × 2 × 2 *k*-points were used. For the 8 × 1 × 1 bulk in Fig. 3D, 1 × 8 × 8 *k*-points were used. The exchange and correlation potentials for both rt-TDDFT and FB-rt-TDDFT were treated in the framework of local spin density approximation and were described by Eq. 2 for the noncollinear spin. Norm-conserving pseudopotential was used in all the calculations. The Hartree and exchange-correlation potentials will be self-consistently updated based on the new charge density matrix in every time step. This is different from conventional analytical perturbation treatments, where the charge density and the corresponding electron potentials are not self-consistently updated.

To quantitatively estimate the average Δθ_{exp}(*T*) under the experimental *m*_{exp}(*T*) in Fig. 3A, we can reasonably assume that the correlation between the nearest spins are the same as in the LLG simulation; however, only the amplitudes of the angle deviations are overestimated by LLG for a given temperature. Thus, to find Δθ_{exp}(*T*) for an experimental reduced magnetic moment *m*_{exp}(*T*) at temperature *T*, one can find the temperature *T*^{′} where *m*_{LLG}(*T*^{′}) = *m*_{exp}(*T*), and then use the corresponding Δθ_{LLG}(*T*^{′}) as the experimental Δθ_{exp}(*T*). Following this procedure, we can compute the whole Δθ_{exp}(*T*) curve, as shown in Fig. 3A.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/6/eaau8000/DC1

Section S1. Notes on rt-TDDFT simulations

Section S2. Notes on FB-rt-TDDFT simulations

Section S3. Atomic LLG model

Fig. S1. Evolution of total DOS.

Fig. S2. Evolution of magnetic moment per atom for Ni bulk under light with linear polarization, left circular polarization, and right circular polarization.

Fig. S3. Change of electron spin angular momentum *S*_{e}, electron orbital angular momentum *L*_{e}, ion orbital angular momentum *L*_{ion}, and their total angular momentum (all in the *z* direction) for Ni dimer under the circular polarization light.

Fig. S4. Magnetic vector potential (*A*) of Gaussian-enveloped linearly polarized laser with 550 nm of wavelength, 22.5 mJ/cm^{2} of fluence, 4 fs of duration, and peak at 2 fs.

Fig. S5. Evolution of normalized spin (*m*) in an 8 × 1 × 1 long system excited by the 4-fs laser of fig. S4.

Fig. S6. Tests of the fixed basis set size on the 8 × 1 × 1 system (80 valence electrons in total).

Fig. S7. The statistical average angle (θ) from the spin direction to the *z* direction in the equilibrant states, as a function of temperature.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**This research used the resources of the National Energy Research Scientific Computing Center (NERSC) and Oak Ridge Leadership Computing Facility (OLCF) that are supported by the Office of Science of the U.S. Department of Energy, with the computational time allocated by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) project NTI009. We thank Dr. Zhi Wang for the very useful help in the coding of the rt-TDDFT program.

**Funding:**This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, under contract no. DE-AC02-05-CH11231 within the Non-Equilibrium Magnetic Materials program (MSMAG).

**Author contributions:**Z.C. and L.-W.W. wrote the plane-wave DFT codes, rt-TDDFT codes and FB-rt-TDDFT codes, performed the simulations, analyzed the data, and wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).