## INTRODUCTION

Over the past decade, topology has taken the center stage in condensed matter physics and materials science, emerging as an organizing principle of the states of matter (*1*). Many topological phases, such as quantum spin Hall insulators, quantum anomalous Hall (QAH) insulators, three-dimensional (3D) topological insulators (TIs), and Weyl semimetals, have been observed (*2*, *3*). Despite tremendous progress, the dominant majority of known topological materials (apart from the QAH insulators) are nonmagnetic materials, whereas studies on magnetic topological materials have been far more limited. In contrast to their nonmagnetic counterparts, magnetic topological materials allow a distinct set of new topological states including the axion insulator, the magnetic Weyl semimetals, the Chern insulators, and the 3D QAH insulators (*2*, *3*). Furthermore, magnetism is a natural way to induce nonnegligible electronic interactions, paving the way for studying the interplay between band topology and correlations. Therefore, topological magnetic materials have emerged as the frontier of the field. Since the realization of the QAH state in the magnetically doped TI thin films Cr_{0.15}(Bi_{0.1}Sb_{0.9})_{1.85}Te_{3} (*4*), recent studies have identified a range of topological magnetic materials such as Fe_{3}Sn_{2}, Co_{3}Sn_{2}S_{2}, Mn_{3}Ge, and Co_{2}MnGa (*5*–*9*). However, the magnetically doped TI thin films are nonstoichiometric systems where disorder and inhomogeneity are unavoidable. Fe_{3}Sn_{2}, Co_{3}Sn_{2}S_{2}, Mn_{3}Ge, and Co_{2}MnGa are all large carrier density metals where substantial topologically trivial bands coexist with the topological bands at the chemical potential (*5*–*9*). To avoid these drawbacks, an intrinsically ferromagnetic (FM) topological material with magnetic ions occupying their own crystallographic sites and with only the topological bands at the charge neutrality energy is strongly desired but has so far remained elusive.

Recently, MnBi_{2}Te_{4}, a van der Waals (vdW) compound composed of the septuple layers (SLs) of [MnBi_{2}Te_{4}], was identified as an intrinsic magnetic topological material with clean band structure (*10*–*32*). Unfortunately, magnetic moments in MnBi_{2}Te_{4} are antiferromagnetically (AFM) coupled across adjacent [MnBi_{2}Te_{4}] planes (*15*, *31*). Therefore, although the quantized anomalous Hall conductance was observed at a record-breaking temperature of 4.5 K, it requires an external magnetic field as large as 12 T to polarize the system into the FM state (*21*, *22*). Can we reduce the interlayer AFM coupling between the adjacent Mn layers to realize intrinsic FM? One material design strategy is to increase the interlayer distance between the adjacent Mn layers. How can we achieve it? Structurally, SL blocks have great compatibility with quintuple (QL) blocks of [Bi_{2}Te_{3}], as suggested by the existence of GeBi_{4}Te_{7}, which has alternating [GeBi_{2}Te_{4}] SL and [Bi_{2}Te_{3}] QL building blocks (*33*). This superior compatibility provides us flexible structural control to reduce the interlayer magnetic coupling by increasing the interlayer distance between the adjacent [MnBi_{2}Te_{4}] layers. Based on our material design strategy, here, we report the discovery of an intrinsic FM TI MnBi_{8}Te_{13}. Despite the interlayer distance between the adjacent [MnBi_{2}Te_{4}] SLs being 44.1 Å, it is striking and surprising that our thermodynamic, transport, and neutron diffraction measurements indicate that MnBi_{8}Te_{13} has long-range FM order below 10.5 K, with the easy axis along the *c* axis. Our first-principles calculations and angle-resolved photoemission spectroscopy (ARPES) measurements further suggest that it is an intrinsic FM axion state. Considering the natural heterostructure nature of MnBi_{8}Te_{13}, our finding provides a superior material realization to explore zero-field QAH effect, quantized topological magnetoelectric effect, and associated phenomena.

## RESULTS

### Crystal structure of MnBi_{8}Te_{13}

Although the existence of MnBi_{8}Te_{13} was mentioned (*34*), its crystal structure was never reported, partially due to the difficulty in growing the MnBi_{8}Te_{13} phase and separating it from the other members in the MnBi* _{2n}*Te

