Highly tunable magnetic coupling in ultrathin topological insulator films due to impurity resonances

We theoretically investigate the exchange interaction between magnetic impurities in ultrathin Bi$_2$Se$_3$ topological insulator films by taking into account the low-energy states produced by the impurities. We find that the locally induced impurity resonances strongly influence the exchange interaction between magnetic moments. In particular, we find a non-collinear alignment being more favorable than the collinear ferromagnetic alignment preferred when impurity states are ignored and only the pristine topological insulator band structure is considered. Moreover, we show that by applying of an electric field perpendicular to the ultrathin film, the exchange interaction can be drastically enhanced. This opens for the possibility of highly tunable magnetism by electric field.


I. INTRODUCTION
Magnetic impurity moments in solids couple indirectly to each other through the electronic carriers of the host material, with the nature of the coupling mechanism determined by the electronic structure of the host. In terms of the spatial extension, magnetic moments couple at distances up to a few nanometers in topological insulators, whereas the corresponding coupling in most semiconducting materials ranges only at most a few Ångströms 1 . This long-range coupling between the magnetic impurities in topological insulators can lead to the existence of a robust magnetic phase, which has also been recently confirmed in experiments, both using angular resolved photoemission spectroscopy (ARPES) and scanning tunneling microscopy (STM) [2][3][4] .
It has been suggested that the indirect magnetic coupling in topological insulators is governed by topology 2 . While the magnetically ordered phase in topologically trivial materials is fragile to decreasing dimensionality, the magnetic phase in topological insulators appears to be strengthened when lowering the dimension. Indeed, it has been observed that the magnetic phase generated by coupled magnetic impurities on the surface of topological insulators is far more robust than the corresponding phase created between bulk impurities 5 . In line with this, magnetic coupling of impurities has also been verified in ultrathin films of topological insulators in experiments measuring the quantum anomalous Hall effect (QAHE) [6][7][8][9][10] . Moreover, in topological insulators thinner than five quintuple layers, the surface states become gapped due to hybridization between the surface states at opposite sides 11 . Thus ultrathin films also provide an opportunity to combine the exotic properties of topological insulators and a finite energy gap into a single unit. In addition, with the band dispersion of ultrathin topological insulator films tunable by the application of an electric field, a tantalizing electric control of magnetism 12 might also be achievable in the ultrathin limit.
Very generally, the indirect magnetic coupling between magnetic impurities is governed by the magnetic susceptibility of the whole system. Thus, in materials with strong spin-orbit interaction, such as topological insulators, the magnetic susceptibility tensor gives three fundamentally different contributions: an isotropic (Heisenberg-like), a symmetrically anisotropic (Ising-like), and an asymmetrically anisotropic (Dzyalosinskii-Moryia-like) contribution [13][14][15] . The isotropic and symmetrically anisotropic contributions lead to collinear configurations of the magnetic moments, whereas the asymmetrically anisotropic contribution favors non-collinear configurations 16 . This asymmetric anisotropy easily leads to the existence of exotic magnetic phases, such as spin-glass, or chiral ferro-and anti-ferromagnetic phases 17,18 .
The nature of the indirect coupling between magnetic impurities is also significantly different in metallic and semiconducting materials. While the interaction in metals is dominated by the excitations around the Fermi level, or Ruderman-Kittel-Kausya-Yosida (RKKY) interactions, the interactions in semiconductors is mainly provided by the interband susceptibility between the valence and conduction band electrons, more related to van Vleck-type magnetism. The latter type based on interband susceptibility has in fact already been observed in topological insulator thin films [19][20][21] .
With the indirect coupling determined by the magnetic susceptibility of the whole system, the impurities themselves can in fact modify the indirect coupling 16,22,23 . In particular, in topological insulators, impurity resonance states easily appear near the Dirac point due to scattering off the non-magnetic part of the impurity potential [24][25][26] . In ultrathin topological insulator films these resonances may even arise inside the band gap, depending on the strength of the impurity scattering potentials and the properties of the topological insulator 24,27 , also recently addressed using ab-initio calculations 19,23,28 . The gapped dispersion in ultrathin films thus makes impurity states even more prominent as they can be the only low-energy states available to mediate the indirect magnetic coupling. Despite the ubiquitousness of low-energy impurity states in topological insulators, most calculations of the magnetic susceptibility in topological insulator systems have been based only on the itinerant electrons of the unperturbed, or pristine, topological insulator surface states 21,29 , with only recent work highlighting the effect of the impurities, and then only in thick topological insulators 16 . This is despite the fact that recent experiments on ultrathin topological insulator films have indicated that the indirect magnetic coupling depends strongly on the nature of the magnetic impurities 1,2,6-9 .
In this work we therefore calculate the influence on the magnetic susceptibility of the impurity induced in-gap resonances in ultrathin topological insulator films, to accurately capture the indirect coupling between surface magnetic impurity moments. We show that the impurity resonances arXiv:2003.05539v1 [cond-mat.mtrl-sci] 11 Mar 2020 strongly enhance both the isotropic Heisenberg and symmetrically anisotropic Ising contributions. We also find a strong energy dependence, which opens up for excellent tuning possibilities of the indirect coupling on-and off-resonance. Most importantly, we find that the asymmetrically anisotropic Dzyalosinskii-Moryia contribution is generated when we include the impurity states. As a direct consequence, a noncollinear collective configuration of the magnetic moments easily becomes favoured, with a finite out-of-plane net magnetization. This is in sharp contrast to the ferromagnetic ground state previously obtained when only considering the pristine ultrathin topological insulator films 21 . Finally, we also show that by simply applying an external electric field, the magnetic coupling becomes extensively tunable, between ferroand anti-ferromagnetic to chiral configurations.
The remaining of the article is organized as follow. In Sec. II, we introduce the model Hamiltonian and the general formalism for calculating the magnetic susceptibility, including the impurity states. In Sec. III we present our results, with a focus on the contribution from impurity resonances and their tunability. Finally, we summarize and offer a few concluding remarks in Sec. IV. Some detailed part of the calculations can be found in the appendices A-C.

