Twisted magnetic topological insulators

Motivated by the discovery of the quantum anomalous Hall effect in Cr-doped \ce{(Bi,Sb)2Te3} thin films, we study the generic states for magnetic topological insulators and explore the physical properties for both magnetism and itinerant electrons. First-principles calculations are exploited to investigate the magnetic interactions between magnetic Co atoms adsorbed on the \ce{Bi2Se3} (111) surface. Due to the absence of inversion symmetry on the surface, there are Dzyaloshinskii-Moriya-like twisted spin interactions between the local moments of Co ions. These nonferromagnetic interactions twist the collinear spin configuration of the ferromagnet and generate various magnetic orders beyond a simple ferromagnet. Among them, the spin spiral state generates alternating counterpropagating modes across each period of spin states, and the skyrmion lattice even supports a chiral mode around the core of each skyrmion. The skyrmion lattice opens a gap at the surface Dirac point, resulting in the anomalous Hall effect. These results may inspire further experimental investigation of magnetic topological insulators.

Motivated by the discovery of the quantum anomalous Hall effect in Cr-doped (Bi, Sb) 2 Te 3 thin films, we study the generic states for magnetic topological insulators and explore the physical properties for both magnetism and itinerant electrons. First-principles calculations are exploited to investigate the magnetic interactions between magnetic Co atoms adsorbed on the Bi 2 Se 3 (111) surface. Due to the absence of inversion symmetry on the surface, there are Dzyaloshinskii-Moriya-like twisted spin interactions between the local moments of Co ions. These nonferromagnetic interactions twist the collinear spin configuration of the ferromagnet and generate various magnetic orders beyond a simple ferromagnet. Among them, the spin spiral state generates alternating counterpropagating modes across each period of spin states, and the skyrmion lattice even supports a chiral mode around the core of each skyrmion. The skyrmion lattice opens a gap at the surface Dirac point, resulting in the anomalous Hall effect. These results may inspire further experimental investigation of magnetic topological insulators.

