Magnetic skyrmion annihilation by quantum mechanical tunneling

Magnetic skyrmions are nano-scale magnetic states that could be used in various spintronics devices. A central issue is the mechanism and rate of various possible annihilation processes and the lifetime of metastable skyrmions. While most studies have focused on classical over-the-barrier mechanism for annihilation, it is also possible that quantum mechanical tunneling through the energy barrier takes place. Calculations of the lifetime of magnetic skyrmions in a two-dimensional lattice are presented and the rate of tunneling compared with the classical annihilation rate. A remarkably strong variation in the onset temperature for tunneling and the lifetime of the skyrmion is found as a function of the values of parameters in the extended Heisenberg Hamiltonian, i.e. the out-of-plane anisotropy, Dzyaloshinskii–Moriya interaction and applied magnetic field. Materials parameters and conditions are identified where the onset of tunneling could be observed on a laboratory time scale. In particular, it is predicted that skyrmion tunneling could be observed in the PdFe/Ir(111) system when an external magnetic field on the order of 6T is applied.


Introduction
Non-collinear localized magnetic states are receiving a great deal of attention, in particular magnetic skyrmions which have been proposed as elements in future spintronics devices [1][2][3]. Along with interesting transport properties, skyrmions exhibit particle-like behavior and carry a topological charge enhancing their stability with respect to a uniform ground state, typically a ferromagnetic phase. A key issue is the lifetime of the skyrmions and its dependence on temperature and applied magnetic field. Various mechanisms for the annihilation of a skyrmion have been characterized by theoretical modeling of atomic scale systems, in particular collapse in the interior of the sample [4][5][6][7] and escape through the boundary of the magnetic domain [6][7][8]. Two skyrmions can also merge into one, as well as the reverse process where a skyrmion is transformed into a pair of identical skyrmions [9].
The activation energy for the various possible transitions can be calculated for a given spin Hamiltonian by finding the minimum energy path (MEP) from the local energy minimum corresponding to the skyrmion state to the final state minimum corresponding to the uniform state. The maximum in energy along the path corresponds to a first order saddle point on the energy surface. The MEP can be found using the geodesic nudged elastic band method [4] and the calculation accelerated by making use of knowledge obtained from previous calculations by focusing only on the region near the maximum [10]. For the classical over-the-barrier mechanism, the pre-exponential in the Arrhenius type rate expression can be estimated using harmonic transition state theory (HTST) for magnetic systems [11,12]. This has, for example, been done for skyrmions in a PdFe overlayer on Ir(111) surface [7,8], a system that has been studied extensively in the laboratory [13,14]. The challenge is to design materials where skyrmions are sufficiently stable at room temperature and still small enough to be used in spintronic devices. Theoretical calculations can help accelerate this development by identifying the dominant annihilation mechanisms and predicting how the stability of skyrmions depends on materials properties.
In most calculations, it is assumed that the system is able to overcome the energy barrier of the transition by thermal activation and this is typically a valid assumption. It is, however, possible that quantum mechanical tunneling brings the system from the metastable skyrmion state to the ground state. Tunneling in systems described by a single magnetic moment, within a macrospin approximation, has been studied extensively [17][18][19][20][21][22][23], in particular in the context of molecular magnets, both experimentally [24][25][26][27] and theoretically [18,22]. For example, the rate of magnetization reversals in Mn 4 monomer and dimer molecular magnets has been calculated as a function of temperature and excellent agreement obtained with the experimentally measured rates [26,27] both the high temperature classical regime and the onset temperature for tunneling, using a Hamiltonian parametrized from various experimental observations [21,22].
Since skyrmion stability is of central importance, it is important to have a way to estimate whether tunneling is an significant annihilation mechanism. An experimental observation of skyrmion tunneling would, furthermore, be an example of what is referred to as macroscopic quantum tunneling and is of considerable interest in the study of quantum phenomena.
Recently, the quantum mechanical nature of skyrmions has received some attention. Roldán-Molina et al [28] calculated the zero-point energy associated with quantum spin fluctuations and found that it can contribute up to 10% of the total Zeeman energy necessary to remove a skyrmion with an applied magnetic field. Diaz and Arovas [29] considered a model where skyrmions are formed by adding a local magnetic field over a circular spot, opposing a uniform magnetic field stabilizing the ferromagnetic phase so as to form a single skyrmion ground state. The rate of skyrmion nucleation by tunneling at zero temperature to this induced ground state was calculated using path integrals and a collective coordinate approximation. Psaroudaki et al [30] generalized the micromagnetic equations of motion to finite temperature using a path integral formalism and predicted a quantum mechanical addition to the effective mass of the skyrmion. Derras-Chouk, Chudnovsky and Garanin [31] derived a Lagrangian describing the coupled dynamics of skyrmion radius and chirality angle and estimated the tunneling rate of skyrmion collapse within an instanton approach. In their calculations the skyrmion is described as a Belavin-Polyakov (BP) soliton [32], an approximation corresponding to small DMI and crystalline anisotropy as compared to the exchange interaction. The BP soliton shape can be a good approximation for some systems but is not of general validity and the widely studied PdFe/Ir(111) system is one example where it does not apply, as illustrated below. In order to search more generally over materials parameters and experimental conditions to evaluate the possibility of observing skyrmion tunneling in the laboratory, it is important to use as general approach as possible to predict both the onset temperature for tunneling as well as the lifetime of the skyrmion at that temperature.
This article presents results of calculations of the onset temperature for quantum mechanical tunneling for a range of materials parameters, with special focus on the PdFe/Ir(111) system. The magnetic moments in the system are taken into account explicitly without introducing collective coordinates a priori. Unlike particle rearrangements where the negative eigenvalue of the Hessian matrix at the first-order saddle point on the energy surface can be used to estimate the onset temperature for tunneling, the onset temperature for magnetic systems, that are governed by the Landau-Lifshitz equation of motion, depends on the whole spectrum of eigenvalues. The method used here represents an extension of the methodology presented earlier for systems described by a single magnetic moment [21,22]. A region in parameter space is identified where the lifetime of the skyrmion is on a laboratory time scale at the onset temperature, thus identifying possible candidate materials for the observation of skyrmion tunneling. The prediction is made that skyrmion tunneling could be observed in PdFe/Ir(111) in the presence of an external magnetic field.