II. MODEL AND METHOD
The low-energy properties of the surface state electrons in the ultrathin topological insulator films can be described by an effective two-dimensional Hamiltonian near the Γ−point describing the two surfaces Here, σ, τ are Pauli matrices in the spin and surface space, respectively, k = (k x , k y ) denotes the two-dimensional wave vector for the surface electrons, and v F is the Fermi velocity. Moreover, V denotes the potential difference between the two surfaces. This potential arises from the combined effect of a substrate and/or an external electric field applied perpendicular to the film. Due to the thikness of the film there is an effective mass hybridization term ∆ that couples the two surfaces of the topological insulator. The model in Eq. (1) leads to the dispersion relation of a gapped Dirac spectrum with gap size 2∆. Here, s = ±1, refers to the conduction and the valence bands, respectively, whereas m = 1, 2 labels the solutions. Further, we model the added magnetic impurities through the Hamiltonian where U and J c S i represent the spin-independent and spindependent part, respectively, of the scattering potential. Here we allow magnetic impurities on sites r i , summing over all impurities, where we restrict ourselves to consider impurities in the top surface only, without limiting the applicability of our results. The only assumption in Eq. (3) is that the impurities behave as classical spins with strength |S i | = S , as appropriate for higher spin impurities. For two local impurity magnetic moments, located at r and r , respectively, the effective indirect coupling, or exchange, Hamiltonian can be written as where χ(r, r ) is the magnetic susceptibility tensor. Using the results in Ref. 16 , we can always write χ(r, r ) as In this expression, G 0 (r, r ; ω) denotes the bare single electron Green's function, i.e. without impurities, the trace runs over the spin degrees of freedom, g(ω) = tr G 0 (k, ω) dk/2, and where the expression in the denominator encode for the contributions from the impurity states. The magnetic susceptibility of the pristine topological insulator, χ 0 , is simply retained by setting U = J c S = 0 in the denominator. We refer to Appendix A for more details on the bare Green's function.
It is convenient to rotate the spin vectors 16,29 in terms of the polar angle ϕ R of the relative distance between the impurities. Then the exchange Hamiltonian takes the form Here H and I refer to the isotropical and symmetrical anisotropic couplings, respectively, whereas D = D (1, −1, 0) denotes the asymmetrical anisotropy. Here H, I, and D can be thought of as Heisenberg-, Ising-, and Dzyaloshinskii-Moriya-like interactions, respectively. Using the equation for the magnetic susceptibility given in Eq. (5) in Eq. (4), these interaction parameters are obtained from the expressions 16 Here G tt and G tt indicate the ↑↑ and ↑↓ components, respectively, of the bare Green's function on the top surface, see Appendix A for more details. Again, the results for a pristine topological insulator, J 0 = H 0 , I 0 , D 0 , is obtained by setting J c S = U = 0 in the denominator of Eqs. (7). In order to provide a physical understanding for the behavior of the exchange parameters in the presence of impurities, we also consider the spin-polarized local density of states (spin-resolved LDOS) ρ ↑,↓ , see Appendix B for a detailed expression derived from the full Green's function.
Before proceeding we also note that the magnetic susceptibility is, in general, a sum of interband and intraband contributions, χ = χ intra + χ inter . At zero temperature and finite occupancy at the Fermi level, such as in a metal, the intraband contribution normally dominates the exchange interactions. Then, the indirect magnetic coupling can be well described in terms of the itinerant electrons of the Fermi surface, giving rise to a (generalized) RKKY interaction. By contrast, in gapped systems, such as in pristine utrathin topological insulator films, the intraband contribution vanishes whenever the chemical potential lies within the band gap, due to the absence of carriers. Hence, only interband contributions are present, although they are also small 29 . The part of the response function which originates from such an interband process between conduction and valence bands is usually referred to as a van Vleck interaction. 19