I. INTRODUCTION
Topological insulators (TIs) have attracted tremendous research interest in the past decade, owing to their fundamental physics and potential applications for next-generation devices [1,2]. Three-dimensional (3D) TIs have a gapped electron structure in the bulk and a metallic surface state whose existence is protected by the nontrivial topology of the bulk bands under time-reversal symmetry (TRS) [3][4][5]. The metallic surface state consists of an odd number of Dirac electrons. A unique nature that makes it different from an ordinary surface state is spin-momentum locking. The spin direction of the Dirac electron is dictated by its momentum direction, and vice versa. It cannot be gapped by nonmagnetic impurities because the TRS is still respected, and the bulk band topology remains nontrivial. In contrast, magnetism from the magnetic impurities would break the TRS. With a ferromagnetic order, the surface Dirac fermion acquires a mass, as demonstrated by experiments [6]. Such gap opening by ferromagnetic impurities can lead to various interesting physical phenomena, such as the half-integer quantum Hall effect [5], the topological magneto-electric effect [5,7], the induced magnetic monopole [8], the quantum anomalous Hall effect [9][10][11], and so on.
On the other hand, the TI surface itself breaks the spatial inversion symmetry between two magnetic adatoms. As a result, the magnetic interaction between the adatoms involves not only the Heisenberg interaction but also the Dzyaloshinskii-Moriya interaction (DMI) [12,13]. The latter is an effect rooted in spin-orbit coupling (SOC) that cannot be ignored on the TI surface, and is nonvanishing in noncentrosymmetric systems in general. TIs such as Bi 2 Se 3 [14] have large SOC, so that their bulk bands are inverted, resulting in their nontrivial band topology. It is conceivable that such a large SOC can bring about significant DMI components in the * These authors contributed equally to this work. † gangchen@hku.hk adatom interactions. In fact, strong DMI has been found in the Co/Pt interface [15][16][17][18][19], where the heavy Pt atoms provide strong SOC to induce DMI between magnetic Co atoms.
With the DMI, the magnetic adatoms on the TI surface could develop much richer magnetic orders, such as the spin density wave, spin spirals, and the skyrmion lattice [20][21][22][23][24][25], than a simple ferromagnet. The DMI between Fe adatoms on Bi 2 Se 3 was calculated to be quite large [26]. Experimental evidence has been found for the existence of magnetic skyrmions at the interface of ferromagnetic Cr 2 Te 3 and TI Bi 2 Te 3 [27]. Besides, a phase diagram that contains noncollinear magnetic structures such as spin spirals and the skyrmion lattice still needs to be established. Conversely, the magnetic structures have an impact on the TI surface states. The low-energy spectrum may no longer be a massless Dirac cone, and the wavefunctions may not be uniformly distributed in space (on a large scale compared with the lattice constant) anymore. It was found that although outof-plane ferromagnetically ordered surface impurities generate a band gap on the surface, the domain wall between up and down magnetization supports one-dimensional gapless chiral modes [28]. Magnetic textures such as domain walls and skyrmions become electrically charged when coupled to a 3D TI surface [29,30]. Magnetic skyrmions can give rise to bound states on TI surfaces [31]. In regard to transport properties, the single magnetic skyrmion was found to induce an anomalous Hall effect on TI surfaces that is different from the conventional topological Hall effect [32]. The skew scattering from the skyrmion is robust against geometric deformation [33]. Most of these analyses focused on Bloch-type skyrmions, and approximations such as hard-wall approximations were employed to model the skyrmion structure, inevitably losing the twisting information of the skyrmion.
In this paper we explore the phase diagram of a magnetic cobalt adatom lattice on a TI surface and investigate the effect of different twisted magnetic structures on the TI surface states. The remaining parts of the paper are organized as follows. In Sec. II, we describe the results of first-principles calculations on the geometry of cobalt atoms adsorbed on the surface of a typical TI, Bi 2 Se 3 . The magnetic properties of Co atoms are also discussed. In Sec. III, we calculate the magnetic interactions between the Co adatoms. Then Ginzburg-Landau theory is used to analyze the magnetic phase diagram of a Co adatom lattice. Magnetic structures of spin spirals and skyrmions are found to be achievable. In Sec. IV, we numerically solve the problem of TI surface Dirac electrons coupled to the twisted magnetic structures. The electronic states under the magnetic background of a spin spiral, a single skyrmion, and a skyrmion lattice are discussed. The paper is concluded in Sec. V.

