Surface Magnetism in Topological Crystalline Insulators

TL;DR

This paper investigates how magnetic impurities induce surface ferromagnetism in topological crystalline insulators, revealing symmetry-dependent magnetic states and analyzing their stability and experimental implications.

Contribution

It introduces a model showing multiple degenerate magnetic ground states influenced by surface symmetry and electron density, with stability analysis of ferromagnetic phases.

Findings

Surface magnetism depends on surface crystalline symmetry.

Magnetic states can have twofold, sixfold, or eightfold degeneracy.

Ferromagnetic states are stable against thermal fluctuations.

Abstract

We study topological crystalline insulators doped with magnetic impurities, in which ferromagnetism at the surface lowers the electronic energy by spontaneous breaking of a crystalline symmetry. The number of energetically equivalent ground states is sensitive to the crystalline symmetry of the surface, as well as the precise density of electrons at the surface. We show that for a SnTe model in the topological state, magnetic states can have twofold symmetry, sixfold symmetry, or eightfold degenerate minima. We compute spin stiffnesses within the model to demonstrate the stability of ferromagnetic states, and consider their ramifications for thermal disordering. Possible experimental consequences of the surface magnetism are discussed.

Equations

$${\vec{S}}^{(a)}=\frac{1}{4}\left(u_a^2{\tilde{\sigma }}_2,u_a^2{\tilde{\sigma }}_1,(u_a^2-v_a^2){\tilde{\sigma }}_3\right).$$

$$H_i\approx E_i+{\alpha }_i(q_2-b_2){\tilde{\sigma }}_1+{\beta }_i(q_1-b_1){\tilde{\sigma }}_2+{\Delta }_i{\tilde{\sigma }}_3,$$

$$\frac{\delta E(Q)-\delta E(0)}{\mathcal{S}}=\frac{1}{2}\sum _{\mu ,\nu =x,y}{\rho }_{\mu ,\nu }{Q}_{\mu }{Q}_{\nu },$$

$$H_{\bar{X}}=v_1{q}_1{\tilde{\sigma }}_1+v_2{q}_2{\tilde{\sigma }}_2+E_{\bar{X}}+{\eta }_1{\tilde{\mu }}_x{\tilde{\sigma }}_1+{\eta }_2{\tilde{\mu }}_y,$$

$$\frac{q_{\pm }^2{\Delta }^2}{4}=[\tilde{D}(h_x+h_y)-\sqrt{2}\tilde{C}(h_x-h_y)q_{\pm }{]}^2$$

$$+2[\tilde{A}{\eta }_2-\tilde{B}{\eta }_1{]}^2(h_x+h_y{)}^2+[\tilde{D}{\eta }_1+u_a^2{\eta }_2/\sqrt{6}{]}^2{h}_z^2,$$

Full text

Introduction – Topological crystalline insulators (TCI’s) are a class of materials in which the energy bands can host non-trivial topology protected by a crystalline symmetry Fu (2011). These systems support surface states Liu et al. (2013) which remain gapless provided the crystal symmetry is unbroken, and are believed to present themselves in (Sn,Pb)Te and related alloys Hsieh et al. (2012); Tanaka et al. (2012); Xu et al. (2012); Dziawa et al. (2012); Okada et al. (2013); Yan et al. (2014). Interesting effects may arise when the symmetry protecting a topological band structure is broken. In topological insulators protected by time-reversal symmetry (TRS), magnetic impurities on a surface break this symmetry and form collective states Liu et al. (2009); Biswas and Balatsky (2010); Garate and Franz (2010); Abanin and Pesin (2011); Liu et al. (2016a), which may be understood in terms of a gap opening in the surface spectrum Efimkin and Galitski (2014).

In contrast, TCI’s are not protected by TRS, so the loss of this symmetry does not by itself energetically favor ordering of magnetic moments Shen and Cha (2014); Assaf et al. (2015). However, a uniform magnetization can undermine one or more relevant crystalline symmetries Zhang et al. ; Fang et al. (2014). Indeed, the most common such symmetry is reflection across a mirror plane, of which there can be several. We show below that spontaneous surface magnetization opens a maximal gap when oriented along axes dictated by the bulk  symmetries of the system. For a generic surface with a single mirror plane, there are two surface Dirac points at different momenta and energies Ando and Fu (2015), and in such cases at low temperature this results in a metallic, Ising-like ferromagnet, with the easy axis determined by the chemical potential $\mu$. Importantly, the number of degenerate low-energy directions is enhanced for surfaces with further symmetries. Rotational symmetries in particular yield multiple mirror planes, and connect distinct surface Dirac cones to one another, yielding a multiplicity of easy axis directions. For sufficiently high symmetry, all the surface Dirac points may be related by symmetry operations, resulting in a fully gapped surface spectrum and a large number of groundstate orientations.

