Topological surfaces of domain wall-decorated antiferromagnetic
topological insulator MnBi_2nTe_3n+1

Introduction. As a representative example of magnetic topological crystalline insulator,zhang2015 the antiferromagnetic (AFM) topological insulator (TI) is characterized by a $\mathbb{Z}_{2}$ invariant protected by a composite anti-unitary symmetry $S=\Theta T_{1/2}$, where $\Theta$ is time-reversal and $T_{1/2}$ is a half lattice translation. Topologically protected gapless states are expected on surfaces respecting the $S$ symmetry. Recently, a family of layered compounds with an intrinsic A-type interlayer AFM order, $\mathrm{MnBi_{2n}Te_{3n+1}}$, have been shown to be AFM TIsotrokov2019 , whose thin films can display quantum anomalous Hall and topological magnetoelectric effects. li2019 ; zhang2019b ; liu2020a ; deng2020a ; ge2020 The surface spectrum measurements have focused on the top surface(otrokov2019, ; vidal2019a, ; hao2019, ; swatek2020a, ) but the topological surface states on symmetry-preserving surfaces have so far evaded detection. Lattice imperfections, such as surfaces and defects, can violate the $S$ symmetry and lead to bandgap or localization of the surface states.ando2015 ; zhang2019b ; li2019 ; deng2020a Of particular interest for a layered AFM TI is the inevitable presence of antiferromagnetic domain walls comprised of a ferromagnetic bilayer, which generally breaks the $S$ symmetry. Then whether the topological surface states remain gapless and/or conductive in an AFM TI with symmetry-breaking domain walls requires elucidation.

In this paper, we investigate the surface spectrum and transport properties based on an electronic model for a domain wall-decorated AFM TI and found the presence of domain walls has an intriguing impact on the topological surface states. Remarkably, the topological surface gap is a function of a single non-statistical variable describing the overall magnetic order, which is termed the Ising moment. The surface spectrum can be gapless even when the requisite $S$ symmetry is broken by domain walls. In the thermodynamic limit, the Ising moment converges to special gapless values when magnetic fluctuation is bounded. Computed Lyapunov exponents of surface transfer matrices show that the topological surface states are generally localized, except with an emergent chiral symmetry a delocalization transition occurs at a single energy. With bounded magnetic fluctuation, the crossover to chiral symmetry is broadened and displays an unconventional critical scaling, providing experimentally accessible transport signatures that in return uncover the nature of antiferromagnetic domains walls.

Surface gap and Ising moment. We will suppose that our model system is topologically nontrivial according to the $S$ symmetry in a perfect AFM configuration, as shown in Fig.1(a). Although the presence of domain walls generally breaks the $S$ symmetry and nullifies the $S$-protected topology, in a supercell a composite symmetry $\bar{S}=\Theta\bar{T}_{1/2}$ may emerge. This is the case if domain walls of opposite magnetization alternate with equal spacings, which can be achieved by inserting magnetic monolayers into a perfect AFM supercell as shown in Fig.1(b). If this insertion process can be seen as a continuous deformation preserving $\bar{S}$ symmetry, the resultant domain wall superlattice is bound to have gapless surfaces. The more relevant and interesting question is what happens to the surface spectrum when domain walls are randomly placed with random magnetizations.

An effective tight-binding model is available for $\mathrm{MnBi_{2n}Te_{3n+1}}$, from regularizing the 4-band $k\cdot p$ model of the topological insulator Bi₂Se₃zhang2009 on a hexagonal lattice (Fig. 1(a)) and adding a layer-staggered Ising-type exchange field describing the AFM order.zhang2020 We generalize this model to include a layer-dependent Ising field in the Hamiltonian

where integer $\ell$ is the magnetic monolayer index and $\bm{q}=(k_{x},k_{y})$ is the in-plane wavevector. Bloch bases in each magnetic monolayer are used, and the field operator $\psi_{\ell}$ and matrices $J_{\bm{q}}$ and $D_{\bm{q}}(m_{l})$ have dimensions of 4. The intralayer exchange interaction is modulated by the Ising order parameter $m_{\ell}=\pm 1$. For instance, $m_{\ell}=(-1)^{\ell}$ for the perfect AFM ordering. Domain walls can be easily introduced by $m_{\ell}$. Descriptions of the parameters, electronic structure and symmetry of this model are provided in the Supplemental Materials (SM).suppl It suffices to mention here that the parameters are chosen to ensure the system is topological in presence of $\bar{S}$ symmetry.suppl