II. GEOMETRIC AND MAGNETIC PROPERTIES OF THE ADATOM
The material Bi 2 Se 3 has a layered structure consisting of Se-Bi-Se-Bi-Se quintuple layers (QLs) stacking along the crystallographic c axis. Each atomic layer forms a triangular lattice, and they stack in the ABCAB sequence within each QL [see Fig. 1(a)]. The two topmost atomic layers of Bi and Se, enclosed by the dashed rectangle in Fig. 1(a), make up a honeycomb lattice when viewed from the z direction, as shown in Fig. 1(b). The equilibrium position of a magnetic cobalt atom adsorbed on the surface is determined by firstprinciples calculations (see details in Appendix A). It is found that among the three typical adsorption sites marked by colored crosses in Fig. 1(b), site A is the most stable one. The height of the Co atom is 0.24Å lower than the top Se layer. Thus the Co adatom can be seen as an interstitial impurity occupying the hollow site buried a little bit into the Bi 2 Se 3 (111) surface.
In the absence of spin-orbit coupling (SOC), the magnetic moment of the Co adatom is 1.0 µ B from our first-principles calculation. To investigate the electron configuration, the band structure of a 3 × 3 supercell is plotted, where the weight of different Co atomic orbitals is color coded as presented in Figs. 2(a)-2(h). It can be seen that the 3d states of the Co atom lie inside the bulk gap of Bi 2 Se 3 . Due to the C 3v point-group symmetry of the adsorption site A, the five 3d orbitals split into three sets in their energies: {d xz , d yz }, {d z 2 }, and {d x 2 −y 2 , d xy }. In contrast to the empty 4s orbitals, most of the 3d orbitals are fully occupied except the {d x 2 −y 2 , d xy } manifold, which has an unpaired electron. From this observation, it can be inferred that the electron configuration of the Co adatom is 4s 0 3d 9 , and the filling is schematically shown in Fig. 2 There is no doping from Co to Bi 2 Se 3 when the SOC in the system is shut off, only that the two 4s electrons are transferred into the 3d shell. This is different from the isolated Co atom, which has the electron configuration 4s 2 3d 7 . The existence of only one unpaired electron explains the 1.0 µ B magnetic moment.
The SOC slightly raises the magnetic moment of the Co adatom to about 1.2 µ B . The increase can be attributed to the doping effect of the Co electrons into the surface state of Bi 2 Se 3 . To see this, the band structures without and with SOC are compared in Fig. 3. The weight of the Co atomic orbitals is color coded. Without SOC, the bands originating from the Co atom are nearly flat, indicating that the Co 3d electrons are well isolated. When SOC is turned on, the Co bands overlap and hybridize with those from Bi 2 Se 3 . The band hybridization leads to a fraction of the Co electrons doping into the Bi 2 Se 3 surface. Because of such doping, the unpaired part of the Co electrons becomes larger considering Hund's rule, which accounts for the larger magnetic moment.
In Fig. 3 one can further see that the surface state of Bi 2 Se 3 is gapped, and only the upper massive Dirac cone is in the bulk band gap. This is because we use a slab of one quintuple layer in our first-principles calculations. The surface states of the top and the bottom layer overlap in the bulk and hybridize, leading to the gap opening [34,35].

