Enhanced Stability of Skyrmions in Two-Dimensional Chiral Magnets with Rashba Spin-Orbit Coupling

Recent developments have led to an explosion of activity on skyrmions in three-dimensional (3D) chiral magnets. Experiments have directly probed these topological spin textures, revealed their nontrivial properties, and led to suggestions for novel applications. However, in 3D the skyrmion crystal phase is observed only in a narrow region of the temperature-field phase diagram. We show here, using a general analysis based on symmetry, that skyrmions are much more readily stabilized in two-dimensional (2D) systems with Rashba spin-orbit coupling. This enhanced stability arises from the competition between field and easy-plane magnetic anisotropy and results in a nontrivial structure in the topological charge density in the core of the skyrmions. We further show that, in a variety of microscopic models for magnetic exchange, the required easy-plane anisotropy naturally arises from the same spin-orbit coupling that is responsible for the chiral Dzyaloshinskii-Moriya interactions. Our results are of particular interest for 2D materials like thin films, surfaces, and oxide interfaces, where broken surface-inversion symmetry and Rashba spin-orbit coupling naturally lead to chiral exchange and easy-plane compass anisotropy. Our theory gives a clear direction for experimental studies of 2D magnetic materials to stabilize skyrmions over a large range of magnetic fields down to T=0.