For an unbiased survey of the effect of domain wall decoration on surface states, a series of supercells made of $N=60$ layers with randomly placed $n$ domain walls are examined. Each configuration has a zero net magnetization, which requires an equal number of $\Uparrow$ and $\Downarrow$ domain walls ($n_{\Uparrow}=n_{\Downarrow}$). The bandgaps of (010) surface (perpendicular to $y$-axis) states at the Dirac point ($\Gamma$) are calculated using the iterative Green’s function methodsancho1985 . Remarkably, the computed surface gap is seen to be a function of a single variable, which we call the Ising moment

where the overbar stands for averaging over all magnetic layers ($\mathcal{I}$ is invariant if calculated in a larger periodic or origin-shifted supercell). This is exemplified in Figs.2(a) and (b) with the numbers of domain walls per supercell $n=4$ and $6$, respectively. The surface states turn out to be gapless at special $\mathcal{I}$ values

And the specific value of $\mathcal{I}^{*}$ depends on the permutation of domain wall magnetization in the supercell.

In particular, a configuration without the $\bar{S}$ symmetry can also be accompanied by gapless (010) surface states. Indeed, each $\mathcal{I}^{*}$ corresponds to multiple domain wall configurations modulo cyclic transformations. This can be understood from how $\mathcal{I}$ changes with domain wall migrations. An elementary domain wall migration is accomplished by a local transposition of a pair of magnetic layers, which can be either (A) $\uparrow\downarrow\,\mapsto\,\downarrow\uparrow$ or (B) $\downarrow\uparrow\,\mapsto\,\uparrow\downarrow$. (A)/(B) changes $\mathcal{I}$ by $\pm 2/N$ and thus $\mathcal{I}$ is unaltered when performing transpositions (A) and (B) on two separate domain walls. As shown in Fig.2(c), this makes either a pair of $\Uparrow$ and $\Downarrow$ travel in the same direction (I$\rightarrow$II), or a pair of like domain walls travel in opposite directions (II$\rightarrow$III); in both cases $\mathcal{I}$ is unchanged.

An pair of $\Uparrow$ and $\Downarrow$ separated only by AFM layers can be brought next to each other by a sequence of local transpositions, and then annihilated by an additional transposition. And all domain walls can be removed by consecutive pair annihilations, ending in the perfect AFM order. Keeping track of the changes in $\mathcal{I}$ in this process of approaching the perfect AFM order ($\mathcal{I}=0.5$) provides us with a formula for $\mathcal{I}$, up to modulo 1,

where $p_{i}-2$ is the number of magnetic layers between the $i$th and $(i+1)$th domain walls, and $\alpha_{i}$ is $n_{\Uparrow}-n_{\Downarrow}$ in from the $(i+1)$th to the $n$th domain walls. The second equality in Eq.(4) is derived from the observation that, given a domain wall permutation, $\mathcal{I}^{*}$ is attained at equal separations (if possible) i.e. $\delta p_{i}\equiv p_{i}-\frac{1}{n}\sum_{i}p_{i}=0$.suppl This result has interesting consequences in the thermodynamic limit, to be returned to shortly.

Surface transfer matrix. We now derive the relation between the surface gap and the Ising moment, with the help of a bond defect model described by the Hamiltonian

where $D_{\bm{q}}$ is independent of $\ell$, signifying a ferromagnetic order. The interlayer hopping $J_{\ell\bm{q}}$ depends on the magnetizations of the two layers it connects in the original domain wall model

where $M_{x}$ stands for reflection about a mirror perpendicular to $x$-axis. In this model, $J_{\ell\bm{q}}$ is uniform ($\mathrm{i}M_{x}J_{\bm{q}}$) except bond defects ($J_{\bm{q}}$) at the original domain walls, as schematized in Fig.3(a). It can be verified through symmetry analysis that the bond defect model Eq.(5) at $k_{x}=0$ can be obtained from the domain wall model Eq.(1) by applying $\mathrm{i}M_{x}$ to every $\downarrow$ layer, and hence the two models are equivalent at $k_{x}=0$ apart from a local gauge transformation.(suppl, ) This bond defect model allows us to analyze the layer-wise transfer matrices of the (010) surface at $k_{x}=0$, and since the surface Dirac point occurs also at $k_{x}=0$, the surface transfer matrices devised can be used to analyze the existence of surface gap.