_{3n+1}family. We managed to grow high-quality MnBi

_{8}Te

_{13}single crystals and solved the crystal structure of MnBi

_{8}Te

_{13}by refining the room-temperature powder x-ray diffraction (XRD) pattern using various structure models. We found that MnBi

_{8}Te

_{13}crystallizes in the

*R*-3

*m*space group with the lattice parameters:

*a*=

*b*= 4.37485(7) Å,

*c*= 132.415(3) Å, α = β = 90°, γ = 120°. The refinement results and structural parameters are summarized in tables S1 and S2. The powder XRD pattern and the Rietveld refinement are shown in Fig. 1A. The low-angle diffraction peaks can be well indexed as the (0 0 L) with lattice parameter

*c*as 132.415(3) Å (see the inset of Fig. 1A), which is longer than the parameter

*c*of 101.825(8) Å in MnBi

_{6}Te

_{10}. The formation of MnBi

_{8}Te

_{13}is also visualized by the scanning transmission electron microscopy (STEM) image, as shown in Fig. 1B. The STEM image shows vdW structures, which is composed of repeating units of one SL block made of seven atomic layers and three consecutive QL blocks made of five atomic layers. The chemical analysis via the wavelength dispersive spectroscopy measurement results in the elemental composition of the sample as Mn:Bi:Te = 0.74(3):8.2(1):13, where the Mn:Bi ratio is slightly smaller than the ideal 1:8. The smaller than ideal Mn:Bi ratio is also observed in other MnBi

_{2n}Te

_{3n+1}because of the Bi substitution on Mn sites (

*16*,

*31*).

The crystal structure is shown in Fig. 1C. It is characterized by the alternating stacking of monolayer of MnTe_{6} octahedra that are well separated by a number of monolayers of BiTe_{6} octahedra running along the *c* axis. The distance between the nearest Mn-Mn interlayers is 44.1 Å, which is much larger than 13.6 Å for MnBi_{2}Te_{4}, 23.8 Å for MnBi_{4}Te_{7}, and 33.9 Å for MnBi_{6}Te_{10}, as shown in the inset of Fig. 1C. The stacking sequence of the MnBi_{2n}Te_{3n+1} series can be rationalized (*31*). Using Bi_{2}Te_{3} shown in Fig. 1C as the starting point, along the *c* axis, the stacking sequence of Bi_{2}Te_{3} is -A-B-C-A-B-C-, where A, B, and C represent the bilayers of BiTe_{6} octahedra whose bottom Te atoms, center Bi atoms, and top Te atoms are on the cell edges, respectively. The “Mn” layer can replace “A” or “B” or “C” bilayers of BiTe_{6} octahedra to make MnBi_{2n}Te_{3n+1}. For example, as shown in Fig. 1C, for MnBi_{2}Te_{4}, the stacking sequence is -Mn(B)-C-Mn(A)-B-Mn(C)-A-; for MnBi_{4}Te_{7}, the stacking sequence is -Mn(B)-C-A-; and for MnBi_{6}Te_{10}, the stacking sequence is -Mn(B)-C-A-B-Mn(C)-A-B-C-Mn(A)-B-C-A-. When it comes to MnBi_{8}Te_{13}, the stacking sequence is Mn(B)-C-A-B-C-Mn(A)-B-C-A-B-Mn(C)-A-B-C-A-, exactly identical to the one we obtained based on our powder x-ray refinement. Therefore, with the stacking rule, we can easily assign the stacking sequence for the yet-to-be-discovered higher *n* members of MnBi_{2n}Te_{3n+1} or design new magnetic TIs with the QL and SL building blocks.

### Long-range ferromagnetism and its strong coupling with charge carriers in MnBi_{8}Te_{13}