III. RESULTS
Using Eqs. (7) we quantify the complete and general indirect exchange coupling in terms of J = H, I, and D between two magnetic impurities on the top surface of an ultrathin topological insulator film. We first present the different components in an unbiased (V = 0) film, and thereafter we discuss the effects of an external electric field. Building on these results, we proceed to study the magnetic ordering of impurity moments and the different magnetic phases in the system.
For the remainder of this work we assume a four quintuple layer thin film of Bi 2 Se 3 , which has an energy gap 2∆ = 70 meV and Fermi velocity v F = 4.48 × 10 5 m/s. We also assume that the chemical potential resides within the band gap, i.e. |µ| < ∆, in order to directly connect which recent experiments 10 . Moreover, we assume an inter-impurity distance of R = 10 Å in the surface plane, unless we explicitly investigate the distance behavior. We also require that 0 ≤ J c S ≤ 1 and 0 ≤ U ≤ 6 (both given in units ofhv F /k c , where k c is the band momentum cutoff) and keep J c S /U ≈ 1/8. We find these values by comparing our spin-resolved LDOS with the results of Refs. [30][31][32][33] . To simplify our plots we express the coupling terms, J, and spin-resolved LDOS in units of (J c /h 2 v 2 F Ω BZ ) 2 and 1/h 2 v 2 F Ω BZ , respectively, where Ω BZ is the area of the first Brillouin zone.

