Boundary conditions for the Néel order parameter in a chiral antiferromagnetic slab
Oleksandr V. Pylypovskyi, Artem V. Tomilo, Denis D. Sheka, Jürgen Fassbender, Denys Makarov
BBoundary conditions for the N´eel order parameter in a chiral antiferromagnetic slab
Oleksandr V. Pylypovskyi,
1, 2, ∗ Artem V. Tomilo, † Denis D. Sheka, ‡ J¨urgen Fassbender, § and Denys Makarov ¶ Helmholtz-Zentrum Dresden-Rossendorf e.V., Institute of Ion Beam Physics and Materials Research, 01328 Dresden, Germany Kyiv Academic University, Kyiv 03142, Ukraine Taras Shevchenko National University of Kyiv, 01601 Kyiv, Ukraine (Dated: January, 22, 2021)Understanding of the interaction of antiferromagnetic solitons including domain walls andskyrmions with boundaries of chiral antiferromagnetic slabs is important for the design of prospec-tive antiferromagnetic spintronic devices. Here, we derive the transition from spin lattice to mi-cromagnetic nonlinear σ -model with the corresponding boundary conditions for a chiral cubic G-type antiferromagnet and analyze the impact of the slab boundaries and antisymmetric exchange(Dzyaloshinskii–Moriya interaction) on the vector order parameter. We apply this model to evaluatemodifications of antiferromagnetic domain walls and skyrmions upon interaction with boundariesfor different strengths of the antisymmetric exchange. Due to the presence of the antisymmetricexchange, both types of antiferromagnetic solitons become broader when approaching the boundaryand transform to a mixed Bloch–N´eel structure. Both textures feel the boundary at the distanceof about 5 magnetic lengths. In this respect, our model provides design rules for antiferromagneticracetracks, which can support bulk-like properties of solitons. I. INTRODUCTION
The requirement for high storage densities and opera-tion speed of devices stimulates the development of anti-ferromagnetic (AFM) spintronics and spin-orbitronics [1–4]. The envisioned devices rely on AFM textures mov-ing in spatially confined channels [5–10] similarly to fer-romagnetic racetracks [11]. One of the most efficientways to control their dynamics are spin-orbit staggeredtorques, which require specific symmetry of antiferromag-nets rendering them chiral with Dzyaloshinskii–Moriyainteraction (DMI) [12–15]. In this respect, the techno-logical progress in design and optimization of AFM race-tracks requires a fundamental understanding of the in-teraction of magnetic solitons with boundaries of a chiralAFM slab.Similarly to ferromagnets [16, 17], sample boundariesin antiferromagnets usually act as an injector of soli-tons [18] and alter the shape of a spatially confined do-main wall in the media with patterned surfaces [19]. Theinhomogeneous DMI of the surface type leads to the sur-face twist of the order parameter in a two-dimensional(2D) antiferromagnet [20]. If a homogeneous DMI ispresent in addition, an enhanced surface magnetizationaccompanies the deviation of the N´eel vector from thecollinear state [20].To understand properties of the ground state and AFMsolitons in confined geometries, it is necessary to make aproper transition from the Heisenberg spin lattice to themicromagnetic model [21]. The behavior of AFM latticescan be described using two alternative approaches. His- ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] torically, the first one was proposed in the seminal worksof Louis N´eel (for review we refer the reader to [22]).This framework utilizes a representation of antiferro-magnetically coupled ferromagnetic sublattices M and M with | M , | = M s being the saturation magnetiza-tion [23]. The corresponding boundary conditions arederived for each of the sublattices [24–27] also allow-ing to take into account demagnetizing fields [24, 25].Within the second approach, the equations of motionfor M , are rewritten in terms of dimensionless vec-tors of N´eel n = ( M − M ) / (2 M s ) and ferromagnetism m = ( M + M ) / (2 M s ). For many practical cases, therelation | m | (cid:28) | n | ≈ n in one dimension (1D) can be done by splitting the spinlattice into dimers [31–34]. A straightforward procedureshows that the continuum Lagrangian contains the topo-logical term proportional to m · ∂ x n , which determinesdifferences between the quantum integer and half-integerspin chains [35], e.g., the Haldane gap. It originates fromthe choice of dimer pairs: the Hamiltonian of a spin chainis not invariant with respect to the sublattice exchange.The latter also leads to the intrinsic magnetization of 1Dtextures [33]. For the case of 2D AFM lattices, the terms m · ∂ i n are not as important as for spin chains [31]. How-ever, the dimerization in 2D bipartite lattices is ambigu-ous since it can be performed along one of two indepen-dent direction ( i = x or i = y ) and may lead to spuriouseffects due to the choice of spin pairs [36, 37]. This issuecan be overcome by splitting e. g., of the square latticeinto tetramers [38, 39]. While still there are only two or- a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n der parameters, the master N´eel vector (director) n andthe slave ferromagnetism vector m , this procedure re-quires two additional auxiliary fields behaving as spatialderivatives of the N´eel vector. This transition allows tobuild a nonlinear σ -model of a chiral 2D antiferromag-net preserving the spatial invariance within the latticeplane [38, 39].Although 1D and 2D cases are fairly well understood,there is no rigorous transition from a spin lattice descrip-tion to a three-dimensional (3D) micromagnetic model.This leads to a gap in our understanding of the impactof boundary conditions on the ground state and antifer-romagnetic solitons in spatially confined chiral AFMs.The boundary of a sample is an additional source of thesymmetry break in AFM lattices. Thus, a unit antifer-romagnetic cell should be properly chosen to correctlydetermine the boundary conditions for the order param-eters.Here, we rigorously derive a transition from the classi-cal spin-lattice Heisenberg Hamiltonian to the nonlinear σ -model for a 3D chiral AFM slab with a simple cubiclattice by means of the octamerization process. We avoidspurious parity-breaking effects and show the influence ofboundary conditions and DMI on the ground state in thelaterally confined sample. We apply this model for twotypes of magnetic solitons, namely, a translational AFMdomain wall and skyrmion, and demonstrate their mod-ification due to the sample boundaries. These texturespossess twists in the order parameter as well as the defor-mation of their shape near the surface, which is analyzedquantitatively.This paper is structured as follows. The spin-latticeHamiltonian and derivation of the corresponding non-linear σ -model with boundary conditions is described inSec. II A. In Sec. II B, this micromagnetic model is ap-plied for the description of the ground state of a chi-ral uniaxial AFM slab. Properties of the domain walland skyrmion in a confined geometry are described inSec. II C and II D, respectively. The main results of thepaper are summarized in III. Further details describingthe transition from the spin lattice to the continuum de-scription are provided in Appendices A and B. Appen-dices C and D contain the description of spin-lattice sim-ulations. The models for the domain wall and skyrmionare discussed in Appendices E and F, respectively. II. RESULTS AND DISCUSSIONA. Nonlinear σ -model We consider a chiral, uniaxial G-type antiferromagnetwith a simple cubic lattice with the lattice constant a .It can be characterized by the following Hamiltonian in- cluding exchange, anisotropy and DMI terms H = J S (cid:88) ρ , ρ (cid:48) µ ρ · µ ρ (cid:48) − K S (cid:88) ρ ( µ z ρ ) + S (cid:88) ρ , ρ (cid:48) d ρ , ρ (cid:48) · [ µ ρ × µ ρ (cid:48) ] . (1a)Here, J > S is the spinlength, µ is the unit magnetic moment in the lattice siteenumerated by vector index ρ = { i, j, k } with ρ (cid:48) runningover all nearest neighbors, K is the constant of uniaxialanisotropy and d ρ , ρ (cid:48) = − d ρ (cid:48) , ρ is the DMI vector[40–42].The dynamics of the magnetic moments is described bythe Landau–Lifshitz equation [43] ∂ t µ ρ = 1 (cid:125) S µ ρ × ∂ H ∂ µ ρ (1b)with (cid:125) being the Planck constant. The characteristiclength scale for the spin lattice is given by the magneticlength (cid:96) = a (cid:112) J / | K | .To develop the continuum counterpart of (1), we di-vide the 3D spin lattice by groups of spin octamers withspins being labeled from A ρ to H ρ within the given oc-tamer, see Fig. 1(a) and Appendix A. This approach is a3D counterpart of the tetramerization scheme used forthe case of a 2D AFM [38, 39]. The magnetic stateof each spin octamer is described by vectors of the to-tal magnetic moment m ρ , the N´eel vector (staggeredmagnetic moment) n ρ and auxiliary fields u ρ i , u ρ i , i = x, y, z , see Fig. 1(c). In the following, we performthe analysis in a long-wave approximation (i.e. spatialand temporal variations of the vector fields are slow) us-ing (cid:15) = (cid:112) | K | / J = a /(cid:96) → K ∼ (cid:15) , the DMI d ρ , ρ (cid:48) and timederivatives are of the order of (cid:15) . The ground state isgiven by | n ρ | = 1. This suggests that the length of themagnetization and auxiliary vectors is of the order of (cid:15) ,see details in Appendix A. Then, the relation betweenthe continuum counterparts of the order parameters andauxiliary fields is given by the linear expansion of Eq. (1b)with respect to (cid:15) m = − (cid:125) J S n × ∂ t n , u i = − a ∂ i n , u i = 0 . (2)Thus, three auxiliary vectors are determined by the spa-tial derivatives of the N´eel vector and other are zero dueto the lattice symmetry. This is similar to the case ofa 2D antiferromagnet [36, 37, 39], where the unit cellconsists of four neighboring spins arranged in square.In the main text, we focus on the chiral AFM slabswith the DMI of the bulk type which is commonly foundin AFM crystals [40, 41]. For completeness, the case ofthe surface DMI is discussed in Appendix B. For the caseof bulk DMI, the DMI vector is d ρ , ρ (cid:48) = d e ρ , ρ (cid:48) with e ρ , ρ (cid:48) being the unit vector in the direction from the spin ρ to ρ (cid:48) , see Fig. 1(d,e). The dynamics of the N´eel vector is FIG. 1.
Chiral G-type antiferromagnet.
Simple cubic crystal lattice with antiferromagnetically coupled spins (a) at theedge of a slab (b). The highlighted octamer ρ = { i, j, k } contains 8 spins, labeled as A , . . . , H (index ρ is omitted). Lessneighbors at the edge as compared to the bulk results in a tilt of spin moments from the anisotropy direction due to theexchange and DMI. (b) AFM slab is schematically colored to show two chiral domain walls, tilted at the edges of the slab.(c) The continuum counterpart of the octamer (a) with the N´eel vector n , magnetization vector m and auxiliary vector fields u x,y,z . (d) The spin octamer in the bulk and its continuum counterpart posses no tilt from the anisotropy direction. The DMIvector is shown for the bulk DMI type. (e) Parametrization of the N´eel vector n in the local spherical reference frame with thepolar and azimuthal angles θ and φ , respectively. governed by the equation, obtained within the harmonicexpansion of Eq. (1b) with respect to (cid:15) n × (cid:20) M s γ Λ ∂ tt n − A ∆ n − Kn z e z + D ∇ × n (cid:21) = 0 (3a)with n ≡ n ( r , t ) being a continuum (micromagnetic)counterpart of n ρ , A = J S / (2 a ) is the exchangestiffness, K = K S / (2 a ) is the anisotropy constantand D = dS /a is the micromagnetic DMI constant.The saturation magnetization of each sublattice is M s = gµ b S/ (2 a ) with g being Land´e factor and µ b is the Bohrmagneton and γ = gµ b / (cid:126) is the gyromagnetic ratio. Thecharacteristic scales of the magnetic field for K > B sf = √ Λ K/M s with Λ = 6 J S /a being the constant of the uniform exchange and spin-flipfield B x = Λ /M s , see Appendix C and D for comparisonwith spin-lattice simulations. The critical DMI value is D c = 4 √ AK/π .The complete formulation of the micromagnetic prob-lem includes boundary conditions for the vector n . Aconventional way to obtain the boundary conditionswithin the model of multiple sublattices is to considerthe difference in torques acting on the boundary spinsand bulk spins [25, 27, 44] and matching the equationsof motion for both of them. The linear in (cid:15) analysis ofthe discrete equations of motion Eq. (1b) for (100), (110)and (111) surfaces provides the following boundary con-ditions n × [2 A ( ˆ ν · ∇ ) n − D ˆ ν × n ] = 0 (3b) with ˆ ν being the surface normal. We note that whenconsidering other crystallographic cuts, Eq. (3b) is notchanged within the linear approximation.In contrast to antiferromagnets, the influence ofboundary conditions is well addressed in ferromagnetism.The behavior of ferromagnets is governed by the Rado–Weertman boundary conditions [45–49] with the DMI re-lated term [50]. The lateral confinement alters the shapeof chiral domain walls [51] and leads to the non-reciprocaldomain wall dynamics [52]. The surface twist of theground state [50, 53] and vortex-like textures at the sur-face [50, 53–57] are observed in chiral ferromagnets. Thisbehavior of antiferromagnetic solitons is not known yetand will be addressed in the following.To study static properties of the ground state and mag-netic solitons, we analyze the micromagnetic energy. Thecontinuum micromagnetic functional of the potential en-ergy for the σ -model, corresponding to the spin-latticeHamiltonian Eq. (1a) reads E = (cid:90) E d r , E = A ( ∂ i n )( ∂ i n ) − Kn z + D n · ∇ × n . (4)As follows from Eq. (2), the magnetization m is aslave variable, m = − M s / ( γ Λ) n × ∂ t n . Unlike spinchains [31, 33], the symmetry breaking term in (4) isabsent. This is a consequence of the possibility to de-rive the micromagnetic model in a spatially-symmetricway in two- [38, 39] and three dimensions. In crystals n θ st e y e x − w + w −
10 0 100 . .
951 Coordinate y/‘ n z D/D c = 0 . D/D c = 0 . D/D c = 0 . . . . D/D c θ s t f o r n , ◦ AnalyticsSimulations(a) (b) (c)
FIG. 2.
Ground state of a chiral antiferromagnetic slab. (a) Schematic representation of the changes of the groundstate in confined AFM sample (width 2 w < ∞ along e y direction). At side faces, the N´eel vector is tilted with the surfacetwist angle θ st . The ground state in the interior of the domain is shown within the transparent part of the sample. (b) Verticalcomponent of the N´eel vector at the center of the sample along the y-axis: for different strength of the bulk DMI D . Symbolsand solid lines correspond to simulations and analytics (Eq. (6)), respectively. (c) Surface twist angle θ st for different strengthof the bulk DMI. The simulations are carried out for a rectangular slab of 120 × ×
160 spins (24 (cid:96) × (cid:96) × (cid:96) ) with themagnetic length (cid:96) = 5 a . with the symmetry lower than the simple cubic one, thesymmetry-breaking term appears as the homogeneousDMI [40, 58]. The dynamic equation (3a) and bound-ary conditions (3b) can be recovered by the variation ofthe Lagrangian (cid:76) = M s γ Λ (cid:90) ( ∂ t n ) d r − E (5)taking into account orthogonality of m and n within (cid:79) ( (cid:15) ).A similar procedure can be applied to other types oflattices to derive a direct correspondence between the pa-rameters of the micromagnetic model and the spin latticeparameters including the boundary conditions. B. Ground state
The ground state of an achiral bipartite antiferromag-net corresponds to the direction of the order parameteralong the easy axis of the anisotropy e z . This is alsotrue for the boundary spins in the absence of DMI. Thepresence of the bulk DMI alters the order parameter uponapproaching the side faces. The order parameter acquiresa tilt at the face surfaces, indicated with a surface twistangle θ st in the schematics in Fig. 2(a). To describe thistwist, we consider a slab of the width 2 w along the y axisassuming the origin in the center of the sample. Usingthe parametrization n = { sin θ cos φ, sin θ sin φ, cos θ } inCartesian reference frame with θ = θ ( y ) = arccos n z and φ = φ ( y ) being polar and azimuthal angles, see Fig. 1(e),the minimum of the energy (4) is reached with n z = tanh arctanh (cid:115) − (cid:18) π DD c (cid:19) + w − | y | (cid:96) ,φ = − π , (6)where the requirement ∂ y θ ( | y | = w ) = D/ (2 A ) comesfrom the boundary conditions (3b). We note that for thecase of the AFM order parameter (director), the stateswith ( θ, φ ) and ( π − θ, φ ± π ) are equivalent. The groundstate is significantly altered by the boundary at the dis-tance about 2 (cid:96) . The order parameter possesses a twist bythe angle θ st = arccos n z ( w ) at the boundary, determinedby the strength of the bulk DMI, see lines in Fig. 2(b).The surface twist angle θ st grows almost linearly with D and can reach 30 ◦ when approaching the critical DMI D c , see line in Fig. 2(c). The surface twist angle is de-termined by the relation D/D c only because the criticalDMI D c holds the exchange and anisotropy scales in thefirst term in (3b), while D governs the twist itself. Theobserved effect is similar to the boundary twists observedin ferromagnets [50, 53].The analytical results shown in Fig. 2(b,c) are con-firmed by spin-lattice simulations, performed for a rect-angular slab containing 120 × ×
160 spins with themagnetic length (cid:96) = 5 a , see Appendix C for details. C. Domain wall
A G-type antiferromagnet supports translational or so-called phase domain states [59], where the wall separatesdomains with the swapped order of sublattices on theatomistic level [Fig. 1(b)]. These domains schematically e x e y e z δφ β into thebulkSurfaceBulk 0 0 . . . D/D c β f o r D W , ◦ z = 0 z = − ‘z = − ‘ . − π − π − π D/D c D o m a i n w a ll ph a s e φ AnalyticsSimulations − − −
101 2∆ /‘ y = 0 z = 0Coordinate x/‘ n z Analytics Simulations . . y = 0surface bulk z Depth | z | /‘ W i d t h ∆ / ‘ DD c = 0 . DD c = 0 . DD c = 0 . FIG. 3.
Chiral domain walls in antiferromagnetic slabs. (a) Schematics of the domain wall (DW) in a confined sampleof a constant width 2 w y < ∞ along e y axis. (b) At the sample’s edge, the domain wall acquires a tilt indicated with the angle β (cid:54) = 0. The tilt vanishes ( β = 0) in the bulk. The panel (b) shows the simulation results for D/D c = 0 . δφ of the domain wall phase φ far from the side faces. (c) Domain wall profile with p = − n z along e x axis. (d) Domain wall width ∆ expands at the surface z = 0 in comparisonwith the bulk. Solid lines and symbols correspond to the analytics (9) and simulations, respectively. (e) Domain wall phase φ = − π/ − δφ [see panel (b)] at the sample’s surface far from edges. (f) Tilt angle β for different DMI. Dashed lines are guidesto the eye. colored in red and blue, are shown in Fig. 3(a). We con-sider the translational domain wall initially located in yz plane with the origin of the reference frame lying at thecenter of the top surface. In the absence of DMI, thedomain wall plane is flat and perpendicular to the sidefaces of the sample being quasi-1D texture with cos θ = p tanh( x/ ∆) and φ = const, where ∆ = (cid:96) is the domainwall width and p = ± φ is not determined. A finite D leads to thepreferred domain wall chirality. Namely, the last term inthe energy density E = A ( ∂ x θ ) + K sin θ + D sin φ∂ x θ forces the stabilization of a Bloch-type domain wall withthe favorable chirality sin φ = C = sign( pD ).A lateral confinement of the domain wall leads to itsdeformation (bent and broadening) and change of its in-ternal structure, see Fig. 3(a,b). We start with the de-scription of the internal structure of the domain wallnear the top surface and far from the side faces of thesample. The domain wall profile can be described us-ing a 2D Ansatz cos θ ( x, z ) = p tanh[ x/ ∆( z )] and φ ( z ) = Cπ/ δφ ( z ) with x ∈ ( −∞ , ∞ ) and z ∈ ( −∞ , δφ ( z ) = φ exp zλ (7)with the parameter λ characterizing the penetrationdepth. The boundary conditions (3b) require ∂ z φ (0) = D/ (2 A ) and ∂ z ∆(0) = 0 at the top surface, while freeboundary conditions at z = −∞ are assumed for bothfunctions. The substitution of the Ansatz in (4) gives thefollowing effective energy density after integration along x axis: E dw = 2 K ∆ + A
12 + π ( ∂ z ∆)
6∆ + πpCD cos (cid:18) λD A e z/λ (cid:19) − D A ∆ e z/λ (2 − e z/λ ) , z ≤ . (8)Here, the first and the second terms correspond to theanisotropy and the part of the exchange energy densitydependent on the derivatives of θ . The third term repre-sents the part of the DMI energy density related to the ∂ x θ . It determines C in the same way as for the Blochdomain wall in the infinite medium. The last term in (8)originates from parts of the exchange and DMI energydensities related to the ∂ z φ . The domain wall width ∆is determined by the following variational equation, in-dependent of C and p : π A ∆ (cid:2) ∂ zz ∆ − ( ∂ z ∆) (cid:3) + 2 A ∆ + D A e z/λ (2 − e z/λ ) = 2 K, z ≤ . (9)with λ being an unknown parameter to be found from theenergy minimization, see Appendix E for details. We findthat the domain wall becomes wider near the top surfaceup to about 10% in comparison with the bulk value, seelines in Fig. 3(d). The function ∆ has a Gaussian-likeshape. Its characteristic half-width z of the ∆ is 1 . (cid:96) for D c and grows with the decrease of D being 2 . (cid:96) at D = 0 . D c . The twist of the domain wall at the topsurface, φ , increases with the the strength of the DMI,see line in Fig. 3(e). The penetration depth of the phase λ behaves in a similar way as z . It equals (cid:96) for D = D c and becomes 2 . (cid:96) for D = 0 . D c . Note, that the increaseof λ and z for smaller values of DMI are accompaniedby a rapid reduction of φ and width ∆(0) at the topsurface.We elaborate the model of the AFM domain wall in aslab using spin-lattice simulations. They show a reason-able quantitative agreement with the analytical predic-tions, see symbols in Fig. 3(c,d,e). The value of ∆ for | z | (cid:38) (cid:96) obtained in simulations is slightly larger than (cid:96) due to effects of discreteness and is reduced with smaller (cid:15) used for numerical investigations. In addition, we numer-ically analyze the domain wall behavior near the sampleedges. The domain wall plane possesses a twist, which isobserved as an “S-shaped” profile at the top surface, seeFig. 3(a,b) for schematics and simulations. We character-ize this distortion by the angle β with respect to the edgenormal within the plane of the top surface. While β = 0in the absence of DMI, a finite D leads to the increase of β up to 20 ◦ , see Fig. 3(f). There is a slow reduction of β to the equilibrium value β eq = 0 far from the sample’sedges with about of 60% of the surface value at the depth z = − (cid:96) . Thus, taking into account the variation of thedomain wall width near the surface and bend at the edgesof the sample, the bulk-like properties of this texture arepreserved for samples significantly thicker and wider than10 (cid:96) . D. Skyrmion
A skyrmion is a chiral texture stabilized by theDMI [31, 60]. In this section, we consider impact of the3D confinement on individual skyrmions in a chiral AFMslab with the focus on the modification of the shape andphase of the skyrmion, see Fig. 4(a). Even for ferromag-nets, a rigorous description of skyrmions of small [61] andlarge [62] radius is a complicated task, which is usuallyaddressed by asymptotic analysis or numerically. Often, models of circular domain walls or numerical integrationare utilized, which allows to explain current-driven dy-namics [9, 39, 63] and excitations [64]. To address a 3Dskyrmion texture, we describe them qualitatively usingan axially symmetric Ansatz θ = θ ( r, z ) and φ = φ ( χ, z ).To highlight the peculiarities of the confined geometry,we consider a semi-infinite slab with z ∈ ( −∞ ,
0] in an-alytics and a sufficiently thick and wide rectangular boxin simulations.We start with the analysis of skyrmions of small radius .Their bulk properties can be addressed with the linearAnsatz [65] θ sk sm ( r ) = π (cid:18) − r R bulksk sm (cid:19) r ≤ R bulksk sm , r > R bulksk sm ,φ ( χ ) = χ + ˜ φ (10)where the phase ˜ φ = const. Here and below we use thedefinition of the skyrmion radius R sk as θ ( R sk ) = π/ R bulksk sm = (cid:96) π | D | D c , sin ˜ φ = C = sign D, (11)see Appendix F for details. This corresponds to a Blochskyrmion with the radius linearly growing with D , seered line in Fig. 4(b). The skyrmion of small radius inthe sample with a sufficiently large lateral size is influ-enced only by the top surface z = 0. To address thisspatial confinement in the vertical dimension, we modifythe Ansatz (10) adding the dependence on the longitu-dinal coordinate z , namely R bulksk sm → R sk sm ( z ) and usingthe phase φ ( χ, z ) = χ + Cπ/ δφ ( z ) with the definitionof δφ ( z ) according to (7). The boundary conditions (3b)lead to ∂ z R sk sm (0) = 0 and ∂ z φ (0) = D/ (2 A ). Sub-stitution of this Ansatz into energy (4) and integrationalong the radial direction allows to obtain a variationalequation for R sk sm ( z ) with the parameter λ similar tothe case of the domain wall (9), see Appendix F for de-tails. The solution of the obtained equation shows thatthe skyrmion possesses a Gaussian-like bottle-neck shape.The skyrmion is narrow in the bulk and becomes widerwhen approaching the surface, see red line in Fig. 4(d).The skyrmion shows a mixed Bloch–N´eel texture at thesurface due to the twist governed by the boundary condi-tions, see red line in Fig. 4(i). Note, that a similar shapedistortion is observed for vortices in easy-plane ferromag-nets with surface anisotropy [66, 67]. Large radius skyrmions in bulk samples can be de-scribed as circular domain walls using the Ansatz [64]cos θ sk lar = tanh r − R bulksk lar ∆ , φ ( χ ) = χ + ˜ φ (12)assuming ∆ (cid:28) R bulksk lar . The energy (4) reaches minimum ‘ e x e y e z un s t a b l e R bu l k s k / ‘ SimulationsSmall R sk Large R sk un s t a b l e . . . . D/D c R s u r f s k / R bu l k s k .
05 1 . − − R sk /R bulksk C oo r d i n a t e z / ‘ AnalyticsSimulations .
05 1 . . /‘ AnalyticsSimulations − − π π DD c = 0 . x/‘ P h a s e φ o f n , r a d topcenterbottom . . π π π Strength of DMI
D/D c P h a s e φ o f n , r a d SimulationsSmall R sk Large R sk π π π − − φ C oo r d i n a t e z / ‘ π π π (a) (b)(c) (d) (e)(f) (g)(h) (i) (j) DD c = 0 . DD c = 0 . DD c = 0 . DD c = 0 . FIG. 4.
Skyrmion in antiferromagnetic slabs. (a) Schematics of a skyrmion in a thick AFM slab. Its radii at the surfaceand in the bulk are marked as R surfsk and R bulksk , respectively. In the simulations we consider a slab consisting of 150 × × (cid:96) = 5 a . (b,c) Skyrmion radius as a function of the DMI strength. “Unstable” marker shows the region where theskyrmion of small radius is unstable. Symbols and solid lines correspond to simulations and analytics (11), (13), respectively.Dashed line is guide to the eye. (d) Skyrmion radius for different axial cross-sections of the sample. Solid lines correspondto the Ansatz of the skyrmion of small radius. (e) Width of the skyrmion of large radius measured at R bulksk . (f) Skyrmionisosurfaces n x,y,z ( r ) = 0 in the center and at the bottom of the slab. While this is a purely Bloch skyrmion in the bulk, the n is tilted when approaching the surface. (g,h) Phase φ at the top surface for different strengths of the DMI at the distance x = R surfsk from the origin. Symbols and solid lines correspond to simulations and analytics, respectively. Dashed line is guide tothe eye. (i,j) Phase φ of n at the distance R bulksk from the origin for skyrmions of small and large radius, respectively. Notationsare the same as in panel (d). with this Ansatz at [68] R bulksk lar ≈ | D | /D c (cid:112) − D /D c , ∆ ≈ | D | D c , sin ˜ φ = C = sign D, (13)see black line in Fig. 4(b). While it is expected that theskyrmion radius should not be significantly influencedby the sample’s boundary, its width ∆ and phase φ arealtered due to confinement. Therefore, the 3D texturecan be described by the replacement ∆ → ∆( z ) and φ ( χ ) → φ ( χ, z ) = χ + Cπ/ δφ ( z ) in (12). We findthat the structure of the circular domain wall stabilized by the DMI is similar to the straight one, considered inSec. II C, see Fig. 4(e,j).To obtain a quantitative description of the skyrmionshape in a confined geometry, we elaborate the aboveanalytics by spin-lattice simulations performed for a slabconsisting of 150 × ×
100 spins with (cid:96) = 5 a . Thestability of the skyrmion is influenced by the discretenessof the system and 3D shape of the texture. We find thatthe skyrmion can be relaxed for D (cid:38) . D c with thesmallest bulk radius R bulksk ≈ . a , see Fig. 4(b). Theskyrmion radius in the bulk grows with D up to 8 . (cid:96) at 0 . D c . The largest size of the skyrmion is limitedby the lateral dimensions of the sample [50]. The ratioof the skyrmion radius at the surface and in the bulk R surfsk /R bulksk found numerically is in agreement with theanalytical model: it is large for small D and reduces to1 when approaching D c , see Fig. 4(c). The skyrmionpossesses a complex structure in the bulk as well as atthe surface. The longitudinal profile shows two maximalradii at the distance of about (cid:96) from the top and bot-tom surfaces, see symbols in Fig. 4(g). The phase of n possesses a radially dependent asymmetric surface twist,which is changed with DMI, see Fig. 4(g,h). While an-alytics quantitatively capture the spatial profile of thephase φ (Fig. 4(i,j)), only qualitative agreement is ob-tained for the bulk and surface skyrmion radii. III. CONCLUSIONS
We derive a nonlinear σ -model with boundary condi-tions for a uniaxial chiral antiferromagnet of G-type witha simple cubic lattice and DMI of surface and bulk types.We establish a correspondence between the spin latticeand micromagnetic parameters relying on the approachwith the N´eel vector order parameter n . The transitionbetween spin lattice and micromagnetic models requiressix auxiliary fields, determined by the spatial derivativesof the N´eel vector. The micromagnetic boundary condi-tions for the N´eel vector match the variational derivationfrom the micromagnetic Lagrangian and are similar tothe Rado–Weertman ones with the DMI term for ferro-magnets [45, 47, 50]. The difference lies in the symme-try of the order parameter: the states of vector–director n and − n are indistinguishable. A procedure describedhere for the case 3D antiferromagnets with a simple cubiclattice can be straightforwardly extended to other typesof lattices.The obtained model is applied to analyze the groundstate and magnetic solitons in a spatially confined sam-ple. In this discussion, we focused on the case when theAFM slabs possesses a bulk DMI. The order parameterin the ground state acquires a chiral surface twist at theboundaries due to the lack of neighboring spins and com-peting exchange and DMI energy terms. Depending onthe DMI strength, the value of the surface twist anglecan reach up to 30 ◦ . The noncollinear textures, suchas domain walls and skyrmions, become modified nearthe boundary with the characteristic penetration depthof about 5 magnetic lengths. The domain wall being lat-erally constrained, possessed an S-shaped bend at thesurface. Both, the domain wall and skyrmion becomeof the mixed, Bloch–N´eel type at the surface. The DMIforces the skyrmions and domain walls to become broadernear the surface. In particular, for skyrmions of small ra-dius, the radius becomes 10% larger when approachingthe face of the sample.The here discussed impact of the confined geometryand DMI on the static magnetic textures provides anestimate for the minimal dimensions of AFM sampleshosting chiral magnetic solitons with bulk-like properties. Furthermore, we note that the change of the size of thetextures when approaching the boundaries is expectedto alter their dynamic properties. In this respect, thepresented model can be applied for perspective designof AFM racetracks and description of AFM textures instructured samples [19]. ACKNOWLEDGMENTS
Authors thank Prof. Patrick Maletinsky, NataschaHedrich, Dr. Kai Wagner and Dr. Brendan J. Shields(University of Basel) for fruitful discussions. This workwas financed in part via the German Research Foun-dation (DFG) grants MA 5144/22-1, MC 9/22-1, MA5144/24-1, Alexander von Humboldt Foundation (Re-search Group Linkage Programme), and by the Ministryof Education and Science of Ukraine (Project 19BF052-01).
Appendix A: Description of the spin lattice
To describe a G-type antiferromagnet, we split thelattice into groups of octamers. Within a single oc-tamer enumerated by the vector index ρ = { i, j, k } ,the spins are labeled by the Latin letters A ρ ,. . . , H ρ ,see Fig. 1(a). Then, the coordinate of each spin is ρ + { α, β, γ } , where α , β and γ running 0, 1. In the fol-lowing, we use η, ζ ∈ { x, y, z } for subscripts and spatialderivatives and η, ζ ∈ { α, β, γ } in exponents. For exam-ple, (cid:80) η ( − η u η = ( − α u x + ( − β u y + ( − γ u z takeseight different values for different α, β, γ . The single in-dex ρ is omitted for simplicity. Then, the unit magneticmoment within an octamer reads µ ρ + { α,β,γ } = m + ( − ξ n + p ( α, β, γ ) , p ( α, β, γ ) = ( − ξ + η u η + ( − η u η (A1)where ξ = α + β + γ , the triple { α, β, γ } enumerate thespin within the given octamer and Einstein summationrule is used. In the following, we apply a multiscale anal-ysis to describe the micromagnetic transition from thespin lattice approach (1) using (cid:15) = (cid:112) | K | / J as a scal-ing parameter.To describe the behavior of the spin system in the con-tinuum limit, the following relations for the neighboringspins along η direction are used: V ρ ± ∆ η = V ( r ) ± (cid:15)(cid:96)∂ η V ( r ) + 2 (cid:15) (cid:96) ∂ η V ( r ) , V = A , . . . , H . (A2)Considering slow spatial and temporal variations of themagnetic moments, we rewrite the equations of mo-tion (1b) using dimensionless time τ = (cid:15) Ω t with Ω = J S/ (cid:126) .The rescaled anisotropy and DMI coefficients are k = (cid:15) and δ = (cid:15)d . This also implies that m and auxiliaryfields u η and u η are of the order of (cid:15) for the N´eel groundstate in the bulk.The linear expansion of (1b) at the site { α, β, γ } reads˙ n = − m × n + 2 (cid:96) n × ( − η ∂ η n + 4 n × [ p + ( − η u η ] , (A3)where overdot means the derivative with respect to τ . The expression (A3) represents eight equations with re-spect to different values of α, β, γ . The solution of (A3)within each octamer is given by Eq. (2). It provides therelations between the primary and auxiliary vector fields,describing each octamer.The harmonic expansion of Eq. (1b) provides equationsof motion for m . For the given { α, β, γ } , they read˙ m + ( − ξ + η ˙ u η + ( − η ˙ u η = − n × ( n z e z ) (cid:124) (cid:123)(cid:122) (cid:125) anisotropy +2 (cid:96)d n × [ ∇ × n ] (cid:124) (cid:123)(cid:122) (cid:125) DMI − (cid:96) n × ∆ n − m × p (cid:124) (cid:123)(cid:122) (cid:125) exchange − m + p ] × (cid:2) p + ( − ξ + η u η (cid:3) − ( − ξ + ζ (cid:96) n × ∂ ζ (cid:104) m + ( − ξ + η + δ [ η,ζ ] u η + ( − η + δ [ η,ζ ] u η (cid:105)(cid:124) (cid:123)(cid:122) (cid:125) exchange , (A4)where δ [ η, ζ ] is the Kronecker delta with respect to sym-bols η and ζ . Summation of (A4) for all possible valuesof α, β, γ within each octamer and excluding m leads toEq. (3a).The boundary conditions can be rigorously obtainedfrom the equations of motion of the boundary spins. Theequations of motion within the continuum limit are thesame in the bulk and at the surfaces, while the spinshave different number of neighbors and experience dif-ferent torques. The boundary conditions arise as thematch between boundary and surface torques. For exam-ple, considering a (111) surface with the normal vector ˆ ν = { , , } / √
3, the boundary spin is H ρ with absentneighbors G ρ + { , , } , E ρ + { , , } and C ρ + { , , } . Thisimplies H × ( G i +1 + E j +1 + C k +1 )+ d ( G i +1 × e x + E j +1 × e y + G k +1 × e z ) = 0 . (A5)Substitution of the expressions for spins (A1) allows toreduce (A5) to (3b).The total energy of the σ -model reads E tot = (cid:90) (cid:18) M s γ ˙ n + E (cid:19) d r , (A6)where E is the potential energy density introduced in (4). Appendix B: DMI of the surface type
The DMI of the surface type can be obtained for theDMI vector d ρ , ρ (cid:48) = d e z × e ρ , ρ (cid:48) . In this case, the energyof the surface DMI reads E surf dm = hD (cid:90) [ n z ( ∇ xy · n ) − ( n · ∇ xy ) n z ] d S, (B1)where h is the sample thickness, D has the same value asfor the DMI of the bulk type and the magnetic texture is assumed to be homogeneous along e z . The derivationof (B1) implies d = 0 for e ρ , ρ (cid:48) (cid:107) e z . This allows to derivethe equation of motion and boundary conditions for theN´eel vector similarly to the case of the bulk DMI: n × (cid:40) M s γ Λ ∂ tt n − A ∆ n − Kn z e z − D [( ∇ · n ) e z − ∇ n z ] (cid:41) = 0 (B2a) n × { A ( ˆ ν · ∇ ) n + D [ n z ˆ ν − ( ˆ ν · n ) ˆ ν ] } (B2b)with all derivatives within the e x,y plane. Appendix C: Spin-lattice simulations
We numerically solve the Landau–Lifshitz equa-tion (1b) with the Gilbert relaxation torque T relax = α g µ ρ × ∂ t µ ρ and α g being the relaxation constant forthe Hamiltonian (1a) using the spin lattice simulatorSLaSi [69]. To analyze the spin-flop and spin-flip be-havior, an additional term H zee = − (cid:80) µ gµ b S µ ρ · B zee with B zee being the external magnetic field is includedto (1a). To model an infinite medium, periodic bound-ary conditions are applied. We use the following pa-rameters: the spin length S = 1, the exchange inte-gral J = 2 . × − J, the constant of a single-ionanisotropy K = 9 . × − J, the Gilbert constant α g = 0 . d is varyingfrom 0 to 2 . × − J. The integration is performedusing the midpoint algorithm at GPU with the time step δt = 0 .
01 ps. The relaxation is performed during 2 ns.Simulations were carried out using the high performanceclusters at the HZDR [70] and TSNUK [71].Figs. 1(b,c) are built taking into account that the geo-metrical width of the sample is 2 w = 24 (cid:96) and the position0 . s p i n - fl o p s p i n - fl i p Field
B/B sf R e du ce d m ag n e t i c m o m e n t .
96 0 .
98 1 1 . . s p i n - fl o p × spins 5 × spins Periodic BC Tilted field B . . . s p i n - fl i p (a) (b)(c) FIG. 5.
Spin-flop and spin-flip transitions. (a) Symbolscorrespond to spin-lattice simulations and line is guide to theeye. Dashed lines correspond to the analytically found valuesof the spin-flop and spin-flip fields B sf and B x . Insets (b,c)show zoom of spin-flip and spin-flop regions in (a). Solidblue and dotted orange lines correspond to the samples withdimensions 50 × ×
20 and 100 × ×
50 spins, respectively.Dashed red and solid green lines correspond to simulationswith periodic boundary conditions (BC). The field is tiltedby 1 ◦ angle from e z in the xz plane. of the N´eel vectors at the boundary in simulations corre-spond to the effective width 2 w eff = 23 . (cid:96) . Appendix D: External magnetic field
The Hamiltonian describing the interaction of the spinlattice with the external magnetic field B zee reads H zee = − gµ b S (cid:88) µ µ ρ · B zee (D1)with the corresponding continuum counterpart E zee = − M s (cid:90) m · B zee d r . (D2)We relaxed the spin lattice exposed to an external mag-netic field using two staggered initial states: along andperpendicularly to the anisotropy axis. The energies ofthe stable states are compared to determine the phasetransition. The spin-flop and spin-flip transitions, forthe case when B zee is applied along the anisotropy axis e z , are shown in Fig. 5. Fig. 5(c) shows the dependencyof the spin-flop field B sf on the boundary conditions insimulations. Smaller samples have smaller B sf , while thesample with periodic boundary conditions, equivalent tothe infinite system, shows the exact agreement with the-ory. We note that the auxiliary fields u x,y,z and u x,y,z do not influence the spin-flop and spin-flip even for thecase of finite (cid:15) for the homogeneous texture. Appendix E: Analysis of the domain wall near thetop surface
To obtain the domain wall shape, we numerically solvethe variational equation for the domain width ∆( z ) (9)using the test value λ = (cid:96) . The obtained function is sub-stituted into the expression of the energy density (8) andintegrated as a function of λ . This allows to determinethe value of the penetration depth λ in the second or-der and substitute it back into (9) to repeat the iterationprocess until convergence. The relative accuracy of 10 − for the domain wall parameters can be obtained within3–5 iterations. The same procedure is used to analyzeskyrmions of small and large radii. Appendix F: Analysis of the skyrmion shape
The energy (4) in the cylindrical reference frame( r, χ, z ) reads E sk = A (cid:2) ( ∂ r θ ) + ( ∂ z θ ) + sin θ ( ∂ z φ ) (cid:3) + K sin θ + D (cid:20) sin 2 θ sin( φ − χ )2 r + sin( φ − χ ) ∂ r θ − sin θ∂ z φ (cid:21) . (F1)To analyze skyrmions of small radius, we substitute theAnsatz (10) into (F1), which leads to the effective energydensity E bulksk sm ≈ . A +2 πA (cid:20) ( R bulksk sm ) (cid:96) − C DD c R bulksk sm (cid:96) (cid:21) . (F2)The condition of the minimum of this expressiongives (11). The 3D Ansatz gives the effective energy den-sity E sk sm ≈ . A + 2 π A ( ∂ z R sk sm ) + 2 πKR − π DR sk sm cos (cid:18) λD A e z/λ (cid:19) − π D A R e z/λ (cid:16) − e z/λ (cid:17) (F3)with the variational equation for R sk sm π A∂ zz R sk sm + D A e z/λ (cid:16) − e z/λ (cid:17) R sk sm = 2 KR sk sm − πD cos (cid:18) λD A e z/λ (cid:19) . (F4)and the boundary conditions ∂ z R sk sm (0) = 0, ∂ z R sk sm ( −∞ ) = 0.The effective energy density of the large radiusskyrmion in the bulk (12) reads E bulksk lar ≈ πA (cid:20) R sk ∆ + ∆ R sk + R sk ∆ (cid:96) − C DD c R sk (cid:96) (cid:21) (F5)1with the minimum reached at (13). Taking into accountthe effect of the surface, the energy density reads E sk lar ≈ πA (cid:40) R bulksk lar (cid:2)
12 + π ( ∂ z ∆) (cid:3)
3∆ + 4 ∆ R bulksk lar (cid:41) + 4 πK ∆ R bulksk lar − π DR bulksk lar cos (cid:18) D A λe z/λ (cid:19) − π D A ∆ R bulksk lar (cid:16) − e z/λ (cid:17) . (F6)This expression leads to the variational equation π A (cid:2) ∂ zz ∆ − ( ∂ z ∆) (cid:3) + 2 A ∆ + D A e z/λ (2 − e z/σ ) = 2 K + 2 A ( R bulksk lar ) . (F7) with ∂ z ∆(0) = 0 and ∂ z ∆( −∞ ) = 0, c.f. (9) for a straightdomain wall.The difference between simulations and analytics is aconsequence of the simplified Ansatz (10) and (12), whichdoes not take into account a fine structure of the radialdependency of φ and asymptotics for θ at the origin andinfinity. For example, taking into account z dependenceof the skyrmion radius in (12) as R bulksk lar → R sk lar ( z ) inaddition to the function ∆( z ), one obtains the boundarycondition ∆(0) ∂ z R sk lar (0) + [ r − R sk lar (0)] ∂ z ∆(0) = 0.This shows that the condition ∂ z ∆(0) is not a strict oneif a fine structure of the soliton near the surface is takeninto account. [1] T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich,Nature Nanotechnology , 231 (2016).[2] O. Gomonay, T. Jungwirth, and J. Sinova, physica sta-tus solidi (RRL) - Rapid Research Letters , 1700022(2017).[3] V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono,and Y. Tserkovnyak, Reviews of Modern Physics ,015005 (2018).[4] H. Yan, Z. Feng, P. Qin, X. Zhou, H. Guo, X. Wang,H. Chen, X. Zhang, H. Wu, C. Jiang, and Z. Liu, Ad-vanced Materials , 1905603 (2020).[5] J. Barker and O. A. Tretiakov, Physical Review Letters , 147203 (2016).[6] O. Gomonay, T. Jungwirth, and J. Sinova, Physical Re-view Letters , 017202 (2016).[7] C. Jin, C. Song, J. Wang, and Q. Liu, Applied PhysicsLetters , 182404 (2016).[8] H. Xia, C. Jin, C. Song, J. Wang, J. Wang, and Q. Liu,Journal of Physics D: Applied Physics , 505005 (2017).[9] L. Shen, J. Xia, G. Zhao, X. Zhang, M. Ezawa, O. A.Tretiakov, X. Liu, and Y. Zhou, Physical Review B ,134448 (2018).[10] L. S´anchez-Tejerina, V. Puliafito, P. Khalili Amiri,M. Carpentieri, and G. Finocchio, Physical Review B , 014433 (2020).[11] S. S. P. Parkin, M. Hayashi, and L. Thomas, Science , 190 (2008).[12] J. ˇZelezn´y, H. Gao, K. V´yborn´y, J. Zemen, J. Maˇsek,A. Manchon, J. Wunderlich, J. Sinova, and T. Jung-wirth, Physical Review Letters , 157201 (2014).[13] X. Zhang, Q. Liu, J.-W. Luo, A. J. Freeman, andA. Zunger, Nature Physics , 387 (2014).[14] A. Manchon, J. ˇZelezn´y, I. M. Miron, T. Jungwirth,J. Sinova, A. Thiaville, K. Garello, and P. Gambardella,Reviews of Modern Physics , 035004 (2019).[15] M. S. Wornle, P. Welter, M. Giraldo, T. Lottermoser,M. Fiebig, P. Gambardella, and C. L. Degen, ArXiv e-prints (2020), 2009.09015v1.[16] W. Jiang, P. Upadhyaya, W. Zhang, G. Yu, M. B.Jungfleisch, F. Y. Fradin, J. E. Pearson, Y. Tserkovnyak,K. L. Wang, O. Heinonen, S. G. E. te Velthuis, and A. Hoffmann, Science , 283 (2015).[17] J. M¨uller, A. Rosch, and M. Garst, New Journal ofPhysics , 065006 (2016).[18] R. Khoshlahni, A. Qaiumzadeh, A. Bergman, andA. Brataas, Physical Review B , 054423 (2019).[19] N. Hedrich, K. Wagner, O. V. Pylypovskyi, B. J. Shields,T. Kosub, D. D. Sheka, D. Makarov, and P. Maletinsky,ArXiv e-prints (2020), 2009.08986v1.[20] M. A. Lund, K. Everschor-Sitte, and K. M. D. Hals,Physical Review B , 180412 (2020).[21] B. A. Ivanov, Low Temperature Physics , 635 (2005).[22] B. Barbara, Comptes Rendus Physique , 631 (2019).[23] E. A. Turov, A. V. Kolchanov, V. V. Menshenin, I. F.Mirsayev, and V. V. Nikolaev, Symmetry and physicalproperties of antiferromagnets (FIZMATLIT, Moscow,2001).[24] R. L. Stamps and R. E. Camley, Journal of AppliedPhysics , 3497 (1984).[25] R. L. Stamps and R. E. Camley, Physical Review B ,1919 (1987).[26] D. Ghader and A. Khater, Scientific Reports , 6290(2019).[27] D. Ghader and A. Khater, Journal of Physics: CondensedMatter , 315801 (2019).[28] I. V. Bar’yakhtar and B. A. Ivanov, Sov. J. Low Temp.Phys. , 361 (1979).[29] A. F. Andreev and V. I. Marchenko, Sov. Phys. Usp. ,21 (1980).[30] H. J. Mikeska, J. Phys. C 13 , 2913 (1980).[31] B. A. Ivanov and A. K. Kolezhuk, Low TemperaturePhysics , 275 (1995).[32] H.-J. Mikeska and A. K. Kolezhuk, in Quantum Mag-netism (Springer Berlin Heidelberg, 2004) pp. 1–83.[33] E. G. Tveten, T. M¨uller, J. Linder, and A. Brataas,Physical Review B , 104408 (2016).[34] O. V. Pylypovskyi, D. Y. Kononenko, K. V. Yershov,U. K. R¨oßler, A. V. Tomilo, J. Fassbender, J. van denBrink, D. Makarov, and D. D. Sheka, Nano Letters ,8157 (2020).[35] I. Affleck, Journal of Physics: Condensed Matter , 3047(1989). [36] N. Papanicolaou, Physical Review B , 15062 (1995).[37] N. Papanicolaou, Physical Review B , 12290 (1997).[38] S. Komineas and N. Papanicolaou, Nonlinearity , 265(1998).[39] S. Komineas and N. Papanicolaou, SciPost Physics (2020), 10.21468/scipostphys.8.6.086.[40] I. Dzyaloshinsky, Journal of Physics and Chemistry ofSolids , 241 (1958).[41] T. Moriya, Physical Review Letters , 228 (1960).[42] H. Yang, A. Thiaville, S. Rohart, A. Fert, andM. Chshiev, Physical Review Letters , 267210 (2015).[43] L. D. Landau and E. M. Lifshitz, Phys. Zs. Sowjet.;Reproduces in in Ukr. J. Phys. 53 25-35 (2008) , 153(1935).[44] W.-M. Huang, T. Hikihara, Y.-C. Lee, and H.-H. Lin,Scientific Reports , 43678 (2017).[45] T. Rado and J. R. Weertman, Journal of Physics andChemistry of Solids , 315 (1959).[46] M. Labrune and J. Miltat, Journal of Magnetism andMagnetic Materials , 231 (1995).[47] A. Hubert and R. Sch¨afer, Magnetic domains: The anal-ysis of magnetic microstructures (Springer Berlin Heidel-berg, Berlin, 2009).[48] V. V. Kruglyak, O. Y. Gorobets, Y. I. Gorobets, andA. N. Kuchko, Journal of Physics: Condensed Matter , 406001 (2014).[49] O. Busel, O. Gorobets, and Y. Gorobets, Journal ofMagnetism and Magnetic Materials , 226 (2018).[50] S. Rohart and A. Thiaville, Physical Review B , 184422(2013).[51] C. B. Muratov, V. V. Slastikov, A. G. Kolesnikov, andO. A. Tretiakov, Physical Review B , 134417 (2017).[52] S. Zhang and O. Tchernyshyov, Physical Review B ,104411 (2018).[53] S. A. Meynell, M. N. Wilson, H. Fritzsche, A. N. Bog-danov, and T. L. Monchesky, Physical Review B (2014), 10.1103/physrevb.90.014406.[54] M. N. Wilson, E. A. Karhu, D. P. Lake, A. S. Quigley,S. Meynell, A. N. Bogdanov, H. Fritzsche, U. K. R¨oßler, and T. L. Monchesky, Physical Review B (2013),10.1103/physrevb.88.214420.[55] Y. M. Luo, C. Zhou, C. Won, and Y. Z. Wu, AIP Ad-vances , 047136 (2014).[56] K. M. D. Hals and K. Everschor-Sitte, Physical ReviewLetters , 127203 (2017).[57] A. Raeliarijaona, R. Nepal, and A. A. Kovalev, Phys.Rev. Materials , 124401 (2018).[58] I. E. Dzialoshinskii, Sov. Phys. JETP , 1259 (1957).[59] S.-W. Cheong, M. Fiebig, W. Wu, L. Chapon, andV. Kiryukhin, npj Quantum Materials , 3 (2020).[60] A. N. Bogdanov, U. K. R¨oßler, M. Wolf, and K.-H.M¨uller, Physical Review B , 214410 (2002).[61] S. Komineas, C. Melcher, and S. Venakides, Nonlinearity , 3395 (2020).[62] S. Komineas, C. Melcher, and S. Venakides, ArXiv e-prints (2019), http://arxiv.org/abs/1910.04818v1.[63] H. Velkov, O. Gomonay, M. Beens, G. Schwiete,A. Brataas, J. Sinova, and R. A. Duine, New Journalof Physics , 075016 (2016).[64] V. P. Kravchuk, O. Gomonay, D. D. Sheka, D. R. Ro-drigues, K. Everschor-Sitte, J. Sinova, J. van den Brink,and Y. Gaididei, Physical Review B , 184429 (2019).[65] A. Bogdanov and A. Hubert, Journal of Magnetism andMagnetic Materials , 255 (1994).[66] O. V. Pylypovskyi, D. D. Sheka, V. P. Kravchuk, andY. Gaididei, Journal of Magnetism and Magnetic Mate-rials , 201 (2014).[67] O. V. Pylypovskyi, D. D. Sheka, V. P. Kravchuk, andY. Gaididei, Low Temperature Physics , 361 (2015),1501.06548.[68] Note, that Ansatz (13) also works for small radiusskyrmions [64].[69] “ SLaSi spin–lattice simulations package”.[70] “High Performance Computing at Helmholtz–ZentrumDresden–Rossendorf,” .[71] “High–performance computing cluster of TarasShevchenko National University of Kyiv,” http://cluster.univ.kiev.ua/eng/http://cluster.univ.kiev.ua/eng/