Subtle metastability of the layered magnetic topological insulator MnBi2Te4 from weak interactions

The metastable layered compound MnBi2Te4 is the first experimentally realized intrinsic antiferromagnetic topological insulator, predicted to host the quantum anomalous Hall effect at high temperatures upon exfoliation to atomically thin layers. While its magnetic ordering and topological properties have generated intensive interest, the mechanism behind its metastability and the ideal crystal synthesis conditions have remained elusive. Here, using a combined first-principles-based approach that considers lattice, charge, and spin degrees of freedom, we investigate the metastability of MnBi2Te4 by calculating the Helmholtz free energy for the reaction Bi2Te3 + MnTe ->MnBi2Te4. We identify a narrow temperature range (767 K to 873 K) in which the compound is stable with respect to the competing binary phases and successfully synthesize high-quality MnBi2Te4 single crystals using the Bi-Te flux method within this range. We also predict the various contributions to the total specific heat, which is consistent with our experimental measurements. Our findings indicate that the degrees of freedom responsible for the van der Waals interaction, magnetic coupling, and nontrivial band topology in layered materials not only enable emergent phenomena but also determine thermodynamic stability. This conclusion lays the foundation for future computational material synthesis of novel layered systems.

Introduction. MnBi 2 Te 4 is an intrinsic antiferromagnetic topological insulator (TI) currently under intensive study [1], in which the interplay between magnetism and topology is expected to cause the quantum anomalous Hall (QAH) effect upon exfoliation to atomically thin layers [1,2]. While bulk MnBi 2 Te 4 is characterized by a Z 2 invariant protected by the combination of time-reversal and half-lattice translational symmetry, the QAH insulator phase is enabled by a nonzero Chern number (Z invariant) in odd-layer thin films [1,3]. MnBi 2 Te 4 is formed by intercalating magnetic MnTe layers, wherein the open d-shell Mn ions form localized magnetic moments, into the quintuple layers of Bi 2 Te 3 , which is topologically protected due to spin-orbit coupling (SOC)induced band inversions [4,5]. Hence, as-synthesized MnBi 2 Te 4 is also a layered material bound by the van der Waals (vdW) interaction.
The experimental synthesis of large, pure, high-quality MnBi 2 Te 4 single crystals suitable for magnetotransport measurements is notoriously difficult. The first reported polycrystalline synthesis of MnBi 2 Te 4 was achieved via congruent melting of stoichiometric melts of raw Bi, Mn, and Te, followed by rapid quenching, then annealing at 808 K [6]. Additionally, by using the binary phases as precursors, single crystals of MnBi 2 Te 4 , stable in a narrow temperature range below 873 K, have been synthesized from a solid-state reaction of Bi 2 Te 3 and α-MnTe followed by water-quenching and annealing at 838.15 K for 10 days [7]. However, the obtained samples likely suffered from antisite defects resulting in nonstoichiometry [7], which could prevent the topological protection of the material due to magnetic disorder. While the pow- * Corresponding author: jsun@tulane.edu der samples were found to be metastable at room temperature [6], bulk single crystals could be cooled to low temperatures without decomposition [7]. More recently, large single crystals of MnBi 2 Te 4 have been obtained using a Bi-Te flux with molar ratio Mn:Bi:Te = 1:10:16 [8]. The self-flux and vertical Bridgman methods have also been used, both requiring extremely precise control of temperatures [2,9]. The above experiments all confirmed that MnBi 2 Te 4 is a metastable phase that can only be synthesized within a narrow temperature range below 873 K. [6][7][8][9]. Without a fundamental understanding of its metastability, the ideal growth conditions for pure, large crystal synthesis will remain elusive, ultimately preventing further experimental study of its unique topological and magnetotransport properties. Furthermore, although there have been a significant number of materials predicted to host similarly interesting topological properties, there are only handful of experimental realizations of such predictions [10][11][12][13]. Hence, the understanding of the metastability of MnBi 2 Te 4 is also highly desirable for the synthesis of other candidate layered topological materials.
Thermodynamic stability of materials can be predicted from chemical reaction free energy based on firstprinciples density functional theory (DFT) calculations. Nevertheless, for complex layered magnetic quantum materials like MnBi 2 Te 4 , such predictions are challenging, due to the increasing importance of various weak interactions (including SOC, magnetic coupling, vibrational anharmonicity, and vdW interactions), most of which are also responsible for their emergent properties, as illustrated in Figure 1. Typically, these materials have different kinds of chemical bonds ranging from strong in-plane covalent bonds to weak interlayer vdW interactions, with interaction strengths across almost 3 orders of magnitude (from ∼1 eV for strong bonds to ∼1 meV for vdW interactions). Both strong chemical bonds and weak interactions are equally important for the thermodynamic stability predictions, demanding a single density functional approximation for simultaneously accurate descriptions. More importantly, the electronic, magnetic, and vibrational thermal excitation energies of these materials at finite temperature can be up to several tens of meV; although this scale is typically negligible for stability predictions of bulk solids [14], it is comparable to the thermal energy required to stabilize MnBi 2 Te 4 . Since DFT is a zero-temperature ground state electronic structure method, post-DFT models have to be used to describe such excitations. Therefore, all these interactions have to be considered and described accurately for the thermodynamic stability of MnBi 2 Te 4 , either at the DFT level or by DFT-based models.
Here, using a set of first-principles-based approaches that consider lattice, charge, and spin degrees of freedom to take into account the various weak interactions, we calculate the temperature-dependent reaction free energy of MnBi 2 Te 4 based on the reaction MnTe + Bi 2 Te 3 → MnBi 2 Te 4 , which identifies a narrow temperature range for the metastability of MnBi 2 Te 4 consistent with experiments. These approaches are then validated by comparing the predicted specific heat capacities of the relevant compounds to experimental results.
After all interactions are taken into account, Figure  2 shows that MnBi 2 Te 4 is unstable with respect to its competing binary phases (MnTe and Bi 2 Te 3 ) for low temperatures, and becomes stable above 767 K, indicated by the calculated reaction free energy (solid blue line). However, since the melting points are ∼853 K and ∼873 K for Bi 2 Te 3 and MnBi 2 Te 4 , respectively [7], the quenching required for single crystal synthesis of MnBi 2 Te 4 without the minor Bi 2 Te 3 phase must occur between ∼853 K and ∼873 K. This narrow temperature window, in addition to the fact that MnBi 2 Te 4 is only more stable than the competing binaries by about 6.8 meV/f.u. in this window, explains why previous experimental attempts at synthesizing large single crystals have been difficult.
Below, we analyze the importance of various weak interactions to the metastability of MnBi 2 Te 4 . We find that the vdW interaction, SOC, magnetic ordering, and lattice vibrations (phonons) play crucial roles in determining the thermodynamic stability of MnBi 2 Te 4 and tend to stabilize MnBi 2 Te 4 with respect to the competing binary phases as temperature increases, while the elec-tronic thermal excitation is negligible.
Van der Waals Interaction. It is well known that layered materials are bound by vdW interactions between layers, while the popular Perdew-Burke-Ernzerhof (PBE) density functional [15] misses most vdW interactions and barely binds layered materials [16,17]. Meanwhile, it has been shown that the strongly-constrained and appropriately-normed (SCAN) density functional [18,19] combined with the revised Vydrov-Van Voorhis (rVV10) nonlocal correlation density functional [20] is effective in accurately describing structural and energetic properties of layered materials [17]. As expected, Figure 2 shows that PBE destabilizes MnBi 2 Te 4 with respect to the competing binaries by about 80 meV/f.u. in comparison with SCAN+rVV10 for the 0 K ground state energies. Since 80 meV/f.u. is significantly larger than the finite temperature contributions to the reaction free energy from the electronic, magnetic, and vibrational degrees of freedom, PBE will falsely predict that MnBi 2 Te 4 cannot be synthesized below its melting temperature, even in a metastable phase. In addition to the vdW interaction, it is likely that SCAN+rVV10 more accurately describes the d electrons of Mn (and hence the magnetic properties of MnTe and MnBi 2 Te 4 ) when compared to PBE, due to the reduced self-interaction errors, a result demonstrated by previous studies on many other transition metal compounds [18,19,[21][22][23][24].
Spin-orbit Coupling. It is usually assumed that SOC is a purely atomic effect and canceled out in a reaction free energy calculation [14]. Figure 2 shows that the inclusion of SOC stabilizes MnBi 2 Te 4 with respect to MnTe and Bi 2 Te 3 by about 10 meV, which is negligible for typical solid reactions but is critical for the metastability of MnBi 2 Te 4 . It is well known that SOC becomes important for heavy atoms, like Bi and Te, and is responsible for the topological properties of Bi 2 Te 3 and MnBi 2 Te 4 , while the open d-shell of Mn provides localized magnetic moments for MnTe and MnBi 2 Te 4 . The synergy of SOC and AFM-ordered localized magnetic moments results in MnBi 2 Te 4 being an intrinsic magnetic TI [25]. Our results further suggest that the presence of the magnetic moments enhances the SOC to lower the total energies more for MnBi 2 Te 4 and MnTe than for Bi 2 Te 3 (see supplementary materials, Figure S1).
Electronic Contribution. Figure 2 shows that the electronic contribution to the reaction free energy is negligible. This is to be expected, since the relevant compounds are semiconductors with band gaps no less than 0.1 eV (see supplementary materials, Figure S2), an energy level much larger than the electronic thermal excitations for the considered temperature range.
Vibrational Contribution (Harmonic Approximation and Anharmonicity). Figure 2 shows that lattice vibrations under the harmonic approximation (HA) destabilize the competing binaries more than the ternary phase. This is because the MnTe binary is much stiffer than the other two compounds, as illustrated by their equations of state (see supplementary materials, Figure S3), and thus has much fewer phonon modes at low frequencies in the computed phonon density of states (see supplementary materials, Figure S4). We also found that the lattice anharmonicity modeled by the Quasi-harmonic Debye (QHA-Debye) model [26,27] and the Slater approximation [27] to the Debye temperature Θ D stabilizes MnBi 2 Te 4 with respect to MnTe and Bi 2 Te 3 . This is probably due to the fact that MnTe is stiffer in its equation of state and experiences less anharmonicity in the temperature range considered here than the other two compounds. The anharmonicity has an effect of softening the phonon modes and therefore likely reduces the phonon frequencies of MnBi 2 Te 4 and Bi 2 Te 3 more than MnTe, favoring thermal population. This effectively stabilizes MnBi 2 Te 4 with respect to MnTe and Bi 2 Te 3 as temperature increases.
Magnetic Contribution. MnBi 2 Te 4 has roughly the same Mn magnetic moment as MnTe, ∼4.3 µ B [28], while the exchange coupling strength J ij is much weaker (see supplementary materials, Figures S5 and S6), likely due to the layered structure. Similar to the lattice vibration effect, the stronger magnetic coupling should result in higher frequency magnons that destabilize MnTe more than MnBi 2 Te 4 as temperature increases, as illustrated in Figure 2.
Specific Heat. Next, we turn to the calculated heat capacity at constant pressure (C p ) of the ternary phase MnBi 2 Te 4 and its competing binary phases Bi 2 Te 3 and MnTe, as shown in Figure 3. We choose to focus on C p because it is closely related to the free energy and is the key quantity obtained from calorimetric measurements that can be used to validate our calculations. Here, we neglect the thermal expansion effect on the electronic and magnetic contributions to C p and thus use C v for these two contributions instead. A complete description of the calculation is included in the computational methods section.
When compared to the most recent experimental data for each phase, our calculations are accurate below 300 K. Since heat capacity is closely related to the free energy, this result reinforces the reliability of our initial stability estimation for MnBi 2 Te 4 . Notable discrepancies between our calculations and experimental results for Bi 2 Te 3 (see Figure 3a) occur above ∼400 K. This is likely due to the formation of point defects at high temperatures, such as vacancies and antisite defects, which has been consistently observed in experiment [7]. While these defects affect C p , entropy, and free energy for Bi 2 Te 3 at high temperatures, similar formation of defects can be expected for MnBi 2 Te 4 [7]. Thus, the effects of defects in the total reaction free energy can mostly be canceled out. For both MnTe (T N ≈ 300 K) and MnBi 2 Te 4 (T N ≈ 25 K), we predict well-defined peaks in C p near the respective antiferromagnetic phase transition temperatures, consistent with experimental data.
Specific heat measurements were also taken up to ∼200 K on a single crystal sample of MnBi 2 Te 4 we synthesized using the Bi-Te flux method. When compared to FIG. 3. Calculated heat capacity from electronic, magnetic, and vibrational contributions, compared with available experimental data, for (a) Bi2Te3 [29][30][31], (b) MnTe [32,33], and (c) MnBi2Te4 [7]. Cv el, Cv mag and Cv HA are contributions from electronic, magnetic, and vibrational (HA) excitations calculated within the constant volume approximation at 0 K equilibrium volume. Cp QHA-Debye is the constant pressure specific heat, involving the thermal expansion contribution.
our calculated results, the low-temperature C p is quite accurate (see inset, Figure 3c). At intermediate temperatures (∼423 K to ∼767 K), thermal vibration is strong enough to excite MnBi 2 Te 4 out of the metastable state and separate it into the binary phases [6], which makes C p measurements nearly impossible.
Summary. Based on first-principles thermodynamics, we have analyzed the metastability of the recent antiferromagnetic TI candidate MnBi 2 Te 4 and have successfully synthesized high-quality single crystals for specific heat measurements. Our DFT-based approach yields specific heat results in good agreement with experimental data, validating our stability estimation of MnBi 2 Te 4 as a function of temperature. We confirm that MnBi 2 Te 4 is a metastable phase which can only be synthesized in a short range of temperatures (767 K to 873 K) with a small reaction free energy less than ∼6.8 meV/f.u., explaining the previous difficulties confronted in experiment. We demonstrate that fundamental weak interactions, including SOC, vdW, magnetic coupling, and lattice vibrations all contribute to this subtle high-temperature metastability. This finding potentially facilitates future computational discoveries of novel stable or metastable 2D magnetic topological materials.
Experimental methods. Single crystals of MnBi 2 Te 4 were synthesized using the Bi-Te flux method [8]. The starting materials of Mn, Bi, and Te powder were mixed with a molar ratio of 1:10:16 and loaded into an Al 2 O 3 crucible, then sealed in a quartz tube under high vacuum. The mixture was heated to 1173.15 K in a muffle furnace and held there for 24 hours to allow homogeneous melting, then slowly cooled down to 863.15 K at a rate of 3 K/h. After removing the excess Bi 2 Te 3 flux through centrifuging, plate-like single crystals were obtained. Their crystal structure was confirmed by X-ray diffraction measurements, and their chemical composition was examined using energy dispersive X-ray spectroscopy. The specific heat of MnBi 2 Te 4 up to 200 K was measured with an adiabatic relaxation technique using a commercial Physical Property Measurement System (PPMS, Quantum Design). This data is plotted in Figure 3(c).
Computational methods. In this work, we carry out DFT [34] calculations using the Vienna Ab-initio Simulation Package (VASP) [35] with the projector-augmented wave (PAW) method [36,37]. The recently developed strongly-constrained and appropriately-normed (SCAN) meta-GGA [18,19] is used for its superior performance in description of different chemical bonds and transition metal compounds [18,19,[21][22][23][24]. A long range vdW correction is combined with SCAN through the rVV10 nonlocal correlation [20], a revised form of VV10, the Vydrov-Van Voorhis non-local correlation functional [38]. The PAW method is employed to treat the core ionelectron interaction. An energy cutoff of 520 eV is used to truncate the plane wave basis. We use Γ-centered meshes with a spacing threshold of 0.15Å −1 for K-space sam-pling. Geometries of the three compounds were allowed to relax until the maximum ionic forces were below a threshold of 1 meVÅ −1 . We use the Phonopy code [39] to calculate the key information of phonon frequency ω q (where q is the wave vector), from harmonic force constants calculated by VASP within the DFPT method. The Perdew-Burke-Ernzerhof GGA [15] is used for this purpose due to its smoother potential energy surfaces for generating forces [40].
The key quantity for thermodynamic stability predictions of solids is the Helmholtz free energy, defined as Here, the P V contribution has been ignored, and the adiabatic approximation has been used, which decouples lattice, charge, and spin degrees of freedom for thermal excitations. P is the pressure, T the temperature, V the volume, and E0 (V ) the total electronic energy at zero temperature, which can be directly calculated using different density functional approximations. F el (V, T ), F vib (V, T ), and F mag (V, T ) are the contributions to the free energy at finite temperature from electronic, lattice, and magnetic degrees of freedom, respectively, which can be modeled based on DFT results. In this study, E0 (V ) was calculated using SCAN+rVV10.
Based on the electronic density of states obtained from DFT calculations (see supplementary materials, Figure S2), the electronic contribution to the free energy, F el (V, T ), can be determined by the finite temperature method [41] according to Fermi-Dirac distribution, following the fixed density of states approximation [42].
The lattice vibrational contribution F vib (V, T ) can be expressed as where ω (q, s) is the phonon frequency associated with wave vector q and band index s, and k B is the Boltzmann constant. ω (q, s) typically depends on V and T due to the anharmonicity of the lattice potential, and it can be calculated or modeled based on DFT at different levels. The HA assumes that the lattice sees a harmonic potential at the equilibrium volume (see Figure  1), and ω(q, s) is calculated based on the frozen-phonon approach, which computes the harmonic force constants from DFT calculations via finite atomic displacements [39]. The quasi-harmonic approximation (QHA) method [39,43] takes a step further to consider the volume dependence for the phonon anharmonicity, while the temperature is assumed to indirectly affect phonon vibrational frequencies through thermal expansion. Typically, the phonon spectra of about ten or more volumes are usually required for a DFT-based QHA simulation, and such calculations are always time-consuming. Here, F vib (V, T ) are instead approximated using the phonon density of states (see supplementary materials, Figure S4) by the Debye model. The Debye temperature is then estimated based on the equation of state from DFT calculations to account for the anharmonicity [26,27]. We fit the static E(V ) curve of each phase (see supplementary materials, Figure S3) to the Vinet equation of state [26,[44][45][46][47]. The total lattice vibrational contribution F vib (V, T ) is treated by combining the HA (as implemented in Phonopy [39]) and the QHA-Debye model (as implemented in Gibbs2 [26,27]) via a scaling factor (a function of the Poisson ratio) to recover the behavior of the HA at low temperature, similar to previous studies [48][49][50]. We chose not to use the conventional QHA method [51], not only because it is much more computationally expensive, but also because it requires more careful validation for layered systems (like GeSe [52]) and even some conventional materials (like Si [53,54]), in addition to the ubiquitous imaginary frequency problems [52,55].
An often-employed approach to capture F mag (V, T ) of systems with localized magnetic moments starts with the effective Heisenberg Hamiltonian mapped from DFT calculations, with first-principlesderived exchange coefficients J ij (see supplementary materials, Figures S5 and S6) that mediate the magnetic exchange between spinsŜ localized at lattice sites i and j. F mag (V, T ) can then be calculated by the eigenenergies ofĤ, similar to the previous expression for F vib (V, T ). However, for most realistic systems, an exact analytical solution is not known. Here, F mag (V, T ) corresponding toĤ is calculated using a rescaled Monte Carlo (rMC) method [56], which maps the classical Monte Carlo (CMC)-obtained thermodynamic quantities to those obtained by quantum Monte Carlo (QMC), scaled by a factor dependent on the spin quantum number S. Although the QMC solution to F mag (V, T ) is desired, it has a limited applicability due to the so-called (negative) sign problem [57,58]. CMC is known to be useful and reliable for magnetic critical temperature predictions and acceptable for high temperature (above critical temperature) heat capacity predictions. At 0 K, however, CMC gives a finite heat capacity. This is incorrect, since the zero-temperature heat capacity contributed by magnons should be zero according to quantum statistics of the magnon excitations. The use of rMC corrects the low temperature part effectively. Explicit coupling terms, e.g., phonon-magnon, phonon-electron, and higher order phonon-phonon interactions are all assumed to be small and neglected. Formation of defects at high temperatures (such as vacancies and antisite defects [7]), especially for Bi 2 Te 3 and MnBi 2 Te 4 , may also contribute to the C p and the free energy. This defect factor, however, is expected to mostly cancel between the two phases and is possibly negligible to the total free energy. Therefore, this factor is not considered in our calculation.