A. Impurity states and their exchange interaction contributions
We start by setting the parameters V = 0 and U = 4, in order to study how the impurities influence the exchange interaction. The plots in Fig. 1(a,b) show the corrections from the impurities, δJ = J − J 0 with J = H, D, I and J 0 = H 0 , I 0 , D 0 , as a function of the chemical potential µ within the energy gap for J c S = U/8 and J c S = U/4. As a reference we also plot the exchange interactions parameters for the pristine topological insulator, J 0 , with black lines. In Fig. 1(c,d) we show the corresponding spin-resolved LDOS as a function of energy ε at the Fermi level for the same values of J c and chemical po- tential µ = 0 as produced by one impurity measured at the position of the other impurity. First, we directly see that the magnetic coupling is independent on the chemical potential, within the gap, in the pristine case, which is expected since this coupling is of van Vleck nature and thus given by interband transitions only. These constant exchange interactions obtained for the impurity-free ultrathin films should be contrasted with the interactions in the presence of impurities, which acquire both significantly different values and a very strong energy dependence, here reflected in the variation as a function of the chemical potential. In fact, first, we observe that the correction δJ to the exchange interaction changes the overall amplitude of the exchange coupling. From this observation we conjecture that the carrier density redistributed from the valence band into the impurity resonances has a substantial influence on the magnetic susceptibility. This is expected since the exchange coupling depends on the density of occupied states, i.e. it is a Fermi sea property, and thus it is natural that the presence of impurity resonances contributes to the overall amplitude for a wide range of chemical potentials. In fact, we see that the isotropic and symmetrically anisotropic corrections δH and δI even obtain opposite signs compared to H 0 and I 0 . The asymmetric anisotropy has, on the other hand, significantly different behavior. Since the band dispersion of a pristine topological insulator film is electron-hole symmetric, the overlap between states in the conduction and valence bands, due to their different helicities, result in a vanishing asymmetric anisotropy, D 0 , for all chemical potential values inside the gap. But, in the presence of the impurity resonances, the electron-hole symmetry is broken and the asymmetric anisotropy becomes finite. Second, there is a very strong energy dependence in a specific range of the chemical potential. More specifically, whenever the chemical potential µ is positioned between the impurity resonances (see e.g. the range −20 µ −15 meV in Fig. 1(a,c)), the amplitudes of the magnetic coupling are strongly enhanced compared to when both impurity resonances lie on the same side of the chemical potential. This strong energy dependence can thus be directly traced back to the emergence of impurity resonances inside the band gap, as seen in Figs. 1(c,d). As a consequence of the finite J c S , the impurity resonances are spin split into two spin-polarized resonances, merging into one non-spin-polarized resonance in the limit J c S → 0 27 .
Considering the properties of the impurities, both the spindependent J c S and spin-independent U parts of the impurity potential can vary between the impurities. Therefore, in Fig. 2 we plot the exchange couplings J = H, I, D, as a function of the spin-independent impurity potential U, again keeping R = 10Å and V = 0 and fixing J c S = U/8, for two different chemical potentials (a) µ = −20 meV and (b) µ = 0 meV. In both cases we observe that the isotropic (symmetric anisotropic) exchange is positive (negative) for low scattering potentials, but then transitions to negative (positive) values before diminishing for large U. The asymmetric anisotropy, on the other hand, vanishes in the absence of a scattering potential, peaks at small values before slowly approaching zero for increasing scattering potentials. However, in Fig. 2(a) this overall smooth dependence on U is interrupted by a sharply defined region with larger values, in the range U ∼ 3 to U ∼ 4. These boundaries exactly mark the energies where, at least one of, the induced impurity resonances coincide with µ, in analogy with the sharp features in Fig. 1(a,b). The substantially increased values of the exchange coupling in this U- range thus originate from that the impurity resonances residing within the band gap. The absence of the corresponding features in Fig. 2(b) is due to the fact that for this range of potentials U, the impurity resonances reside in the valence bands only and thus their effect is not present for the choice of in-gap value of µ.