To illustrate this physics, we present detailed calculations for the (111) surface of (Sn,Pb)Te Li et al. (2013); Yan et al. (2014); Taskin et al. (2014); Shen et al. (2014), using a known model Hamiltonian Littlewood et al. (2010); Hsieh et al. (2012). The (111) surface states are characterized in this system by four surface Dirac points, one at the $\bar{\Gamma}$ point and one at each of three $\bar{M}$ points Hsieh et al. (2012) [Fig. 1(a).] When the system is doped by substitutional isolectronic magnetic impurities, with $\mu$ adjusted to the bulk gap, the absence of free carriers in the volume suggests that bulk magnetism will not occur. However, it can be stabilized on the surface when it opens a gap in the surface spectrum. Moreover, the magnitude of the various gaps are sensitive to the direction of the magnetization [Fig. 1(b)]. Since the gap centers are at different energies, favored magnetization directions are determined by $\mu$ [Figs. 1 (c) and (d)]. When in the vicinity of the $\bar{\Gamma}$ Dirac point energy, there is a single easy axis, yielding an Ising ferromagnet. With $\mu$ near the $\bar{M}$ Dirac point energies, because these points are connected by a three-fold rotational symmetry, there is a six-fold degenerate set of groundstate orientations. Remarkably, the sensitivity to $\mu$ implies that the low-energy orientations can be controlled externally by a gate.

The (111) surface is an example of how symmetry leads to a multiplicity of magnetic groundstates. This also occurs at a (001) surface, where four surface Dirac cones are supported by two distinct mirror planes. As explained below, the high symmetry leads to a potential 8-fold degeneracy of groundstate magnetization directions. Moreover, the energetic coincidence of all the surface Dirac points allows for the surface states to be fully gapped, yielding insulating behavior for a range of $\mu$.

The existence of the degenerate magnetic groundstates should be detectable via the behavior of their domain wall (DW) excitations, which proliferate at thermal disordering transitions, or can be frozen in when the system is zero-field cooled. Furthermore, because DW’s connect regions with different Chern numbers, they necessarily support bound conducting states Jackiw and Rebbi (1976). Their energetics can also behave rather differently depending upon the placement of $\mu$ relative to the surface bands. We discuss possible effects of the DW’s below.

Bulk Hamiltonian and Surface States – Our analysis employs a tight-binding Hamiltonian $H_{bulk}$ for materials in the (Sn,Pb)Te class, which is a rocksalt structure (fcc lattice). $H_{bulk}$ involves twelve orbitals: for each spin there are $p_{x},\,p_{y},\,p_{z}$ states on each of two sublattices, labeled $a$ and $b$, with on-site energies $m_{a,b}$ (see Supplementary Material Sup ). The model represents a direct gap semiconductor with smallest gaps at the $L$ points [${\bf k}={\bf k}_{1},{\bf k}_{2},{\bf k}_{3},{\bf k}_{4}\equiv(\frac{\pi}{2},\frac{\pi}{2},\frac{\pi}{2})$, $(-\frac{\pi}{2},\frac{\pi}{2},\frac{\pi}{2})$, $(\frac{\pi}{2},-\frac{\pi}{2},\frac{\pi}{2})$, $(\frac{\pi}{2},\frac{\pi}{2},-\frac{\pi}{2})$ in units of the inverse nearest neighbor separation]. For $\bf k$ precisely at an $L$-point states have well-defined sublattice index (with on-site energies $m_{a}$ and $m_{b}$). Adjusting $m_{b}-m_{a}$ to an appropriate value $m_{0}$ brings $a$ and $b$ states into energetic coincidence, forming the basis of a Dirac point at the Fermi energy.

