Recurrent network classifier for ultrafast skyrmion dynamics

TL;DR

This paper presents a recurrent neural network trained via supervised learning to classify ultrafast skyrmion dynamics in two-dimensional magnetic systems, enabling accurate recognition of different magnetization processes.

Contribution

The study introduces a novel application of recurrent neural networks for classifying ultrafast skyrmion dynamics, including unseen data and finite temperature effects.

Findings

High classification accuracy across parameter ranges

Effective on unseen data including finite temperature processes

Potential for autonomous control in skyrmion-based devices

Abstract

By using the supervised learning we train a recurrent neural network to recognize and classify ultrafast magnetization processes realized in two-dimensional nanosystems with Dzyaloshinskii-Moriya interaction. Our focus is on the different types of skyrmion dynamics driven by ultrafast magnetic pulses. Each process is represented as a sequence of the sorted magnetization vectors inputted into the network. The trained network can perform an accurate classification of the skyrmionic processes at zero temperature in wide ranges of parameters that are the magnetic pulse width and damping factor. The network performance is also demonstrated on different types of unseen data including finite temperature processes. Our approach can be easily adapted for creating an autonomous control system on skyrmion dynamics for experiments or device operations.

Equations

$$\begin{array}{r}H=-\sum _{i\ne j}{J}_{ij}{\mathbf{S}}_i{\mathbf{S}}_j-\sum _{i\ne j}{\mathbf{D}}_{ij}[{\mathbf{S}}_i\times {\mathbf{S}}_j]-\sum _{i}K({\mathbf{S}}_i^z{)}^2.\end{array}$$

$$\begin{array}{r}\frac{d{\text{S}}_i}{dt}=-\frac{\gamma }{1+{\alpha }^2}{\text{S}}_i\times [-\frac{\partial H}{\partial {\text{S}}_i}+b_i(t)]-\\ -\frac{\gamma }{|{\text{S}}_i|}\frac{\alpha }{1+{\alpha }^2}{\text{S}}_i\times ({\text{S}}_i\times [-\frac{\partial H}{\partial {\text{S}}_i}+b_i(t)]),\end{array}$$

$$\begin{array}{r}{\text{B}}_p(t)=B_{0}exp\left(-\frac{(t-t_p{)}^2}{2t_w^2}\right){\text{e}}_B,\end{array}$$

$$h_j^t=\sigma (y_j^t)=\frac{1}{1+e^{-y_j^t}},$$

$$y_j^t=\sum _{k=1}^N{x}_k^t{V}_{kj}+\sum _{m=1}^{N_h}{h}_m^{t-1}{U}_{mj},$$

$$o_i=softmax(z_i)=\frac{e^{z_i}}{{\sum }_{n=1}^{N_o}{e}^{z_n}},$$

$$z_i=\sum _{j=1}^{N_h}{h}_j^t{W}_{ji},$$

$$E=-\sum _{i=1}^{N_o}{t}_i\log o_i,$$

$$W_{ji}=W_{ji}+\frac{\delta E}{\delta W_{ij}},$$

$$V_{kj}=V_{kj}+\sum _{p=0}^{N_t-1}\frac{\delta L_j^{t-p}}{\delta V_{kj}},$$

$$U_{mj}=U_{mj}+\sum _{p=0}^{N_t-1}\frac{\delta L_j^{t-p}}{\delta U_{jm}},$$

$$\frac{\delta E_i}{\delta W_{ij}}=(t_i-o_i)h_j,$$

$$L_n^{t-p}=\underset{p}{\underbrace{\sum _{j,m,l,...,n=1}^{N_h}}}\sum _{i=1}^{N_o}{E}_i{W}_{ij}\underset{p}{\underbrace{(U_{jm}{U}_{ml}...U_{cn})}},$$

$$\frac{\delta L_j^t}{\delta V_{kj}}=\sum _i(W_{ji}^2\frac{\delta E_i}{\delta W_{ij}})(1-h_j^t)x_k^t,$$