III. MAGNETIC INTERACTIONS OF ADATOMS AND MAGNETIC PHASE DIAGRAM OF THE ADATOM LATTICE
To calculate the magnetic interactions between two Co adatoms on neighboring A sites, we use a 3 × 3 supercell as shown in Fig. 1(c). The spin-spin interactions between them can be described by the Hamiltonian (1) The system has a mirror symmetry σ v with respect to the plane which perpendicularly bisects the line joining atoms 1 and 2.
With such a symmetry constraint, the tensor ↔ J takes the general form   in which J ii (i = x, y, z) are Heisenberg interactions, D x and D z are DMIs, and Γ xz is the off-diagonal pseudodipolar interaction.
The results of the parameters in ↔ J are listed in Table I (see methods in Appendix B). For convenience, we have rescaled the interaction parameters for a normalized spin |S| = 1. It can be seen that the symmetry requirement of Γ xy = Γ yz = D y = 0 is confirmed by first-principles calculations. It is clear that the Co adatoms are mainly coupled by the Heisenberg ferromagnetic interaction. Besides, there are the DMI and pseudodipolar interaction. Both of them are relativistic SOC effects. The former is made nonzero by the broken inversion symmetry of the TI surface. Among them, the DMI can lead to more interesting magnetic structures than ferromagnetism. In the following, we assume a high coverage of Co adatoms located on the A sites and arranged in a triangular lattice as shown in Fig. 1(d).
To discuss the magnetic ground state, we derive the continuum model from the microscopic magnetic interaction model and make use of the Ginzburg-Landau theory. In the long wavelength limit, one can expand the spin distribution to the second order, and obtain the zero-temperature free energy under external magnetic field B = (0, 0, B) up to quadratic order of S(r) where the coefficients are related to the microscopic parameters as Here, a = 4.11Å is the adatom lattice constant, and a contribution from single-ion anisotropy A has been incorporated inÃ. The first-principles calculation finds a tiny A = −0.04 meV. Note that the low-energy continuum model, Eq. (4), has the same form as in Ref. [36] except for the latticedependent coefficients in Eqs. (5)- (8). This is because in the Ginzburg-Landau theory the underlying microscopic differences are coarse grained and it relies solely on the C 3v pointgroup symmetry. The magnetic structures of the free energy with the form of Eq. (4) have been discussed in detail [36]. Here, we review the main conclusions. The free energy can be split into two parts, F = F 0 + F , in which F 0 contains theJ 1 ,J 2 ,J 3 ,D, and B terms. One can make the approximationJ 1 ≈J 2 ≈J 3 and treat the difference F as a perturbation according to their values in Table I. The ground state of F 0 is a spin spiral S(r) = φ 1 sin(q · r) cos θ, φ 1 sin(q · r) sin θ, φ 1 cos(q · r) + φ 0 (9) with the propagation vector q = Q (cos θ, sin θ, 0), Q = |D/J|. The additional ferromagnetic component φ 0 comes from the partial spin polarization under the external magnetic field and satisfies the global normalization constraint [37] |S (r) The contour degeneracy of arbitrary θ for the free energy is quite similar to the case in spiral spin liquids [38]. The inclusion of F leads to a hexagonal warping in the free energy landscape and six discrete propagation directions favored by the ground-state spin spiral. As more terms beyond quadratic order are considered, spin spirals with different directions of propagation vectors can interact with each other, resulting in the superposition of several spin spirals. The leading term of such a correction is a quartic one ∆F ∝ d 2 r |S (r)| 4 . It makes possible the ground state being the skyrmion lattice (SkX) where The SkX is composed of three spin spirals whose propagation vectors have a relative angle of 120 • with one another. Their amplitudes are subject to the constraint 3 i=1 |φ i | 2 + φ 2 0 = 1. After comparing the free energy including the quartic term between the spin spiral, Eq. (9), and the SkX, Eq. (11), it was found that when the external magnetic field strength lies in the range the SkX will have lower energy than the spin spiral [36]. Here, J is the average ofJ 1 ,J 2 , andJ 3 . We draw the phase diagram on the B-D 2 x /|J| plane in Fig. 4(a), in which B is converted to units of milliteslas and J is the average of J 1 , J 2 , and J 3 .
The calculated material parameters are indicated by the vertical black line. The magnetic field strength corresponding to the SkX phase is about tens of milliteslas, which is quite accessible in the laboratory. The wavelength of the spin spiral is λ = 2π/Q = 364 nm. A typical spin distribution of the SkX phase is presented in Fig. 4(b). The distance between two adjacent skyrmion centers is 2λ/ √ 3 = 421 nm. The spin spiral and the SkX are both of the Néel type, with the spins rotating in the plane spanned by the wave vector and the surface normal. Note that although the external magnetic field influences the energetics among different magnetic structures and adds some ferromagnetic spin components into the spin spiral and SkX, it does not alter their periods.

IV. DIRAC ELECTRONS COUPLE TO MAGNETIC TEXTURES
In the previous section, we have investigated the twisting effect of Dirac electrons on the magnetic interactions between adsorbed Co atoms. It gives rise to versatile magnetic orders as shown in Fig. 4(a). Once the spin spiral or SkX is formed on the adatom lattice, the Dirac electrons could interact with the large-scale magnetic texture S(r) conversely via the Kondo (or Hund's) coupling. In general, this coupling tends to align the spin of the electrons with the orientation of local moment, an inverse process of the spin transfer torque effect. However, further complexities would inevitably emerge when the spinmomentum locking is incorporated for the Dirac electrons. Such a process can be modeled by the following Hamiltonian: The first term is a Rashba-type Hamiltonian describing the surface state of Bi 2 Se 3 [1], and the second term is Kondo (or Hund's) coupling between the electron spin and the spin of the magnetic adatom. Because the ferromagnetic exchange coupling J ex is typically much larger than the Zeeman energy arising from the external magnetic field, we neglect the latter hereinafter for the sake of simplicity. In this section, we explore the fate of Dirac electrons under the fixed back ground spin texture S(r) of a spin spiral, a single skyrmion, and a skyrmion lattice.