When $m_{b}-m_{a}=m_{0}+m$ with $m<0$, there is a band inversion and associated nontrivial band topology Hsieh et al. (2012), protected in this system by mirror symmetries, so that surfaces respecting any of them support gapless states Liu et al. (2013). Low energy forms of these may be constructed Liu et al. (2010); Silvestrov et al. (2012); Brey and Fertig (2014), as we describe for the specific case of the (111) surface in the Supplementary Material Sup . For this surface, there are Dirac points residing at $\bar{\Gamma}$ and each of the three $\bar{M}$ points – as illustrated in Fig. 1 (a) – with energies $E_{\bar{\Gamma}}$ and $E_{\bar{M}}$ respectively. Note that in this situation the system has a three-fold rotational symmetry, which maps the two degenerate states at $E_{\bar{\Gamma}}$ onto one another, forming a two-dimensional representation of this rotation group. States at the three ${\bar{M}}$ points form a six dimensional representation. For each Dirac point, approximate explicit forms of the the wavefunctions may be constructed Sup , which can be written as eigenstates of a mirror operator $\tilde{\sigma}_{1}$ with eigenvalues $\pm 1$. Projection of $H_{bulk}$ onto these surface states allows us to construct effective Hamiltonians in the vicinity of each Dirac point.

Magnetic Impurities and Surface Hamiltonians – It has long been known that metals in the (Sn/Pb)Te class Inoue et al. (1975, 1977, 1979); Story et al. (1986); Karczewski et al. (1992); Geist et al. (1996, 1997); Prinz et al. (1999); Łusakowski et al. (2002), may be doped with magnetic ions which in some circumstances order ferromagnetically at low temperature. In these systems the magnetic ions enter substitutionally for Sn/Pb atoms, and the coupling of the magnetic moments with the conduction electrons can be understood rather well using an $s-d$ model Dietl et al. (1994), $H_{sd}=J\sum_{i}\vec{S}({\bf r}_{i})\cdot\vec{s}_{i}$, where $\vec{s}_{i}$ represents an impurity spin at location ${\bf r}_{i}$ and $\vec{S}({\bf r}_{i})$ is the conduction electron spin density Liu et al. (2016b). We consider the situation where the chemical potential is in a gap of the bulk spectrum, so that free carriers are not present and bulk magnetic ordering is not expected. In the TCI state, however, surface electrons couple the magnetic moments of the substitutional impurities near the surface, and may lead to ferromagnetism Rosenberg and Franz (2012); Lasia and Brey (2012). We model this by assuming magnetic impurities are present in the system, on one sublattice, near the surface.

The explicit surface wavefunctions allow us to project the electron spin operators onto surface states for the ${\bar{\Gamma}}$ and $\bar{M}$ points. As discussed in the Supplementary Material Sup , the spin operators on a single (say, the $a$) sublattice for either the $\bar{\Gamma}$ or an $\bar{M}$ point may be written

[TABLE]

In this expression, $\tilde{\sigma}_{3}$ is $2\times 2$ matrix whose eigenvectors yield the combinations of $\tilde{\sigma}_{1}$ eigenstates with well-defined eigenvalues under $2\pi/3$ rotations around a bulk $\Gamma-L$ direction Sup , and $\tilde{\sigma}_{2}=-i\tilde{\sigma}_{3}\tilde{\sigma}_{1}$. The quantities $u_{a}$, $v_{a}$ are real coefficients involving the tight-binding parameters Sup . We assume the impurity spins ferromagnetically order and treat the Hamiltonian in mean-field theory; the linear stability of the state against formation of a spin-density wave can then be checked. Projecting $H_{sd}$ onto the subspace of surface states for a Dirac point using Eq. 1 leads to an effective Hamiltonian of the form

[TABLE]

