CC BY
29
THERMALLY INDUCED TRANSITIONS AND MINIMUM ENERGY PATHS FOR MAGNETIC SYSTEMS
E. V. Alexandrov ITMO University, St. Petersburg, Russia aloraman@live.com
PACS 75.10.Hk, 82.20.Kh, 02.60.-x DOI 10.17586/2220-8054-2015-6-1-57-62
Thermally induced magnetic transitions are rare events as compared with vibrations of individual magnetic moments. Timescales for these processes differ by 10 orders of magnitude or more. Therefore, the standard Monte-Carlo simulation is not suitable for the theoretical description of such phenomena. However, a statistical approach based on transition state theory is applicable for calculations of the transition rates. It presupposes finding the minimum energy path (MEP) between stable magnetic states on the multidimensional energy surface of the system. A modification of the Nudged Elastic Band (NEB) method for finding the energy barriers between states is suggested. A barrier on the energy surface corresponds to the difference between maximum energy along the MEP (highest saddle point) and the initial state minimum. The NEB procedure is implemented for spin rotations in Cartesian representation with geometric constraint on the magnitude of the magnetic moment. In this case, the effective magnetic forces are restricted to the tangent plane of the magnetic momentum vector. Keywords: Potential energy surface, minimum energy path, nudged elastic band method, numerical optimization, quick min-mode. Received: 31 December 2014
1. Introduction
In recent years, there has been large interest in the study of magnetization reversal processes under the action of thermal fluctuations and external perturbations. The stability of magnetic states is of great importance for research and development of miniaturized magnetic storage devices [1]. For estimation of the rates of activated transitions it is essential to know the response of the system to random perturbations. Direct Monte-Carlo simulation of temperature induced magnetization reversal is effectively impossible due to the time-scale difference between high-frequency vibrations of individual moments and low-frequency magnetic transitions between magnetic states of the whole system. However, a statistical approach known as 'transition state theory' [2,3] can be used in this case. This approach describes transitions through description of behavior of a system in the least observable 'transition state'. Stable states of the system correspond to the minima on the multidimensional energy surface of the system considered as a function of all parameters determining the magnetic configuration. For magnetic systems, these parameters have to determine all spin vectors. These could be Cartesian coordinates of the magnetic moments. However, due to the large difference in the time scale of longitudinal and transverse relaxation of the magnetic moment, the length can be taken to be fixed and thereby decrease the number of independent variables. Transition rates depend on the height of energy barriers, which correspond to saddle points on minimum energy paths (MEP), as well as on the shape of the energy surface at the minimum and saddle points [2]. A popular and efficient approach used for calculation of MEP for atomic rearrangements is the nudged elastic band (NEB) method [3]. A transition path is represented there by a chain of images (replicas) of the system between two local minima on the potential energy surface; each image
corresponds to specific configuration of the system with known potential energy. The images are connected with artificial springs to keep them well separated and prevent them from falling into the minima. The images give a discrete representation of the path on the energy surface. An iterative algorithm moves the images under the action of a force that contains the transverse component of the real force and parallel component of the spring force. Detailed description of this procedure, taking into account the constraint on the modulus of magnetic moments, is given below.
2. Method
2.1. Minimization technique
The system is defined by the Hamiltonian which gives the energy of the system for any configuration of the spin vectors. It can be either an analytical function of spin vectors, or numerical procedure, such as DFT calculations [4]. A constraint on the length of the spin vectors confining them to the surface of a sphere or a self-consistent determination of the modulus of magnetic moments require specific measures to be incorporated into the implementation of the NEB method. Common constrained-optimization techniques, such as Lagrange multipliers method, are not suitable in this case due to the computational cost which tends to be high for non-analytical Hamiltonians [3,5]. In addition, in most cases pure analytical methods for finding extrema are unsuitable due to the computational difficulty for evaluating the second derivatives of the Hamiltonian. Thereby, the implementation method itself should preferably incorporate these constraints. One way to do this is to use a coordinate system that is bounded to surface of the sphere, such as spherical or stereographical coordinates [5]. This approach is quite successful; however, it suffers from topological artifacts of a sphere, such as coordinate periodicity and problems with pole regions. As an alternative to using sphere-specific coordinates and unconstrained transformations, we propose to use Cartesian coordinates and sphere-specific constrained transformations, e.g. rotations. Using this approach in the optimization procedure allows one to overcome most of the topological artifacts of a sphere and prevent violation of native constraints.
We use a minimization technique known as 'quick-min' or 'velocity projection optimization' [4,5], which is essentially a gradient descent method, optimized with accumulation of velocity in direction of minimization movement. An iteration step can be formulated in four stages as follows:
where m and At are optimization parameters and x is Heaviside function. The method is inspired by the velocity Verlet integration of motion equation. The native constraint requires the force to be restricted to the tangent plane of the current spin vector. If so, movement on tangent plane can be calculated. An infinitesimal movement in the tangent plane is equal to movement on the surface of the sphere; thereby we can transform that shift into rotation in
movement plane. We can construct rotation matrix in axis-angle representation from n - normal to movement plane and B - rotation angle:
R (n, 0) =
( nX + cos 0 (1 — nX) nxny (1 — cos 9) — nz sin 0 nxnz (1 — cos 0) + ny sin 0^
nxny (1 — cos 0) + nz sin 0 n2 + cos 0 (1 — n^) nynz (1 — cos 0) — nx sin 0
\nxnz (1 — cos 0) — ny sin 0 ny nz (1 — cos 0) + nx sin 0 n^ + cos 0 (1 — n2)
(1)
The minimization step is performed by using matrix (1) in following form:
S-+1 = R (n-, B-) ■ Sj.
In principle, instead of rotation matrices, another type of object can be used - unit quaternions, also known as versors. From n, B, quaternion describing the same rotation is defined as follows:
0
0
q (s,v)=^cos2,n sin2J . (2)
If quaternion is defined as in (2), then minimization step can be written in following
form:
S-+1 = qi -Sj ■ qi 1,
where Si is converted to quaternion by adding unit scalar part to it.
Accurately implemented, both rotation matrix approach and quaternion approach have similar computational cost. Quaternions have major computational benefit: the composition of the rotations can be calculated by simple multiplication of corresponding quaternions, which is much faster than multiplication of rotation matrices, but the difference was not significant in the present applications. However, it seems worth looking into yet another spin representation for this task: each spin vector as unit vector of a surface of a sphere can be constructed by rotation from a selected start vector. Therefore, it is possible to use rotations as coordinates. It may be of large interest to implement the procedure described above solely using quaternion terms.
2.2. Calculation of minimum energy path
The NEB method is implemented as minimization of all intermediate images of the path, using the specially modified force. In addition to the true force, each image is affected by forces of artificial springs that connect it to its neighbors. This force for each spring is determined from length of arc segment between two spin vectors in different images. To prevent the 'corner-cutting' effect [3], nudging direction is introduced using 'upwind tangent" approach [6]. In this approach, nudging direction is chosen as direction to highest energy neighbor. After that, forces are projected dependent on that direction:
ft = Ft ± + FE ||.
With this modified force, the iterative minimization procedure moves the images to the MEP. The saddle point, which corresponds to maximum energy on MEP, can be estimated from two neighboring images with highest energy. In addition, 'climbing-image' (CI) modification [7] can be used: The highest energy image marked as 'climbing', is free from spring forces and parallel component of the true force is inverted to move that image higher along MEP to the saddle point.
3. Test cases
Two simple systems were chosen as test cases. Both represent a single atom with two degrees of freedom to make it easier to visualize. In that case potential energy surface can be shown as a two-dimensional surface in three-dimensional space. In these particular cases, spherical coordinates were used for plotting: 0 - azimuth angle and 9 - elevation angle. The potential energy surfaces, plotted as function of these angles, are shown in Fig. 1 and Fig. 2. In both cases there are two minima corresponding to 0 = 9 = 0 and 0 = n, 9 = 0. The initial path corresponds to the chain shown with empty circles and is chosen to be the same in both cases. Images of a system in that state were intentionally distributed non-equidistantly through predetermined randomized procedure to show that NEB procedure is not constricted to a limited set of initial path configurations. The final MEP, obtained using the algorithm described above, corresponds to the chain shown with solid circles. The saddle point is denoted by a cross mark.
Fig. 1. Potential energy surface for single spin with spin Hamiltonian (3). Calculation of MEP without Climbing Image. Saddle point is directly between minima. Maximums are in pole regions (9 = ±n/2)
part:
The first system is a single spin with Hamiltonian containing only internal anisotropy
H = ST
A 0
V°
0 \
0
1 0 0 -9
5.
(3)
For this system, the MEP is expected to be a straight line between the local minima. In Fig. 1, this case is displayed, with even number of intermediary images and without CI. The saddle point is estimated to be between the two images with highest energy on opposite slopes. Due to symmetry, the point in the center is the saddle point. Also, if an odd number of intermediate images were used, then the middle image would be expected to get to saddle point even without using CI. Note that in Fig. 1, the distances between adjacent images in the final configuration of the path are equal.
The second system is single spin with Gaussian Hamiltonian:
Fig. 2. Potential energy surface for single spin with Gaussian Hamiltonian (4). Calculation of MEP with Climbing Image. Saddle point is exactly in climbed image. There are maximums in pole regions (0 = ±n/2), and there is maximum directly between minima (^ = n/2, 0 = 0), which results in curved form of MEP
P
/-1\
0 0
w
H = Wi exp
ai
0.35
I -5, -5, 1, 3.5, 4, 4 }
1
0 0
0
-1 0
0
1 0
0 0
-1
0
0 1
(4)
where ai is the angle between Pi and spin vector.
For this system, the MEP is non-trivial; it is a curved line since now there is a maximum in the energy at the center. In Fig. 2, the results of a calculation which was carried out with CI turned on. The saddle point is estimated to be exactly at the climbing image. The usage of CI results in splitting the band in two parts, separated by saddle point, with a different distance between adjacent images on the two sides. Due to symmetry, if either number of images were odd or CI were turned off, then the images would have even distribution of distances between them. Note that in Fig. 1, there are two possible MEPs, and in Fig. 2, there are four. This illustrates that the NEB method finds the local MEP closest to the initial path, not necessarily the global MEP. The final configuration depends in that way on the initial path.
These cases show that the NEB implementation described above is in accordance with theoretical description [3] and behaves in the expected way. There are no problems specific to pole regions or periodicity of coordinates, which is good for application to magnetic systems. In the future, this implementation should be compared with other implementations in terms of speed of convergence. Also, it may be interesting to try a pure quaternion implementation.
Acknowledgements
This work was partially supported by RFBR grant No. 14-02-00102 and Nordic-Russian Training Network for Magnetic Nanotechnology (NCM-RU10121). Helpful discussions with Pavel Bessarab, Valery Uzdin and Hannes Jonsson are gratelfully acknowledged.
References
[1] Loth S., Baumann S., et al. Bistability in atomic-scale antiferromagnets. Science, 335, P. 196-198 (2012).
[2] Bessarab P.F., Uzdin V.M., Jonsson H. Harmonic transition-state theory of thermal spin transitions. Phys. Rev. B, 85, P. 184409 (2012).
[3] Jonsson H., Mills G., Jacobsen K.W. Nudged elastic band method for finding minimum energy paths of transitions. Proceedings of the Conference 'International School of Physics. Classical and quantum dynamics in condensed phase simulations', Lerici, Italy, July 1997, P. 385-404 (1997).
[4] Henkelman G., Johannesson G., Jonsson H. Methods for finding saddle points and minimum energy paths, in Progress on Theoretical Chemistry and Physics. Kluwer Academic Publishers, Dordrecht, the Netherlands, 2000, P. 269-300 (2000).
[5] Bessarab P.F., Uzdin V.M., Jonsson H. Potential energy surfaces and rates of spin transitions. Zeitschrift fur Physikalische Chemie, 227, P. 1543 (2013).
[6] Henkelman G., Jonsson H. Improved tangent estimate in the NEB method for finding minimum energy paths and saddle points. J. Chem. Phys., 113, P. 9978 (2000).
[7] Henkelman G., Uberuaga B., Jonsson H. A Climbing-Image NEB Method for finding saddle points and minimum energy paths. J. Chem. Phys., 113, P. 9901 (2000).
