Topological magnetic textures in magnetic topological insulators

The surfaces of intrinsic magnetic topological insulators (TIs) host magnetic moments exchange-coupled to Dirac electrons. We study the magnetic phases arising from tuning the electron density using variational and exact diagonalization approaches. In the dilute limit, we find that magnetic skrymions are formed which bind to electrons leading to a skyrmion Wigner crystal phase while at higher densities spin spirals accompanied by chiral 1d channels of electrons are formed. The binding of electrons to textures raises the possibility of manipulating textures with electrostatic gating. We determine the phase diagram capturing the competition of intrinsic spin-spin interactions and carrier density and comment on the possible application to experiments in magnetic TIs and spintronic devices such as skyrmion-based memory.

Both momentum space and real space topology can play crucial roles in magnetic topological insulators, where Dirac electrons interact with topological spin textures associated with magnetic moments 19,20 . Recently the stoichiometric compound MnBi 2 Te 4 , a TI containing a periodic sublattice of Mn 2+ ions with magnetic moment 5µ B , was synthesized and studied for the first time 21,22 , opening the new field of intrinsically magnetic TIs. The compound comprises alternating quintuple layers of the topological insulator Bi 2 Te 3 and ferromagnetic (FM) MnTe bilayers. A single septuple layer (SL) thin film exhibits ferromagnetic order. The exchange coupling between FM-ordered magnetic moments and Dirac surface electrons on the adjacent TI layer can open up a gap at the Dirac point and give rise to nontrivial topology in momentum space, characterized by quantized Chern number 23 . The resulting Chern insulator and zero-field quantum anomalous Hall state have recently been observed in transport measurements on exfoliated few-layer Mn-Bi-Te flakes.
Unlike TIs with randomly doped magnetic impurities, intrinsic magnetic TIs are stoichiometric compounds containing a lattice of magnetic atoms coupled to topological electrons. This feature not only promises a magnetic ordering and quantum anomalous Hall effect at temperatures as high as 50 K 24,25 , but also opens the exciting possibility of magnetic control of topological elec-tronic properties and electrical control of magnetic order. In particular, magnetic textures such as spirals and skyrmions can arise at surfaces and interfaces from the chiral Dzyaloshinskii-Moriya (DM) interaction 26-28 due to broken inversion symmetry. Magnetic TI surfaces provide a platform of this kind. This motivates us to consider the effect of real-space magnetic structures on topological Dirac electrons, which may enable the manipulating magnetic domains and textures by electric currents and electrostatic gating.
Motivated by recent advances in magnetic TIs and skyrmion physics, we study the magnetic phenomena arising from doping Dirac electrons. We focus separately on the cases of (a) a single added electron, (b) low carrier density and (c) high carrier density. We predict the formation of skyrmion textures with localized Dirac modes for (a) and (b) and the formation of stripes with 1d channels of chiral modes for (c). Because the ground state magnetic order parameter depends on carrier density, we obtain electrically tunable skyrmion and stripe phases whose periods vary with density.
As a potential application, we envision that intrinsic magnetic TIs can open the path for a memory device storing information using magnetic skyrmions in a low-carrier density material. The use of skyrmions for information storage has been reported in magnetic metals, where the exchange interaction couples itinerant electrons and localized spins. However, due to the large carrier density of metals, a large current density is required to drive skyrmions in writing and reading out data, which results in Joule heating and considerable power consumption. This drawback can be alleviated by working with skyrmions in a low-carrier density, bulk insulating material such as a magnetic TI. As we shall show, a skyrmion in a magnetic TI carries tightly bound electric charge. Under the right conditions, these charged skrymions are the only charge carrier at low doping. Therefore, a relatively small current is sufficient to drive the skyrmion motion.
Model. For the Hamiltonian of the 2d surface of a magnetic TI describing Dirac electrons coupled to a dense array of N classical spins S = (S x , S y , S z ) we take where the first two terms, and are a Rashba Hamiltonian for massless Dirac electrons and an exchange interaction, respectively. The last term consists of intrinsic spin interactions, including a ferromagnetic exchange interaction, DM interaction, Zeeman field, and out-of-plane anisotropy, namely where z is the out-of-plane direction. We have chosen a normalization so that S 2 r = 1. The form of the DM interaction above is appropriate for a lattice with C nv symmetry, favoring the formation of spirals rotating along the wavevector (Néel type), as opposed to orthogonal to the wavevector (Bloch type).
The phase diagram of the magnetic moments without coupling to electrons has been studied both theoretically 17,29 and experimentally 8,30-32 . With vanishing Zeeman field and anisotropy, the ground state of H S is a periodically modulated spiral texture with wavevector q ∼ D/A and Néel wall-like rotation. With increasing anisotropy, the spiral degenerates into a system of domain walls which is energetically favorable over the uniform state when D 2 /AK is above a threshold. Turning on an out-of-plane Zeeman field penalizes the large areas of antiparallel spins, resulting in a phase transition to a skyrmion crystal (SkX) above a critical field and to a uniform state above a higher critical field. As the temperature is increased the spiral or SkX order is destroyed.
We now consider the exchange coupling between localized spins and topological surface electrons. As we shall show below, the interplay between magnetic and electronic degrees of freedom leads to (1) a carrier density dependent magnetic phase diagram; and (2) chiral electronic states. Furthermore the exchange coupling is capable of driving the formation of charged skyrmions at small electron doping and Coulomb repulsion is capable of leading to a skyrmion Wigner crystal (SWX).
The remainder of this paper is organized as follows. In Sec. II we study the formation of skyrmion textures from doping a single electron and the formation of a SWX for a very dilute density of electrons. In Sec. III we study the stripe phase which forms at finite density using exact diagonalization and discuss the 1d channels of chiral modes bound to the stripes. We then present a phase diagram of competing stripe orders when H eS and H S are both important, and in Sec. IV we conclude.
Left, skyrmion profiles θ(r) where θ = 0 corresponds to +ẑ, for the shape parameters shown (c.f. Eq. (5)). Right, the cutoff kUV below which a skyrmion of the given radius becomes favorable upon doping a single electron. (Distances measured in lattice spacings.)