For Eq.(5), the FM bulk without bond defects has a pair of gapless (010) surface modes, for it corresponds to the perfect AFM state before the gauge transformation. The surface Dirac cone is unfolded in the FM Brillouin zone, as shown in Fig.3(b). Thus, the in-gap surface spectra describe a two-channel ballistic conductor with the transfer matrix at an energy

where $T(\xi)=\mathrm{diag}[e^{\mathrm{i}\xi},e^{-\mathrm{i}\xi}]$, and $k_{z}$ is the wavevector of the surface mode moving in $-z$ direction.

We define a defect zone comprised of $2L$ layers centered at a bond defect, depicted in Fig.3(a). $L$ is large enough so that the evanescent wave escaping the defect zone is negligible. A superlattice of defect zones supports gapless surface states since it corresponds to an AFM configuration with $\bar{S}$ symmetry before the gauge transformation. Its surface spectra are shown in Fig.3(c), from which we can write down the surface transfer matrix of a defect zone

where $\theta$ is the wavevector in the supercell Brillouin zone of the mode propagating in $-z$ direction at a given energy. The matrix $W(v)=[[1,u]^{T},[u^{*},1]^{T}]$ accounts for the scattering of the FM surface modes as an electron enters the bond defect. suppl

Consider a supercell of $N$ (even) layers with $n$ (even) bond defects. Between the $i$th and $(i+1)$th defect zones there are $d_{i}$ ferromagnetic layers. The surface transfer matrix of the supercell is then

where $\delta\phi_{i}$ is difference of $\phi_{i}\equiv d_{i}k_{z}$ with its average $\bar{\phi}=\bar{d}k_{z}$, and in the second line we have introduced $T_{b}(\theta^{\prime},u^{\prime})=T(\bar{\phi})T_{b}(\theta,u)$. If the supercell indeed corresponds to an AFM supercell with zero magnetization, then for the energy at which $\theta^{\prime}=0$, the trace of $T_{n}$ to second order is found to be suppl

where $\mathcal{I}-\mathcal{I}^{*}$ is given by Eq.(4) using $\delta\phi_{i}=k_{z}\delta p_{i}$.

As revealed by Eq. (9) the trace invariant of the surface transfer matrix, like the surface gap, is a function of $\mathcal{I}$ . This actually is not a coincidence, since $\det T_{n}=1$, $T_{n}$ supports propagating modes when $|\operatorname{tr}T_{n}|<2$. Eq.(9) thus indicates a surface gap at $\mathcal{I}\neq\mathcal{I}^{*}$ when $\operatorname{tr}T_{n}>2$. And at $\mathcal{I}=\mathcal{I}^{*}$, $\operatorname{tr}T_{n}=2$ implies that the eigenvalues of $T_{n}$, $\lambda_{1}=\lambda_{2}=1$, corresponding to a surface Dirac point. Consequently, Eq.(9) confirms the empirical relation in Eq. (4), regarding the existence and value of $\mathcal{I}^{*}$ when a domain wall-decorated AFM TI possesses gapless (010) surface states. Moreover, the Dirac point of $\mathcal{I}=\mathcal{I}^{*}$ supercells occurs at $\theta^{\prime}=0$ implies that the Dirac point energy depends only on the number density of domain walls, and not on their arrangement.

Delocalization in thermodynamic limit. In the thermodynamic limit where $n\rightarrow\infty$ and $N/n\rightarrow$ const., the localization of (010) surface at $k_{x}=0$ is characterized by the Lyapunov exponent (Goldsheid89, ; Kramer93, ) of $T_{n}$:

which is related to the dimensionless conductance through $g\sim\operatorname{sech}^{2}n\gamma$.(Pichard86, ) We calculate the Lyapunov exponent Geist90 according to Eq. (8) with $T^{d_{i}}_{0}(k_{z})$ calculated as $\sigma^{d_{i}}_{z}T(\phi_{i})$. On the premise that the domain walls are placed randomly owing to weak interactions, the nearest-neighbor domain wall separations independently follow an identical exponential distribution. Accordingly, $\phi_{i}$ is sampled as a continuous variable via $\phi_{i}=\bar{\phi}x$ with $x$ drawn according to the probability density $p(x)=\exp(-x),x\in[0,\infty)$. As only whether $d_{i}$ is even or odd enters into $\sigma_{z}^{d_{i}}$, they are sampled as a binary sequence, fulfilling the condition $n_{\Uparrow}=n_{\Downarrow}$. Concerning the domain wall magnetizations, two ensembles (I and II, to be described) have been examined.

Ensemble I represents the non-interacting limit, where the domain wall magnetizations are uncorrelated so that $\Uparrow$ and $\Downarrow$ appear entirely by chance. As plotted in Fig.4(a), $\gamma$ in ensemble I is computed as a function of $\theta$ and $\bar{\phi}$, with $u=u_{0}/\cos\theta$ to account for increasing back-scattering by a bond defect with larger $\theta$. Since the low-energy surface modes of the perfect FM bulk and defect zone superlattice show linear dispersion $E=v_{0}k_{z}+\omega_{0}=v_{d}\theta+\omega_{d}$, $\bar{\phi}$ and $\theta$ are also related linearly, i.e. $\theta=\kappa\bar{\phi}+\theta_{0}$. $\kappa\bar{d}$ and $\theta_{0}$ are fixed by the model parameter specification which determines the values of Fermi velocities $v_{0}$, $v_{d}$ and Dirac points energy $\omega_{0}$, $\omega_{d}$. Typically, $\gamma$ values are seen to be finite over the energy range examined, indicating a generic localization of the surface states of a domain wall-decorated AFM TI when the domain wall magnetizations are uncorrelated.

The configurations in ensemble I generally correspond to $\mathcal{I}\neq\mathcal{I}^{*}$. Although $n_{\Uparrow}-n_{\Downarrow}=0$ over the entire sample, $|n_{\Uparrow}-n_{\Downarrow}|$ is of order $O(\sqrt{L})$ over a segment with $L$ domain walls due to statistical fluctuation. This means that typically $|\alpha_{i}|\sim O(\sqrt{n})$, and the fluctuation in Ising moment $\mathcal{I}-\mathcal{I}^{*}\sim(\overline{\alpha^{2}}/n)^{1/2}\sim O(1)$, with $\overline{\alpha^{2}}=\frac{1}{n}\sum_{i}\alpha^{2}_{i}$. Since the excess magnetic moments are carried by domain walls, $\alpha_{i}$ provides a measure of the macroscopic magnetization fluctuation. Consequently, $\mathcal{I}$ fails to converge in the thermodynamic limit, due to the unbounded magnetization fluctuation $|\alpha_{i}|$, signaling macroscopic breaking of time-reversal symmetry.

It follows immediately that if the fluctuation $|\alpha_{i}|$ is bounded in the thermodynamic limit, then $\mathcal{I}\rightarrow\mathcal{I}^{*}$ and the surface states are expected to be gapless. This might correspond to the physical situation for two considerations, namely, the mediated AFM exchange interactions between domain walls that favors cancellation of magnetic moments, and the magnetic dipole interaction that suppresses magnetization on macroscopic scales(Ashcroft76, ). For a demonstration, domain wall sequences in ensemble II are generated using only three kinds of segments, ”$\Uparrow\Downarrow$”, ”$\Uparrow\Uparrow\Downarrow\Downarrow$” , ”$\Uparrow\Uparrow\Uparrow\Downarrow\Downarrow\Downarrow$”, and their cyclic permutations. In this case $|\alpha_{i}|\leq 3$ is bounded. Numerical results in Fig.4(b) show an oval-shaped region with vanishing $\gamma$ near $\theta=\overline{\phi}=0$ tilted along $\theta^{\prime}=0$, which corresponds to the surface Dirac point of these $\mathcal{I}=\mathcal{I}^{*}$ configurations. The vanishing $\gamma$ (i.e. diverging localization length) suggests a delocalization transition at the Dirac point in ensemble II when the model line $\theta=\kappa\bar{\phi}+\theta_{0}$ crosses the oval region.