As a comparison, the physical properties of MnBi_{6}Te_{10} where two [Bi_{2}Te_{3}] QLs are sandwiched between the adjacent [MnBi_{2}Te_{4}] SLs are also presented. The inset of Fig. 2 (A and B) presents the temperature-dependent specific heat data of MnBi_{8}Te_{13} and MnBi_{6}Te_{10}, respectively. A specific heat anomaly associated with the magnetic phase transition is observed, which determines the ordering temperature as 10.5 K for MnBi_{8}Te_{13} and 11.0 K for MnBi_{6}Te_{10}. The magnetic properties of MnBi_{8}Te_{13} are shown in Fig. 2 (A, C, and E), whereas those of MnBi_{6}Te_{10} are presented in Fig. 2 (B, D, and F). The data indicate that, with the *c* axis as the easy axis, MnBi_{8}Te_{13} is FM below 10.5 K, while MnBi_{6}Te_{10} is AFM below 11.0 K. Figure 2 (A and B) presents the zero-field-cooled (ZFC) and FC magnetic susceptibility data, χ* ^{c}* (

*H*∥

*c*) and χ

*(*

^{ab}*H*∥

*ab*), measured at 0.1 kOe for MnBi

_{8}Te

_{13}and MnBi

_{6}Te

_{10}, respectively. Sharp contrast can be seen. A large bifurcation of ZFC and FC data of χ

*appears below 10.5 K in MnBi*

^{c}_{8}Te

_{13}, where upon cooling the ZFC data decrease but the FC data increase, suggesting the formation of FM domains. However, for MnBi

_{6}Te

_{10}, we observed a sharp cusp feature centering at 11.0 K in χ

*, similar to the ones in AFM MnBi*

^{c}_{2}Te

_{4}and MnBi

_{4}Te

_{7}(

*17*,

*35*) but with a small bifurcation of ZFC and FC data below 9 K. Furthermore, at 2 K, the magnitude of the FC χ

*in MnBi*

^{c}_{8}Te

_{13}is orders larger than that in AFM MnBi

_{2}Te

_{4}, MnBi

_{4}Te

_{7}, and MnBi

_{6}Te

_{10}(

*17*,

*35*). This strongly suggests different types of ground states in these two materials, with MnBi

_{8}Te

_{13}being FM and MnBi

_{6}Te

_{10}being AFM. Our conclusion is further confirmed by the hysteresis loop of isothermal magnetization curves

*M*(

^{c}*H*) (

*H*∥

*c*) shown in Fig. 2 (C and D) for MnBi

_{8}Te

_{13}and MnBi

_{6}Te

_{10}, respectively. At 2 K, unlike the

*M*(

^{c}*H*) data where multiple-step features are observed because of the spin-flop transition in MnBi

_{2}Te

_{4}(

*15*,

*17*) and spin-flip transition in MnBi

_{4}Te

_{7}(

*35*) and MnBi

_{6}Te

_{10},

*M*(

^{c}*H*) in MnBi

_{8}Te

_{13}shows a typical hysteresis loop for FM materials with coercivity of

*H*= 0.75 kOe and saturation remanence of

_{c}*M*= 3.1 μ

_{r}*/Mn. Upon warming,*

_{B}*H*decreases as the hysteresis loop shrinks in area due to the enhanced thermal fluctuations. Figure 2E shows the

_{c}*M*(

^{ab}*H*) and

*M*(

^{c}*H*) of MnBi

_{8}Te

_{13}up to 7 T, where a field of 10.4 kOe is required to force spins to saturate in the

*ab*plane. This value is 11.3 kOe for MnBi

_{6}Te

_{10}(Fig. 2F). Both sets of data indicate Ising anisotropy, with the

*c*axis as the easy axis in these two compounds. Ising anisotropy has also been suggested in MnBi

_{2}Te

_{4}by the inelastic neutron scattering measurement (

*28*). The saturation moment measured at 7 T for MnBi

_{8}Te

_{13}is

*M*= 3.5(

_{s}*1*) μ

*/Mn, whereas it is 3.4(*

_{B}*1*) μ

*/Mn for MnBi*

_{B}_{6}Te

_{10}. In both cases,

*M*is smaller than 4.6 μ

_{s}*/Mn, the obtained value from our density functional theory (DFT) calculation, but similar to 3.6 μ*

_{B}*/Mn in MnBi*

_{B}_{2}Te

_{4}(

*17*) and 3.5 μ

*/Mn in MnBi*

_{B}_{4}Te