A. The spin spiral case
We first consider a spin spiral propagating along the x direction, whose spin configuration is S (r) = (sin Qx, 0, cos Qx) .
Because of the periodicity along the x and y directions, the electron wave function ψ(r) can be expanded by plane waves with a cutoff N, By taking the coupling parameter J ex = 1 in units of v F Q, we solve this eigenvalue problem in N = 5. The Fermi energy is determined by requiring half filling of the energy bands, and is tuned to E = 0. Then without the coupling to the spin spiral, the Fermi energy lies at the Dirac point of the surface state. From the band structures shown in Figs. 5(a) and 5(b), one can see that after coupling to the spin spiral, the TI surface remains a semimetal, with the Fermi energy still lying at a Dirac cone. The spatial distributions of two wave functions with k = (0, ±k 0 )Q on the negative branch of the Dirac cone are presented in Fig. 5(c). Here, k 0 = 0.1. By comparing with the spin configuration of the spiral structure, it can be seen that the wave functions localize around the domain walls with S z = 0. The wave function at the domain wall with S x >0 (S x < 0) has a velocity v y > 0 (v y < 0). The spatial distribution and velocity of the states near the Fermi energy are reminiscent of the chiral edge states of the quantum Hall effect.
For a small Q, the magnetic structure varies rather smoothly, so that an electron may only feel the environment surrounding it. In the S z 0 areas, the electron behaves as if its band structure has a Dirac mass whose sign is determined by the sign of S z . In the transition region where S z changes its sign, the absolute value of Dirac mass decreases to zero and then increases. The existence of gapless states at the domain walls can be understood qualitatively in this way.
We further comment on the possibility of thermal fluctuations in the spin spiral state. As the energy scale of the exchange is much lower than the electron bandwidth, the thermal fluctuations would affect the magnetism first. It is known that, the thermal fluctuations often convert the coplanar spin spiral state into a collinear spin density wave state before the system enters the paramagnetic state. If the collinear spin density wave state modulates with the S z component, the chiral edge mode associated with the domain walls persists until the system becomes paramagnetic.
It is illuminating to make a connection with the early work on magnetic domains on the surface of topological insulators. There, on the domain wall between two neighboring ferromagnetic domains, there exists a conducting chiral edge mode. This is because the neighboring ferromagnetic domains with opposite magnetic orders perpendicular to the surface have opposite Chern numbers [1]. This domain chiral edge mode has been realized experimentally on the surface of magnetic topological insulators [39,40], and with a more controlled method with the strong permanent magnet placed atop the GaAs/AlGaAs quantum Hall heterojunction [41]. The spin modulation of the spin spiral state with one spiral period can be approximately viewed as one domain wall, and there will naturally be a chiral edge mode associated with it.

