Ultrafast demagnetization in bulk vs thin films: an ab-initio study

We report on {\it ab-initio} simulations of the quantum dynamics of electronic charge and spin when subjected to intense laser pulses. By performing separate calculations for a Ni thin film and bulk Ni, we conclude that surface effects have a dramatic influence on amplifying the laser induced demagnetization. We show that the reason for this amplification is due to increased spin-currents on the surface of the thin film. This enhancement is a direct consequence of the broken symmetry originating from the surface formation. We find that the underlying physics of demagnetization, during the early femtoseconds, for both bulk and thin film is dominated by spin-flips induced by spin-orbit coupling. After the first $\sim 40$ fs this changes in that the dominant cause of demagnetization is the flow of spin-currents, which leads to stronger demagnetization in the film compared to that of the bulk.

Most of the theoretical work to-date deals with the effect of laser pulses on the bulk of the magnetic material. Research that deals with interface/surface are model calculations [26,27] focusing only on a single aspect namely the diffusion of electrons across the interface. However, realistic devices contain surfaces, interfaces and bulk regions and a detailed understanding of the behaviour of these regions separately, under the influence of a laser pulse, is crucial [41,42] to unraveling the complete physics of demagnetization. Experimentally it is a very challenging task to distinguish between contributions coming from various regions of the sample, but theoretically, by considering separate calculations for thin film and bulk of the same material, these effects can be disentangled. In the present work we use Ni thin film and bulk calculations to show that the underlying physics responsible for ultra-fast demagnetization in early femtoseconds is the same in both; spin-orbit induced spin-flips. Interestingly, we find that symmetry breaking originating from the formation of a surface (or interface) greatly enhances the demagnetization process. By explicitly treating spin and charge-currents in our simulations we find that the main reason behind this enhanced demagnetization in thin films is the presence of large spin-currents generated in these broken symmetry systems.
In order to do this theoretical work we use the ab-initio method of time-dependent density functional theory (TDDFT) [43][44][45][46], which is, in principle, an exact theory for studying light-matter interaction and makes no assumptions about the form of the electron dynamics. The study of dynamics of spins using TDDFT requires an extension where the magnetization density [47] is treated as an unconstrained vector field. Such an extension was recently performed [28], facilitating the present work. The calculations presented here are purely electronic in nature and include contributions like (a) spinorbit induced spin-filps (b) restricted set of magnon excitations (by means of a super-cell calculations) and (c) spin-diffusion (or spin currents) to the process of demagnetization. Processes such as Elliott-Yafet scattering and electron-phonon (or lattice) induced spin-relaxation are ignored. A comparison with experiments then allows one to quantify the contribution of purely electronic processes to the physics of demagnetization and time scales at which other processes like Elliott-Yafet scattering become significant. In the present work we find that for times-scales below 120fs purely electronic processes dominate the physics of demagnetization.

Theoretical Aspects
Within TDDFT [43] a Kohn-Sham(KS) Hamiltonian is used for time evolving the electronic wave-function. The scalar-relativistic KS Hamiltonian used in the present work reads:Ĥ where c is the speed of light,σ are the Pauli spin operators and A ext (t) is the vector potential representing an applied laser field. In the present work it is assumed that the wavelength of the applied laser pulse is much longer than the length of a unit cell. This assumption allows for the so-called dipole approximation and it implies that the spatial dependence of the vector potential A ext can be disregarded. The final term of Eq. (1) is the spin-orbit coupling term which is included in the most general form and thus automatically includes any derived forms, e.g. Rashba-Dresselhaus coupling.
The KS effective potential v s (r, t) = v ext (r, t) + v H (r, t) + v xc (r, t) is decomposed into the external potential v ext , the classical electrostatic Hartree potential v H and the exchange-correlation (XC) potential v xc . Similarly, the KS magnetic field is written as B s (r, t) = B ext (t) + B xc (r, t) where B ext (t) is the magnetic field of the applied laser pulse plus possibly an additional magnetic field and B xc (r, t) is the XC magnetic field.
In the present work we use the adiabatic local density approximation [48,49] (ALSDA) for the XC functional. In order to analyze the contribution of each term in this Hamiltonian (1) to the dynamics of the magnetization M(t) = m(r, t)d 3 r = σn(r, t) d 3 r, wheren is the density operator and m(r, t) is the magnetization density, we start by calculating the dynamics of the magnetization density using Ehrenfest's theorem: SubstitutingĤ s from Eq. 1 into Eq. 2 leads to Evaluating the commutators results in can be evaluated by integrating Eq. 4; the integral over −∇ · ← → J (r, t) vanishes due to Gauss's law and the integral over 1 4c 2 [∇n(r, t) × ∇v S (r, t)] vanishes upon integrating by parts. For ALSDA B xc (r, t) × m(r, t) = 0 and in the absence of any external magnetic field (i.e. B ext =0), the dynamics of the magnetization is given by: Here each spin-current density describes the flow of the respective spin-component. The z-component of Eq. 5 reads:

Results
In order to distinguish the behaviour of spins on the surface (or interface) from those in the bulk of a sample subjected to an intense laser pulse, we study two cases in the present work; laser induced spin-dynamics in (a) bulk Ni and (b) a free standing film of Ni. Bulk Ni is easy to simulate in an electronic structure code which uses periodic boundary conditions (Elk code is used for all simulations [50]), but the computation of a thin film requires special care; in the present work we have used a 5 atomic layer thick film with a 5 layer thick vacuum. The z-axis points in plane of the film and y-axis out of the plane. The ground state of this thin film is ferromagnetic with the magnetization pointing along the z-axis (in-plane of the film) (with a layer resolved moment of 0.71, 0.68 and 0.66 µ B starting from the layer adjacent to the vacuum). The pump laser pulse is then applied perpendicular onto this surface. At t = 0 we begin the simulation from the ground-state for both systems (bulk and film) and the results for the relative moment Comparison of the averaged layer-resolved moment to the experimental data of Ref. [11]. The parameters of the pulse were also taken from this reference.
M z (t)/M z (t = 0) as a function of time are shown in Fig. 1. In both cases, under the influence of the laser pulse of fluence 8.05 mJ/cm 2 , the system is first optically excited (i.e. electrons are excited to energetically higher lying states), followed by a global loss in the magnetic moment. The two cases differ in the amount of demagnetization; the bulk shows a small loss of moment (∼ 8% during the simulation) while for the thin film this loss is much larger (∼ 20%). Furthermore the calculations done in the film geometry can be compared to realistic experiments, where the laser pulse is applied perpendicular to the surface of the sample and such a comparison is made in Fig. 1c. The results show that we only reach an agreement with the XMCD data of Ref. [11] for the first ∼120fs beyond which the theoretical work saturates and experimental data continues to demagnetize. The reason behind this disagreement between theory and experiment is the fact that in the present work we have included electronic only contribution (like spinflips, spin-current and magnons) to the process of demagnetization and electron-lattice, electron-phonon, other Elliott-Yafet like mechanisms and radiation losses are ignored. These results thus point to two important findings (a) how significant and at what time scales is the contribution of the electronic processes to the total (experimental) demagnetization and (b) the time scales at which dissipative processes like Elliott-Yafet mechanism start to be significant. We find that, in the present case, around 120fs such processes start to contribute and become dominant for longer times. At this point, it is natural to ask if the underlying mechanism for ultrafast demagnetization in bulk differs from a thin film. To answer this we note that in a purely electronic simulation there are two distinct spin excitation processes that can lead to a loss in the moment: processes that lead to local moment loss like magnons or non-collective canting of spins between atoms and processes that lead to global moment loss like spin-flips (or Stoner like excitations). We find that in both cases, for the first ∼120fs, the observed loss in magnetization along the z-direction is not accompanied by an increase in magnetization in the x or y directions, indicating that the long-range non-collinearity plays very little role on these time scales and the dynamics of Fig. 1 is dominated by spin-flip processes for both bulk and thin films. Furthermore, the term in the Hamiltonian (Eq. 1) responsible for this magnetization dynamics is the spin-orbit coupling and setting this term to zero leads to no global demagnetization.
Despite the similarity in the underlying physics, the amount of demagnetization is markedly different in the two cases (film and bulk). The question then arises: what leads to this strong difference? From Eq. 5 it is clear that the spin-current tensor, ↔ J , is the sole quantity responsible for demagnetization. Hence to understand the difference between the thin film and the bulk demagnetization, in Fig. 2, we plot the x-spin-component of the spin-current, j x , which can be defined following Eq. 4 as j x (r) = σ x ⊗ 1 2 {n(r),p + 1 c A ext } = σ x ⊗ĵ . This quantity can be interpreted as the vector field showing the flow of spins orientated along the x-direction. As an aid to visualizing these currents, we can create streamlines by following the direction of the vector at each point with the strength of this flow shown in colour. Fig. 2 thus shows how the spins pointing in the x-direction are flowing (circulating) around the y-axis. From these results it is clear that for both the bulk and the film the spincurrent loops clockwise about the y-axis, as is required by Eq. 6 to cause a change in the moment. However, the magnitude of the spin-current in the case of the thin film is much larger than for the bulk and this enhanced spin-current then leads to larger demagnetization in the film. The reason behind this we find to be the broken symmetry due to surface/interface formation, which allows for large surface currents due to the presence of localized electronic wave-function (a natural consequence of the lower symmetry). The flow of this surface current towards the center of the film then causes large spin-currents every where in the film resulting in a significant demagnetization.
To gather further insight into different behaviour of thin films and bulk it is instructive to compare the layer resolved magnetization of the film to that of the bulk. These results, for an ultra-short pulse with FWHM=17fs, are plotted in Fig. 3 for the bulk, the top layer of the film and the very central layer of the film. From this data it is clear that during the first ∼30 fs the demagnetization in the bulk and the central layer of the film are similar while the top layer of the film demagnetizes faster. Beyond 40fs the layers of the film continue to demagnetize while the bulk saturates. These results raise two interesting questions; why the demagnetization is almost the same for the bulk and the film in the first 30fs and why do the layers of the film continue to demagnetize whereas the bulk saturates beyond 40fs? One of the major reasons for these effects is that below ∼30fs the spin-currents have almost the same magnitude for the bulk and the film, subsequently the broken symmetry allows for stronger currents in the film, while for the case of the bulk these currents stay small. In case of the film these spin-currents flow from one layer to another and continue to flow beyond the first 40fs, leading to further demagnetization of the layers of the film. This explains the observed physics of demagnetization in the two cases. It would be interesting to know the thickness at which the very central layer of the film start to behave like the bulk of Ni. But TDDFT, a state-of-the-art method to deal with light-matter interaction, comes at a very high computational cost. As a result, treating a film greater than a few atomic-layer thickness is not feasible with the resources available to us.
Finally, it is crucial to mention that in Fig. 3 we have used the predictive power of TDDFT to study the influence of ultra-short laser pulses on the spin dynamics. Presently there are no experiments reported which utilize such pulses, however the parameters are within the range of available technology. The effect of this ultra-short pulse is faster and much enhanced demagnetization; for a film 51% of the moment is lost while for bulk this value is 17%. The corresponding values for the long pulse (see Fig. 1) are 20% for the film and 8% for the bulk.