Methods
Within HTST for over-the-barrier transitions in magnetic systems, the mechanism and rate of thermally induced transitions are characterized by the first order saddle point on the energy surface, s † with energy E † , representing the highest point along the MEP connecting the initial state and the final state [11,12]. At a first order saddle point, the Hessian matrix has one negative eigenvalue. HTST can be used to estimate the lifetime of a metastable skyrmion state when the dominant annihilation mechanism is over-the-barrier transitions as is the case at high enough temperature. As the temperature is lowered, the rate of such transitions drops and eventually the temperature independent quantum mechanical tunneling becomes the  Instanton theory can be used to define and evaluate T c [33][34][35]. It corresponds to the highest temperature where a periodic solution of the equation of motion in imaginary time exists within an infinitesimal vicinity of the first-order saddle point on the energy surface. Such a path is referred to as an instanton and from its period, β, a corresponding temperature can be found as T = /k B β. The appearance of such a path as temperature is lowered can be deduced from an analysis of the Euclidean (imaginary-time) action. For a system with N spins with magnitude μ and orientation defined by a set of unit vectors s i , i = 1, . . . , N, the action is [36] where s = {s 1 , s 2 , . . . s N }, s k (τ ) is a closed trajectory, s k (−β/2) = s k (β/2), and H( s) is the Hamiltonian. The first term in this equation is the Berry phase and A k is referred to as Berry connection [37,38]. It is related to the area of on a sphere for each spin bounded by the trajectory and can be expressed as [36] A k ∂ τ s = 2 arctan ∂ τ s · (s × s 0 ) 2 + 2s · s 0 + ∂ τ s · s 0 where s 0 is some reference direction.
To find an instanton in the vicinity of s † , the action in equation (1) is expanded up to second order