B. The single-skyrmion case
Being an example of a topologically nontrivial magnetic texture, a skyrmion could influence the flow of electron spins and even nucleate itinerant electrons through the spin transfer torque effect. This is quite different from the spin spiral case. Before considering the SkX, we first clarify the physics of Dirac electrons within a single skyrmion in order to get some insight. To achieve this, the effective Hamiltonian in Eq. (13) is applied to an open system with the disk geometry of radius R as shown in Fig. 6(a). Note that only the spin momentum σ is involved in the exchange interaction and the possible contribution from the orbital momentum L = r × p due to the nucleation is not considered at this stage. As analyzed in the previous section, the skyrmion on the surface of Co/Bi 2 Se 3 is of the Néel type and hence can generally be written in the polar coordinates r = (r, θ) as so that at the skyrmion core r = 0 we have S(r) = [0, 0, 1] and at the skyrmion boundary r = R we have S(r) = [0, 0, −1]. It can be verified that the z component of the total angular momentum J z = −ı ∂ θ + /2σ z commutes with the effective Hamiltonian in Eq. (13) and thus provides a good quantum number j z . The wave function has the general form with an integer , satisfying J z Φ (r, θ) = j z Φ (r, θ) with j z = + 1/2. The Hamiltonian in Eq. (13) can be further separated into angular and radial parts, and the latter reads where The radial Hamiltonian can be solved in each subspace of fixed angular momentum quantum number by expanding the radial wave functions u (r) and v (r) in the Fourier-Bessel series with a cutoff N [42], The orthonormal basis is defined by where the parameter j ,n is the nth zero of the -order Bessel function of the first kind J (x). The Dirichlet boundary condition has been assumed implicitly. This approximation reduces the radial Hamiltonian to a 2N × 2N matrix eigenvalue problem where ψ T = (u 1 , . . . , u N , v 1 , . . . , v N ). The matrix elements of the Hamiltonian are given by (C , ) n,n = J ex 2 R 0 cos π R r ϕ ,n (r)ϕ ,n (r)rdr, The Fermi velocity of a Bi 2 Se 3 surface Dirac electron is about v F = 329 meV nm. For numerical convenience, a characteristic length r 0 is defined such that the energy unit v F /r 0 = 100 meV, which gives r 0 ≈ 3.3nm. Based on the analysis of density functional theory (DFT) and Ginzburg-Landau theory, the typical radius of a single skyrmion is estimated to be R = 60 in units of r 0 . With the coupling J ex = 40 meV and a cutoff N = 100 on the Fourier-Bessel basis, the low-energy spectrum is computed and depicted in  Fig. 6(b). In the total angular momentum space, the Dirac electrons open an energy gap immediately, and a chiral edge model emerges concurrently inside. This can be understood through the skyrmion topology. The swirling structure of a skyrmion is characterized by the topological charge, which counts the times that S(r) wraps a virtual sphere. For a single Néel-type skyrmion with a unit topological number shown in Fig. 6(a), the winding of spins divides the disk into two regions bounded by r = R/2. Both the net magnetization and its resulting Dirac mass have opposite signs inside and outside the boundary. Therefore there must be a closing and reopening of the bulk energy gap when crossing from one region to another, which produces edge modes near the boundary. In Fig. 6(c), the probability density |Φ | 2 for one of such edge states with near-zero energy is shown and compared with the skyrmion spin texture. In addition, the degeneracy of states with ± j z is lifted by the coupling due to the spin-momentum locking, and only the one with smaller energy can be localized on the boundary. This endows the edge mode with a chirality.