Discussion
By performing ab-initio simulations for laser pulse excited ferromagnetic thin films, we have demonstrated that thin films show enhanced spin-orbit mediated demagnetization compared to a bulk. We have further demonstrated that this enhancement in the demagnetization near the surface of a film is due to the broken symmetry which increases the rotating spin currents in the system. These calculations show the importance of treating the spin-orbit coupling as well as charge-and spin-currents at the same footing.
From our calculations one can conclude that for a typical laser pump pulse (as currently used in experiments), demagnetization (caused by purely electronic processes) is strong within the first few atomic layers of a material but will then decrease to the smaller bulk value as we get deeper into the sample and it becomes more bulk-like. For longer duration pulses, such as that used in Fig. 1, for the first ∼120fs purely electronic processes dominate the physics of demagnetization. Beyond first 120fs, these electronic mechanisms will exist in addition to dissipative mechanisms such as Elliott-Yafet, and it will require more ingenious experiments to disentangle and distinguish them. This early time purely electronic regime is of great importance for optimal control and device production at ultrashort timescales as it allows for the possibility of coherent control of the electrons.

Computational Details
Considering future coherent control of spins by light, in the present work we explicitly concentrate on the electronic degrees of freedom during the early femtoseconds. Dissipative processes that induce decoherence such as radiation and phonons are not included (nuclei are kept fixed during the simulation). At this point it is important to mention that collective magnetic excitations (e.g. magnons) are also purely electronic in nature and have been included in the present work through usage of a large super-cell. However, we found that magnons only have a minor contribution to the spin-dynamics in the time scales studied in the present work [28].
State-of-the art full potential linearized augmented plane wave (LAPW) method implemented within the Elk code [50] is used in the present work. The core electrons (with Eigenvalues below 95eV) are treated using the radial Dirac equation while higher lying electrons are treated using the scalar relativistic Hamiltonian in the presence of the spin-orbit coupling. To obtain the 2-component Pauli spinor states, the Hamiltonian containing only the scalar potential is diagonalized in the LAPW basis: this is the first-variational step. The scalar states thus obtained are then used as a basis to set up a second-variational Hamiltonian with spinor degrees of freedom [51]. This is more efficient than simply using spinor LAPW functions, but care must be taken to ensure that there is a sufficient number of first-variational eigenstates for convergence of the second-variational problem. In the present work 300 empty states per k-point are used.
For Ni a face centered cubic crystal structure with lattice spacing of 3.52Å is used. A k-point grid of 8 × 8 × 8 is used for the bulk and 1 × 8 × 8 for the film calculations. A full geometry optimization for the Ni-film was performed, however, we found that for intense laser pulses, such as used in present work, geometry optimization does not change results significantly. The maximum augmented plane-wave cutoff for the orbitals was 3.5 au −1 corresponding to an energy cutoff of approximately 160 eV, and for the density and potential a value of 12 au −1 was used. An angular momentum cutoff of 8 for orbitals and 7 for densities and potential was used. A time step of 0.002 fs was used for time propagation [52]. The initial state for all TDDFT calculations was the DFT ground-state, which is at a temperature of 0K. We use the adiabatic local density approximation [48,49] for the XC functional. The film calculations required 900000 CPU hours, running on the HYDRA supercomputer at RZG Garching for 1 month.

Acknowledgements
TM and SS would like to thank QUTIF-SPP for funding. PE and SS would like to acknowledge funding from DFG through SFB762 project.