Demonstration of Optimal Control of Laser Induced Spin-Orbit Mediated Ultrafast Demagnetization

Laser induced ultrafast demagnetization is the process whereby the magnetic moment of a ferromagnetic material is seen to drop significantly on a timescale of $10-100$s of femtoseconds due to the application of a strong laser pulse. If this phenomenon can be harnessed for future technology, it offers the possibility for devices operating at speeds several orders of magnitude faster than at present. A key component to successful transfer of such a process to technology is the controllability of the process, i.e. that it can be tuned in order to overcome the practical and physical limitations imposed on the system. In this paper, we demonstrate that the spin-orbit mediated form of ultrafast demagnetization recently investigated [arXiv:1406.6607] by ab-initio time-dependent density functional theory (TDDFT) can be controlled. To do so we use quantum optimal control theory (OCT) to couple our TDDDT simulations to the optimization machinery of OCT. We show that a laser pulse can be found which maximizes the loss of moment within a given time interval while subject to several practical and physical constraints. Furthermore we also include a constraint on the fluence of the laser pulses and find the optimal pulse that combines significant demagnetization with a desire for less powerful pulses. These calculations demonstrate optimal control is possible for spin-orbit mediated ultrafast demagnetization and lay the foundation for future optimizations/simulations which can incorporate even more constraints.


INTRODUCTION
Control of quantum dynamics using tailored laser pulses is a long standing goal of modern physics [1][2][3][4][5][6][7] as it opens up a whole new world of possibilities for future technologies. Faster, smaller, and more efficient devices could be constructed if we could master control over the charge and spin dynamics of electrons on the nanoscale [8]. However precisely at these very short length and time scales, quantum effects are strong, which makes it difficult to exert this control. With the advent [9] of laser pulse shapers that can tailor the laser field to a given shape, there was now a tool that could be used for control of quantum dynamics. The challenge is finding the shape of the laser pulse that produces the desired dynamics.
Optimal control theory (OCT) is a method developed [10,11] in both Mathematics and Engineering to solve the problem of finding a particular control variable that gives a desired outcome. In our case, we will search for the electric field (t) of a laser pulse to control the properties of our system. In general OCT works by creating a target functional of the control field calculated from simulation of the system. Then any constraints on the system are incorporated using penalty functionals, before extremizing the total functional to find the optimal field. OCT can be extended to the realm of quantum mechanics by constructing the target functional using observables given by the time-dependent schrödinger equation (TDSE).
For more than a handful of electrons, propagating the TDSE is a computationally intractable problem due to the coulomb interaction between electrons and an alternative approach must be used. Time-dependent density functional theory (TDDFT) is one such approach, which works by mapping the problem to a non-interacting system [12], referred to as the Kohn-Sham (KS) system. This system is defined such that propagating electrons in this system will reproduce the same time dependent density (the probability to find an electron at any given point) as propagating in the exact system using the TDSE. As the KS system is non-interacting, the problem is now computationally tractable. Although, in principle, this mapping is exact, in practice an approximation must be used. The most commonly used approximation, adiabatic local density approximation (ALDA), has been shown to successfully predict absorption spectra of a large range of atoms, molecules, and solids [13][14][15]. Thus TDDFT is an outstanding candidate to couple to OCT [16] in order to predict laser pulses for control of quantum dynamics, and has been used successfully for control of charge transfer [17], HHG [20], strong-field ionization [18,19], bond-breaking [21], among others.
Laser-induced ultrafast demagnetization was first observed in the mid 1990s, whereby a strong femtosecond laser pulse caused a significant loss of the magnetic moment of a thin film of Ni in a time less than 1ps [22]. Since then, this phenomena has been the subject of much experimental [23][24][25][26][27][28][29][30][31][32][33] and theoretical [34][35][36][37] endeavor and several mechanisms have been proposed to explain the demagnetization. In Ref. 37, ab-initio TDDFT simulations were performed to investigate the demagnetization and found that when spin-orbit interaction was included in the system Hamiltonian, a loss of moment was ob-arXiv:1506.02332v1 [cond-mat.mtrl-sci] 8 Jun 2015 served for very short (5fs), very intense (1×10 15 W/cm2) laser pulses. It is this system we wish to control by varying the intensity and frequency of the laser pulse, subject to several practical constraints, in order to maximize the total loss of moment. To do so we utilize the framework developed in Refs. 16,21, and 38 which combines OCT with quantum simulations of spin dynamics.