$$\frac{\delta L_j^t}{\delta U_{mj}}=\sum _i(W_{ji}^2\frac{\delta E_i}{\delta W_{ij}})(1-h_j^t)h_k^{t-1},$$

$$\frac{\delta L_n^{t-p}}{\delta V_{kn}}=\sum _q^N\sum _r^{N_h}(U_{nr}^2\frac{\delta L_r^{t-p+1}}{\delta V_{qr}}{x}_q^{t-p+1})h_n^{t-p}(1-h_n^{t-p})x_k^{t-p},$$

$$\frac{\delta L_n^{t-p}}{\delta U_{ln}}=\sum _{q,r}^{N_h}(U_{nr}^2\frac{\delta L_r^{t-p+1}}{\delta U_{qr}}{h}_q^{t-p+1})h_n^{t-p}(1-h_n^{t-p})h_l^{t-p-1}.$$

Full text

I Introduction

Understanding of relaxation processes in magnetic materials and systems due to the ultrafast pulses realized with different experimental techniques plays a crucial role for creating next-generation magnetic storage, spintronics and quantum technologies trotzky ; ultrafast1 ; ultrafast2 ; ultrafast3 . Nowadays pico- and femtosecond magnetization dynamics of a system is accessible with X-ray free electron lasers XFEL and Lorentz microscopy Lorentz . At these time scales it becomes possible to trace the behaviour of magnetic systems frame by frame or spectrum by spectrum. The amount of such a data grows rapidly and one of the important problems is classification of the dynamical processes. The main specific here is that each experimental frame is a part of a particular excitation process and can not be analyzed separately, which makes this problem to be similar to classification of video content video1 .

In material science a special focus is on the magnetic skyrmions that are topologically protected multi-spiral magnetic excitations formed in two- and three-dimensional materials. Each skyrmion is a complex object with non-collinear ordering of the atomic magnetic moments stabilized by the Dzyaloshinskii-Moriya interaction. In the static case the magnetic state recognition and classification problems can be effectively solved neither with standard methods of machine learning such as support vector machine, k-means and others or with neural network approach Iakovlev1 ; Iakovlev2 . The latter is very useful for analysis of the transitional areas between different phases, for instance, between skyrmion and spin spiral or between skyrmion and ferromagnetic states. Nevertheless, the elementary functions (erase and write) with skyrmions for next-generation electronics are entirely related to dynamics of skyrmions in external electric or magnetic fields Nagaosa . A practical implementation of these elementary functions requires an autonomous control on the system’s dynamics, which motivates us to develop a neural network approach for recognition and classification of dynamical skyrmionic processes.

In our work we follow the idea ”classify a movie rather than snapshot” Evert and propose a recurrent neural network classifier for ultrafast skyrmion processes of picosecond scale in two-dimensional materials. A simplified scheme of our approach is presented in Fig. 1. A stream of magnetic process frames is inputted into a recurrent neural network (RNN) frame by frame. We show that the process classification can be performed by tracing only the dynamics of the $z$ components of spins. The trained network provides robust results independent on the frame time frequency and pulse details during data collection. Our study paves the way to construct precise process diagrams for magnetic systems in wide range of external parameters.

II Model, Method and Network

In this work we use the following spin Hamiltonian to describe two-dimensional materials hosting skyrmions

[TABLE]

Here $J_{ij}$ and ${\bf D}_{ij}$ are the isotropic exchange interaction and Dzyaloshinskii-Moriya vector, respectively. ${\bf S}_{i}$ is a unit vector along the direction of the $i$th spin, $K$ is the strength of the uniaxial anisotropy in the $z$ direction. We take into account the only interactions between the nearest neighbours. The summation for inter-spin couplings runs twice over every pair. The Hamiltonian is defined on the 100$\times$100 square lattice. In our simulations we used the following parameters: $J=0.03676$ mRy, $|{\bf D}|=0.008824$ mRy and $K=0.00735$ mRy. Dzyaloshinskii-Moriya interaction vector for each pair of spins is parallel to the corresponding inter-site radius vector. Such a Hamiltonian was previously used in Ref. Blugel, to simulate switching of magnetic skyrmions by picosecond magnetic field pulses via transient topological states.