B. Electrical tunability of magnetism
Having seen a strong dependence for the exchange coupling on the chemical potential, we next turn to the possibility to easily tune this behavior by applying an external electric field perpendicular to the plane of the film. To explore the signature of such potential difference in the magnetic exchange interaction, we present H, I, D in terms of the parameter V for different values of the chemical potential in Fig. 3(a,b), while in Fig. 3(c,d) we present the relevant spinresolved LDOS extracted at energy ε = µ and plotted as a function of V. Clearly, we see how all the magnetic coupling terms are significantly higher for Vs between two distinct values. In-between these two V-values the exchange coupling is in fact very large, for instance in (b) the Heisenberg coupling increases by ∼ 32 times with respect to the unbiased case. For different chemical potentials the region of enhanced exchange coupling shifts, but it still exists equally prominently. The existence of a region with giant exchange couplings and its behavior with chemical potential and bias is explained by looking at the spin-resolved impurity resonance positions in Fig. 3(c,d). As been shown before in the Ref. 27 , the application of an electric potential alters the position of impurities resonance peaks inside the gap, which then also moves the corresponding region with giant exchange couplings. Thus we find a huge electric tunability, with an extreme sensitivity, of the exchange interactions in an ultrathin topological insulator film. It is worth mentioning here that qualitatively, all features exposed in Figs. 1-3 are also valid for other impurity distances.

C. Orientation of magnetic moments
Having shown how two have a highly unusual magnetic impurities mutual interaction and also with a large tunability, we next calculate the spin configuration of two magnetic moments. We continue to assume classical spins, which means that the Hamiltonian (6) can be rewritten as H ex =|S | 2 (H + I) cos θ 1 cos θ 2 + H cos ∆ϕ + cosφ 1 cosφ 2 sin θ 1 sin θ 2 + D sin θ 1 cos θ 2 cosφ 1 − sin θ 2 cos θ 1 cosφ 2 , (8) where θ 1,2 and ϕ 1,2 are the polar and azimuthal angles of the spin vectors, S 1,2 , respectively. Here, the azimuthal angles are considered with respect to ϕ R , withφ 1(2) = ϕ 1(2) − ϕ R and we also define ∆ϕ = ϕ 2 − ϕ 1 . Following straightforward calculations presented in Appendix C, we find that the minimum energy of two magnetic impurities coupled to each other is given by tan ∆θ = D/(H + J), where ∆θ = θ 2 − θ 1 and ϕ 1 = ϕ 2 = ϕ R . The misalignment between the impurities is thus described in terms of the phase ∆θ. This phase is finite whenever the asymmetric anisotropy is finite (D 0), but vanishes in its absence. By introducing new spin variablesS (θ) with |S | = |S|, azimuthal angle ϕ = ϕ R , and polar angle θ, the effective spin Hamiltonian can for this arrangement be written as H =S (0) ·S (∆θ) 14 .
To further investigate the spin configurations, we plot in Fig. 4 the relative polar angle ∆θ between two magnetic mo- ments, as a function of both the inter-impurity distance R and electric field V for both µ = −20 meV (a) and µ = 0 (b). In both cases the configuration tends towards becoming ferromagnetic, i.e. ∆θ = 0, when the distance between the two impurities diminishes, R → 0, for all values of V. At finite distances, however, the phase diagrams display a wide range of different configurations, spanning from ferromagnetic through non-collinear to anti-ferromagnetic configurations. In particular, at both chemical potentials distinctive regimes exist where the moments align toward an anti-ferromagnetic-like configuration, ∆θ > π/2, for all distances R > 10 Å: in (a) −10 < V < 0 meV and in (b) 35 < V < 50 meV. We trace these regimes directly trace back to the entry and exit of the chemical potential between the impurity resonances, see e.g. Fig. 3(a,b). Beyond this electric field regime we find that most of the phase space is consists of clearly non-collinear configurations, where the relative angle is both far from 0 and π, with the exact configuration determined by R. This is due to the large influence from the asymmetric anisotropy, which renders the collinear cases less favorable compared to a noncollinear arrangement.
Extending the results of Fig. 4 to a multi-impurity setup, we conclude that impurities favour pair configurations that can be represented by the angle ∆θ in the ρz-plane, whereρ defines the in-plane direction between impurities. Since the impurities are located in different directions,ρ, with respect to each other, the only common axis of all magnetic moments is along the z-axis. Hence, the resulting phase may most likely be ascribed a non-collinear ferromagnetic nature, but with generally a z-axis component, i.e. an out-of-plane common component.