Q[ s]
where s = s † + δ s and Q † = βE † . In addition to the boundary conditions, δs k (−β/2) = δs k (β/2) = 0, a normalization constraint needs to be added for each spin, δs k · s k = 0. At the saddle point δQ = 0. The second order variation of the action in the vicinity of the saddle point can be written as where∇ 2 H is the Hessian matrix evaluated at s † and the tilde specifies a restriction to the tangent space [9]. The linearized Landau-Lifshitz equation in imaginary time can be obtained by setting δ 2 Q = 0 and taking the cross product with s † from the left This can be written in matrix form as The matrix M has N pairs of complex conjugate eigenvalues, λ k . One of the pairs, chosen here to correspond to k = 1, is real valued but the other N − 1 pairs are purely imaginary λ k = ±iη k A discussion of the properties of the matrix M can be found in reference [39,40]. The general solution of equation (6) has the following form where the u k are the eigenvectors and λ k the eigenvalues of M, and c k are the expansion coefficients. The real pair of eigenvalues, λ 1 , corresponds to the periodic solution with eigenfrequency ω 1 = iλ 1 /μ and this is the instanton. The period of motion along this trajectory relates to the frequency as 1/β = |ω 1 |/2π and gives the onset temperature for tunneling as Equation (8) has the same form as the onset temperature in particle systems [35], but there is a significant difference. In the case of particle rearrangements, ω 1 = √ −λ 1 where λ 1 is the negative eigenvalue of the Hessian matrix at the first-order saddle point on the energy surface, whereas for magnetic systems, governed by the Landau-Lifshitz equation of motion, ω 1 depends on the determinant of the Hessian matrix, i.e. all the eigenvalues. For a single spin described by spherical polar coordinates θ and φ [21] where E † θθ , E † φφ , and E † θφ are second derivatives of the energy function. For a multi-spin system, it is better to work with Cartesian coordinates, as has been done in the present case, because of the problems that can arise if any of the spins points in the vicinity of the two poles.
Once the onset temperature for tunneling has been found, the transition rate at that temperature can be estimated using HTST [11]. The quantum mechanical effects on the transition rate can, however, be large already at T c . For example, in atomic rearrangements involving chemical reactions or diffusion events, the transition rate has in some cases been found to be a couple of orders of magnitude larger than the classical rate at T c [41]. Significant quantum effects can, therefore, be evident well above T c .

Model
The calculations presented here are based on an extended Heisenberg Hamiltonian for the system where D ij is the Dzyaloshinskii-Moriya vector lying in the plane of the lattice parallel to the vector pointing between two nearest neighbor sites i and j, thereby supporting Bloch type skyrmions, J is the exchange coupling parameter, K the out of plane anisotropy constant, and B the uniform external magnetic field applied perpendicular to the lattice plane. The sum ij includes distinct nearest neighbor pairs. The system consists of 2500 spins on a triangular lattice with lattice spacing α = 1 and periodic boundary conditions. This model can, in particular, represent well the PdFe/Ir(111) system. Parameter values obtained from density functional theory calculations [16] have been used in calculations of skyrmion lifetime using HTST and found to give results that are consistent with experimental measurements of PdFe/Ir(111) [8]. Figure 1 shows the calculated onset temperature, T c , for the tunneling of the skyrmion to the ferromagnetic state as a function of scaled Dzyaloshinskii-Moriya parameter, D/J, and scaled anisotropy parameter, K/J, when the applied magnetic field is B = 0.73B D where B D = D 2 /(μJ) is the critical field [42]. The value obtained for T c varies from 1 K to 4 K over the range of parameter values used. The shift in annihilation mechanism from over-the-barrier to tunneling could be observed by measuring the lifetime of skyrmions as a function of temperature. It is manifested by a break between a temperature dependent lifetime above T c to temperature independent lifetime below T c . The Arrhenius graph displaying the logarithm of the lifetime vs. inverse temperature shows an intersection between a straight line with a slope determined by the activation energy above T c to a line with zero slope below T c . In order for this change in mechanism to be detectable, the skyrmion lifetime at T c needs to be on laboratory time scale, i.e. on the order of seconds or minutes. It is, therefore, important to also estimate the lifetime of the skyrmion at T c and this is done here using HTST. Figure 2 shows the calculated lifetime as a function of scaled parameter values in the Hamiltonian. The size of the skyrmion is also indicated by the color scale. The stability of skyrmions is to a large extent related to their size [43], as can be seen from figure 2. The smaller the skyrmion, the shorter the lifetime. The lifetime changes remarkably strongly over this small range of parameter values. Since the lifetime is also a strong function of temperature, the window for observing the change in mechanism is limited to a narrow range of parameter values. But, this range does exist and it should be possible to find materials where such measurements can be performed.