BACKGROUND AND METHODS
We begin by briefly reviewing TDDFT and OCT, a more thorough discussion can be found here [16]. Then we review the results of Ref. 37 that showed ultrafast demagnetization in bulk ferromagnets (Fe,Co,Ni) for short, strong, laser pulses.

TDDFT
The electronic density is defined as where N is the total number of electrons, r is the spacial coordinate, t is the time, and Ψ is the wavefunction of the TDSE: for Hamiltonian:Ĥ =T +V ext +V ee composed of the kinetic energy,T , the electron-electron interaction,V ee , and the external potential,V ext , which includes both the electron-nuclear interaction and the electric fields of any laser pulses. We use atomic units throughout unless otherwise stated. TDDFT is founded upon the Runge-Gross theorem [12] which proves a 1 − 1 correspondence between the time-dependent density and the time-dependent external potential (up to a timedependent constant) for any electron-electron interaction. Hence all observables of the system are, in principle, unique functionals of the density. In particular, a non-interacting KS system [39] can be defined with a unique KS potential that reproduces the time-dependent density of the interacting system and thus predicts all observables of the true system without requiring the costly propagation of Eq. 2. The TDKS equation is: with the total density given by The KS potential, v S (r, t), consists of three pieces: where v H (r, t) is the usual Hartree potential of the instantaneous density, and v XC (r, t) is the exchange correlation (XC) potential and is a functional of the density at all previous times, the interacting initial state, and the non-interacting initial KS state. In practice, it must be approximated, with the most common approximation being the ALDA: which uses just the instantaneous density inputed into the ground-state DFT LDA XC functional and e unif XC (n) is the XC energy density of the uniform electron gas. The initial KS state is typically the ground-state found from a DFT calculation.
From this starting point, TDDFT has been extended to include non-collinear magnetism and magnetic fields [40]. For this case, we have a non-interacting Pauli KS Hamiltonian [41] which is used to propagate 2 component spinors, from which the density and magnetization density exactly replicate those of the interacting system. This is the formulation we will use for our simulations. The magnetization density operator may be written as: wheren(r) is the density operator and in the twocomponent spinors propagated in our calculations, S = g/2 σ where {σ x , σ y , σ z } are the familiar Pauli spin matrices and g is the electronic gyromagnetic ratio. For periodic boundary conditions, the total moment is then where Ω is a single unit cell. The KS Hamiltonian for our simulations is: (10) wherep is the momentum operator,Ŝ is the vector spin operator, and c is the speed of light. The laser pulse electric field is written as a vector potential, A ext (t) in the velocity gauge as it allows Bloch's theorem to be utilized. The KS magnetic field is written as is the magnetic field of the applied electromagnetic field and B XC (r, t) is the XC magnetic field. The ALDA can be extended to B XC using the LDA rotation method of Kübler [42]. The final term of Eq. (10) is the spin-orbit coupling (SOC) term, which can be thought of as the interaction between the spin of an electron and the effective magnetic field caused by relativistic motion thought a scalar potential. In a centrosymmetric potential, this term reduced to the well knownL ·Ŝ coupling. Propagation with Hamiltonian Eq. (10) is implemented in ELK [43], an all-electron electronic structure code, which was also used for all ground state and time-dependent calculations.

