Chiral magnetism: a geometric perspective

We recast the model of a chiral ferromagnet with the Dzyaloshinski-Moriya interaction as a Heisenberg model with nontrivial spin parallel transport. The Dzyaloshinskii-Moriya vectors serve as a background SO(3) gauge field. A combination of analytical and numerical arguments suggests that the ground state of this gauged Heisenberg model in 2 spatial dimensions is a hexagonal skyrmion crystal in a wide range of applied magnetic fields.

Chiral magnetic order, exemplified by helicoidal and more complex periodic structures, has a long history in the field of magnetism. These periodic spatial modulations arise from a competition of the Heisenberg exchange and of a weaker Dzyaloshinskii-Moriya (DM) interaction induced by the relativistic spin-orbit coupling [1,2]. In a continuum model, where the local spin or magnetization is represented by a 3-component vector field m(x) [3] of unit length, the competing energies have densities The DM vectors d i determine the spatial period of the magnetic order.
In recent years, chiral magnetism has received renewed interest in connection with the experimental discovery [4] of the skyrmion crystal, a magnetic analog of the Abrikosov vortex lattice [5] predicted by Bogdanov and Yablonskii [6]. Both superconducting vortices and magnetic skyrmions are examples of topological solitons. Even the simplest theories allowing for such soliton lattices have a strongly nonlinear character, which makes finding analytical solutions a highly nontrivial problem [7]. Whereas Abrikosov found an exact solution for his vortex lattice in the Ginzburg-Landau theory of superconductivity, no such feat has been accomplished for a skyrmion crystal in models of chiral magnetism to the best of our knowledge.
In this Letter, we take a different perspective on chiral magnetism. Instead of treating its emergence as a result of competing spin interactions (1), we view it as an outcome of a modified geometry of spin parallel transport [8]. This approach is similar in spirit to the geometrization of gravity in theory of relativity [9].
The formal basis for this geometric perspective is an extension of the Heisenberg model (1a), whose energy is invariant under global rotations of the spin reference frame, to a gauged version of the same model, whose energy is invariant under local frame rotations. The analog in relativity theory is the extension from special relativity with Lorentz transformations to general relativity with arbitrary coordinate transformations.
The transition from global to local SO(3) rotations requires the introduction of an SO(3) gauge field (the spin connection) A i (x) and of the covariant derivative [10] The condition D i m = 0 specifies the rule of spin parallel transport: a spin moved from x to x + dx rotates by the angle A i dx i . Parallel transport has a nonvanishing curvature if a trip around a loop changes the spin orientation. A spin moved around the boundary of an infinitesimal parallelogram with sides dx and dx rotates by the angle F ij dx i dx j , where the curvature is the magnetic field of the SO(3) gauge theory. Much like the regular Heisenberg model (1a) imposes an energy penalty for any inhomogeneity in the spin field, ∂ i m = 0, the gauged version does so if the spin field fails to follow the rules of parallel transport, D i m = 0: To the zeroth order in A i , the gauged theory (4) reduces to the pure Heisenberg model (1a). The first-order term −∂ i m · (A i × m) yields the spin-orbit coupling (1b), provided that we set A i = d i . Thus the gauged Heisenberg model reproduces the basic chiral model (1). It is remarkable that the Dzyaloshinskii-Moriya interaction is now encoded in the spin geometry, wherein the DM vectors serve as the spin connection. At the second, and final, order in A i , the gauged theory acquires a local anisotropy term, The addition of this term is not necessarily a problem. For example, in a magnet with a cubic symmetry such as MnSi [11], this correction yields an m-independent constant and the gauged model (4) is equivalent to the chiral one (1). Furthermore, Shekhtman et al. [12] pointed out that, at least in simple microscopic models, the DM term (1b) is accompanied by the matching anisotropy (5 There are further reasons to study the gauged Heisenberg model (4). For example, the global SO(3) symmetry of the pure Heisenberg exchange (1a) ensures local conservation of the spin current j i = −m × ∂ i m. Effects of spin-orbit coupling, such as the DM term (1b), break the global symmetry of spin rotations and therefore violate the conservation of spin current. The gauged Heisenberg model (4) retains this symmetry in a local form, which still allows to define a conserved spin current as a functional derivative of the action with respect to the gauge field A i to obtain j i = −m × D i m [13].
Despite its inner beauty, the gauged Heisenberg model has not been widely used. Aside from situations with zero curvature [10,13], it has been applied to the DM interaction on a lattice [12] and most recently in the continuum [14]. The geometric approach to the magnetism of conduction electrons, wherein an SU(2) gauge field acts on the spinor wavefunction, has been advocated by Fröhlich and Studer [15] and others [16][17][18].
To demonstrate the utility of the gauged Heisenberg model, we have studied it in 2 spatial dimensions. Our analytical arguments and numerical simulations strongly suggest that the ground state of the model is a hexagonal skyrmion crystal. The skyrmion lattice is a notoriously fickle magnetic phase of matter that usually occupies a small portion of a phase diagram and requires finely tuned temperature and applied magnetic field [4].
In the gauged Heisenberg model, the skyrmion crystal appears to be robust, staying stable in a wide range of magnetic fields around zero.
The choice of the background gauge fields A i = d i is restricted by the symmetry of the magnetic material [6]. For symmetry classes C nv , D n (n = 3, 4, 6) and D 2d , the directions of the DM vectors are fixed and the only choice is the overall magnitude κ, Table I. Here e 1 and e 2 are unit vectors in the plane of the film and e 3 = e 1 × e 2 .
Bogomolny states. To gain theoretical understanding of this model, we employ the method of Belavin and Polyakov [19], who found a special class of field configurations m(x) minimizing the energy of the regular Heisenberg model. Known in the broader context of field theory as Bogomolny solitons [20,21], these field configurations obey a first-order differential equation and saturate a lower bound for the energy, E ≥ E min , given by a topological number. For the regular Heisenberg model, the Bogomolny equations and energy bounds are The topological charge Q is the skyrmion number of the physical region Ω. Thus skyrmions turned out to be elementary excitations of the Heisenberg model in d = 2 spatial dimensions with the energy 4π. In what follows, we will stick with the upper signs in Eqs. (6) and their gauged counterparts (8). The lower signs correspond to time-reversed situations.
The Bogomolny solutions are most efficiently represented by switching to complex coordinates and fields, Here The omitted term ∓ ∂Ω dx i A i ·m comes from the boundary ∂Ω and can be ignored in the thermodynamic limit. The explicit dependence of the energy bound on m(x), and not just on the topological invariant Q, creates a problem. The Bogomolny soultions are not even local energy minima: δE/δm(x) = ±F 12 = 0. Thus they are not even stationary states of the gauged model (4). We can still find a good use for the Bogomolny solutions if we extend the model by adding a magnetic field normal to the film plane, h = he 3 , This problem has been recently explored by several authors [14,[22][23][24][25][26]. When the applied field matches the SO(3) magnetic field, h = ±F 12 , the Bogomolny states are stationary and the lower bound is topological [22], with the same boundary term omitted. Bogomolny equations (8a) for the upper sign are listed in Table I. Their solutions, expressed in terms of the complex field ψ (for symmetry class D 2d ) orχ (for symmetry classes C nv and D n ), are superpositions of an arbitrary meromorphic function w(z) and a term linear inz [22]. In what follows, we shall focus on the symmetry class D n , for which the Bogomolny states have the form False vacuum. Bogomolny solutions (11) include just one uniform state,χ = ∞, or m = +e 3 . It has Q = 0 and E = 0, so it is appropriate to call it the vacuum. The field h = F 12 = +κ 2 e 3 is critical in the following sense: the vacuum is locally stable for h > κ 2 and unstable for h < κ 2 . However, even at the critical field the vacuum is not the state of lowest energy; as we show below, the true ground state is most likely an antiskyrmion crystal.
High-energy skyrmion crystal. Bogomolny solutions (11) also include skyrmion crystals. The function w(z) = C/z yields a Q = +1 skyrmion with m = +e 3 at the origin. A skyrmion lattice is given by a meromorphic function w(z) with periodically arranged simple poles. The Weierstrass ζ function [27,28] fits the bill. For a square lattice with the spatial period a, the ζ function has periods 2ω 1 = a and 2ω 2 = ia. For a hexagonal lattice, 2ω 1 = a and 2ω 2 = ae 2πi/3 . The ζ function is not strictly periodic but quasiperiodic: ζ(z + 2ω n ) = ζ(z) + 2η n , n = 1, 2. To make a periodic function, it suffices to add terms z andz in the right proportions [28]. The function yields square and hexagonal skyrmion crystals with energy E = 4π per unit cell of area S ∝ a 2 . Antiskyrmion. To construct a state with the energy below zero, we may try the following strategy. If adding a skyrmion raises the energy by 4π, perhaps we should try adding antiskyrmions [29]. For w(z) = 0 in Eq. (11), a Bogomolny state with Q = −1. Eq. (10), with the upper sign, suggests that the energy E = −4π is lower than that of the vacuum. This result is encouraging. However, we encounter two obstacles on our way to success. First, the Bogomolny bound (10) omits the boundary term − ∂Ω dx i A i · m [14] that adds 8π to the energy so that overall the energy E = 4π is positive. This is not a serious problem. In the thermodynamic limit, the bulk energy scales with the area of our two-dimensional sample, whereas the boundary energy scales with its perimeter and can therefore be neglected for large systems [30]. Therefore, if we could construct a Bogomolny state with a finite density of antiskyrmions, in the thermodynamic limit the energy density would be −4π per unit cell.
We now encounter a second problem. The lowest possible skyrmion number for a Bogomolny state (11), Q = −1, is achieved by setting w(z) to a constant. If w(z) is a polynomial of degree N > 1 then Q = N . (For N = 1 it could be −1 or +1, depending on the relative amplitude of the z andz terms.) We simply cannot construct a Bogomolny state with a skyrmion charge Q < −1. This problem is also not fatal. If Bogomolny states do not work, we can look beyond them.
Well-separated antiskyrmions. The Bogomolny state with Q = −1 (13) contains a single antiskyrmion with a characteristic radius 2/κ. At distances |z| 2/κ, ψ → 0 and m returns to the vacuum value +e 3 . Take a state with N antiskyrmions, If the antiskyrmions are separated by distances well beyond 2/κ we expect their interactions to be weak and the energy to be E ≈ −4πN in the thermodynamic limit. For two antiskyrmions with separation a 2/κ, we find that the energy behaves asymptotically as [30] E ∼ 8π + 512π (κa) 2 ln (Cκa), with C a numerical constant. Again, the leading term is positive because of the boundary contribution. The second term represents a long-range repulsive interaction between antiskyrmions. Antiskyrmion crystal. It is easy to go from a function with N single poles (14) to the Weierstrass ζ function with a periodic lattice of poles. In analogy with Eq. (12) for a skyrmion crystal, we write down an Ansatz for a (square or hexagonal) crystal of antiskyrmions: The trial parameter here is the lattice constant a. In the limit a → ∞, we recover the Bogomolny state with one antiskyrmion (13). At large antiskyrmion separations a 2/κ, the energy per unit cell asymptotically approaches −4π. By analogy with the Q = −2 case (15), we expect the leading correction to be (κa) −2 ln (κa). In the thermodynamic limit, the appropriate intensive variables are the (macroscopic) skyrmion and energy densities, ρ = Q/S and U = E/S: To corroborate these theoretical results, we ran Monte Carlo simulations for a lattice version of the gauged Heisenberg model [17], It reproduces the continuum theory (9) with an appropriate choice of the exchange constant J and the "gyromagnetic ratio" γ. For a hexagonal lattice of spins, J = 1/ √ 3 and γ = √ 3/2. R ij is an SO(3) matrix with the angle of rotation κ and the axis parallel to the link ij for the D n symmetry class. We worked on lattices with up to 30 × 30 sites and periodic boundary conditions. The critical field for the vacuum state m = +e 3 in the lattice model is h = 4e 3 sin 2 (κ/2) ∼ κ 2 e 3 as κ → 0. For κ = π/6, used in our simulations, the lattice and continuum critical fields differ by about 2%. This gives a rough estimate for the expected discrepancy between the lattice and continuum models. Monte Carlo simulations were run at a low temperature T = 0.01. For direct comparison with theory for the ground state, we subtracted the thermal energy of spin waves equal to T per spin.
To check the local stability and accuracy of the antiskyrmion crystal Ansatz (16), we used this trial state as a starting point in the simulations. Topological stability of the skyrmion number enabled us to work at fixed skyrmion densities. The Ansatz turned out to be quite accurate when antiskyrmions are well separated, |ρ| κ 2 . The energy density of the final state (open circles in Fig. 1) agrees well with the theory predictions (filled circles) in the range ρ = −0.02κ 2 . . .−0.01κ 2 , which includes the optimal density ρ 0 ≈ −0.0172κ 2 . With the local stability of the antiskyrmion crystal confirmed, we have also searched for alternative ground states by starting simulations with a random high-energy (T = ∞) state and quenching it to T = 0.01. The final states (open triangles in Fig. 1) usually turned out to be the same hexagonal antiskyrmion crystals as those obtained by starting from Ansatz (16). Fig. 2 shows two such final states. In other cases, the magnet was trapped in metastable states of higher energy.
It is thus reasonable to conclude that the ground state of the gauged Heisenberg ferromagnet (9) in the critical field h = κ 2 e 3 is the hexagonal antiskyrmion crystal with the lattice constant a 0 ≈ 8.19κ −1 . A similar conclusion was reached recently by Ross et al. [25]. The hexagonal antiskyrmion crystal is also a viable candidate for the ground state of the chiral model (9) in a broad range of fields, from 0, where it coexists with the hexagonal skyrmion crystal, and up to h ≈ +1.3κ 2 e 3 , where the vacuum state m = +e 3 takes over. Ansatz (16) can be readily extended to other regular lattices, allowing a straightforward way to compute the elastic properties of the antiskyrmion crystal and to determine the spectrum of low-frequency spin waves [31]. These results and further details will be described in a forthcoming paper [30].