Since we are interested in simulation of dynamical processes, to solve Eq. (1) we used the Landau-Lifshitz-Gilbert (LLG) equation as it was implemented in the Uppsala Atomistic Spin Dynamics (UppASD) package UppASD1 ; UppASD2

[TABLE]

where $\gamma$ is the gyromagnetic ratio, $\alpha$ is the damping parameter and $b_{i}(t)$ is a stochastic magnetic field with a Gaussian distribution arising from the thermal fluctuations. In principle, one can train a deep neural network to simulate the magnetization dynamics of the system, as it was done in Ref. Kovacs, .

To realize different processes related to skyrmion dynamics we used a time-dependent magnetic field pulse defined by a Gaussian distribution as it was proposed in Ref. Blugel,

[TABLE]

where $B_{0}$ is the amplitude of the magnetic field, $t_{w}$ is the Gaussian width and $t_{p}$ is the time position of the pulse maximum. The real-space orientation of the magnetic pulse $\textbf{e}_{B}$ is described by the polar and azimuthal angles $\theta$ and $\varphi$. We used $B_{0}=2$ T and $\varphi=0^{\circ}$. Other parameters were varied depending on the type of the simulation, as described below.

Recurrent Neural Network— The magnetic configurations calculated with LLG approach at different times of the relaxation process are used as input for the recurrent neural network. We have found that it is of crucial importance to perform a data preprocessing that is ordering of the magnetization vector containing only z components of spins as shown in Fig. 2. Namely, such a sorting provides a very accurate separation of the magnetic phases in the static case as discussed in Ref. Iakovlev2, . It also gives us opportunity to use a moderate number of the hidden layer neurons that is 512 for process classification. Without sorting one needs to substantially increase the number of the hidden layer neurons up to 5000 to get reasonable results.

To train the recurrent neural network we use magnetic processes generated with the $\theta=50^{\circ}$ magnetic pulses. For the training set we generated 400 configurations for each of breathing, switching and collapse processes from the internal areas of the corresponding diagram. A complete description of the technical details concerning the neural network structure and training is presented in Appendices.

III Process diagram

Fig. 3 gives the switching diagram classifying the skyrmion dynamics in the system depending on the $\theta=40^{\circ}$ magnetic pulse width and damping factor. To construct such a diagram we used a 44$\times$40 grid (1760 points in total) on the $t_{w}$/$\alpha$ plane. The process corresponding to each point of the diagram was classified manually. In accordance with the results of Ref. Blugel, there are four different types of the skyrmion switching processes: breathing, collapse, switching and mixed. Examples of the collapse, switching and breathing processes are presented in Fig. 4. The breathing phase is related to a gradual change of the skyrmion size over time. As a result the initial and final profile frames are close to each other. It is not the case for switching and collapse processes which are characterized by a global change in the system spins orientation. One can see that these processes are related to the evolution of single skyrmion in the system and clearly distinguishable on the level of the sequences of the profiles. Following Ref. Blugel, the processes for which the final state composed of multiple skyrmions and domain walls are classified as mixed state. The latter will be discussed below.

Neural network classification. The same Fig. 3 contains processes classification performed by the trained recurrent neural network. One can see that the RNN provides very accurate classification of the breathing and collapse relaxation processes. On the other hand the neural network classifier shifts the boundary between mixed and switching processes to smaller values of the damping parameter. Within the manual classification discussed above a switching process will be associated with the mixed phase either if the system hosts more than one skyrmion or the final configuration is characterized by a combination of the skyrmion and domain wall. Some examples of the processes that we define as mixed ones are represented in Fig. 5.