_{7}(

*35*). The reduced Mn saturation moment in MnBi

_{8}Te

_{13}is likely due to the Mn site disorder as suggested in MnBi

_{2}Te

_{4}and MnBi

_{4}Te

_{7}(

*16*,

*34*). The fitted Weiss temperatures of MnBi

_{8}Te

_{13}are

= 12.8 K and

${\mathrm{\Theta}}_{w}^{c}$= 12.1 K, suggesting FM exchange interactions. The fitted effective moments are

${\mathrm{\mu}}_{\text{eff}}^{\mathrm{ab}}=5.4{\mathrm{\mu}}_{B}$/Mn and

${\mathrm{\mu}}_{\text{eff}}^{c}=5.1{\mathrm{\mu}}_{B}$/Mn, indicating weak single-ion anisotropy. Although the effective moment is smaller than 5.9 μ* _{B}*/Mn for a Mn

^{2+}ion, it is similar to MnBi

_{2}Te

_{4}, MnBi

_{4}Te

_{7}, and MnBi

_{6}Te

_{10}(

*17*,

*35*).

Figure 2G presents the temperature-dependent in-plane (ρ* _{xx}*) and out-of-plane resistivity (ρ

*) of MnBi*

_{zz}_{8}Te

_{13}. Above 20 K, upon cooling, both ρ

*and ρ*

_{xx}*decrease semilinearly. Upon further cooling, when approaching*

_{zz}*T*

_{c}, because of the stronger scattering from enhanced spin fluctuations, both values increase (

*36*,

*37*). Then, a sharp drop appears in both curves below 10.5 K due to the loss of spin scattering. The drop in ρ

*below*

_{zz}*T*

_{c}is distinct from that in MnBi

_{2}Te

_{4}, MnBi

_{4}Te

_{7}, and MnBi

_{6}Te

_{10}(see Fig. 2H) (

*17*,

*35*), where the antiparallelly aligned Mn moments enhance resistivity via the spin-flip scattering. Therefore, the drop in ρ

*again suggests FM ordering in MnBi*

_{zz}_{8}Te

_{13}because the parallelly aligned Mn moments along the

*c*axis eliminate such spin-flip scattering and thus reduce electric resistivity.

To confirm the long-range FM phase transition in MnBi_{8}Te_{13}, we have performed single-crystal neutron diffraction experiments at various temperatures. The quality of the instrument resolution limited magnetic Bragg peaks indicates that the stacking faults in this crystal are minimal. Figure 3A shows the Bragg reflection (1 0 1) at 4.5 and 15 K. It is evident that the peak intensity enhances substantially from 15 to 4.5 K, indicating the development of long-range magnetic ordering with the propagation vector **k** = 0. Starting from the nuclear space group *R-3m* with Mn atoms occupying the Wyckoff position (0, 0, 0) and **k** = 0, through symmetry analysis (*38*), we obtained the *k*-maximal magnetic subgroup as *R-3m*′ whose symmetry only allows the presence of an FM arrangement with the Mn spins along the *c* axis. Its magnetic structure is illustrated in Fig. 1C. The order parameter of MnBi_{8}Te_{13} was measured from 4.5 to 15 K. As shown in Fig. 3B, a power law was used to fit the intensity of magnetic reflection (1 0 1) as a function of temperature. The fit yields the critical temperature 9.8(3) K and a critical exponent β = 0.4(1). The former is consistent with our thermodynamic measurements. The latter is considerably larger than the value of 0.125 expected for the Ising model in two dimensions, excluding this possibility (*39*). More data at lower temperatures could better define the critical behavior.

Strong coupling between charge carriers and magnetism is observed in MnBi_{8}Te_{13} through the magnetotransport measurements, as shown in Fig. 4. ρ* _{xx}*(

*H*), ρ

*(*

_{zz}*H*), and ρ

*(*

_{xy}*H*) follow the same hysteresis as that in

*M*(

*H*). Using

*n*=

*H*/

*e*ρ

*, our 50-K data correspond to an electron carrier density of 1.66 × 10*

_{xy}^{20}cm

^{−3}, similar to that of MnBi

_{2}Te

_{4}(

*26*). Figure 4F shows the ρ