where $i$ denotes either $\bar{\Gamma}$ or one of the $\bar{M}$ points, and $q_{1,2}$ represent wavevector components along the surface. (Note the relationships between $(q_{1},q_{2})$ and $(q_{x},q_{y},q_{z})$ depend on the specific Dirac point $i$.) As expected on general symmetry grounds, $\alpha_{\bar{\Gamma}}=\beta_{\bar{\Gamma}}$, but $\alpha_{\bar{M}}\neq\beta_{\bar{M}}$. The offsets $b_{1}$ and $b_{2}$ are proportional to components of the impurity magnetization perpendicular to ${\bf k}_{i}$, while $\Delta_{i}$ is proportional to the component along it. The resulting spectra, $\varepsilon_{i}=E_{i}\pm\sqrt{\alpha_{i}^{2}(q_{1}-b_{1})^{2}+\beta_{i}^{2}(q_{2}-b_{2})^{2}+\Delta_{i}^{2}},$ provides the important observation that when the moments align along a $\Gamma-L$ direction and $\mu\sim E_{i}$, a gap opens in the corresponding surface spectrum that lowers its contribution to the total electron energy Efimkin and Galitski (2014).

Numerical Studies – To test this idea we have numerically computed the electronic energy of a TCI slab with open (111) surfaces, using the tight-binding model $H_{bulk}$ Sup , and adding an effective magnetic field $\vec{b}$ near the surface on only the $a$ sublattice, in such a way that their coupling to the two states associated with the surface Dirac cones is the same. Our tight-binding parameters are adapted from Ref. Fulga et al., 2016, and we have verified the presence of four surface Dirac points. (Some technical details are provided in the Supplementary Material Sup .)

Initially we consider a slab with primitive unit cell presenting only a single site of one sublattice on the surface, and introduce a surface magnetization as described above. While this represents a relatively large density of impurities (relative to experiment), it captures the correct qualitative physics, and allows us to study a wide enough slab that the surfaces are effectively decoupled. As expected from the above discussion, among the four Dirac cones the one with largest magnetization projection along its corresponding $\Gamma$ - $L$ direction develops the largest gap. Figs. 1 (c) and (d) illustrate the resulting total electronic energy, for $\mu\sim E_{\bar{M}}$ in (c) and $\mu\sim E_{\bar{\Gamma}}$ in (d). The energy is minimized in the former case for $\vec{s}$ along a $\Gamma-L$ direction associated with an $\bar{M}$ point, while in the latter minimization occurs for $\vec{s}$ along ${\bf k}_{1}$. This leads to two degenerate minima for $\mu\sim E_{\bar{\Gamma}}$ and six for $\mu\sim E_{\bar{M}}$. (Only half of these can be seen in the figures.) Analogous results are found when $\mu$ rather than particle number is fixed. As expected, the number of minimal energy states reflects the symmetry of the surface.

To further substantiate this, we also studied a more dilute magnetic moment model, in which the impurities are present only for $2/9$ of the atoms of one sublattice near the surface. Fig. 2 shows the difference in Gibbs free energy of the system ($\langle H_{bulk}\rangle-\mu N$ with $N$ the number of electrons) when the magnetic moments are oriented in the $(111)$ direction and in the $(11\bar{1})$ direction, as a function of chemical potential. The results again demonstrate that energetically favored directions are determined by $\mu$. We have also used this geometry to verify that orienting the two surface magnetic moments in different directions always raises the energy of the system, supporting our assumed ferromagnetic ordering, and that different placements of the impurities in the unit cell on one of the sublattices has little effect Sup . The latter suggests that disorder in the impurity location has little impact on our results Liu et al. (2009). Finally we note that density of states resonances Wehling et al. (2014) that appear from vacancies and impurities in, for example, graphene, do not appear to be present for the kind of disorder we consider.

The ferromagnetic ordering may be further substantiated by considering what happens to the electronic energy when the effective field is allowed to vary spatially with some wavevector $\vec{Q}$ along an average direction, $b_{1,2}({\bf r})=b_{1,2}^{(0)}+\delta b_{1,2}\cos\left({\bf Q}\cdot{\bf r}\right)$, $\Delta({\bf r})\equiv\Delta^{(0)}+\delta\Delta\cos\left({\bf Q}\cdot{\bf r}\right)$. To compute this we adopt as our basic Hamiltonian Eq. 2, assuming for simplicity $\alpha_{i}=\beta_{i}\equiv\alpha$, use the directions associated with $q_{1}$ and $q_{2}$ to define $x$ and $y$ directions on the surface, and compute the change in energy to second order in $\delta b_{1,2}$, $\delta\Delta$, and $Q$. Interestingly, for the valleys in which there are Fermi surfaces, this turns out to be independent of $Q$, as is the case for graphene Brey et al. (2007). Thus the spin stiffness (quadratic dependence of the energy correction on $Q$) comes from any valley(s) for which the Fermi energy passes through a gap. After an involved calculation (see Supplementary Material Sup ), one finds a correction of the form