Since the neural network does not distinguish the skyrmion configurations with different topological charges (Fig. 5 top panels), the switching processes accompanied with creation of two or more skyrmions are classified as switching ones, which is not the case for manual classification introduced in Ref. Blugel, . Thus, the neural network has refined the process diagram boundaries. However, the network is able to identify correctly configurations which are far from the skyrmion ones. They are mostly located at small damping parameter and between clear breathing and switching phases. Typical example is given in Fig. 5 bottom panels. For such processes the values of output neurons are less than the threshold which is equal to $0.8$.

A narrow parameters area of the process diagram between switching phase defined manually and mixed phase defined by the RNN is of technological importance since one can create two and more skyrmions from single skyrmion by magnetic pulses. Previously, the duplication mechanism of magnetic skyrmions was found and analyzed in Ref. Blugel2, . We found that there are many small regions in phase diagram Fig. 3 in which it is possible to observe the creation of two skyrmions. The biggest one is located at the parameters $\alpha\in(0.09;0.12)$ and $t_{w}\in(16;32)ps$.

One of the main advantages of the RNN is the ability of classification of the whole time sequence of the magnetization frames instead of individual ones. To examine whether it works in our case, we calculated the values of the output neurons at each time step and analyzed the following characteristics for different processes: output neuron values, components of the total magnetization and skyrmion number. As can be seen from Fig. 6 (top panels) at the beginning all the configurations are classified as the breathing state because we used $softmax$ function for the output layer and the sum of all outputs has to be equal to unity. Besides, initial state has the negative magnetization which is one of the key features of the breathing process. All changes start at the fourth time step, which corresponds to 45 ps in the considered case. In other words, the RNN responds to the magnetic pulse. It is also clearly seen that the values of the output neurons change slowly step by step, which means that network relies not only on the current magnetic configuration.

Since we used the only $z$ components of spins the question arises whether the RNN with single hidden layer performs classification on the basis of the total $M_{z}$. We address it by plotting the magnetization and topological charge behaviour for the same processes. One can see from Fig. 6 (middle panels) that the switching and collapse have almost the same $M_{z}$ profiles while the network outputs are significantly different. Indeed, for collapse process we have a large value on the switching neuron at time steps 5 and 6 before vanishing the switched skyrmion. However, we can traced out the similar behaviour of topological charge which was not given to RNN. It gives us the opportunity to conclude that the network learns some topological parameters of the magnetic configurations.

IV Network tests

Variation of the time scale. Since we have trained the RNN on the processes with the fixed number of the time steps and center of the Gaussian magnetic pulse another interesting question whether we can apply it to classify processes obtained at different parameters. The answer is yes we can! We found that it does not matter how many time steps we have in our process and when we act on the system with the magnetic pulse. To illustrate this we plotted the network outputs for the same processes we discussed in previous paragraph, but with the parameters $t=90$ ps and $\Delta t=5$ ps instead of $t=105$ ps and $\Delta t=15$ ps. Since we used the same $t_{p}=40$ ps in both cases maximum of the magnetic pulse moved from the fourth to the eighth time step. As can be seen from Fig. 7 all the dependencies have the same behaviour as in Fig. 6 (top panels). It gives us the opportunity to use our approach for classification of the experimental ultra-fast processes.

Temperature effect. The results we present above were obtained at the zero temperature. We found that the trained network can perform a robust and accurate classification of the magnetic processes up to 1.5 K. The example of this is presented in Figs. 8 and 9. Taking into account that the transition to the paramagnetic state is at about 6 K, we conclude that one can use the zero-temperature trained network for analysis of the low-temperature processes. It follows from Fig. 9 that within chosen range of magnetic pulse width values the process diagram changes its shape at the finite temperature. The transition from breathing to switching phase shifts towards smaller values of $t_{w}$. At high values of the pulse width ($t_{w}>$ 40 ps) we observe the formation of the collapse phase instead of switching one at T = 0.

V Conclusions