*(*

_{xx}*H*) curves up to 9 T at various temperatures. The overall “W” shape was observed at low temperatures. The W-shaped

*MR*(

*H*) was previously observed in MnBi

_{4}Te

_{7}

^{38}, where it was suggested to be a combination of nonmagnetic parabolic MR contribution and negative MR originating from FM fluctuations. Unlike MnBi

_{4}Te

_{7}where the W-shaped behavior was observed far above its transition temperature 13 K, in MnBi

_{8}Te

_{13}, the MR quickly changes from the deepest W shape at

*T*

_{c}into a parabolic shape at just a few degrees above it. This may suggest weak FM fluctuations above

*T*

_{c}in MnBi

_{8}Te

_{13}.

### FM axion state revealed by DFT calculation

To investigate the band topology, in the FM configuration with the spin oriented along the *c* axis, we calculated the electronic band structure of MnBi_{8}Te_{13}. Using the experimentally determined lattice parameters and atomic positions, our calculation shows only a continuous bulk energy gap (note S3). In contrast, using the experimental lattice parameters with relaxed atomic positions, our calculation indicates a 170-meV global energy bandgap (Fig. 5B). By comparing our DFT calculation with the experimental ARPES data, the one with relaxed atomic positions describes the material well.

To highlight the spin splitting in the presence of FM ordering, we present the <*S _{z}*> resolved band structure in Fig. 5C. The band structure projected on the Bi

*p*and Te

*p*orbitals shows that the bands near the Fermi level mostly originate from the Bi

*p*and Te

*p*orbitals. As shown in Fig. 5D, there are clear band inversions between the Bi

*p*and Te

_{z}*p*states. The Bi

_{z}*p*orbitals reach deep into the valence bands, indicating multiple possible band inversions that originate from the different [Bi

_{z}_{2}Te

_{3}] QL and [MnBi

_{2}Te

_{4}] SL of MnBi

_{8}Te

_{13}.

The presence of clear band inversions around the Γ point hints toward a topological phase. To unravel the exact topology of this system, we first compute the Chern number in the *k _{z} = 0* and

*k*π planes. In both planes, the Chern number is found to be zero. Next, we compute the parity-based higher-order

_{z}=*Z*

_{4}invariant, which is given by

(1)

Here, ξ* _{n}* (Γ

*) is the parity of the*

_{i}*n*th band at the

*i*th time reversal invariant momentum (TRIM) point Γ

*and*

_{i}*n*

_{occ}is the number of occupied bands. The

*Z*

_{4}invariant is well defined for an inversion symmetric system, even in the absence of time reversal symmetry (

*40*–

*42*). The odd values of

*Z*

_{4}(1, 3) indicate a Weyl semimetal phase, while

*Z*

_{4}= 2 implies an insulator phase with a quantized topological magnetoelectric effect (axion coupling θ

*=*π) (

*43*). A detailed list of the number of occupied bands with even (

*n*

^{+}_{occ}) and odd (

*n*

^{−}_{occ}) parity eigenvalues at the eight TRIM points is shown in table S3. Based on this, the computed

*Z*

_{4}invariant is found to be 2, which demonstrates that MnBi

_{8}Te

_{13}is an intrinsic FM axion insulator.

It is noted that, in calculations with or without atomic relaxation, the characteristics of band inversion and the topology of whole system remain the same, i.e., an FM axion state. In addition, we have investigated the evolution of the band structures by changing lattice constants *a*, *b* (in-plane), and *c* (out-of-plane) simultaneously. Within the range from −3 to +3% in MnBi_{8}Te_{13}, we did not find any crossing point or the additional band inversion feature. Our calculation reveals that the axion phase in MnBi_{8}Te_{13} is quite stable.

### Surface state revealed by ARPES and DFT

To investigate whether surface states appear in MnBi_{8}Te_{13}, we have performed small-spot (20 μm× 50 μm) ARPES scanned across the surfaces of MnBi_{8}Te_{13}. According to the crystal structure of MnBi_{8}Te_{13}, four different terminations are expected during the cleave, as visualized in the cartoon pictures in Fig. 6. Our ARPES data and DFT calculations reveal distinguishing surface states for the four different terminations of MnBi_{8}Te_{13}, which are summarized in Fig. 6.