II. SKYRMION FORMATION
Eq. (4) has a uniform ground state S r = ±ẑ when the direct exchange interaction dominates, and the coupled electrons form two bands at energies ± J 2 + v 2 F |k| 2 . At charge neutrality, the lower band is filled and the addition of another electron incurs an energy cost of 2J unless a nontrivial texture forms. Eq. (4) is known to support skyrmions as excitations, so we investigate the electron spectrum in the presence of a skyrmion. We plot this spectrum in Fig. 1a. We observe mid-gap bound states, indicating the possibility that the system will spontaneously form a skyrmion to accommodate the extra electron and a skyrmion Wigner crystal for a very dilute density of electrons (Fig. 1b).
Let us work in the strong J limit, neglecting energy costs associated to H S . As a variational approach, we consider a 50×50 square lattice hosting a skyrmion texture and add a single electron to the charge-neutral state. Then we compare the total energy to that of the uniform state with a single added electron. We take a tightbinding approximation for H e suitable for the lattice, which we describe in App. B; it has single Dirac cone at k = 0. We choose a momentum cutoff |k| < k U V to indicate the region of linear dispersion. We parameterize the skyrmion profiles, shown in Fig. 2, using the following ansantz: To our knowledge this is a new parameterization, which was chosen for its compact description of a fourparameter family of skyrmion profiles: γ adjusts in-plane spin orientation (we set γ = π for a Néel-type skyrmion), b sets the radius, α − β is the vorticity, and α + β defines the sharpness of the domain wall. For example, a hard-wall magnetic bubble of radius b is recovered by taking α = β → ∞. Indeed, any smooth magnetic texture can be captured by some smooth function w(x, y). Skyrmion textures resulting from holomorphic w, a special case, were discussed in Ref. [33] and a one-parameter family of skyrmions 34 is recovered by taking α = 0, β = −1, γ = π/2. Higher-winding skyrmions, multiple skyrmions, or even a skyrmion crystal can be achieved by taking appropriate sums of the above ansatz. Such a sum may be preferable to the typical sinusoidal ansatz when the core size of the skyrmions is uncorrelated with their spacing, as in the skyrmion Wigner crystal we discuss shortly.
In Fig. 2 we show that the system prefers to form a Néel skyrmion (with vorticity 1) to accommodate the extra electron assuming a cutoff k U V a consistent with a Dirac cone dispersion, where a is the lattice spacing. A Bloch skyrmion would be favored for a non-Rashba spin-orbit coupling. The added electron binds to the skyrmion with profiles we numerically determined in Fig.  1. The radius b of the skyrmion satisfies a b 1/k U V . The cutoff dependence indicates that in real materials for which the Dirac cone is a low-energy approximation, the full dispersion may be relevant in determining skyrmion size. For instance from ARPES data 21 for MnBi 2 Te 4 the surface states' Dirac cone remains a good low energy approximation up to k U V ∼ .05Å −1 , beyond which the bulk bands are important. Since a ≈ 4Å our numerics would yield a skyrmion texture of b/a ∼ 3, an extended object comprising ∼ 10 Mn atoms. We note that although using the family of textures in Eq. (5) shows that in principle a skyrmion forms we have not optimized over all skyrmion shapes. Moreover, for real materials the direct spin-spin interactions, including the ferromagnetic exchange and DM interaction, may also help to stabilize skyrmions and affect their size.
Let us approximate the spectrum of Dirac modes in the presence of a skyrmion. Working in polar coordinates, we consider the skyrmion of radius b with δ b. This is an idealized approximation to the skyrmions shown in Fig. 1. Dirac bound states are localized near b, with a half integer angular momenta m satisfying |m|v F /b J. We solve the continuum Dirac equation in detail in Appendix A. We find mid-gap energies where f is a bounded, positive function of the skyrmion size whose explicit form is given in Eq. (A16). Thus small in-plane component δ > 0 results in a downward spectral shift, breaking particle-hole symmetry. This provides another way to understand why skyrmion formation is favorable.
The binding of Dirac modes to skyrmions was previously studied in Refs. [19] and [20]. Ref. [19] studied the general relations between magnetic textures and their induced electric charge, as well as their motion in an electric field (see also [35] and [36]). Ref. [20] studied the simplified case of a step-function out-of-plane skyrmion profile, and found that the charged skyrmion is only energetically favorable with an external magnetic field due to particle-hole symmetry at zero field. However, realistic skyrmion profiles are smooth and contain an in-plane region, which breaks particle-hole symmetry. Our results show that for realistic profiles with broken particle-hole symmetry, charged skyrmions may be favorable even at zero external field.
The binding of charges to skyrmions is also seen in the quantum Hall ferromagnet, for instance at the filling factor ν = 1. In this state, the Coulomb repulsion favors the spontaneous polarization of electron spins and doping the ν = 1 ferromagnet leads to electrically-charged skyrmion textures [37][38][39] . The quantum anomalous Hall state in intrinsic magnetic topological insulators differs from the quantum Hall ferromagnet in several key aspects. First, its ferromagnetism comes from the ordering of local magnetic moments (with ∼ 5µ B Bohr magneton per Mn atom in the case of MnBi 2 Te 4 ), rather than the spins of lowdensity charge carriers. Second, spin-momentum locking in the Dirac surface states of magnetic topological insulator leads to a strongly anisotropic magnetic response to the exchange field: an out-of-plane field opens a gap at charge neutrality, while an in-plane field does not. Therefore, the Chern-Simons effective theory for quantum Hall ferromagnets cannot be applied to describe the skyrmion physics in magnetic topological insulators.
Very dilute limit. Since the large-S local magnetic moments are treated as classical, the skrymions and Dirac electrons bound to them are immobile. When the density of added electrons is sufficiently low, the Coulomb repulsion is expected to drive a Wigner crystal phase of charged skyrmions. There are two distinct length scales in this phase: the lattice constants of the skrymion Wigner crystal (SWX), set by the density of doped electrons, and the skyrmion size, set by the exchange coupling and the spin susceptibility of Dirac electrons. The first length scale can far exceed the second. We emphasize that the SWX proposed here is driven by Wigner crystallization of charged skrymions. This mechanism is similar to the SWX phase observed in the lightly doped ν = 1 quantum Hall ferromagnet [40][41][42] , but different from the skrymion crystal phase in helimagnets that is driven by DM interaction.

