Journal of Magnetism and Magnetic Materials 549 (2022) 168972 Contents lists available at ScienceDirect Journal of Magnetism and Magnetic Materials journal homepage: www.elsevier.com/locate/jmmm Research articles Exploiting symmetries in skyrmionic micromagnetic simulations: Cylindrical and radial meshes Josep Castell-Queralt, Leonardo González-Gómez, Nuria Del-Valle, Carles Navau ∗ Departament de Física, Universitat Autònoma de Barcelona, Edifici C (campus UAB), 08193 Bellaterra, Barcelona, Catalonia, Spain ARTICLE INFO ABSTRACT Keywords: Micromagnetic simulations are an essential tool in the theoretical study of magnetic skyrmions. When confined Micromagnetism into nanometric samples, these structures can serve as bits of information among other possible applications. Skyrmions Accurate simulations are one of the major sources of theoretical results. In the case of confined skyrmions, it Simulation is known that the boundaries play a critical role in their stabilization. However, most of the micromagnetic simulations are done using a finite-difference method with quadrilateral meshes, that do not exactly fit the boundaries. The use of this mesh can introduce a significant numerical error that can completely change the results of the simulations. We present here two different finite-difference meshes to study skyrmions in confined disks, taking advantage of the symmetry of that geometry. A two-dimensional cylindrical mesh for non-symmetric scenarios but geometrically symmetric (boundary conditions) that reduces the propagation of the numerical error, and the particular case of a one dimensional mesh for axisymmetric scenarios where the computation time is hugely reduced. 1. Introduction Within the micromagnetic model, the ferromagnet is described by a continuous magnetization field, 𝐌(𝐫, 𝑡) (𝐫 is the position and 𝑡 the Ferromagnetic materials are at the core of a large and diverse time). This magnetization varies in time due to an effective magnetic number of devices. Besides their description at macroscopic scales, at field, 𝐇ef f (𝐫, 𝑡), which is not only produced by external sources but the submicrometer scale, the magnetization dynamics of these materials also by internal interactions and, thus, depends on the magnetiza- is usually modeled using micromagnetics. This model accounts for tion function itself. The magnetization dynamics is described by the the dynamics of the local magnetization function, that represents the Landau–Lifshitz–Gilbert (LLG) equation [20], average of atomic interacting magnetic moments [1–4]. d𝐌 𝛼𝛾 The micromagnetic model has been (and is currently being) used (1 + 𝛼 2 ) = −𝛾𝐌 × 𝐇ef f − 𝐌 × 𝐌 × 𝐇ef f , (1) to simulate systems ranging from nanomagnetic data storage [5–7], d𝑡 𝑀𝑠 spintronic logic devices [8,9], or sensors [10,11]. Devices based on where 𝛼 is the Gilbert damping parameter, 𝛾 is the gyromagnetic ratio, skyrmions and other skyrmionic structures are not an exception. Many and |𝐌| = 𝑀𝑠 is the saturation magnetization, which is assumed to be important advances in their study have been done after micromag- independent of position and time. netic calculations (for example, [12–15]). In the case of skyrmions We shall show that in a circular disk and starting from a cylindri- in thin films, they are usually stabilized thanks to the interfacial cally symmetric equilibrium magnetization distribution, like a skyrmion, Dzyaloshinskii–Moriya interaction (iDMI) appearing at the contact sur- the time dependent 𝐌(𝐫, 𝑡) evolving according to the LLG equation face with a heavy-metal substrate [16]. iDMI is added to the exchange, should preserve this symmetry (if the external field has also that same the anisotropy, the dipolar interaction and, eventually, to the interac- symmetry). However, when solving numerically the LLG equation with tion with external fields (or, in general, to other external agents). More- quadrilateral meshes (QM) but having rounded boundaries results in over, if the skyrmion is confined, the presence of the iDMI, imposes extra restrictions to the equilibrium magnetization distribution [17]. inevitable numerical errors. These errors, usually larger at the bound- The use of skyrmions in confined geometries has been proposed as aries, can affect the whole sample. Refining the meshes to retain the elements for magnetic memories in bit patterned recording media or symmetry for short enough calculation times is a valid strategy but at for artificial neural network computation [18,19], among others. the price of increasing the calculation time. ∗ Corresponding author. E-mail address: carles.navau@uab.cat (C. Navau). https://doi.org/10.1016/j.jmmm.2021.168972 Received 6 October 2021; Received in revised form 17 December 2021; Accepted 21 December 2021 Available online 7 January 2022 0304-8853/© 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). J. Castell-Queralt et al. Journal of Magnetism and Magnetic Materials 549 (2022) 168972 √ Here we propose the use of cylindrical meshes (CM) in skyrmionic dimensions with respect to the exchange length, 𝑙ex = 2𝐴∕𝜇0 𝑀𝑠2 , and micromagnetic simulations in two-dimensional disks in order to reduce the time with respect to (𝛾𝑀𝑠 )−1 . Thus, we define the numerical errors triggered by the usual use of rectangular dis- 𝐇𝛽 cretizations, without compromising (even increasing) the performance 𝐡𝛽 = , (3) 𝑀𝑠 of the calculation in terms of computation time and accuracy. In the 𝐌 case of skyrmions confined in disks, the boundary conditions are deter- 𝐦 = , (4) 𝑀𝑠 mined by the iDMI, which is also the interaction responsible of their 𝑡 stabilization [17,18]. In that case, larger numerical errors at borders 𝜏 = , (5) (𝛾𝑀𝑠 ) could be critical for the validity of the simulation. 2𝐾 We shall proceed as follows. In Section 2, we present the general 𝜅 = , (6) 𝜇0 𝑀𝑠2 geometry, boundaries and effective fields used in the calculations and 𝑙 𝐷 in the problem definition. After some general comments on the error 𝜉 = ex , (7) 𝐴 propagation due to quadrilateral meshes with cylindrical boundary where 𝛽 stands for ext, 𝐴, 𝐷, 𝐾, or ef f . The normalized effective field conditions (Section 3), we introduce (Section 4) the two-dimensional is rewritten as (note that the ∇ operator is now also a normalized CM for the cases where the boundaries are axisymmetric but 𝐌(𝐫, 𝑡) is operator, which means that the derivatives are done with respect to not. The use of this mesh reduces the numerical error and, more im- normalized magnitudes) portantly, avoids the appearance of erroneously computed metastable [ ] magnetic structures due to the 𝜋∕4-symmetry of a QM in a circular 𝐡ef f = ∇2 𝐦 + 𝜅𝑚𝑧 𝐳̂ + 𝜉 (∇ ⋅ 𝐦)̂𝐳 − ∇𝑚𝑧 + 𝐡ext , (8) boundary, as we will show in Section 5. Next, in Section 6, we adapt the and the normalized LLG equation becomes study to the case of rings with a central hole, where the roundness of d𝐦 the boundaries have more influence. As a particularly interesting case, (1 + 𝛼 2 ) = −𝐦 × 𝐡ef f − 𝛼𝐦 × 𝐦 × 𝐡ef f . (9) d𝜏 we shall also introduce, in Section 7 a variation of the CM to a one- When considering a circular dot geometry (with radius 𝑅) and iDMI, dimensional mesh with cylindrical symmetry for the cases where both, we solve Eq. (9) with boundary condition [17] 𝐌 and the boundaries are axisymmetric. We finish with the conclusions (Section 8). 𝜕𝐦 || 𝜉 = (̂𝐳 × 𝝆) ̂ × 𝐦, (10) 𝜕𝜌′ ||𝑅′ 2 2. Geometry, effective fields, and boundaries where 𝜌′ = 𝜌∕𝑙ex and 𝑅′ = 𝑅∕𝑙ex . For all the examples in this work we have used the following parameters: 𝑙𝑒𝑥 = 8.42 nm, (𝛾𝑀𝑠 )−1 = 7.84 ps, 𝜉 = 1.85, 𝜅 = 0.3, and 𝛼 = 0.3, which are in the range of the common We consider a planar ferromagnetic disk of radius 𝑅 and thickness micromagnetic parameters for skyrmionic systems [12,22]. From now, 𝑡𝐹 𝑀 . We assume that 𝑡𝐹 𝑀 ≪ 𝑅 so that the magnetization and any and for the rest of the paper, although we will not explicitly write it, other magnetic quantity are uniform across the thickness. The problem all the magnitudes are considered as normalized magnitudes, including is thus reduced to a two-dimensional problem. The disk is located on coordinates. To simplify notation, we will skip the primes in 𝜌′ and 𝑅′ . the 𝑥𝑦-plane with the center at the origin of coordinates. To fix the notation, we will use both standard Cartesian coordinates (𝑥, 𝑦, 𝑧) and 3. General considerations on the symmetry of the solutions cylindrical coordinates (𝜌, 𝜑, 𝑧), with the components of the vectors in these coordinate systems indicated by the corresponding subscript. Given Eqs. (8), (9), and (10) in a circular dot, one sees that if We consider here that the effective fields acting on the ferromagnet both the initial magnetization distribution and the external field are are the exchange field 𝐇𝐴 , the anisotropy field 𝐇𝐾 , the iDMI field 𝐇𝐷 , axisymmetric, 𝐦(𝐫, 𝜏 = 0) = 𝐦(𝜌, 𝜏 = 0), and 𝐡ext (𝐫, 𝜏 = 0) = 𝐡ext (𝜌, 𝜏 = and an externally applied field 𝐇ext . Physically, we consider that the 0), then the initial effective field also shares the same symmetry: ferromagnet is on top of a heavy metal that produces iDMI which not 𝐡ef f (𝐫, 𝜏 = 0) = 𝐡ef f (𝜌, 𝜏 = 0). If both initial 𝐦 and 𝐡ef f are axisymmetric, only affects the effective field but also the boundary conditions. We also it follows that 𝐦 × 𝐡ef f and 𝐦 × (𝐦 × 𝐡ef f ) are also axisymmetric. Hence, consider uniaxial anisotropy with the easy axis perpendicular to the according to Eq. (9), the initial variation over time of the magnetization should also be axisymmetric: 𝑑𝐦 (𝐫, 𝜏 = 0) = 𝑑𝐦 (𝜌, 𝜏 = 0). Finally, if both disk plane (in the 𝑧-direction). Since we consider ultra-thin disks, the 𝑑𝜏 𝑑𝜏 the magnetization and its time derivative have radial symmetry, then demagnetization field is included into the anisotropy field by renormal- this symmetry should be conserved over time: 𝐦(𝐫, 𝜏) = 𝐦(𝜌, 𝜏). izing the anisotropy constant [21]. This assumption is valid as long as In summary, in a ferromagnetic ultra-thin dot, if we initially have the sample is an ultrathin film, since in this case the 𝑧 component of the an axisymmetric magnetic structure and we apply an axisymmetric demagnetizing field is approximately proportional to the 𝑧 component external field, the structure can evolve in time, but it must conserve of the local magnetization, justifying the renormalization. Although the the symmetry all the time. numerical treatment we present would be also valid for other sets of However, when solving the LLG using finite differences in a QM we effective fields that could appear considering bulk materials or non- shall show that, even though initially both the initial magnetization and uniaxial anisotropies, we focus on this particular case because it is the the applied field are axially symmetric, one could obtain structures that usual way for studying skyrmionic structures in thin films. Extensions do not conserve the axial symmetry. As an example, consider a ferro- considering other types of effective fields would be straightforward. magnetic dot of radius 𝑅, with an initial magnetization distribution as Indeed, the total effective field is expressed as follows: ( ) 𝜉𝜌 𝐇ef f = 𝐇𝐴 + 𝐇𝐾 + 𝐇𝐷 + 𝐇ext = 𝑚𝜌 = sin , (11) 2 2𝐴 2𝐾 2𝐷 [ ] = ∇2 𝐌 + 𝑀𝑧 𝐳̂ + (∇ ⋅ 𝐌)̂𝐳 − ∇𝑀𝑧 + 𝐇ext , (2) 𝜇0 𝑀𝑠2 𝜇0 𝑀𝑠2 𝜇0 𝑀𝑠2 𝑚𝜑 = 0, (12) ( ) where the terms in the second line are expressed in the same order 𝜉𝜌 𝑚𝑧 = cos . (13) as in the first one. 𝐴 is the exchange constant, 𝐾 the effective uni- 2 axial anisotropy constant, 𝐷 the iDM constant, and 𝜇0 the vacuum This initial magnetization satisfies the boundary conditions, Eq. (10). permeability. We can compute the effective fields (𝐡ext = 0), obtaining In order to work with non-dimensional magnitudes, we normalize ( 2 ) ( ) ( ) 𝜉 1 𝜉𝜌 𝜉 𝜉𝜌 the magnetization and all the fields with respect to 𝑀𝑠 , the length ℎef f ,𝜌 = − sin + cos , (14) 4 𝜌2 2 2 2 2 J. Castell-Queralt et al. Journal of Magnetism and Magnetic Materials 549 (2022) 168972 Fig. 1. Difference between the effective field components multiplied by 103 [(a): radial, (b): angular, and (c): axial] when calculated numerically using a QM (square- shaped cells of side 𝑅∕100) and when evaluated with the analytical expressions of Eqs. (14)–(16). Fig. 2. Sketch of the proposed cylindrical mesh (right) and comparison with the quadrilateral mesh. The main advantage is that the boundary conditions can be much more accurately established in the cylindrical mesh for rounded geometries. However, ℎef f ,𝜑 = 0, (15) ( ) ( ) the magnitudes do not need to be axisymmetric since angular finite differences are 𝜉 𝜉𝜌 𝜉2 𝜉𝜌 allowed within numerical precision. The small green dots are the grid points of our ℎef f ,𝑧 = sin + cos . (16) system. 2𝜌 2 4 2 In Fig. 1 we show the difference between the numerically calculated effective field using a QM and the effective field analytically obtained from Eqs. (14)–(16), considering the same magnetization distributions 1. Set the number of radial crowns, 𝑁 + 1. 𝑅 given in Eqs. (11)–(13). It is clear that there is a slight difference 2. Define 𝛥𝜌 = 𝑁+1∕2 . This will determine the radial step. The between both but the interesting part is that the symmetry of that radial crowns where we evaluate the functions are indicated by numerical difference is not cylindrical but it has a 𝜋∕4-rotation symme- subscript 𝑖 = 0, 1, … , 𝑁. The radial positions are 𝜌𝑖 = (𝑖 + 1∕2)𝛥𝜌. try. That slight difference will produce, following the numerical time At the extremes, 𝜌0 = 12 𝛥𝜌, and 𝜌𝑁 = 𝑅. integration of Eq. (9), a (numerically) non-symmetric magnetization 3. The angular step, 𝛥𝜑𝑖 , depends on 𝑖, which means that the num- distribution. At each time step this numerical error propagates with the ber of points that form each radial crown is different. Set 𝛥𝜑0 = possible appearance of magnetic structures with a rotation symmetry of 𝜋 , which is the maximum angle that satisfies the inequality in 𝜋∕4 or 𝜋∕2, the same symmetries than the error distribution and the QM 2 the next point have. 4. Compute 𝛥𝜑𝑖 , for 0 < 𝑖 ≤ 𝑁, according to The numerical error produced by the QM is higher at the bound- aries, due to the fact that the boundary cannot be exactly fitted. This • If 𝛿𝛥𝜌 > 𝜌𝑖 𝛥𝜑𝑖−1 , then 𝛥𝜑𝑖 = 𝛥𝜑𝑖−1 . error is specially relevant when iDM interaction is present since it • If 𝛿𝛥𝜌 ≤ 𝜌𝑖 𝛥𝜑𝑖−1 , then 𝛥𝜑𝑖 = 12 𝛥𝜑𝑖−1 . affects not only to the effective field but also directly to the boundary conditions. In materials with larger 𝜉 (large iDM constant) we have Where we set 𝛿 = 32 , so that the error in the radial derivative is steeper space derivatives at the boundaries [Eq. (10)] and thus, since of the same order as the error in the angular derivative. The 𝑖th the numerical error using finite differences is proportional to the spatial radial crown will have 𝐿𝑖 angular points uniformly distributed derivatives, larger numerical errors with the QM mesh. along the crown: 𝐿𝑖 = 2𝜋∕𝛥𝜑𝑖 . 4. The cylindrical mesh, CM Consider an arbitrary function evaluated at a given point, 𝑓 (𝜌, 𝜑). Its value at the considered mesh points (𝜌𝑖 , 𝜑𝑗 ) = ((𝑖 + 1∕2)𝛥𝜌, 𝑗𝛥𝜑𝑖 ), In order to precisely simulate skyrmionic structures confined in dots, for 𝑖 = 0, 1, … , 𝑁 and 𝑗 = 0, 1, … , 𝐿𝑖 − 1, would be 𝑓𝑖,𝑗 . Notice that we propose a new type of mesh which, without being fully axially when 𝑖 = 0, 𝜌 ≠ 0. This is done to avoid the cylindrical coordinates symmetric, takes advantage of the cylindrical symmetry. indeterminacy at 𝜌 = 0. To evaluate the effective fields [Eqs. (17)–(20)] It is important to note that, while the introduced CM is axisymmet- using the presented mesh, one also needs the discrete derivatives of ric, neither the applied field nor the magnetization distribution need to the (components of) magnetization function. In Fig. 3 one can see the be cylindrically symmetric. stencils used for evaluating the radial derivatives. The detailed expres- The effective fields, expressed in standard cylindrical coordinates sions for these discrete radial and angular derivatives are presented in (𝜌, 𝜑, 𝑧), can be expressed as: Appendix A. ( ) 𝜕 2 𝑚𝜌 1 𝜕𝑚𝜌 𝑚𝜌 2 1 𝜕 𝑚𝜌 2 𝜕𝑚𝜑 𝐡𝐴 = + − + − 𝝆̂ + 𝜕𝜌2 𝜌 𝜕𝜌 𝜌2 𝜌2 𝜕𝜑2 𝜌2 𝜕𝜑 5. Comparative results between quadrilateral and cylindrical ( ) meshes 𝜕 2 𝑚𝜑 1 𝜕𝑚𝜑 𝑚𝜑 2 1 𝜕 𝑚𝜑 2 𝜕𝑚𝜌 + + − + − ̂ + 𝝋 𝜕𝜌2 𝜌 𝜕𝜌 𝜌2 𝜌2 𝜕𝜑2 𝜌2 𝜕𝜑 We have performed the same micromagnetic simulation using two ( 2 2 ) 𝜕 𝑚𝑧 1 𝜕𝑚𝑧 1 𝜕 𝑚𝑧 different meshes: one with the QM and the other with the presented + + + 𝐳̂ , (17) 𝜕𝜌 2 𝜌 𝜕𝜌 𝜌 𝜕𝜑2 2 CM. The physical process followed is the following: we set 𝐦(𝜏 = 0) = 𝐳̂ , [ ( ) ] 𝜕𝑚 1 𝜕𝑚𝜑 1 𝜕𝑚𝜌 𝜕𝑚𝜑 and ℎext (𝜏 = 0) = 0.4. Then we let the system evolve as the external field 𝐡𝐷 = 𝜉 − 𝑧 𝝆̂ + 𝝋̂ + 𝑚𝜌 + 𝜌 + 𝐳̂ , (18) varies 𝜕𝜌 𝜌 𝜕𝜑 𝜌 𝜕𝜌 𝜕𝜑 ( from)0.4 to −0.4 following a lineal decay with time, ℎext (𝜏) = 0.4 1 − 2 𝑇𝜏 , setting 𝑇 =6400, the total simulation time. We choose a 𝐡𝐾 = 𝜅𝑚𝑧 𝐳̂ , (19) radius 𝑅 = 7.5 and 𝑁 = 100, so that 𝛥𝜌 = 𝛥 = 𝑅∕100.5, being 𝛥 the side 𝐡ext = ℎext 𝐳̂ . (20) of the squared cell in the QM. The time step is set to 𝛥𝜏 = 5 ⋅ 10−4 . The The proposed mesh is sketched in Fig. 2, compared with the stan- total number of cells are very similar in the QM and in the CM (about dard QM. Note that the CM has uniform radial discretization step 10000 cells, more details in Section 7 and Fig. 7). but adaptative angular discretization step (depending on the radial We show in Fig. 4 some results for the magnetization obtained distance). The mesh is constructed using the following algorithm. Given during a total simulation time equivalent to 50 ns, comparing the two a disk of radius 𝑅: meshes at different times. 3 J. Castell-Queralt et al. Journal of Magnetism and Magnetic Materials 549 (2022) 168972 Fig. 3. Sketch for the finite difference of the radial derivatives of a function using the cylindrical mesh. Blue points are where the finite difference is being evaluated, red points are the ones used in the stencil and purple points are the ones where we apply the boundary conditions. In (a) All the possible stencils are sketched, and the relation between the 𝑗 index of different radial crowns is shown. In (b), (c), and (d) we show the different stencils we can find for 0 < 𝑖 ≤ 𝑁. In (e) we show the points used to apply the boundary conditions. As seen, using the QM we can find magnetic structures with 𝜋∕2 or 𝜋∕4 symmetry, while with the CM the axial symmetry is preserved, as should be. It is worth to note that in the QM the 𝜋∕2 or 𝜋∕4 symmetry appears from the beginning of the calculation. As time goes by, this numerical error becomes more prominent and propagates to the whole sample, triggering the nucleation of a magnetic structure of 𝜋∕4 symmetry. Our conclusion is that the cylindrical mesh is more convenient since it offers coherent physical results. Moreover, the computational complexity of the algorithm is similar for both methods (there is approximately the same amount of cells in both methods for the same numerical precision). Although with a finer mesh the numerical error with the QM could be reduced, the computation time needed to have similar numerical error than with the CM would be excessive for a typical simulation. When reducing the number of cells, reducing the dot size, or/and when there is a relatively large iDMI (affecting considerably the bound- ary conditions), the use of the presented cylindrical mesh becomes more efficient in terms of accuracy/time. 6. Cylindrical mesh for ferromagnetic rings Another important advantage of the CM is that it allows to simulate dots with small circular holes, of radius 𝑅0 , at their center. If one uses a QM, a cell size much smaller than 𝑅0 is required, yielding prohibitive computational time. With the CM one can simulate these ferromagnetic rings without increasing the computational time since the inner boundary is perfectly fitted regardless the cell size. Since the 𝜌 = 0 point is not in the sample, now the algorithm to generate the mesh, Fig. 5, can be optimized varying slightly the general algorithm presented for non-holed dots. Given a ring with internal and Fig. 4. 𝑧-component of the magnetization distribution for an ultra-thin ferromagnetic external radii 𝑅0 and 𝑅, respectively: disk calculated using the cylindrical mesh (I: top-half) and the quadrilateral mesh (II: bottom-half). In each half, each plot corresponds to the different times (a): 0 ns (b): 1. Fix the number of radial crowns, 𝑁 + 1. 𝑅−𝑅 18.75 ns (c): 25 ns (d): 37.5 ns (e): 46.25 ns, and (f): 50 ns. The central bar is the 2. Define 𝛥𝜌 = 𝑁 0 . This will determine the radial step. The radial colorbar for the 𝑚𝑧 values. crowns are indicated with subscript 𝑖. 𝑖 = 0, 1, … , 𝑁. The radial position of the crowns where the magnitudes are evaluated are 𝜌𝑖 = 𝑅0 + 𝑖𝛥𝜌. We have 𝜌0 = 𝑅0 and 𝜌𝑁 = 𝑅. 4 J. Castell-Queralt et al. Journal of Magnetism and Magnetic Materials 549 (2022) 168972 Fig. 5. Sketch of the proposed cylindrical mesh (right) and comparison with the quadrilateral mesh (left), in the case of a ring. 3. The angular step, 𝛥𝜑𝑖 will depend on 𝑖, which means that the number of cells that form each crown will be different. Set 𝛥𝜑0 = 22𝜋𝑆 , being 𝑆 is the minimum integer that would satisfy the inequality 32 𝛥𝜌 ≥ 𝜌0 𝛥𝜑0 . 4. Compute 𝛥𝜑𝑖 for 𝑖 > 0 using the following relation also with 𝛿 = 32 • If 𝛿𝛥𝜌 ≥ 𝜌𝑖 𝛥𝜑𝑖−1 , then 𝛥𝜑𝑖 = 𝛥𝜑𝑖−1 . • If 𝛿𝛥𝜌 < 𝜌𝑖 𝛥𝜑𝑖−1 , then 𝛥𝜑𝑖 = 12 𝛥𝜑𝑖−1 . The 𝑖th radial crown has 𝐿𝑖 angular points: 𝐿𝑖 = 2𝜋∕𝛥𝜑𝑖 The discrete version of the radial and angular derivatives in this case is detailed in Appendix B. We repeat the simulation of the previous section but using a ferro- magnetic ring with 𝑅 = 6 and 𝑅0 = 1. The physical process simulated Fig. 6. 𝑧-component of the magnetization distribution for a magnetic disk calculated and the numerical details are exactly the same as in Section 5. The using the cylindrical mesh (I: top-half) and the quadrilateral mesh (II: bottom-half). In results of the simulations are shown in Fig. 6. Again the differences each half, each plot correspond to the different times (a): 0 ns (b): 21.25 ns (c): 25 ns using the QM and the CM are critical. In particular, we observe that (d): 26.25 ns (e): 28.75 ns, and (f): 50 ns . The central bar is the colorbar for the 𝑚𝑧 the presence of the hole affects substantially the magnetic structure. In values. this case the numerical errors in the QM appearing at the central hole are even larger than those at the external border. cells that each method (QM, CM and RM) needs to simulate a given 7. A particular case: radial mesh ferromagnetic dot. To compare, we plot the total number of cells 𝑁𝑡 as a function of the number of radial divisions 𝑁𝜌 , using 𝛥 = 𝛥𝜌 for From the previous simulations one may wonder why we have the QM. For a fixed 𝛥𝜌 the horizontal axis is proportional to the radius computed the angular derivatives if we are sure that the axial symmetry of the dot. We observe that the amount of cells is always comparable must be conserved. Actually, if the initial magnetization distribution between QM and CM while in the RM the cells used are much less, and the applied field are both cylindrically symmetric, angular deriva- specially as the radius of the dot increases (with 𝛥𝜌 fixed). Actually, tives are zero. In this case, the system is actually one dimensional (only for the RM 𝑁𝑡 ∼ 𝑅 while for the QM and the CM 𝑁𝑡 ∼ 𝑅2 . As 𝜌 dependence) and we could set 𝐦 = 0 in the cylindrical mesh, resulting for the computation times, the number of operations needed for each in only radial divisions. We call it radial mesh (RM). cylindrical is, in average, only slightly larger to those needed for a 2 In this case 𝜕𝐦 𝜕𝜑 = 0, 𝜕𝜕 2𝐦 𝜑 = 0, and the effective fields are quadrilateral cell. Actually, in all the simulations the running times for [( ) ( ) ] a simulation using the CM is less that 5% larger than that for the QM, 𝜕 2 𝑚𝜌 𝑚𝜌 𝜕 2 𝑚𝜑 𝑚𝜑 𝜕 2 𝑚𝑧 with similar number of cells. Note that the RM and the CM do not allow 𝐡𝐴 = − 𝝆̂ + − ̂ + 𝝋 𝐳̂ , (21) for Fast Fourier Transform algorithms (which are used sometimes to 𝜕𝜌2 𝜌2 𝜕𝜌2 𝜌2 𝜕𝜌2 [ ] compute demagnetizing fields) since they are not regular meshes. 1 𝜕(𝜌𝑚𝜌 ) 𝜕𝑚𝑧 𝐡𝐷 = 𝜉 𝐳̂ − 𝝆̂ , (22) 𝜌 𝜕𝜌 𝜕𝜌 8. Conclusions 𝐡𝐾 = 𝜅 𝑚𝑧 𝐳̂ , (23) We have presented some studies on magnetic skyrmionic structures, The involved functions only depend on 𝜌. There is no angular when confined in circular dots or rings. The specially designed CM takes dependence, nor angular variations. Thus, the evaluation points would advantage of the axisymmetric boundaries, even though the whole be defined by a single radial index. We use finite differences as we did system may not necessarily be axisymmetric. We have seen, by imple- for the CM, the exact expressions are presented in Appendix C. menting the CM in a typical skyrmionic problem in thin disks, that the For the boundaries, 𝑖 = 0, 𝑁; we use the boundary condition, use of standard QM can yield numerical error propagation that could | Eq. (10), except for the dot at 𝑖 = 0 where we use 𝜕𝐦 | 𝜕𝜌 |𝜌=0 = 0, which invalidate the obtained results. The presented mesh drastically reduces must be satisfied if 𝐦 is axisymmetric. this error, without compromising the computational performance. Naturally, with this mesh we recover the results obtained in the A real sample would never be perfectly circular, and one would simulation of the previous section done with the CM, with a substan- expect that the real magnetization distribution was not perfectly sym- tially less computation time. In Fig. 7, we compare the amount of metric. Nevertheless, in simulations, if the symmetry is broken not by 5 J. Castell-Queralt et al. Journal of Magnetism and Magnetic Materials 549 (2022) 168972 When 𝑖 = 0, the points 𝑓𝑖−1,𝑗 do not exist, so we use forward finite differences, (Fig. 3-d), 𝜕𝑓0,𝑗 𝑓2,𝑗 − 𝑓0,𝑗 ≃ , (31) 𝜕𝜌 2𝛥𝜌 𝜕 2 𝑓0,𝑗 𝑓2,𝑗 − 2𝑓1,𝑗 + 𝑓0,𝑗 ≃ . (32) 𝜕𝜌2 (𝛥𝜌)2 When 𝑖 = 𝑁, the boundary condition of the system, Eq. (10), represented by the purple stencil in (Fig. 3-e), will be used to directly compute the different components of the magnetization, according to: 𝑚𝜌,𝑁,𝑗 − 𝑚𝜌,𝑁−1,𝑗 𝜉 = 𝑚𝑧,𝑁,𝑗 , (33) 𝛥𝜌 2 𝑚𝜑,𝑁,𝑗 − 𝑚𝜑,𝑁−1,𝑗 = 0, (34) 𝛥𝜌 Fig. 7. Number of computational cells as a function of the number of radial divisions 𝑚𝑧,𝑁,𝑗 − 𝑚𝑧,𝑁−1,𝑗 𝜉 for the different meshes used in this work:QM (blue), CM (red), and RM (violet). = − 𝑚𝜌,𝑁,𝑗 . (35) 𝛥𝜌 2 The angular derivatives are evaluated taking into account that 𝑗 is controlled notches or defects but as the result of numerical error, this a periodic index, which means that 𝑓𝑖,𝐿𝑖 = 𝑓𝑖,0 and 𝑓𝑖,−1 = 𝑓𝑖,𝐿𝑖 −1 . Then, leads to misleading results. We have shown here that these errors can 𝜕𝑓𝑖,𝑗 𝑓𝑖,𝑗+1 − 𝑓𝑖,𝑗−1 be greatly avoided. ≃ , (36) 𝜕𝜑 2𝜌𝑖 𝛥𝜑𝑖 𝜕 2 𝑓𝑖,𝑗 𝑓𝑖,𝑗+1 − 2𝑓𝑖,𝑗 + 𝑓𝑖,𝑗−1 Declaration of competing interest ≃ . (37) 𝜕𝜑2 (𝜌𝑖 𝛥𝜑𝑖 )2 The authors declare that they have no known competing finan- cial interests or personal relationships that could have appeared to Appendix B. Discrete derivatives for the CM mesh in rings influence the work reported in this paper. In the case of a ring with the CM mesh, the numerical radial Acknowledgments and angular derivatives for 0 < 𝑖 < 𝑁 are the same Eqs. (31)–(32) (Appendix A). At the boundaries (𝑖 = 0 and 𝑖 = 𝑁) we use the boundary We acknowledge Catalan project, Spain 2017-SGR-105 and Spanish conditions Eq. (10). For 𝑖 = 𝑁 we use Eq. (33)–(35) (Appendix A), for project PID2019-104670GB-I00 of Agencia Estatal de Investigación/Fondo 𝑖 = 0 the magnetization is given by Europeo de Desarrollo Regional (UE) for financial support. J. C. Q. 𝑚𝜌,1,𝑗 − 𝑚𝜌,0,𝑗 𝜉 = − 𝑚𝑧,0,𝑗 , (38) acknowledges a grant (FPU17/01970) from Ministerio de Ciencia, 𝛥𝜌 2 Innovación y Universidades (Spanish Government). 𝑚𝜑,1,𝑗 − 𝑚𝜑,0,𝑗 = 0, (39) 𝛥𝜌 Appendix A. Discrete derivatives for the CM mesh 𝑚𝑧,1,𝑗 − 𝑚𝑧,0,𝑗 𝜉 = 𝑚𝜌,0,𝑗 . (40) 𝛥𝜌 2 In the case of a circular dot, using the CM mesh, the radial deriva- tives, when 0 < 𝑖 < 𝑁, are evaluated as (see the sketches in Fig. 3), Appendix C. Discrete derivatives for the RM + − 𝜕𝑓𝑖,𝑗 𝑓𝑖,𝑗 − 𝑓𝑖,𝑗 ≃ , (24) Using the same notation as in the CM, the expression used for the 𝜕𝜌 2𝛥𝜌 first and second radial derivatives are, + − 𝜕 2 𝑓𝑖,𝑗 𝑓𝑖,𝑗 − 2𝑓𝑖,𝑗 + 𝑓𝑖,𝑗 𝜕𝑓𝑖 𝑓𝑖+1 − 𝑓𝑖−1 ≃ , (25) ≃ (41) 𝜕𝜌2 (𝛥𝜌)2 𝜕𝜌 2𝛥𝜌 + − depends on the number of angular where the value of 𝑓𝑖,𝑗 and 𝑓𝑖,𝑗 𝜕 2 𝑓𝑖 𝑓𝑖+1 − 2𝑓𝑖 + 𝑓𝑖−1 divisions of the crowns 𝑖 + 1 and 𝑖 − 1, respectively, according to (in ≃ (42) 𝜕𝜌2 (𝛥𝜌)2 the following 𝑗∕2 is the integer division of 𝑗): References • If 𝛥𝜑𝑖+1 = 𝛥𝜑𝑖 , then + [1] H. Kronmuller, M. Fahnle, Micromagnetism and the Microstructure of 𝑓𝑖,𝑗 = 𝑓𝑖+1,𝑗 , (26) Ferromagnetic Solids, Cambridge University Press, 2003. • If 𝛥𝜑𝑖 = 𝛥𝜑𝑖−1 , then [2] A. Aharoni, Introduction to the Theory of Ferromagnetism, Oxford University Press, 2000. − [3] J. Fidler, T. Schrefl, Micromagnetic modelling - the current state of the art, 𝑓𝑖,𝑗 = 𝑓𝑖−1,𝑗 . (27) 2000, http://dx.doi.org/10.1088/0022-3727/33/15/201, https://iopscience.iop. org/article/10.1088/0022-3727/33/15/201/meta. • If 2𝛥𝜑𝑖+1 = 𝛥𝜑𝑖 , then [4] J. Leliaert, J. Mulkers, Tomorrow’s micromagnetic simulations, J. Appl. Phys. 125 + (18) (2019) 180901, http://dx.doi.org/10.1063/1.5093730, http://aip.scitation. 𝑓𝑖,𝑗 = 𝑓𝑖+1,2𝑗 . (28) org/doi/10.1063/1.5093730. [5] A. Wachowiak, J. Wiebe, M. Bode, O. Pietzsch, M. Morgenstern, R. Wiesendanger, • If 2𝛥𝜑𝑖 = 𝛥𝜑𝑖−1 and 𝑗 is even, (Fig. 3-b), then Direct observation of internal spin structure of magnetic vortex cores, Science − 298 (5593) (2002) 577–580, http://dx.doi.org/10.1126/science.1075302, https: 𝑓𝑖,𝑗 = 𝑓𝑖−1,𝑗∕2 . (29) //science.sciencemag.org/content/298/5593/577.abstract. [6] S. Mangin, D. Ravelosona, J.A. Katine, M.J. Carey, B.D. Terris, E.E. Fuller- • If 2𝛥𝜑𝑖 = 𝛥𝜑𝑖−1 and 𝑗 is odd, (Fig. 3-c), then ton, Current-induced magnetization reversal in nanopillars with perpendicular anisotropy, Nature Mater. 5 (3) (2006) 210–215, http://dx.doi.org/10.1038/ − 1 𝑓𝑖,𝑗 = (𝑓 + 𝑓𝑖−1,𝑗∕2+1 ). (30) nmat1595, https://www.nature.com/articles/nmat1595. 2 𝑖−1,𝑗∕2 6 J. Castell-Queralt et al. Journal of Magnetism and Magnetic Materials 549 (2022) 168972 [7] B. Van Waeyenberge, A. Puzic, H. Stoll, K.W. Chou, T. Tyliszczak, R. Hertel, [15] J. Sampaio, V. Cros, S. Rohart, A. Thiaville, A. Fert, Nucleation, stability and M. Fähnle, H. Brückl, K. Rott, G. Reiss, I. Neudecker, D. Weiss, C.H. Back, current-induced motion of isolated magnetic skyrmions in nanostructures, Nature G. Schütz, Magnetic vortex core reversal by excitation with short bursts of an Nanotechnol. 8 (3) (2013) 839–844, http://dx.doi.org/10.1038/nnano.2013.210. alternating field, Nature 444 (7118) (2006) 461–464, http://dx.doi.org/10.1038/ [16] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, nature05240, https://www.nature.com/articles/nature05240. G. Bihlmayer, S. Blügel, Spontaneous atomic-scale magnetic skyrmion lattice in [8] D.E. Nikonov, I.A. Young, Benchmarking spintronic logic devices based on two dimensions, Nat. Phys. 7 (9) (2011) 713–718, http://dx.doi.org/10.1038/ magnetoelectric oxides, 2014, http://dx.doi.org/10.1557/jmr.2014.243, https: nphys2045, http://www.nature.com/articles/nphys2045. //link.springer.com/article/10.1557/jmr.2014.243. [17] S. Rohart, A. Thiaville, Skyrmion confinement in ultrathin film nanostructures [9] S. Luo, M. Song, X. Li, Y. Zhang, J. Hong, X. Yang, X. Zou, N. Xu, in the presence of dzyaloshinskii-moriya interaction, Phys. Rev. B 88 (18) L. You, Reconfigurable skyrmion logic gates, Nano Lett. 18 (2) (2018) 1180–1184, http://dx.doi.org/10.1021/acs.nanolett.7b04722, https://pubs.acs. (2013) 184422, http://dx.doi.org/10.1103/PhysRevB.88.184422, http://link.aps. org/doi/abs/10.1021/acs.nanolett.7b04722. org/doi/10.1103/PhysRevB.88.184422. [10] D.E. Heim, J. Tsang, V.S. Speriosu, B.A. Gurney, M.L. Williams, R.E. Fontana, [18] D. Cortés-Ortuño, N. Romming, M. Beg, K. Von Bergmann, A. Kubetzka, O. Design and operation of spin valve sensors, IEEE Trans. Magn. 30 (2) (1994) Hovorka, H. Fangohr, R. Wiesendanger, Nanoscale magnetic skyrmions and 316–321, http://dx.doi.org/10.1109/20.312279. target states in confined geometries, Phys. Rev. B 99 (2019) 214408, http: [11] D. Suess, A. Bachleitner-Hofmann, A. Satz, H. Weitensfelder, C. Vogler, F. //dx.doi.org/10.1103/PhysRevB.99.214408. Bruckner, C. Abert, K. Prügl, J. Zimmer, C. Huber, S. Luber, W. Raberg, [19] X. Zhang, W. Cai, X. Zhang, Z. Wang, Z. Li, Y. Zhang, K. Cao, N. Lei, W. Kang, T. Schrefl, H. Brückl, Topologically protected vortex structures for low-noise Y. Zhang, H. Yu, Y. Zhou, W. Zhao, Skyrmions in magnetic tunnel junctions, magnetic sensors with high linear range, Nat. Electron. 1 (6) (2018) 362–370, ACS Appl. Mater. Interfaces 10 (19) (2018) 16887–16892, http://dx.doi.org/10. http://dx.doi.org/10.1038/s41928-018-0084-2. 1021/ACSAMI.8B03812. [12] R. Tomasello, E. Martinez, R. Zivieri, L. Torres, M. Carpentieri, G. Finoc- [20] T.L. Gilbert, A phenomenological theory of damping in ferromagnetic materials, chio, A strategy for the design of skyrmion racetrack memories, Sci. Rep. IEEE Trans. Magn. 40 (6) (2004) 3443–3449, http://dx.doi.org/10.1109/TMAG. 4 (2014) 6784, http://dx.doi.org/10.1038/srep06784, http://www.nature.com/ 2004.836740. srep/2014/141029/srep06784/full/srep06784.html. [21] K.J. Harte, Theory of magnetization ripple in ferromagnetic films, J. Appl. [13] J. Iwasaki, W. Koshibae, N. Nagaosa, Colossal spin transfer torque effect on Phys. 39 (3) (1968) 1503–1524, http://dx.doi.org/10.1063/1.1656388, http: skyrmion along the edge, Nano Lett. 14 (8) (2014) 4432–4437, http://dx.doi. org/10.1021/nl501379k, http://pubs.acs.org/doi/abs/10.1021/nl501379. //aip.scitation.org/doi/10.1063/1.1656388. [14] A. Thiaville, Y. Nakatani, J. Miltat, Y. Suzuki, Micromagnetic understanding [22] J. Castell-Queralt, L. González-Gómez, N. Del-Valle, A. Sanchez, C. Navau, of current-driven domain wall motion in patterned nanowires, Europhys. Accelerating, guiding, and compressing skyrmions by defect rails, Nanoscale 11 Lett. (EPL) 69 (6) (2005) 990–996, http://dx.doi.org/10.1209/epl/i2004- (26) (2019) 12589–12594, http://dx.doi.org/10.1039/C9NR02171J, http://xlink. 10452-6, http://stacks.iop.org/0295-5075/69/i=6/a=990?key=crossref. rsc.org/?DOI=C9NR02171J. 08cfdeef6135f63254ff9d17d22c1198. 7