Discussions. The computed delocalization transition near the Dirac point is closely related that in a 1D random hopping model due to an emergent chiral symmetry at half filling.balents1997 ; steiner1999 ; evers2008 A configuration in ensemble II is comprised of segments with zero magnetization, whose transfer matrices $\{M_{i}\}$ satisfy an effective time-reversal symmetry $\sigma_{1}M_{i}\sigma_{1}=M_{i}^{*},\;\forall i$.mello1991 Consequently, the random matrices $\{M_{i}\}$ conform to the orthogonal symmetry classevers2008 and equivalently describe the solution of a 1D stochastic Dirac equationsuppl ; comtet2010

which is invariant under $\sigma_{1}\mathcal{K}$ ($\mathcal{K}$ for complex conjugation). Here, the energy-dependent potential $V_{E}(x)$ and mass $m_{E}(x)$ arise from stochastic pointer scatterers described by $M_{i}$. This model is known to possess delocalized solutions when $V_{E}(x)=0$ and $\langle m_{E}(x)\rangle=0$.balents1997 ; mathur1997 ; steiner1999 ; evers2008 Indeed, the oval region of delocalization along $\theta^{\prime}=0$ for ensemble II (FIG.4(b)) is close to the delocalization criticality of the Dirac equation: on the one hand, $\theta^{\prime}=0$ puts the energy at the Dirac point in ensemble II, which corresponds to the Dirac point of Eq.(11) (e.g. when $m_{E}=0$) appearing at $\langle V_{E}(x)\rangle=0$; on the other hand, $\mathcal{I}-\mathcal{I}^{*}=0$ implies $\langle m_{E}(x)\rangle=0$ at the energy where $\theta^{\prime}=0$, rendering the solution of Eq.(11) also gapless.

The mapping to Eq.(11) is also feasible for ensemble I. suppl But since $\mathcal{I}-\mathcal{I}^{*}\neq 0$ generally indicates $\langle m_{E}(x)\rangle\neq 0$, ensemble I typically shows localization, except at $\theta=k_{z}=0$ when the surface spectra become chiral symmetric ($\omega_{0}=\omega_{d}$) and the total transfer matrix is an identity matrix. For the chiral symmetric models with $\theta_{0}=0$, the dip in the $\gamma$-$\bar{\phi}$ plot for ensemble I displays a standard critical scaling $\gamma\propto 1/\log|E-E_{c}|$ (see Fig.4(c) and note $\bar{\phi}\sim E$), as is also expected from a Eq.(11) that describes a chiral symmetric Dirac Hamiltonian (i.e. when $V_{E}(x)=E_{c}-E$).balents1997 ; evers2008 ; ramola2014 A caveat is that at $E=E_{c}$, $\gamma$ may not be interpreted as the inverse localization length owing to large fluctuations,Texier10 ; ramola2014 ; balents1997 ; mathur1997 ; evers2008 as is also seen in our numerical results in FIG.4(c). In stark contrast, Fig.4(d) shows the delocalization transition in ensemble II exhibits an unusual algebraic scaling $\gamma\propto|E-E_{c}|^{\nu}$ ($\nu\sim 1.76$), and a substantial energy range with vanishing $\gamma$. Moreover, the sample-to-sample fluctuation in $\gamma$ is suppressed toward the critical point, pointing to an intriguing scenario that a sample adopting ensemble II turns out to be a “good conductor” at sufficiently low energies on its topological surfaces.

Experimentally, these results indicate a few intriguing phenomena to be studied on the surfaces of layered AFM TIs, concerning especially the $\mathrm{MnBi_{2n}Te_{3n+1}}$ family (zhang2015, ; otrokov2019, ; zhang2019b, ; deng2020a, ; ando2015, ; li2019, ). We emphasize that our results apply generally to any topological surface preserving a $S$-symmetry and a $M\Theta$-symmetry with mirror $M$ parallel to $z$.suppl Spectroscopic measurements may be employed to monitor the bandgap on any topological surfaces satisfying the above condition, and it is clearly interesting to find out whether such surfaces are gapped. Transport measurement will reveal the critical behavior of a material realization of the crossover from orthogonal to chiral orthogonal ensemble, where the unusual scaling and conductance fluctuation also offer valuable information on the magnetic correlation among AFM domain walls in the bulk material.