Optimal Control Theory
The central equation of OCT is the target functional G[u]: where u is the control field and Ψ[u] contains the information on how the system responds to the control field. In Quantum OCT (QOCT), Ψ[u] is then the wavefunction, which is a functional of the control field via the TDSE and from which any system observables to be controlled may be calculated. The target functional is generally separated into two pieces, J 1 [u] which contains information on the desired dynamics and J 2 [u] which is a penalty function in order to satisfy any constraints on the system or control field. The magnitude of the penalty functional is determined by how strongly a constraint must be satisfied.
Once a relevant target functional has been constructed, the goal of OCT is to extremize it and thus find the optimal control field u to best satisfy the balance between desired dynamics and the constraints. There are many choices for the algorithm to perform this optimization, some are general, such as the Nelder-Mead [44] or NEWOAU [45] algorithms, while some are developed for specific types of problem, e.g. in QOCT the ZBR scheme [46] adds a time dependent auxiliary wavefunction, which is also propagated in time and the overlap with the true wavefunction used to construct the control field.
For our system, we wish to maximize the loss of magnetic moment in a given time interval [0,T] while including practical and physical constraints on the type of laser pulse. Thus (t), the electric field of the laser pulse is the control field and is the target functional to be minimized, i.e. if we choose the initial magnetization M z (0) of the ferromagnet to be along the z-axis, then minimizing M z (T )/M z (0) will maximize the loss of moment. The constraints on the electric field are that the pulses satisfy Maxwell's equations (details below) and only certain frequencies are used to construct the pulse. The second constraint is of practical nature, as experimentally, pulses containing arbitrary frequencies cannot be constructed and often access to a single frequency (or multiples thereof) is only available. From Maxwell's equations, the following constraints on the electric field must be physically satisfied: Following Ref. 21, we can satisfy all these constraints by writing the electric field (t) = Nω n=1˜ n cos(ω n t) (15) where N ω is the number of frequencies to be used and n are the coefficients to be optimized. It can be seen that this choice automatically satisfies the constraints from Maxwell's equations. The frequencies used for our demonstration are Ultrafast Demagnetization In Fig. 1, the results of TDDFT simulations [37] for several laser pulses with differing peak intensities but the same profile are shown. A loss of the total magnetic moment was observed in all cases, with the fraction of moment lost dependent on the field intensity. It was also shown in Ref. 37 that the fraction loss depends on the frequency of the applied laser pulses. Hence the system was a strong candidate for optimal control. The purpose of this paper is test that hypothesis.
It should be pointed out that for less-intense longerduration pulses and for more realistic system geometries, the required peak intensity and fluence can be reduced by several orders of magnitude [47], however the underlying mechanism of demagnetization is the same. Thus we can demonstrate control of this process by focusing on the short strong laser pulses.

RESULTS
To demonstrate optimal control of the ultrafast demagnetization in Ni, we will attempt to maximize the loss of moment after T = 600 au = 14.4 fs. We will optimize a pulse of the form given by Eq. (15) using N ω = 4 different laser frequencies defined by Eq. (16). All calculations are performed with a time step of 0.1 au and 8x8x8 k-points. For the optimization we choose to use the gradient-free Nelder-Mead simplex algorithm. To initialize the Nelder-Mead algorithm, N ω + 1 starting points are required, it is instructive to examine these before moving to the results of the optimization.

Random Initial Pulses
To initialize our calculation, we construct 5 different pulses where the coefficients of Eq. (15) are chosen at random in a suitable range. This range is chosen such that the peak intensity is similar to the demagnetization pulses observed in Ref. 37. The pulses may be seen in the upper panel of Fig. 2, note that this is the vector potential which can be calculated from the electric field via The dynamics of M z (t) are shown in the lower panel and it can be seen that all pulses display demagnetization.
If we look at the final time, the average loss of moment is approximately 20%. If the optimal control is successful, then this percentage loss should be significantly increased.