IV. CONCLUDING REMARKS
In summary, we have investigated effects of impurity resonances on the magnetic exchange coupling between the magnetic moments located on the surface of ultrathin topological insulator films. We find that the contribution from the impurity resonances, typically, is of the same order, as the bare contribution originating from the unperturbed surface states, but become much larger under certain resonance conditions. We further analyze the importance of the impurity resonances on the magnetic interactions in terms of the isotropic, and symmetric and anti-symmetric anisotropy components. For a pristine surface, the first two components are finite whereas the last one vanishes identically. We find that the contribution from the impurity states on the symmetric anisotropy and isotropic components, which both leads to collinear alignment of the magnetic moments moments, are within the same order of magnitude as the corresponding contributions from the intrinsic electronic structure, however, with opposite signs. Overall, this has a tendency to lead to a weakened collinear coupling. Most importantly, the non-collinear asymmetric anisotropic interaction, which is zero for pristine films, acquires a large contribution from the impurity states and imply that the collective ground state of the magnetic impurities should be strongly non-collinear. Furthermore, we show that the applications of an electric field perpendicular to the ultrathin film can be used as a mechanism to shift the energy of the impurity resonances, which opens up the possibility to electrically tune the properties of the magnetic interactions. In effect, this mechanism should provide a tool for tuning between isotropic and anisotropic interactions, something which clearly has a great impact on the magnetic state.
Based on our findings we conclude that calculations of the magnetic exchange interactions, without considering effects originating from the impurities themselves are oversimplified 34 . Hence, there is a great risk of losing interesting and important features in the system. In particular, the exchange interaction based on the pristine topological insulator surface states misses the anti-symmetric anisotropic component, which inevitable leads to the prediction of a noncollinear ferromagnetic phase. In fact, the easy-axis of many magnetic impurities on the surface of topological insulators is in-plane 28 , where calculations using the pristine system give in-plane ferromagnetic phase, while some experiments have already shown a perpendicular magnetic phase 35,36 . However, taking into account the influence from the impurity states but merely including the isotropic and symmetric anisotropic components of the exchange coupling (H + I) does not provide a sufficient description, as it leads to an anti-ferromagnetic phase. Here, we have shown the importance of the antisymmetric anisotropy and the necessity to include it in calculations of the magnetic phase of ultrathin topological insulators films. It is this latter component term that leads to a chiral ferromagnetic phase. Such chiral magnetic phase will clearly affect previous theoretical studies on QAHE experiments 37 . The QAHE has previously been assumed to be proportional to the net magnetization in the system. However, it has more recently been shown that the QAHE persists also in chiral ferromagnet and chiral antiferromagnetic systems with zero net magnetization 38 , extending the effect to large parts of the phase diagram uncovered in this work. In this Appendix we provide the components of the bare Green's function, i.e. for the pristine ultrathin topological insulator film without any impurities, which is used in the main text in Eq. (5) and also the on-site bare Green's function g(ω). The Matsubara Green's function in reciprocal space is given by G 0 (k; ω) = [iω + µ − H 0 ] −1 , which can be transformed into real-space by taking the Fourier transformation as G 0 (R; ω) = 1 Ω BZ dk e i k·R G 0 (k; ω).
After some straightforward calculations, the general matrix form of the bare Green's function is a 4 × 4 matrix reading where ϕ R = arctan (R y /R x ) denotes the polar angle of the vector R = r − r between the two impurities. As we assume the impurities to be located on the top surface of the ultrathin topological insulator film, we need the upper right block of this matrix, in which G tt and G tt represent the ↑↑ and ↑↓ spin configurations of the bare Green's function, respectively, given by whereω = iω + µ andR = i R (V − s √ω 2 − ∆ 2 ) /hv F , and K 0,1 are the modified Bessel functions of the second kind. From these expressions we also find the analytic expression of the on-site Green's function (A4)