[TABLE]

where the coefficients $\rho_{\mu,\nu}$ are all second order in the deviations $\delta b_{1,2}$, $\delta\Delta$, and the eigenvalues of the $2\times 2$ matrix it represents are positive Not . This demonstrates that if the effective field from the surface magnetization has a spatial oscillation, the resulting energy increases with increasing oscillation wavevector, as should be for a ferromagnetically aligned groundstate. Note that the stiffnesses $\rho_{\mu,\nu}$ all diverge as $1/\Delta^{(0)}$ Sup , a property with interesting consequences to which we will return.

(001) Surface – The (001) surface of SnTe supports gapless states that differ qualitatively from those of the (111) surface in that they involve valley admixed states Liu et al. (2013) in the vicinity of an $\bar{X}$ point of the surface Brillouin zone. The effective Hamiltonian for such surface states may be written in the form

[TABLE]

where $q_{1}=(q_{x}+q_{y})/\sqrt{2}$, $q_{2}=(q_{x}-q_{y})/\sqrt{2}$, $\tilde{\mu}_{x,y}$ are a Pauli matrices acting on a two-fold valley space, and $\eta_{1}$, $\eta_{2}$ are phenomenological valley-mixing parameters Liu et al. (2013). This expression may be arrived at by explicit projection of two valleys (e.g., ${\bf k}_{1}$ and ${\bf k}_{4}$) onto (001) surface states using the approach described above. Dirac points for this Hamiltonian lie at the momenta ${q}_{1}=\pm(\eta_{1}^{2}+\eta_{2}^{2})^{1/2}\equiv q_{\pm}$, $q_{2}=0$. Explicit forms for the two zero energy surface states at each of these points may be obtained using the same approach as for the (111) surface, and effective electron spin operators on the $j=a$ sublattice derived which couple to magnetic impurities. Assuming these are ferromagnetically aligned, in mean-field theory they add a term of the form $\vec{h}\cdot\vec{S}^{(a)}$, where $\vec{h}$ is the average magnetization, and induce a gap $\Delta$ at the Dirac points satisfying

[TABLE]