The PdFe/Ir(111) system
A great deal of experimental effort has focused on skyrmions in the PdFe/Ir(111) system [13,15]. Parameter values obtained from density functional theory [16] are found to give results that are consistent with the experimental observations [8]. Figure 1 shows calculated results using those parameter values when the applied field is B = 3T. Under such conditions, the onset temperature for tunneling is 4 K but the lifetime is much too long at that temperature for experimental measurements of the lifetime. By increasing the applied magnetic field, the size of the skyrmion is reduced and the lifetime thereby shortened, as illustrated in figure 3. The variation of T c with respect to the strength of the anisotropy and the DMI is shown in more detail in figure 4. With an applied field of B = 6.4T, the lifetime is brought down to a couple of minutes and the onset temperature for tunneling is still not too low, about 1.5 K. The lifetime is an extremely strong function of the field. In a field of 6.5T, it is predicted to be 10 s. Our calculations, therefore, predict that tunneling of skyrmions could be observed in this system on a laboratory time scale by applying a strong enough magnetic field. Other materials may of course be better suited for the observation of skyrmion tunneling, but the PdFe/Ir(111) system is at least one possible candidate.
The DMI and crystalline anisotropy are relatively strong as compared to the exchange interaction in PdFe/Ir(111), as shown in figure 5. Therefore, they cannot be considered as small perturbations for this system. Also, even when the size of the skyrmion is reduced by applying a field of B = 8T, the exchange interaction has not reached the BP limit of 4π √ 3J. The BP soliton is, therefore, not a good approximation Figure 5. Various contributions to the energy of a skyrmion in PdFe/Ir(111) with respect to the energy of the ferromagnetic state as a function of applied magnetic field. Even for a field of B = 8T, the DMI has a larger magnitude than the exchange interaction and the total energy has not reached the Belavin-Polyakov limit of 4 √ 3πJ (indicated by a dashed line).
for this system and the more general methodology presented here without pre-determined shape function is needed in order to obtain accurate results.

Conclusions
The calculations presented here make a prediction regarding conditions and materials parameters that make it possible to observe a change in annihilation mechanism of skyrmion collapse from classical over-the-barrier to quantum mechanical tunneling, an interesting example of macroscopic tunneling. It is important to assess to what extent the tunneling mechanism needs to be taken into account when estimating the lifetime of skyrmions for practical applications. The method presented here includes all spins in the system without assuming a preconceived effective shape of the skyrmion or a predefined collective reaction coordinate. The onset temperature is estimated using an instanton approach and the lifetime of the skyrmion at that temperature is estimated from the harmonic approximation of transition state theory. This approach has previously been shown to work well in calculations of tunneling of a macrospin representing a molecular magnet but is extended here to systems of multiple spins.