Figure 6 (A to D) shows the isoenergy surfaces at the Fermi level, which uniquely fingerprints each of the four possible terminations of MnBi_{8}Te_{13}. We identify a unique circular Fermi surface for the SL termination, which has been observed in other MnBi_{2n}Te_{3n+1} (*n* = 1, 2, 3) (*44*–*46*). We find that the QL* _{n}* terminations show a dominant sixfold symmetric Fermi surface with decreasing cross-sectional area from QL

_{1}to QL

_{3}. Our assignment of the isoenergy surfaces to the respective terminations for MnBi

_{8}Te

_{13}is consistent with previous measurements on the simpler MnBi

_{2}Te

_{4}, MnBi

_{4}Te

_{7}, and MnBi

_{6}Te

_{10}compounds (

*24*,

*44*–

*46*). Figure 6 (E to H) presents the experimental ARPES

*E-k*maps along the

*M*→ Γ →

*M*high symmetry direction. Different types of surface states are observed for each of the four terminations. For the SL termination, a gapless or nearly gapless surface state appears (Fig. 6E), consistent with the gapless surface state observed on the SL termination in previous ARPES measurements of MnBi

_{2}Te

_{4}, MnBi

_{4}Te

_{7}, and MnBi

_{6}Te

_{1}(

*24*,

*44*–

*46*). The QL

_{3}termination, which is unique to MnBi

_{8}Te

_{13}, shows mass renormalization near the Fermi level (Fig. 6H). Furthermore, a large surface gap of ~105 meV centering around the charge neutrality point of −0.35 eV is revealed in the QL

_{1}termination and is highlighted by the arrow (Fig. 6F). We find that the spectra taken below and above the transition temperature are very similar in terms of the size of the gap, except for thermal broadening of the electrons (note S2). This has been observed in various magnetic TIs, including MnBi

_{2}Te

_{4}, MnBi

_{4}Te

_{7}, and MnBi

_{6}Te

_{10}, and the origin of it is under debate (

*24*,

*47*). This aspect will also be discussed more in note S2.

To confirm the topological nature of MnBi_{8}Te_{13}, we calculated the surface spectral weight throughout the (001) surface Brillouin zone (BZ) using the semi-infinite Green’s function approach for SL, QL_{1}, QL_{2}, and QL_{3} terminations correspondingly, as shown in Fig. 6 (I to L). Each state is plotted with a color corresponding to the integrated charge density of the state within the topmost QL or SL. In addition, we shift the Fermi level to match the experimentally observed Fermi level.

Excellent agreement between the experimental ARPES data and DFT calculation is achieved for the QL* _{n}* terminations as shown in Fig. 6. Our calculation suggests that the sizable surface bandgap presented in the QL

_{1}termination around −0.35 eV is a hybridization gap induced from the hybridization effect between the topmost QL and the nearest-neighboring SL (note S4). Although the calculated bandgap value is smaller than that from ARPES, it is a well-known problem in that generalized gradient approximation (GGA) generally underestimates the bandgap in semiconductors and insulators (

*48*,

*49*). The hybridization gap seems universal in MnBi

_{2n}Te

_{3n+1,}which has been suggested for MnBi

_{6}Te

_{10}(

*50*) and MnBi

_{4}Te

_{7}(

*51*). However, until now, it is unclear if the hybridization gap exists and supports QAH in the 2D limit. To settle down this issue, we have performed DFT calculation on a seven-layered finite-sized slab model with the arrangement of vacuum-[QL-SL-QL-QL-QL-SL-QL]-vacuum, i.e., a structure that has the QL

_{1}termination on both surfaces. The main features of the slab results shown in Fig. 5E are consistent with our semi-infinite Green’s function approach. We then calculate the Chern number of this slab model based on the Wilson loop method (

*52*). Our result shows a nontrivial Chern number (

*C*= 1) in this hybridization gap (Fig. 5F), demonstrating that MnBi

_{8}Te

_{13}is a QAH insulator in its 2D limit if the Fermi level is gated to the middle of the large hybridization gap.