where $\tilde{A}=(u_{a}^{2}-v_{a}^{2})/12\sqrt{6}+4u_{a}v_{a}/\sqrt{12}$, $\tilde{B}=(u_{a}^{2}+2v_{a}^{2}-4\sqrt{2}u_{a}v_{a}$, $\tilde{C}=u_{a}^{2}/2-v_{a}^{2}$, and $\tilde{D}=(v_{a}^{2}-u_{a}^{2})/\sqrt{12}$. For fixed magnitude of $\vec{h}$, it is easy to see that there are two oppositely oriented directions that maximize $\Delta$; moreover, these directions are different for $q_{+}$ and $q_{-}$. Together with the four-fold symmetry that guarantees equivalent behavior for the coupled $\bf{k}_{2}-\bf{k}_{3}$ valleys, we conclude that there are eight minimal energy directions for $\vec{h}$ for the (001) surface when the Fermi energy is in the vicinity of $E_{\bar{X}}$. Furthermore, because all the Dirac cones are centered at this energy, we expect the gaps will generally overlap, so that the chemical potential may pass through all of them simultaneously. The spontaneously magnetized (001) surface then allows for insulating electronic behavior, in contrast to the (111) surface which for dilute magnetic impurities remains metallic.

Discussion and Speculations – We next consider some physical consequences of the surface magnetism discussed above, focusing on temperature ranges where the impurity magnetic moments may be treated classically. As mentioned above, the sensitivity of the magnetization directions for energy minima to $\mu$ should allow it to be controlled via a gate potential, which in principle would be observable in direct magnetization measurements. Another basic observation is that the gap openings induce a Berry’s curvature in the surface bands, which generically induces an anomalous Hall effect. For the SnTe system with (111) surfaces we do not expect it to be quantized Zhang et al. ; Fang et al. (2014), since the chemical potential typically cannot pass through a gap for all the surface Dirac species at the same time. By contrast, for (001) surfaces where the Dirac cones points are initially all at the same energy, the surface magnetization may allow all the induced gaps to overlap. With chemical potential in this gap, one does expect a quantized anomalous Hall effect Zhang et al. ; Fang et al. (2014).

It is interesting to consider possible consequences of $\rho_{\mu\nu}\sim 1/\Delta^{(0)}$ as discussed above. In particular we expect that the multiple minima presented in Fig. 1 imply that there should be DW excitations in the system, with energy per unit length scaling as $\sqrt{\Delta^{(0)}\rho_{0}}$, with $\rho_{0}$ an appropriate average of $\rho_{\mu\nu}$’s. This remains finite even as $\Delta^{(0)}$ vanishes, as should happen at high enough temperature. The divergence of the spin stiffness, $\rho_{0}\sim 1/\Delta^{(0)}$, reflects the fact that as the gap vanishes, the quantity $\delta E({\bf Q})$ is no longer analytic in ${\bf Q}$, and in particular rises linearly with $Q$ in the long wavelength limit Brey et al. (2007), suggesting a non-local interaction among spin gradients. Presuming $\mu$ ends up at the Dirac point as $\Delta^{(0)}$ vanishes, the simplest model of the system is a clock model with long-range interactions, which in the Ising case would approach the transition with mean-field exponents Paulos et al. (2016). For other values of $\mu$ it is possible that DW’s with finite energy per unit length can be stabilized in this circumstance, but if so would not have the simplest structure Rajaraman (1989). In the case of the (111) surface, for $\mu\sim E_{\bar{\Gamma}}$ the system presumably will undergo a second order phase transition in the Ising universality class. By contrast, for $\mu\sim E_{\bar{M}}$ case with six different minima the system could be represented by a six-state clock model. For the (001) surface and $\mu\sim E_{\bar{X}}$, an eight-state clock model is relevant. In both these cases, presuming $\mu$ does not evolve precisely to a Dirac point as the magnetization disorders and the gaps close, the phase transitions should be in the Kosterlitz-Thouless universality class José et al. (1977).

We also note that DW excitations in this system may accumulate charge, both due to mid-gap states Jackiw and Rebbi (1976); Schaakel (2008), and from the surface valleys which have Fermi surfaces allowing low-energy scattering. The existence of mid-gap states is a necessary by-product of the topology of the gapped Dirac points described by Hamiltonians of the form in Eq. 2. The essential effect of the DW on the electrons is that $\Delta_{i}$ changes sign as one moves through it. This means the Chern number evolves from $-1/2$ on one side of the DW to $1/2$ on the other, introducing gapless states bound to the DW interior Fang et al. (2014). At the critical temperature $T_{c}$ where the transition occurs, one expects DW’s to proliferate, opening a channel for conduction which is absent below $T_{c}$. This could lead to singular behavior (e.g., a cusp) in the conductivity at the transition Jungwirth and MacDonald (2001); Dhochak et al. (2015), and should also have a signature when the surface is probed via tunneling. A further possibility to probe the physics is by looking for differences in conductivity between field-cooling and zero-field cooling of the system through its critical temperature. The latter leads to nucleation of groundstate domains with random orientation, and DW’s between them, which cannot relax on the time scale of an experiment. Thus one expects stronger surface conduction from a zero field-cooled sample Assaf et al. (2015); Ueda et al. (2015); Tian et al. (2016). Finally, the presence of charged DW’s on the surface might be detected directly via coupling to electro-magnetic waves Ma et al. (2015), whose scattering should be sensitive to the proliferation of DW’s.

In summary, the surface of a magnetically-doped TCI hosts magnetic ordering in the topological state even when the bulk is disordered. The unique electronic structure of a TCI surface leads to a richer set of possible ordered states than would be expected from time-reversal symmetry protected topological insulators, and implies a number of unusual physical behaviors.

Acknowledgements – The authors thank Fernando de Juan for helful comments, and U. Nitzsche for technical assistance. This work was supported by the NSF through Grant Nos. DMR-1506263 and DMR-1506460, by the US-Israel Binational Science Foundation, and by MEyC-Spain under grant FIS2015-64654-P. Computations were carried out on the ITF/IFW and IU Karst clusters. HAF also thanks the Aspen Center for Physics where some of this work was performed.