Skyrmions are topological objects that first arose in the study of hadrons in high energy physics, but in recent years they have been discussed in connection with a wide range of condensed matter systems including quantum Hall effect, spinor Bose condensates and especially chiral magnets [1][2][3] . There has been tremendous progress in establishing exotic skyrmion crystal (SkX) phases in a variety of magnetic materials that lack inversion symmetry, ranging from metallic helimagnets like MnSi 1,4 to insulating multiferroics 5 using both neutrons 4 and Lorentz transmission electron microscopy 6 . Skyrmions also lead to unusual transport properties in metals like the topological Hall effect [7][8][9] , and may be related to the observed non-Fermi liquid behavior [10][11][12] .
Spin-orbit coupling (SOC) in magnetic systems without inversion symmetry gives rise to the chiral Dzyaloshinskii-Moriya (DM) 13,14 interaction D·(S i ×S j ). This competes with the usual S i ·S j exchange to produce spatially modulated states like spirals and SkX.
The 2D case is particularly interesting. Even in materials that break bulk inversion, thin films show enhanced stability 15,16 of skyrmion phases, persisting down to lower temperatures. Inversion is necessarily broken in 2D systems on a substrate or at an interface, and this too may lead to textures arising from DM interactions. Spin-polarized STM 17 has observed such textures on magnetic monolayers deposited on non-magnetic metals with large SOC. Very recently, we have proposed 18 chiral magnetism at oxide interfaces. In these systems the 2D electron gas at the interface between two insulating oxides has a large gate-tunable Rashba SOC 20 that leads to a tunable DM exchange leading to spirals 18 and SkX phases at finite temperatures 19 .
Motivated by this, we investigate 2D chiral magnets with a free energy dictated by general symmetry considerations, given SOC and broken inversion in the zdirection. Our results are summarized in the T = 0 phase diagram in Fig. 1 as a function of perpendicular magnetic field H and anisotropy A. The easy-axis anisotropy (A < 0) has been studied before 2,16 , but the the easy-plane regime (A > 0) has not. It is precisely here that we find an unexpectedly large SkX phase. Skyrmions not only gain DM energy, but are also an excellent compromise between the field and anisotropy A > 0. Moreover, we show that the skyrmion has a nontrivial structure in the spatial variation of its topological charge density (see Fig. 2) in this easy-plane region.
Can such easy-plane anisotropy of the required strength arise naturally in real materials? We present a microscopic analysis of three exchange mechanisms -superexchange in Mott insulators, and double exchange and Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction in metals -and show that the same SOC that gives rise to the DM interaction D also leads to an easy-plane compass anisotropy A c . The compass term is usually ignored since it is higher order in SOC than DM. We show, however, that its contribution to the energy is comparable to that of DM, with A c |J|/D 2 1/2 for all three mechanisms, where J is the exchange coupling. This striking fact seems not to have been clearly recognized earlier, possibly because these microscopic mechanisms have been discussed in widely different contexts using different notation and normalizations. We also discuss how additional single-ion anisotropies enter the analysis. Our microscopic considerations should serve as a guide for material parameters of 2D chiral magnets such that a large SkX region can be probed experimentally.
Ginzburg-Landau Theory: The continuum freeenergy functional F [m] = d 2 rF(m) for the local magnetization m(r) of a 2D system is given by The isotropic term F iso = F 0 (m) + (J/2) α (∇m α ) 2 consists of F 0 , which determines the magnitude of m and a stiffness J that controls its gradient (α = x, y, z). For our T = 0 calculation we replace F 0 with the constraint m 2 (r) = 1. Broken z-inversion and SOC lead to the DM term . We can rewrite this as m·(∇×m) with a π/2rotation of m about the z axis. SOC also leads to the The A c > 0 "compass" terms give rise to easy-plane anisotropy, while the singleion A s term can be either easy-axis (A s < 0) or easy-plane (A s > 0). We define length in units of lattice spacing a so that J, D, A c and A s all have dimensions of energy. While the form of the free energy (1) follows entirely from symmetry, the microscopic analysis -presented in the second half of this paper -gives insight into the relative strengths of the various terms. The origin of both the DM and compass terms lies in Rashba SOC whose strength λ t, the hopping, in the materials of interest. Thus we obtain a hierarchy of scales with the exchange J D ∼ J(λ/t) A c ∼ J(λ/t) 2 . Naively one might expect the compass term to be unimportant, however its contribution to the energy O(A c ) is comparable to that of the DM term O(D 2 /J). Note that while the DM term is linear in the wave-vector q of a spin configuration, its energy must be O(q 2 ). Thus compass anisotropy, usually ignored in the literature, must be taken into account whenever the DM term is important.
We will show below that, for a wide variety of exchange mechanisms independent of whether the system is a metal or an insulator, the ratio A c J/D 2 1/2. We will also discuss the origin and strength of the single-ion A s term. Here we only note that the effective anisotropy in model (1) is governed by A = A c + A s , which is easy-axis for A < 0 and easy-plane for A > 0.
Phase Diagram: We begin by examining the T = 0 phase diagram of (1) for a fixed D J as function of magnetic field H = Hẑ and the dimensionless anisotropy parameter AJ/D 2 , which we explore by varying A s with A c J/D 2 = 1/2.
We extend the simple spiral above to incorporate more general 1D modulation described by m spiral (r) = sin[θ(Q 0 .r)]Q 0 + cos[θ(Q 0 .r)]ẑ, where θ varies only alonĝ Q 0 , chosen to bex without loss of generality. In contrast to the linear variation in the simplest ansatz, here θ(x) is an arbitrary function with m(x + R) = m(x) where R is the period. We numerically minimize (1) with the variational parameters θ(x) and R (see Supplement). This more general 1D periodic modulation stabilizes spiral solution relative to FM beyond |A|J/D 2 = 1 to 1.25 at H = 0; see Fig. 1 We next turn to H = 0. For the A > 0 FM state, the easy-plane anisotropy competes with the field alonĝ z so that the magnetization points at an angle θ tilt = cos Skyrmions: At finite magnetic field, we also need to consider the SkX phase 2 in addition to FM and spiral. A skyrmion 1 is a"hedgehog"-like spin-texture with a quantized topological charge or chirality q = (4π) −1 d 2 rm · (∂ xm × ∂ ym ), which is restricted to be an integer. For example, the q = −1 skyrmion in Fig. 1(a) is a smooth spin configuration with the topological constraint that the central spin points down while all the spins at the boundary point up.
The SkX state is a periodic array of skyrmions. It is often described by multiple-Q spiral condensation 23,24 .
Here we use the 'single-cell approximation' of Bogdanov 2 , where we impose the topological constraint for the centre and boundary spins within a unit cell. We then find the optimal skyrmion configuration within a single unit-cell, whose size R is also determined variationally. We discuss in the text the results from a 'circular-cell' approximation 2 . This leads to an effectively 1D problem (along the radial direction) that is computationally much simpler than the full 2D numerical minimization of the energy (1) using the conjugate-gradient method described in the Supplementary Information. We find the results of the two methods are very similar. In the circular-cell approximation, we replace the unit cell of the 2D crystal by a circular cell of radius R and take a skyrmion configuration m skyrmion (r) = sin θ(r)r + cos θ(r)ẑ with the topological constraint θ(0) = π and θ(R) = 0. We minimize the energy (1) with θ(r) and the cell radius R as variational parameters. To construct the crystal, we make hexagonal packing of the optimal circular cells and recalculate the energy by filling the space between the circles by up spins.
To get a preliminary idea about the stability of SkX, we simplify further. The linear ansatz 2 θ(r) = π(1−r/R) with the single parameter R, the skyrmion size, has the great virtue of being essentially analytically tractable. It leads to the dashed lines in the anisotropy-field phase diagram of Fig. 1(b) and gives us the first glimpse of the large SkX phase for easy-plane anisotropy.
We obtain the SkX-spiral and SkX-FM phase boundaries shown in Fig. 1(b) by numerical energy minimization using the more general form of eq. (2) and discretizing θ(r) on a 1D grid. We see that this confirms the qualitative observations from the linear approximation and yields an even larger SkX phase on the easy-plane side. Our 2D square cell calculations essentially reproduce the same phase diagram (see Supplementary Information).
Easy-plane vs. easy-axis anisotropy: Our results for the phase diagram in the easy-axis region (A < 0) agree well with previous studies 2, 16 . One might have thought that the perpendicular field H and easy-axis anisotropy would both be favorable for a skyrmion, all of whose spins are pointing up far from the center, but then the FM state is even more favorable.
The remarkable result in Fig. 1(b) is that the SkX phase is much more robust for easy-plane anisotropy (A > 0). We can physically understand this as follows. The skyrmion obviously gains energy from the DM term due to the twist in the spin configuration, but it is also the best compromise between field alongẑ and the easyplane anisotropy. This is why the large SkX region in the phase diagram is more or less oriented around H = 2A, the dashed line in Fig. 1  Here the core has a large 'transition' region (yellow-orange) from down (centre) to up (boundary) in m leading to an unusual two-peak structure for |2πrχ|.
The internal structure of skyrmion within a unit cell gives us further insight into the stability of the SkX phase. In Fig. 2 we plot m z (r) and the (angular averaged) topological charge density |2πrχ(r)|, where χ(r) = [m·∂ x m×∂ y m]/4π. It is conventional to define the 'core radius' of a skyrmion from the maximum of |dm z /dr|, which for the rotationally symmetric ansatz (2) is given by dm z /dr = |2πrχ(r)|.
For the easy-axis case the skyrmion core shows a conventional structure with a single peak in |2πrχ| in Fig. 2(a,c). In contrast, easy-plane anisotropy can lead to a non-trivial core with a double peak in |2πrχ(r)|; see Fig. 2(b,d). As the spins twist from down at the center (θ(0) = π) to all up (θ(R) = 0) at the boundary, it is energetically favorable to have an extended region where θ(r) θ tilt (defined above), the best compromise between the field and easy-plane anisotropy. As a result, |2πrχ| shows a two-peak structure with the topological charge split into two spatially separated parts.
Phase transitions: We next describe the nature of the various phase transitions in Fig. 1(b) within our variational framework. All the phase boundaries between the spiral state and either FM or SkX are first order transitions with a crossing of energy levels. On the other hand the SkX to easy-axis FM transition in the regime H > 2A is second-order with the optimal SkX unit cell size diverging at the transition; (see Supplementary Information). The SkX to tilted FM transition for H < 2A is first order, as the SkX cell size remains finite across it. There is a tricritical point where the SkX phase boundary intersects the line H = 2A. Another interesting feature of Fig. 1(b) are the reentrant transitions from FM → SkX → FM for AJ/D 2 > ∼ 1. Microscopic Analysis: Our results for the enhanced stability of skyrmions for easy-plane anisotropy followed from the phenomenological free energy (1) for 2D chiral magnets. We next present a detailed microscopic derivation of eq. (1) which shows that the parameter regime of interest arises naturally for three very different exchange mechanisms in the presence of SOC.
Moriya's original paper 14 considered superexchange with SOC, and this was further elaborated in a way relevant to our analysis in ref. 25 . The RKKY interaction with SOC was first discussed for spin-glasses 26 and the relation between DM and anisotropy was analyzed 27 in the context of quantum dots. Double exchange ferromagnets with SOC were analyzed in our recent work 18 . In all these cases, it was found by explicit calculation that A c |J|/D 2 = 1/2 (in the notation of this paper).
We sketch here a "unified" way of thinking about these very different problems. Consider the microscopic Hamil- Here t is the nearest-neighbor hopping on a 2D square lattice with sites r i , λ is the SOC, σ are Pauli matrices andd ij =ẑ × r ij /|r ij | with r ij = r i − r j . Our analysis can be easily generalized to further-neighbor hopping and arbitrary lattices.
The interaction H int can be chosen to model several different situations. (i) Hubbard repulsion H int = U i n i↑ n i↓ with U t at half-filling, gives rise to antiferromagnetic (AF) superexchange with SOC. (ii) Coupling of conduction electrons with a lattice of localized spins S i via H int = −J H i s i ·S i leads to Zener doubleexchange with SOC, where s i = (1/2) αβ c † iα σ αβ c iβ and the Hund's coupling J H t. (iii) The H int of (ii) with a Kondo coupling |J K | t leads to an RKKY interaction between moments mediated by electrons with SOC.
The common feature of (i,ii,iii) is that, in each case, the effective Hamiltonian can be derived by considering pairwise interaction between spins. We discuss (i) and (ii) below. The RKKY case (iii) is discussed in the Supplementary Information. To focus on the lowenergy sector, we consider a two-site problem with nearest neighbor sites i and j and rewrite H 0 as H 0 = −t αβ (c † iα [e iϑσ.dij ] αβ c jβ + h.c.) witht = √ t 2 + λ 2 and tan ϑ = λ/t. Next we gauge away the SOC with SU (2) rotations on the fermionic operators at the two sites, via a iα = [e −i(ϑ/2)σ.dij ] αβ c iβ and a jβ = [e i(ϑ/2)σ.dij ] αβ c jβ .
Expressed in terms of the transformed fermions H 0 looks like the kinetic energy without SOC, while the form of local H int is left unchanged, provided the spins are also suitably transformed. The usual analysis for superexchange for case (i) then leads 25 to the result H SE = J AF <ij> s i · R(2ϑd ij )s j where J AF = 4t 2 /U , and R(2ϑd) is the orthogonal matrix corresponding to a rotation by angle 2ϑ aboutd. In case (ii), the same R enters the double-exchange result for classical spins, where J F = κt with κ a constant that depends on the density of itinerant electrons.
At low-temperatures, the effective spin model for both cases (i) and (ii) can be written in a common form (after expanding the square-root in case (ii) and a sublattice rotation in case(i)). We get whereμ =x,ŷ. Here J =J cos 2ϑ withJ = J AF for super-exchange andJ = J F for double-exchange. The SOC-induced terms are the DM term with D = J sin 2ϑ and the compass anisotropy A c =J(1 − cos 2ϑ). Since tan ϑ = λ/t 1, we get the microscopic result A c J/D 2 1/2.
It is straightforward to derive the continuum free energy (1) from the lattice model (4). The only term in (1) that does not come from (4) is the the phenomenological anisotropy A s (m z ) 2 arising from single-ion or dipolar shape anisotropy 21 . For moments with S < 2, the single-ion anisotropy vanishes 22 . In some cases, a simple estimate of dipolar anisotropy is much smaller than the compass term 18 . For larger-S systems, the singleion anisotropy is non-zero and can even be varied using strain 16 . We cannot, however, ignore compass anisotropy since its contribution to the energy is comparable to the DM term, as already emphasized.
Conclusions: We have shown the enhanced stability of skyrmions when the effective anisotropy parameter (A c + A s ) > 0 is easy-plane. The compass term A c is intrinsically easy-plane and thus our results suggest that experiments should look for systems with suitable single ion anisotropies A s , or ways to tune it using strain, e.g., so as to enhance the SkX region. Theoretically, it would be interesting to study in the future the finite temperature phase diagram for easy-plane anisotropy and electronic properties, such as the anomalous Hall effect and possible non-Fermi liquid behavior of itinerant electrons coupled to the SkX spin texture in this regime.