As a consequence of strong exchange fields from the Mn magnetic layer in the SL termination, our calculation on the SL termination results in a parabolic large gapped surface band dispersion in the bulk energy gap (Fig. 6I). This is in sharp contrast with the gapless Dirac surface state revealed by the ARPES (Fig. 6E). For the SL termination, the deviation between ARPES and DFT calculation as well as the gapless Dirac surface state are universal in the MnBi_{2n}Te_{3n+1} family (*n* = 1, 2, 3) (*24*, *44*–*46*). The unexpected gapless surface state is argued to be caused by the surface spin reconstruction when the magnetic [MnBi_{2}Te_{4}] layer is at the vacuum-sample interface (*24*, *45*, *46*). Further experiments in identifying the spin reconstruction on the SL termination are urged to settle down this issue.

## METHODS

### Sample growth and characterization

We grew single crystals of MnBi_{2n}Te_{3n+1} (*n* = 3 and 4) using self-flux (*17*, *35*). Around 20 g of Mn, Bi, and Te elements with the molar ratio of MnTe:Bi_{2}Te_{3} as 15:85 was loaded into a 5-ml alumina crucible and then sealed inside a quartz tube under one-third atmosphere of Ar. The sample ampule was heated up to 900°C, stayed for 5 hours, and then quickly air-quenched to room temperature. It was then placed inside a furnace, which was preheated to 595°C. The ampule stayed at 595°C for 5 hours and then slowly cooled down to the decanting temperature in 48 hours. We then let the ampule stay at the decanting temperature for about 24 hours before we separated the single crystals from the liquid flux using a centrifuge. By varying the decanting temperature, we could obtain different phases of MnBi_{2n}Te_{3n+1} (note S1). We found out that the growth window for MnBi_{4}Te_{7} and MnBi_{6}Te_{10} is around 2°, while the one for MnBi_{8}Te_{13} is limited to be only 1° barely above the melting temperature of the flux. The method to determine the decanting temperature for each phase is described in note S1 in detail. Bi_{2}Te_{3} is the inevitable side product. We also noticed that with increasing *n* numbers, the chance of the intergrowth between Bi_{2}Te_{3} and MnBi_{2n}Te_{3n+1} (*n* = 3, 4) got enhanced. Therefore, extra care was paid in screening out the right piece. XRD at low angles for both the top and bottom (0 0 *l*) surfaces as well as powder XRD were performed on a PANalytical Empyrean diffractometer (Cu Ka radiation). Structural determination based on powder XRD was done using FullProf Suite software (*53*). STEM was measured on a sample prepared by focused-ion beam, and images are made on a Cs-corrected transmission electron microscope (Titan, Thermo Fisher Scientific, FEI) at 300 kV. Images are Bragg-filtered. Electric resistivity was measured in a Quantum Design (QD) DynaCool Physical Properties Measurement System. All samples were shaped into thin rectangular bars, and the four- and six-probe configurations were used for electrical resistivity and Hall resistivity, respectively. The magnetoresistivity was symmetrized using ρ* _{xx}*(B) = (ρ

*(B) + ρ*

_{xx}*(−B))/2, and the Hall resistivity was antisymmetrized using ρ*

_{xx}*(B) = (ρ*

_{xy}*(B)−ρ*

_{xy}*(−B))/2. The magnetization was measured using the QD Magnetic Properties Measurement System. The piece used for the magnetization measurement was later ground into fine powder, whose powder XRD pattern showed no impurity and was used for structural determination.*

_{xy}### Single-crystal neutron diffraction

Single-crystal neutron diffraction experiments were performed at the HB-3A DEMAND single-crystal neutron diffractometer equipped with a 2D scintillation Anger camera with 0.6-mm spatial resolution at the High Flux Isotope Reactor at Oak Ridge National Laboratory in the temperature range of 4.5 to 15 K (*54*). Neutron wavelength of 1.551 Å was selected by using a bent perfect Si-220 monochromator. Rocking curve scans, which cover the full peak profile of the Bragg reflection (1 0 1), were carried out at 4.5 and 15 K, respectively. The peak (1 0 1) was selected for revealing the magnetic ordering parameter from 4.5 to 15 K.

### First-principles calculations