III. STRIPE PHASES
In this section we consider a finite electron density above charge neutrality. We work with the same tightbinding model (App. B) and the strong J limit so that electrons actively dictate the ground state. We find the ferromagnetic state is unstable to a state with spatially modulated magnetization above a critical electron den- sity, as shown in Fig. 3a. Above the transition the preferred wavevector is ∼ 2k F . Let us comment on the band structure of Dirac electrons in the presence of a spiral phase. We expect a network of 1d modes localized along the regions where spins are in-plane, as in Fig. 3c. To see this, let us focus on a single such region, which we take along the y axis, and approximate the spin texture in the vicinity as S(x, y) = (1, 0, qx).
The normalization S 2 = 1 is maintained to linear order in x. Solving the corresponding continuum equation we obtain the bound state with linear dispersion E = −sgn(q)(v F k y + J). This state is Gaussian-localized in the x direction over a length ∼ v F /Jq. It follows that in the presence of a spiral texture S = (cos qx, 0, sin qx) well-localized, approximately degenerate chiral modes form when J/v F q, in which case Eq. (8) is a valid approximation.
In Fig. 3b we plot the full spectrum of Dirac electrons in a Néel spiral texture to show that the linear branch of chiral modes descends into the gap from the conduction band, giving another picture for why a spiral state becomes favorable. Moreover, a spiral state in the appropriate limit would exhibit highly anisotropic conductance; the conductance would be much greater in the y direction due to the 1d channels of bound states. Examples of such bound states are shown in Fig. 3c.
Since the exchange coupling J is small compared to the bandwidth of TI surface states (on the order of eV), the leading-order effect of J on the spin degrees of freedom is to modify the spin-spin interaction through the RKKY mechanism. The effective Hamiltonian for the array of spins, after integrating out the electrons, is given by: where χ ab tot (q) = J 2 χ ab eS (q) + χ ab S (q).
The spin susceptibility χ tot (q) is a 3×3 Hermitian matrix whose largest eigenvalue and corresponding eigenvector describe the wavevector and polarization of the groundstate spin texture. The first contribution is derived from second order perturbation theory in Ja/v F while the second term can be read off from H S when the Zeeman field is absent. Their explicit forms are given in Eq.s (13) and (18). In the strong J limit we drop χ S and focus only on χ eS . The latter can can be written as U †χ U with U a 3 × 3 unitary effecting a π/2 rotation aboutẑ and (13) with ξ ks = s|v F |k − E F and f ks = [1 + e βξ ks ] −1 the Fermi distribution 43 . The F a are given by where e iθ k = (k x + ik y )/k. In Fig. 4 we plot the dominant eigenvalue χ eS (q) of the spin susceptibility varying the electron doping k F (taking v F = 1). A square lattice with lattice constant a is used as a UV completion. The plot has a peak at q = 0 for small doping which trades dominance with a peak at q = 2k F at some critical k F . This corresponds to the transition from a uniform state to a 2k F spiral, in accord with Fig. 3a. The corresponding eigenvector describes a spiral texture with Néel wall-like rotation. For a non-Rashba spin-orbit coupling a Bloch wall-like rotation is . Left and right insets display the inter-and intraband spin susceptibilities, respectively. The 2kF peak is entirely due to the latter and the UV dependence is captured entirely in the former.
favored. Fig. 4 is also consistent with a skyrmion state comprising a sum of multiple magnitude-2k F wavevectors, because χ eS is radially symmetric up to isotropybreaking terms originating from the underlying lattice. We refer to Eq. (13) with s 1 = s 2 as the intraband susceptibility and with s 1 = −s 2 as the interband susceptibility. The interband susceptibility is a UV divergent contribution (cut off by the lattice spacing) which occurs even when the chemical potential is zero and depends weakly on q. The intraband susceptibility is independent of UV cutoff and dominated by Fermi surface contributions. These are both plotted in Fig. 4

