On quantifying the topological charge in micromagnetics using a lattice-based approach

An implementation of a lattice-based approach for computing the topological skyrmion charge is provided for the open source micromagnetics code MuMax3. Its accuracy with respect to an existing method based on finite difference derivatives is compared for three different test cases. The lattice-based approach is found to be more robust for finite-temperature dynamics and for nucleation and annihilation processes in extended systems.


Introduction
The topological charge or skyrmion number associated with an O(3) field, m(r) = 1, is given by This quantity is used to characterize the topology of spin textures such as vortices and skyrmions in two-dimensional systems (see, e.g., Ref. [1]), where m represents the orientation of the magnetic moments. When m(r) is projected onto the unit sphere, Q measures the number of times the moments wrap around the surface of this sphere. For vortices and merons, Q = ±1/2, while for skyrmions, Q = ±1. Higher-order halfand full-integer charges are also possible. In numerical micromagnetism, a common approach involves discretising m(r, t) using the method of finite differences [2][3][4]. The underlying assumption is that cell-to-cell variations in m are sufficiently small such that the exchange energy, approximated to lowest order as (∇m) 2 , remains meaningful. Issues can arise under certain conditions, such as in the nucleation and annihilation of vortices and skyrmions, or in the stochastic dynamics with random fields, where large spatial variations in m can occur which reduce the accuracy of the finite-difference approximations of Eq. 1 and result in nonphysical values of Q. Consider an isolated skyrmion. Fig. 1(a) shows the equilibrium profile computed with the mumax3 code [4] and the parameters in Ref. [5]. The corresponding map of m onto the unit sphere is given in Fig. 1(b), where dots represent the orientations of m and the lines indicate bonds between nearest-neighbour finite difference cells [6,7]. The entirety of the sphere is covered by this mesh, which remains intact everywhere and reflects the fact that the spin texture in Fig. 1(a) possesses a nontrivial topology. Eq. (1) gives Q = −0.99978290 for this configuration, which is acceptably close to the theoretical value of Q = −1.
Consider now the effect of disorder, e.g., due to thermal fluctuations, where each moment is deviated away randomly from its equilibrium orientation in Fig. 1(a), as shown in Fig. 1(c). The corresponding map onto the unit sphere for this disordered case is presented in Fig. 1(d). While the mesh appears distorted, it retains the same topology as the case in Fig. 1(d) and therefore possess an identical charge. However Eq. (1) gives Q = −0.97115153 in this case, which reflects a loss in accuracy of the finite difference derivatives.
In this note, we discuss a lattice-based approach for computing Q that does not require rely on spatial derivatives. We discuss two different implementations of this scheme for finite difference micromagnetics and provide three examples against which these implementations are tested.

Lattice-based implementation for finite difference schemes
We follow the approach of Berg and Lüscher [8], which has been employed in atomistic spin dynamics and Monte Carlo simulations [9,10]. Consider the four moments in Fig. 2(a), each of which represent the average magnetization orientation in a finite difference cell. We treat these moments as lattice spins and set aside all aspects related to the interactions between them. Fig. 2(a) represents one unit cell of this lattice. The topological charge is given by the sum over the ensemble of elementary signed triangles q ijk on the unit sphere, where tan which is invariant under a cyclic permutation of the indices ijk. Fig. 2(a) shows two of such signed triangles that make up the unit cell, q 124 (grey) and q 234 (white). Fig. 2(b) represents another definition that is equally valid. ijk in Eq. (2) indicates that the summation is restricted to unique triangles as shown in Figs. 2(a) or 2(b). Fig. 2(c) illustrates a variation of this scheme that allows a local charge density analogous to Eq. 1 to be defined at a site (i, j), which is commensurate with the coordinates of the finite difference cells in which m i,j is defined. The method involves averaging over the two unit cells comprising the four triangles spanned by (i, j) with its nearest-neighbour spins, (i + 1, j), (i, j + 1), (i − 1, j), and (i, j − 1). This approach uses both of the conventions in Figs. 2(a) and 2(b), and takes the average of the two, thereby assigning a weight of 1/2 to each triangle q ijk . For finite-sized systems, the same averaging procedure cannot be applied at curved boundary edges because not all signed triangles are present, i.e., only one of the definitions, Fig. 2(a) or Fig. 2(b), produces the necessary orientation to cover the three spins that comprise the boundary, e.g., the blue triangles in the top left and bottom right of Fig. 2(d). In such cases, we assign a weight of 1 to the isolated signed triangle.
We provide an implementation of this method for mumax3 through the extension ext topologicalchargelattice which is publicly available in mumax3. 10 [11]. This extension provides a local charge density at site (i, j) in units of m −2 , analogous to the quantity provided by the finite-difference implementation of Eq. 1 through the extension ext topologicalcharge, which is obtained by dividing the q ijk by the surface area of the unit cell.  Gilbert damping of α = 0.3. These parameters model a three monolayer-thick Co film with a Curie temperature of 550 K [12]. Dipolar interactions are neglected for simplicity. The evolution of Q(t) over 100 ns is presented in Fig. 3 for four different temperatures, where an adaptive time-step integration method is used to solve the stochastic Landau-Lifshitz equation [13]. Q is computed at 5 ps intervals using finite difference derivatives [Eq. (1)] as implemented by the existing ext topologicalcharge extension in mumax3, and with the lattice-based implementation in ext topologicalchargelattice. Large fluctuations are seen in the Q(t) computed with finite difference derivatives, whose distribution spreads as the temperature increases as shown by the histograms in  Fig. 3. These represent thermally-driven nucleation and annihilation of meron and skyrmion states, respectively.