C. The skyrmion lattice case
To investigate the effect arising from the pure spin texture of the skyrmion lattice, we neglect the possible ferromagnetic component in the spin configuration. Although the SkX phase has a lower energy than the ferromagnetic and spin spiral phase only in the application of the external magnetic field, there are energy barriers among them, meaning that even though the magnetic field is turned off, the SkX state can still be metastable, with no ferromagnetic component.
The SkX solution of the Ginzburg-Landau free energy is given by Eq. (11). After turning off the ferromagnetic component φ 0 , the spin distribution of the SkX under investigation is which corresponds to the case of φ 1 = 1 and θ 1 = 0 in Eq. (11). The wave vector Q = |D/J| implies a SkX lattice constant of 421 nm. It has the symmetry of a triangular lattice as shown in Fig. 4(b). Therefore Bloch's theorem can be employed to expand the wave function as where k lies in the first Brillouin zone and b 1 , b 2 are the two primitive vectors of the reciprocal lattice. As in the spin spiral case, Hund's coupling is also chosen to be J ex = 1 in units of v F Q. We use the cutoff |n 1 | and |n 2 | ≤ 5 to solve the eigenvalue problem, and the band structure is plotted in Fig. 7(a).
In contrast to the spin spiral case, when the Dirac electrons couple to a skyrmion lattice, an energy gap opens at the Fermi level. Such a gap opening goes beyond the first order perturbation, because the second term of Eq. (13) does not couple the two degenerate states at the unperturbed Dirac point when S z (r) = 0. Hence the energy gap is small compared with the value of J ex . The spatial distributions of the wave functions of the conduction and valence states at the Γ point are presented in Fig. 7(b). It can be seen that the conduction states concentrate around the skyrmion cores, while the valence states concentrate around the boundaries between skyrmions. It is more comfortable for the Dirac electrons to settle in the regions where S z ≈ 0 as expected. The wavefunction distributions are similar to the case of topologically trivial electrons coupled to antiferromagnetic skyrmions [43], in which, when the electrons concentrate at the skyrmion centers, they are confined in a ringlike configuration. The localized electrons charge the skyrmions and have certain overlaps with each other. As a result, one can construct a tight-binding model through the effective Wannier orbitals on a charged SkX and investigate its low-energy properties. This physics has been discussed in a very recent work [44]. The coupling between Dirac electrons and the magnetic structure breaks the time-reversal symmetry. Because a gap opens at the Dirac point, the Berry curvatures in general can be nonzero. We calculate the Berry curvatures of the valence and conduction bands around the Γ point and present the results in Fig. 7(c). It is found that in this region, most of the Berry curvatures distribute near the Γ point and have opposite signs for the conduction and valence bands. Because Bi 2 Se 3 is usually n doped, there are some electrons at the bottom of the conduction band. The nonvanishing Berry curvature of these electrons will bring about the anomalous Hall effect (AHE) on a SkX [45]. It should be emphasized that the AHE in this system is not quite like the "conventional" one, in which the Hall resistivity is empirically proportional to the magnetization, which is zero in our case.
Phenomenologically, the AHE is similar to the topological Hall effect (THE) [46] commonly found on skyrmion lattices. In the THE, the electron experiences an emergent magnetic field that is essentially the real-space Berry phase of the electron hopping on the magnetic texture due to the scalar spin chirality, with the electron spin aligned with that of the magnetic moment. In our case, the electron is described by a Rashba Hamiltonian with a strong spin-momentum locking, and already has the momentum-space Berry curvature. The spin-orbit coupling prevents the electron spin from aligning with the local magnetic moment. Still, the electron spin rotates slightly during the hopping process, and a Berry phase is accumulated for a closed loop in the real space. So this is an example where both real-space Berry curvature and momentum-space Berry curvature are present. The Berry curvature distribution for both bands is depicted in Fig. 7(c).

V. DISCUSSION AND CONCLUSIONS
In our first-principles calculation of the magnetic interactions between the Co adatoms, Bi 2 Se 3 is assumed to be intrinsic, i.e., the Fermi surface of the surface electrons is right at the Dirac point. In this case, the magnetic interaction mainly comes from direct exchange and superexchange mechanisms. In reality, Bi 2 Se 3 crystals grown in laboratories are n doped due to the Se vacancies and antisite defects [47,48]. Nevertheless, the Fermi level can be tuned by further doping [49,50]. If the surface electrons have finite Fermi surface, the Ruderman-Kittel-Kasuya-Yosida (RKKY) mechanism will contribute to the magnetic interactions [28]. The RKKY interaction was found to be always ferromagnetic when the Fermi level lies near the surface Dirac point [28]. Also, the RKKY interaction has a DMI component when considering the spin-orbit coupling [51]. Therefore, even for the case of finite doping, the magnetic interaction tensor (2) remains valid, and the Heisenberg part retains its ferromagnetic nature. We expect only quantitative change in the phase diagram of Fig. 4.
The discussion of the topological surface electrons coupling to the magnetic structures can be generalized to other topological materials, especially stoichiometric magnetic topological insulators such as MnBi 2 Te 4 [52,53], as mentioned in Ref. [54]. These platforms are easier to fabricate in laboratories.
In conclusion, we calculate the magnetic interactions between cobalt adatoms on the topological insulator Bi 2 Se 3 (111) surface. The Heisenberg part of the interaction is ferromagnetic. Because of the broken inversion symmetry and strong spin-orbit coupling on the surface, there is also Dzyaloshinskii-Moriya interaction, which twists the spins of Co atoms from a perfectly parallel alignment. We use the Ginzburg-Landau theory to establish the phase diagram of a Co adatom lattice. With the aid of a small external magnetic field, a spin spiral and a skyrmion lattice can be stabilized, besides the ferromagnetic phase. The topological surface state under the influence of a spin spiral, a single skyrmion, and a skyrmion lattice is numerically solved. Chiral conducting modes are found on the domain walls with zero out-ofplane magnetic moment in a spin spiral. Similar chiral modes are also found on the boundary of a single skyrmion. For a skyrmion lattice, a gap opens at the surface Dirac point, leading to the anomalous Hall effect.
Note added. Recently, we become aware of two recent works by Paul and Fu [54] and by Divic et al. [44], in which some similar results are obtained.

