Optimal control of laser-induced spin–orbit mediated ultrafast demagnetization

The electronic density is defined as

$$\begin{eqnarray}n({\bf{r}},t) & = & N\displaystyle \int {\rm{d}}{{\bf{r}}}_{2}\ldots {\rm{d}}{{\bf{r}}}_{N}\ {{\rm{\Psi }}}^{*}({\bf{r}},{{\bf{r}}}_{2},\ldots ,{{\bf{r}}}_{N},t)\\ & & \times {\rm{\Psi }}({\bf{r}},{{\bf{r}}}_{2},\ldots ,{{\bf{r}}}_{N},t),\end{eqnarray} \tag{ 1 }$$

where N is the total number of electrons, ${\bf{r}}$ is the spacial coordinate, t is the time, and Ψ is the wavefunction of the TDSE:

$$\begin{eqnarray}&&{\rm{i}}\displaystyle \frac{\partial }{\partial t}{\rm{\Psi }}=\hat{H}{\rm{\Psi }}\end{eqnarray} \tag{ 2 }$$

for Hamiltonian:

$$\begin{eqnarray}&&\hat{H}=\hat{T}+{\hat{V}}_{\mathrm{ext}}+{\hat{V}}_{\mathrm{ee}}\end{eqnarray} \tag{ 3 }$$