I. SUPPLEMENTARY INFORMATION
Here we give details of the calculations reported in the main text. We discuss (1) the variational calculation of the phase diagram, (2) results for skyrmion length scales, and (3) the microscopic analysis of RKKY interactions wit SOC.

1) Variational calculation of the phase diagram:
We consider the FM, spiral and SkX phases in turn. We use A = A c + A s as the effective anisotropy, and omit additive constants in the energy, which are common to all phases.

Spiral:
For the spiral solution we take the general 1D periodic modulation m(x) = sin θ(x)x + cos θ(x)ẑ, and minimize the energy, where ∂ x θ = (∂θ/∂x). We use conjugate gradient minimization with respect to the size R and the function θ(x) which is discretized on a 1D grid. We use the periodic boundary condition θ(R) = θ(0) + 2πn where n is an integer. This form allows for a spiral solution with a net magnetization m z in the presence of a perpendicular magnetic field. For analytical calculation in zero field, one can take a more restrictive (linear) variational ansatz θ(x) = 2π(x/R). In this case the energy of the spiral can be easily evaluated by minimizing with respect to R. This gives the spiral pitch R = R sp = 2π(J/D) and the energy F sp = −D 2 /2J + A/2.

Skyrmion crystal:
We have discussed in the main paper the method used to construct a hexagonal SkX solution using the circular cell approximation 1 with rotationally symmetric form of eq. (2).
As discussed in the main paper, to qualitatively understand the stability of SkX over FM and spiral states, one can use a simple linear ansatz θ(r) = π(1−r/R) and minimize the energy by choosing an optimal R. This leads to the solution R sk ≈ πJ/D for the optimal skyrmion cell size, with the energy given by  Fig. 1(b) in the main paper. The symbols and parameters used are exactly the same as described in the caption for Fig. 1(b). Note that the two calculations, although quite different in their computational complexity, nevertheless lead to essentially identical results for the overall phase diagram.
Here Ci(x) = − ∞ x dt cos t/t is the cosine integral and γ is the Euler constant. The result for F sk makes it clear that SkX gains energy from both DM and Zeeman terms.
For the more general θ(r) variation within the circular cell approximation, we need to numerically minimize We need to find the optimal cell size R and optimal values of θ(r), which we discretize on a 1D grid in the radial direction. We have carried out 1D conjugate gradient minimization using Mathematica on a laptop, using grids of up to 250 points.  Fig. 2, however the nontrivial structure of the skyrmion core in the easy-axis case is qualitatively similar to that in the circular cell calculations. a full 2D minimization by discretizing the GL functional (1) over a square grid. For the 2D calculation, we used up to 100 × 100 grids with polar and azimuthal angles (θ(r), φ(r)) of m(r) at each grid point as variational parameters. The 2D conjugate gradient calculations are done using a Numerical Recipes 2 subroutine in C on a local cluster of computers. This 2D minimization is much more computationally intensive than the 1D calculation for the circular cell approximation.
The 2D square cell result shown in Fig. 3 for the phase diagram is essentially the same as that obtained from the circular cell calculation; see Fig. 1(b) in the main text. We show in Fig. 4 the internal structure of the skyrmion as calculated from the full 2D square-cell minimization. This figure should be compared with the results from a circular cell calculation in Fig. 2 of the main paper. Note that the parameters used here are slightly different from those used in Fig. 2, however the nontrivial structure of the skyrmion core in the easy-axis case -the two-peak structure in the topological charge density |2πrχ(r)| -is qualitatively similar to that in the circular cell calculations.
FIG. 5: Skyrmion length scales: Plots of the Hdependence of the skyrmion cell radius R sk (denoted by R for simplicity in main text) and the core radii defined by the location of the maxima of |dmz/dr|. For the ansatz of eq. (2), dmz/dr = |2πrχ(r)|. (a) In the easy-axis region, both R sk and the (inner) core radius Rin are finite at the first order spiral-to-SkX phase boundary, but R sk diverges while Rin remains finite at the second-order SkX-to-FM transition. (b) In the easy-plane region, there are two core radii corresponding to the two maxima in |2πrχ(r)|. These inner and outer core radii Rin and Rout, and the cell radius R sk , all remain finite at the two first order phase transitions out of the SkX phase.
we show the skyrmion cell size R sk and core radiii as a function of field for (a) easy-axis anisotropy with AJ/D 2 = −0.5 and (b) easy-plane anisotropy with AJ/D 2 = 1.35. As described in the main paper, and shown in Fig. 2, there is only one length scale associated with skyrmion core size for the easy axis case, where as two-length scales appear for the easy-plane side, near the re-entrant region of the SkX phase diagram. We also show the skyrmion cell radius in this plot and we see that this is the divergent length scale at the second order phase boundary between the SkX and out-of-plane FM in the easy-axis case in Figs. 5(a).

3) RKKY interaction in the presence of SOC:
Here we discuss the case of RKKY 3 interaction between local moments embedded in a metallic host described by eq (3). In this case, the magnetic exchanges 27 , namely the isotropic, DM and compass, between two moments at r 1 and r 2 turn out to be J 12 =J(r 12 ) cos 2ϑ 12 , D 12 =J(r 12 ) sin 2ϑ 12 and A 12 =J(r 12 )(1 − cos 2ϑ 12 ). Here r 12 = r 2 − r 1 , ϑ 12 = k R r 12 with k R ≡ (λ/ta) andJ(r) −(J 2 K a 2 /4π 2 t) sin (2k F r)/r 2 . This result is obtained 27 for k F r 12 1, where k F is the Fermi wavevector, and by approximating 2D tight-binding energy dispersion by a parabolic band, as appropriate for lowdensity of conduction electrons. Evidently, for λ t and k −1 F r 12 k −1 R , the ratio A 12 J 12 /D 2 12 1/2 is maintained. We consider a set of moments that are regularly distributed on a square lattice with a spacing a such that the ratio AJ/D 2 1/2 for nearest-neighbor exchanges. If we neglect longer-range part of the RKKY, then we obtain the effective spin Hamiltonian (4) of the main paper.