ACKNOWLEDGMENTS
We thank Prof.
Yayu Wang for comments on the manuscript.
The first-principles calculations for this work were performed on TianHe-2.
We are thankful for the support from the National Supercomputing Center in Guangzhou (NSCC-GZ). In this appendix we describe the model geometry and technical details of the first-principles calculations carried out in this paper.
For the QL-stacked Bi 2 Se 3 , we consider one QL slab with experimental lattice parameters [55], and build a slab model of a 20-Å-thick vacuum layer to investigate the adsorption of cobalt atoms on the surface and magnetic interactions between them. The first-principles calculations based on density functional theory (DFT) [56,57] are performed with plane wave basis sets and pseudopotential method, as implemented in the quantum espresso package [58,59]. Perdew-Zunger parametrization of the local density approximation [60] is employed for the exchange-correlation functional. The projector augmented-wave [61] pseudopotentials in the pslibrary [62, 63] (version 1.0.0) are adopted. The energy cutoff of the plane wave basis set is chosen to be 70 Ry. The position of the Co adatom is relaxed by the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm until the Hellmann-Feynman force is less than 0.001 Ry/bohr. The Brillouin zone is sampled by a 3 × 3 × 1 grid in structure relaxations and self-consistent charge density calculations. In magnetic interaction calculations, a 5 × 5 × 1 grid is employed. Spin-orbit coupling (SOC) is included in all the energy calculations unless explicitly stated otherwise.
Three typical adsorption positions are calculated, and their energies are compared. As marked by colored crosses in Fig. 1(b), they are above a hollow (site A), above a Bi atom (site B), and above a Se atom (site C). The equilibrium adsorption positions of these sites have different heights relative to the surface Se atomic layer. Specifically, the site A is 0.24Å below the Se layer, while the site B and the site C are 0.65 and 2.12Å above the Se layer, respectively. The adsorption energies on these sites are further computed using a 3 × 3 supercell. It is found that the most stable adsorption site is site A, with energy 0.71 and 3.14 eV lower than the energies of sites B and C, respectively.

Appendix B: The determination of spin interaction parameters
To calculate the magnetic interactions between two Co adatoms on neighboring A sites, we use a 3 × 3 supercell as shown in Fig. 1(c). In addition, we use the energy-mapping method [64,65] to determine the interaction parameters in Eq. (2) from DFT calculations. For example, we would like to calculate the J xy component. Then two magnetic configurations are constructed: (1) S 1 = (S , 0, 0), S 2 = (0, S , 0); and (2) S 1 = (S , 0, 0), S 2 = (0, −S , 0). We constrain the magnitude and direction of the magnetic moment of the two Co atoms and calculate the total energy of the system, E 1 and E 2 , by DFT. By the spin-spin interaction model in Eq. (1), where E 0 is the energy of other parts of the system. Then J xy can be obtained from E 1 and E 2 , A similar method can be used to calculate J yx ; then Γ xy = J xy + J yx /2 (which should be zero by symmetry) and D z = J xy − J yx /2 can be obtained.