composed of the kinetic energy, $\hat{T}$, the electron–electron interaction, ${\hat{V}}_{\mathrm{ee}}$, and the external potential, ${\hat{V}}_{\mathrm{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 [14] which proves a 1−1 correspondence between the time-dependent density and the time-dependent external potential (up to a time-dependent 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 [42] 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 equation (2). The TDKS equation is:

$$\begin{eqnarray}&&{\rm{i}}\displaystyle \frac{\partial }{\partial t}{\phi }_{j}({\bf{r}},t)=\left[-\displaystyle \frac{{\rm{\nabla }}}{2}+{v}_{{\rm{S}}}({\bf{r}},t)\right]{\phi }_{j}({\bf{r}},t)\end{eqnarray} \tag{ 4 }$$

with the total density given by

$$\begin{eqnarray}&&n({\bf{r}},t)=\displaystyle \sum _{j\quad =\quad 1}^{N}| {\phi }_{j}({\bf{r}},t){| }^{2}\end{eqnarray} \tag{ 5 }$$

The KS potential, ${v}_{{\rm{S}}}({\bf{r}},t)$, consists of three pieces:

$$\begin{eqnarray}&&{v}_{{\rm{S}}}({\bf{r}},t)={v}_{\mathrm{ext}}({\bf{r}},t)+{v}_{{\rm{H}}}({\bf{r}},t)+{v}_{\mathrm{XC}}({\bf{r}},t),\end{eqnarray} \tag{ 6 }$$

where ${v}_{{\rm{H}}}({\bf{r}},t)$ is the usual Hartree potential of the instantaneous density, and ${v}_{\mathrm{XC}}({\bf{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 adiabatic local density approximation (ALDA):

$$\begin{eqnarray}&&{v}_{\mathrm{XC}}[n]({\bf{r}},t)={v}_{\mathrm{XC}}^{\mathrm{LDA}}[n({\bf{r}},t)]={\left.\displaystyle \frac{{\rm{d}}{e}_{\mathrm{XC}}^{\mathrm{unif}}}{{\rm{d}}n}\right|}_{n=n({\bf{r}},t)}\end{eqnarray} \tag{ 7 }$$

which uses just the instantaneous density inputed into the ground-state DFT LDA XC functional and ${e}_{\mathrm{XC}}^{\mathrm{unif}}(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 [43]. For this case, we have a non-interacting Pauli KS Hamiltonian [44] which is used to propagate two 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:

$$\begin{eqnarray}&&\hat{{\bf{m}}}({\bf{r}})=\hat{{\bf{M}}}\;\hat{n}({\bf{r}}),\end{eqnarray} \tag{ 8 }$$

where $\hat{n}({\bf{r}})$ is the density operator and in the two-component spinors propagated in our calculations, $\hat{{\bf{M}}}=-g{\mu }_{B}\hat{{\bf{S}}}$ and ${\bf{S}}={\boldsymbol{\sigma }}/2$ where $\{{\sigma }_{x},{\sigma }_{y},{\sigma }_{z}\}$ are the familiar Pauli spin matrices and g is the electronic gyromagnetic ratio (g = 2 for our purposes). For periodic boundary conditions, the total moment is then

$$\begin{eqnarray}&&{\bf{M}}(t)={\displaystyle \int }_{{\rm{\Omega }}}{{\rm{d}}}^{3}r\ {\bf{m}}({\bf{r}},t),\end{eqnarray} \tag{ 9 }$$

where Ω is a single unit cell. The KS Hamiltonian for our simulations is:

$$\begin{eqnarray}{\hat{H}}_{{\rm{S}}}(t) & = & \displaystyle \frac{1}{2}{\left(\hat{{\bf{p}}}+\displaystyle \frac{1}{c}{{\bf{A}}}_{\mathrm{ext}}(t)\right)}^{2}+{v}_{{\rm{S}}}(\hat{{\bf{r}}},t)\\ & & +\displaystyle \frac{1}{c}\hat{{\bf{S}}}\cdot {{\bf{B}}}_{{\rm{S}}}(\hat{{\bf{r}}},t)+\displaystyle \frac{1}{2{c}^{2}}\hat{{\bf{S}}}\cdot ({\rm{\nabla }}{v}_{{\rm{S}}}(\hat{{\bf{r}}},t)\times \hat{{\bf{p}}}),\end{eqnarray} \tag{ 10 }$$

where $\hat{{\bf{p}}}$ is the momentum operator, $\hat{{\bf{S}}}$ is the vector spin operator, and c is the speed of light. The laser pulse electric field is written as a vector potential, ${{\bf{A}}}_{\mathrm{ext}}(t)$ in the velocity gauge as it allows Bloch's theorem to be utilized. The KS magnetic field is written as ${{\bf{B}}}_{{\rm{S}}}(\hat{{\bf{r}}},t)={{\bf{B}}}_{\mathrm{ext}}(t)+{{\bf{B}}}_{\mathrm{XC}}(\hat{{\bf{r}}},t)$, where ${{\bf{B}}}_{\mathrm{ext}}(t)$ is the magnetic field of the applied electromagnetic field and ${{\bf{B}}}_{\mathrm{XC}}(\hat{{\bf{r}}},t)$ is the XC magnetic field. The ALDA can be extended to ${{\bf{B}}}_{\mathrm{XC}}$ using the LDA rotation method of Kübler [45]. The final term of equation (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 reduces to the well known $\hat{{\bf{L}}}\cdot \hat{{\bf{S}}}$ coupling. Propagation with Hamiltonian equation (10) is implemented in ELK [46], an all-electron electronic structure code, which was also used for all ground state and time-dependent calculations.

The central quantity of OCT is the target functional $G[u]$:

$$\begin{eqnarray}&&G[u]=G[{\rm{\Psi }}[u],u]={J}_{1}[u]+{J}_{2}[u],\end{eqnarray} \tag{ 11 }$$

where u is the control field and ${\rm{\Psi }}[u]$ contains the information on how the system responds to the control field. In quantum OCT (QOCT), ${\rm{\Psi }}[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 [47], NEWOAU [48], or conjugate gradient [49] algorithms, while some are developed for specific types of problem, e.g. in QOCT the ZBR scheme [50] 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 $\epsilon (t)$, the electric field of the laser pulse is the control field and

$$\begin{eqnarray}&&{J}_{1}[\epsilon ]=\langle {\rm{\Psi }}[\epsilon ](T)| {\hat{{\bf{M}}}}_{z}| {\rm{\Psi }}[\epsilon ](T)\rangle ={M}_{z}(T)\end{eqnarray} \tag{ 12 }$$

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:

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{T}{\rm{d}}t\epsilon (t)=0,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\epsilon (0)=0=\epsilon (T).\end{eqnarray} \tag{ 14 }$$

Following [23], we can satisfy all these constraints by writing the electric field

$$\begin{eqnarray}&&\epsilon (t)=\displaystyle \sum _{n=1}^{{N}_{\omega }}{\tilde{\epsilon }}_{n}\mathrm{sin}({\omega }_{n}t){\mathrm{sin}}^{2}(\displaystyle \frac{\pi t}{T}),\end{eqnarray} \tag{ 15 }$$

where N_ω is the number of frequencies to be used and ${\tilde{\epsilon }}_{n}$ are the coefficients to be optimized. It can be seen that this choice automatically satisfies the constraints from Maxwell's equations if ${\omega }_{n}T$ is a multiple of $2\pi$. A convenient choice for the frequencies is

$$\begin{eqnarray}&&{\omega }_{n}=\displaystyle \frac{2\pi n}{T}\end{eqnarray} \tag{ 16 }$$

which will be used in the first of our optimizations. We emphasize that in a realistic control problem these frequencies are dictated by the experimental setup and not by the optimal control formalism. In our demonstration, they where simply chosen for convenience.

To initialize our calculations³ , we construct five different pulses where the coefficients of equation (15) are chosen at random in a suitable range. These pulses may be seen in the upper panel of figure 1. 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 7%. If the optimal control is successful, then this percentage loss should be significantly increased.

Zoom In Zoom Out Reset image size

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 figure 2, 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 not only applicable to this problem but also is crucial from obtaining higher demagnetization.

Zoom In Zoom Out Reset image size

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 38 iterations, the best pulse the optimal control procedure has found causes a 20% loss of moment. This is at least twice as much 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 figure 3 we show the electric field of this best pulse and also the magnetization dynamics, compared to an initial pulse. Examining the pulse shape compared to the initial pulses of figure 1, there is no obvious reason why one leads to a larger demagnetization. This is the power of QOCT to find such pulses. In [40], it was found that the demagnetization process has a highly nonlinear dependence on the pulse intensity of fluence. The optimal pulse shown in figure 3 has a significantly higher fluence compared to the initial pulse which combined with the nonlinearity of the problem leads to a larger loss of the magnetic moment. It was also found in [40] that demagnetization can occur for a period after the laser pulse. In this case we can see in figure 3 that if we simulate for a longer time, we will likely see further demagnetization with our pulses. However, the optimal control algorithm has no knowledge of this and simply maximizes the loss of moment at time T, regardless of what might happen after. It suggests for future optimizations, a better target functional might be to look at a time slightly later than when the pulse finishes or to average the moment during a period after the pulse.

Zoom In Zoom Out Reset image size

For practical reasons, the fluence of the applied pulse should be constrained. The primary reason is simply efficiency—achieving the same dynamics as a higher energy pulse but using a lower 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 include a fluence constraint into our calculations.

If we add to equation (11), the constraint

$$\begin{eqnarray}&&{J}_{2}[\epsilon ]=\alpha {\displaystyle \int }_{0}^{T}{\rm{d}}t\ {\epsilon }^{2}(t)\end{eqnarray} \tag{ 17 }$$

which is proportional to the laser fluence. The free parameter α determines how strong the constraint is, for this calculation we choose $\alpha =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 this example we choose a time period of T = 16.21 fs and frequencies of ${\omega }_{1}=0.075\;{\rm{au}}$, ${\omega }_{2}=\displaystyle \frac{9}{8}{\omega }_{0}$, ${\omega }_{3}=\displaystyle \frac{10}{8}{\omega }_{0}$ and ${\omega }_{4}=\displaystyle \frac{11}{8}{\omega }_{0}$. These also satisfy the constraints of equations (13) and (14) and approximately span the optical frequency range.

In figure 4, we show the value of the total target functional, equation (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 equation (11) by either increasing the demagnetization or decreasing the fluence. We again initialize the search using random coefficients. As the optimization is computationally expensive⁵ , we stopped the optimization after 26 iterations, although this was sufficient to see the trend and demonstrate optimal control.

Zoom In Zoom Out Reset image size

To see the power of this optimization, in figure 5 we plot the electric fields and the dynamics of M_z(t) for two different pulses. These correspond to the best point of figure 4 and a reference pulse corresponding to the G = 1.191 point. This point was chosen to demonstrate the need for optimal control as the total fluence is significantly more (8872 mJ cm⁻²) than the best point found (3665 mJ cm⁻²), yet yielded a worse final moment (33% loss compared to 50%). Thus we can use OCT to find better pulses that cost far less in terms of energy. Even comparing to the best random initial guess, the algorithm was able to reduce the fluence by 15% while achieving the same amount of demagnetization. As was detailed in reference [40], the fluence required to see demagnetization of this type may be decreased by several orders of magnitude when longer duration pulses are used. So while the fluences in this demonstration are quite large, we are confident the results will be similar if applied to longer duration pulses. Furthermore, in subsequent work [51] to [40], we have simulated a pulse based on experimental parameters (intensity of 3.8 × 1011 W cm⁻², fluence 8.05 mJ cm⁻², and a FWHM of 40 fs) and still see a 20% loss in moment due to spin–orbit effects. The underlying physics in this case is the same as for the short pulses. Thus we anticipate that the controllability found in this work will also be true for lower energy laser pulses.

Zoom In Zoom Out Reset image size