(inset).
To see how the skrymion Wigner crystal (SWX) phase fits into the picture presented in this section, we remark that the presence of two length scales in the SWX phase at low doping leads to multiple peaks in the Fourier transform of its magnetic structure, and the largest peak is at wavevector q = 0. Thus the closest single-wavevector approximation is the uniform state, consistent with Fig. 3a and Fig. 4 in the low density regime. Higher order effects in the perturbative expansion in the exchange coupling J, higher harmonics in q and the inclusion of Coulomb repulsion are needed to fully capture the SWX phase.
Intermediate regime. We have previously worked at strong coupling, neglecting the effects of intrinsic magnetism. Now we allow the couplings of H S and H eS to be of the same order. The dominant eigenvalue of the total spin susceptibility (Eq. (12)) determines the ground state properties of the system, where χ eS was discussed in the previous section and χ S can be read off from H S in the case of zero Zeeman field. Since a square lattice was used for χ eS we use the same for χ S , although the only UV sensitivity will come from the interband susceptibility of χ eS , which affects the transition from q = 0 to q = 2k F . Thus we have with χ S (q) given by In terms of a rescaled coupling the dominant eigenvalue of χ S is peaked at q = 0 for d < 1 and at for d ≥ 1 if q is taken along an axis of symmetry. The ground state of the system is a spin texture with wavevector q = 0, 2k F , or q 0 , determined by the peak of the dominant eigenvalue of χ tot . This competition yields the phase diagram shown in Fig. 5, plotted in the (d, k F ) plane with A = K = 0.1, J = 0.5. This phase diagram exhibits several interesting features. Increased electron doping expands the q = q 0 phase, thereby acting as an "effective DM interaction," an effect noted in recent works 44,45 . The line q 0 = 2k F broadens into the RKKY phase with increasing exchange coupling J. When J vanishes, this phase vanishes and the q = 0 and q 0 phases are separated by the transition at d = 1. In the opposite limit J A, K, the q = q 0 phase vanishes and we recover the transition in Fig. 4 as a function of k F . We note that for weak exchange coupling J one expects these phases to be stripe orders (single-q) rather than multiple-q because for J = 0 the stripe phase is known to be the ground state at zero external field.