Maximize Demagnetization
From the initial pulses, the Nelder-Mead algorithm then calculates a new set of coefficients from a simple set of rules and then tests how this affects the target functional by performing a TDDFT simulation with the laser field given by these coefficients. It then iterates this procedure and traverses the multidimensional parameter space, searching for the optimal set of coefficients.
In Fig. 3, we plot the ratio of the final moment after time T to the initial moment for each of the iterations. Although individual iterations can worsen the loss percentage, there is a clear downward trend as better and better pulses are found during the search, indicating that the optimal control is working. Each set of coefficients is a point in the parameter space, at each iteration, the Nelder-Mead algorithm reflects the worst point through the center of mass of the other points. Depending on whether this new point improves upon the next worst point, the algorithm can expand or contract in this direction, otherwise it can reduce all points towards the best point. This explains why individual points may worsen the ratio M z (T )/M z (0).
If we look at the result after approximately 30 iterations, the best pulse the optimal control procedure has found causes a 40% loss of moment. This is twice as good as the random initial pulses used to start the algorithm and is a clear demonstration that the moment can be successfully controlled using OCT. In Fig. 4 we show the electric field of this best pulse and also the magnetization dynamics, compared to the initial pulses. Examining the pulse shape compared to the initial pulses, there is no obvious reason why one leads to a larger demagnetization. This is the power of QOCT to find such pulses. We also observe that the magnetization dynamics is a highly nonlinear process, in particular for these short pulses, again demonstrating the need for QOCT in order to control the moment.

Fluence Constraint
For many practical reasons, the fluence of the applied pulse should be constrained. The simplest reason is simply efficiency, i.e. using a a pulse with lower energy to achieve the same dynamics as a higher energy pulse. Other reasons include surface damage to the material due to high fluence pulses, heating of the sample (and problems associated with cooling it), or physical restrictions on the laser itself preventing production of high fluence pulses. All of these present significant problems to future technological application, hence we wish to demonstrate how a fluence constraint can be incorporated into our calculations.
If we add to Eq. (11), the constraint which is proportional to the laser fluence. The free parameter α determines how strong the constraint is, for this calculation we choose α = 0.05. This parameter was based on examining the results of the previous optimization and choosing α to favor a lower fluence while still maintaining significant demagnetization in the set. In Fig. 5, we show the value of the total target functional, Eq. (11), for each iteration of the optimization algorithm. Unlike the previous case, we cannot attach a physical meaning to the target functional. so the actual value is not significant, only the trend. Furthermore, when choosing α, it was clear that the parameter space is a more complicated environment than the previous case, as a pulse could have the same value of Eq. (11) by either increasing the demagnetization or decreasing the fluence. Due to computational constraints, we stopped the optimization after 12 iteration, although this was sufficient to see the trend and demonstrate optimal control. We initialize the algorithm using the 5 best pulses found by calculating the target functional for all the previous pulses. These are not shown in Fig. 5, but were part of the optimization search. By using these to start the optimization, we save a large amount of computational time as opposed to using random pulses (although in general this may not be the case if the constraints significantly change the parameter landscape).
To see the power of this optimization, in Fig. 6 we plot the electric fields and the dynamics of M z (t) for two different pulses. These correspond to the best point of Fig.  5 and a reference pulse corresponding to the G = 1.1021 point. While this point is not the worst point, it was chosen as the final moment is similar to the best point, −0.537 and −0.555 respectively. Thus, the magnetic moments at time T are very close, however the fluences are very different. If we insert the pulses shown in Fig. 6 into the fluence formula given in Eq. 18, we find values of 6.437 a.u for the reference pulse and 2.437 a.u. for the best pulse. Therefore by using OCT we have reduced the required fluence by ≈ 60%.

CONCLUSIONS
In summary we have successfully demonstrated that optimal control of spin-orbit mediated ultrafast demagnetization is possible. For a short time interval, we showed that the loss of moment can be at least doubled (compared to randomly chosen typical pulses) for a system where the available laser frequencies (used to tailor the laser pulse) are constrained. Furthermore we extended the control problem to include a constraint on the laser fluence and demonstrated that QOCT could successfully find a pulse that balances the fluence and demagnetization requirements. Compared to a reference pulse, this optimal pulse produces almost identical magnetization dynamics, while reducing the fluence by over a factor of 2. Control of the system is of upmost importance for future technological application (for example, in spintronics), where the desired dynamics and constraints are dictated by practical concerns. Any physical phenomenon must be robust to these concerns, and as we have demonstrated, this form of ultrafast demagnetization meets this criteria. Simulation and QOCT of more complicated scenarios, such as longer pulse durations or further constraints on the fluence, intensity, and robustness of the demagnetization, can be build upon this foundation.