The bulk band structures of MnBi_{8}Te_{13} were computed using the projector augmented wave method as implemented in the VASP package (*55*–*57*) within the GGA (*58*) and GGA plus Hubbard *U* (GGA + *U*) (*59*) scheme. On-site *U* = 5.0 eV was used for Mn *d* orbitals. The spin-orbit coupling was included self-consistently in the calculations of electronic structures with a Monkhorst-Pack *k*-point mesh 5 × 5 × 5. The experimental lattice parameters were used. The atomic positions were relaxed until the residual forces were less than 0.01 eV/Å. We used Mn *d* orbitals, Bi *p* orbitals, and Te *p* orbitals to construct Wannier functions, without performing the procedure for maximizing localization (*60*).

The surface spectral weight throughout the (001) surface BZ was calculated using the semi-infinite Green’s function approach for four distinct surface terminations. We set up four different models to calculate the surface state to compare with the ARPES data, for example, for the SL termination, the semi-infinite Green’s function model was set as vacuum-[SL-QL-QL-QL]-[SL-QL-QL-QL]-infinite; to get that of the QL_{1} termination, the model was set to be vacuum-[QL-SL-QL-QL]-[QL-SL-QL-QL]-infinite. This method was valid because the thickness of the QL or SL is about 10 Å, which roughly equals the typical escape depth of photoelectrons.

### ARPES measurements

ARPES measurements on single crystals of MnBi_{8}Te_{13} were carried out at the Stanford Synchrotron Research Laboratory (SSRL) beamline 5-2 with photon energies of 26 eV with linear horizontal polarization and a 20 μm× 50 μm beam spot. Single-crystal samples were top-posted on the (001) surface and cleaved in situ in an ultrahigh vacuum better than 4×10^{−11} torr and a temperature of 15 K.

**Acknowledgments: ****Funding:** Work at UCLA was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0011978. Work at CU Boulder was supported by NSF-DMR 1534734. Work at ORNL was supported by the U.S. DOE BES Early Career Award KC0402010 under contract DEAC05-00OR22725. This research used resources at the High Flux Isotope Reactor, DOE Office of Science User Facilities operated by the ORNL. We acknowledge M. Hashimoto and D.-H. Lu for help with the ARPES measurements. Use of the Stanford Synchrotron Radiation Lightsource, SLAC National Accelerator Laboratory, is supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under contract no. DE-AC02-76SF00515. B.G. wants to acknowledge CSIR for the senior research fellowship. B.G. and A.A. thank CC-IITK for providing the HPC facility. T.-R.C. was supported by the Young Scholar Fellowship Program from the Ministry of Science and Technology (MOST) in Taiwan, under a MOST grant for the Columbus Program MOST108-2636-M-006-002, National Cheng Kung University, Taiwan, and National Center for Theoretical Sciences, Taiwan. This work was supported partially by the MOST, Taiwan, under grant MOST107-2627-E-006-001. This research was supported, in part, by Higher Education Sprout Project, Ministry of Education to the Headquarters of University Advancement at National Cheng Kung University (NCKU). TEM characterization and analysis were supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences (BES), Division of Materials Science and Engineering (DMSE) through Early Career Research Program under award KC0203020:67037. The work at Northeastern University was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences under grant no. DE-SC0019275, and benefited from Northeastern University’s Advanced Scientific Computation Center and the National Energy Research Scientific Computing Center through DOE grant no. DE-AC02-05CH11231. **Author contributions:** N.N. conceived the idea and organized the research. C.H., N.N., S.M., and J.L. grew the bulk single crystals and carried out x-ray, thermodynamic, and transport measurements and data analysis. H.C. and L.D. carried out structure determination and neutron diffraction measurements. K.N.G., H.Li, A.G.L., and D.D. carried out the ARPES measurements and data analysis. B.G., S.-W.L., H.-J.T., C.-Y.H., P.V.S.R., B.S., A.A., A.B., T.-R.C., S.-Y.X., and H.L. performed the first-principles calculations. M.S. and D.L. performed TEM measurements. All authors wrote the manuscript. **Competing interests:** The authors declare that they have no competing interests. **Data and materials availability:** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.