The neural network approach we propose is general nature and can be used for classification of different dynamical processes realized in real magnetic materials and systems. For instance, the first problem closely related to that we consider here is the refinement of the process boundaries for vortex dynamics reported in Ref. vortex, . It is also not limited by two-dimensional case, the profile procedure we used for preprocessing magnetization data gives reliable classification results in the static case for three-dimensional magnets with Dzyaloshinskii-Moriya interaction Iakovlev2 . The generalization power of the network can be increased by adding more hidden neurons layers. In this case one would expect that the network will correctly classify processes at high temperatures close to the paramagnetic phase.

VI Acknowledgements

We would like to thank Yaroslav Kvashnin and Anders Bergman for fruitful discussions and technical assistance with UppASD package. This work was supported by the Russian Science Foundation Grant 18-12-00185.

Appendix A Network details

The recurrent neural network consists of three layers: input, hidden and output. The hidden layer is self cycled, which means that the states of the hidden neurons pass forward to the next steps as shown in Fig. 10. There are 3 weight matrices, $V$ is for connecting input layer and hidden layer neurons, $W$ is for connecting the hidden layer and output layer neurons, and $U$ is for connecting hidden layer neurons at different time steps. Thus, the states of the hidden layer neurons depend on the current input magnetization data and previous state of itself. The input layer is characterized by 10000 neurons, as the number of spins on the lattice. For the hidden layer we have found out that the minimal number of the neurons is 512. If the hidden layer contains less neurons, RNN is not able to reach acceptable error rate during the training procedure. The output layer has 3 neurons, by the number of the processes we want to classify.

Appendix B Machine learning details

For the training set we generated 400 configurations for each of breathing, switching and collapse processes from the internal areas of the corresponding diagram Fig. 11. Such a small amount of data is due to the fact that we used the LLG equation (2) to simulate the magnetization dynamics. Results of such simulations does not changes from time to time at fixed parameters in contradistinction to the Monte Carlo ones.

As an input for our recurrent neural network (Fig. 10) we used the sorted vectors containing $S^{z}$ components of a two-dimensional magnetic configurations. The hidden layer neurons activate by means of the sigmoid function,

[TABLE]

[TABLE]

where $x_{k}=S_{i}^{z}$ is the value of $i$th input neuron, $V_{kj}$ is the weight between the $k$th input neuron and $j$th hidden neuron, $U_{mj}$ — weight between the $m$th hidden neuron at time step $t-1$ and $j$th hidden neuron at time step $t$, $N=L\times L$ — number of the input neurons, $N_{h}$ is the number of the hidden neurons.

In turn, for output neurons we use the softmax function that is given by

[TABLE]

[TABLE]

where $N_{o}$ is the number of the output neurons, $W_{ji}$ — weight between the $j$th hidden neuron and $i$th output neuron.

During the learning process, we randomly choose 10% of training set for cross-validation to avoid overfitting and define the stopping point where error is less than the required value. To evaluate the error we used cross entropy function given by the following expression:

[TABLE]

where $t_{i}$ represents the ground truth labels.

The network optimization was fulfilled through back-propagation through time algorithm bptt . All weight matrices were initialized by means of the ratio of 1 divided by a number of neurons in the previous layer to avoid extremely large sigmoid arguments: $V_{kj}=\frac{1}{10000}$, $W_{ji}=\frac{1}{512}$, $U_{mj}=\frac{1}{512}$. We used the following expressions for new weights in order to avoid getting stuck in a local minima

[TABLE]

[TABLE]

[TABLE]

where $N_{t}$ is the total number of time steps in our process, $p$ shows the index of the previous time step with respect to the final one. In these expressions $\frac{\delta E}{\delta W}$, $L_{j}^{t-p}$, $\frac{\delta L}{\delta V}$ and $\frac{\delta L}{\delta U}$ are given by the equations below. All the additional indices here we introduced to show the summation rules. For simplicity the equations for $\delta L^{t-p}_{n}$ were written in recurrent form.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]