IV. DISCUSSION
We have investigated the phases of an array of magnetic moments coupled to Dirac electrons upon tuning the electron density, motivated by intrinsic magnetic TIs subject to electrostatic gating.
We found at very dilute densities electrons bind to magnetic skyrmions and we conjecture an SWX results; at higher densities the SWX gives way to spin spirals bound to chiral electron channels. The DM interaction resulting from broken inversion symmetry is essential for the formation of skrymions and stripes. Tuning the DM interaction and electron density reveals a phase diagram of stripe orders captured in Fig. 5. The interplay of real-space magnetic structures and topological Dirac electrons suggests the manipulation of magnetic domains and textures by electric currents and electrostatic gating is possible with low dissipation, an attractive prospect for skyrmion-based information devices. An interesting future direction is the study of the effects of skyrmion and spiral phases on transport phenomena and quantum oscillations. In a skyrmion crystal phase the Berry flux attached to each skyrmion should lead to topological DOS oscillations which impact all physical observables 46 . The emergent orbital magnetic field of the nontrivial spin textures leads to a topological contribution in the Hall resistivity ρ xy (B), observed in Mn 2 CoAl thin films 47 , Mn-doped Bi 2 Te 3 quintuple layers 48 , and correlated oxide thin films 49,50 at low temperatures. However, anomalous features in Hall resistivity could originate from surface and bulk ferromagnetism instead of skyrmion physics 51 . We leave the precise form of the contributions from topological magnetism to physical observables to future work. We hope that the magnetic phase transitions predicted in this paper may spark further interest in magnetic topological insulators.
Note added.
Recently two related works on spin textures in magnetic topological insulators have appeared 52,53 .
Yuan and Hiroki Isobe for participation and contribution in the early phase of this project. This work is supported by DOE Office of Basic Energy Sciences, Divi-sion of Materials Sciences and Engineering under Award DESC0018945. L.F. is partly supported by a Simons Investigator award from the Simons Foundation.