Soliton pair generation in a ferromagnetic track
We turn our attention to nucleation of skyrmion-antiskyrmion pairs due to spin-transfer torques [14,15]. The geometry comprises a 1000 × 125 × 0.6 nm film with the same magnetic parameters as in Sec. 3.1, except for D = 0.1 mJ/m 2 and α = 0.05. A nucleation zone is defined by a 25-nm diameter circular region of the track, in which the uniaxial anisotropy is oriented along y instead with K u = 0.5 MJ/m 3 . A conventional current flows along −x with a density of 25 TA/m 2 and a spin polarisation of P = 1, with nonadiabatic terms being neglected. The spin-transfer torques, combined with the nonuniform effective fields seen at the nucleation zone, result in skyrmion-antiskyrmion pairs being shed from this site, which then undergo Kelvin motion and propagate along the x direction before separating and annihilating. Figure 4 presents Q(t) and snapshots of the micromagnetic state. Four different cell sizes are considered to test the relative accuracy of Eq. (1) with respect to Eqs. (2) and (3). For the smallest [ Fig. 4(a)], there is good agreement between the two methods where only a handful of points with noninteger Q are obtained with Eq. (1), which occur at the transitions involving the nucleation and annihilation of (anti)skyrmions. As the cell size is increased [ Figs. 4(b)-(d)], a greater number of noninteger Q is obtained with Eq. (1), with smooth variations observed in Figs. 4(c) and 4(d). Meanwhile, the latticebased approach provides clear plateaus in Q close to integer values for all cases, which suggests that the smooth variations in noninteger Q are related to the loss in accuracy of the finite difference derivatives. The nucleation events differ between the four cases because the circular nucleation zone is discretised differently. Fig. 4(e) shows snapshots of the micromagnetic state at different instances where nucleation, Kelvin motion of skyrmion-antiskyrmion pairs, and (anti)skyrmion annihilation can be seen.

Isolated skyrmion in confined structures at finite temperatures
In systems with DMI, boundary edges result in a tilt in the background magnetization away from the z-axis, even in a nominally uniformly-magnetized state, as a result of chiral boundary conditions [16,17]. The tilt orientation is determined by the sign of D and the resulting Q can take on noninteger values. We consider a disc, 200 nm in diameter and 0.6 nm in thickness, which is discretised with 256 × 256 × 1 finite difference cells (all other magnetic parameters are identical to those in Sec. 3.1). A disc with an isolated skyrmion at T = 0 K is found to give Q of −0.9304853 using ext topologicalcharge and −0.9468352 using ext topologicalchargelattice. Deviations from Q = −1 represent the contribution from the edge magnetization tilts.  Q(t) for this disc is shown in Fig. 5 for four different temperatures. In contrast to Fig. 3, both the derivative-and lattice-based methods give fluctuations in Q, albeit to a lesser extent for the latter. Based on the results above, the variations in Q seen with the lattice-based method in Fig. 5 can be attributed to the thermal fluctuations of the edge magnetization states. Boundary edges also facilitate annihilation of the isolated skyrmion, which can be seen at T = 300 and 400 K [Figs. 5(e,f) and 5(g,h), respectively], as evidenced by a sharp transition in the time-averaged curves toward Q = 0. Minor oscillations in these time-averaged curves also appear, which result from partially-reversed states at the boundaries that occur during the annihilation process. This example shows that deviations from noninteger (and non half-integer) values of Q can be expected in confined structures when nucleation and annihilation of topological charges take place, in the presence of thermal fluctuations with chiral boundary conditions, or both.

Conclusion
Spurious variations in the topological charge due to inaccuracies in finite-difference derivatives can be mitigated by using a lattice-based approach, for which we provide an implementation for the mumax3 micromagnetics code. While the results do not necessarily call into question the validity of published work (since the topological charge is often used as a proxy for magnetization gradients), they do highlight the care with which noninteger values of Q(t) should be interpreted, particularly when processes such as nucleation, annihilation, and thermal fluctuations are at play.