A first-principles investigation of Topological Phase Transition in face-centred cubic LiMgBi

TL;DR

This study uses first-principles calculations to reveal a topological phase transition in LiMgBi under volume expansion, transforming it from a band insulator to a strong topological insulator with potential applications in electronics and spintronics.

Contribution

It is the first to identify and characterize a topological phase transition in face-centered cubic LiMgBi driven by volume expansion using first-principles methods.

Findings

LiMgBi undergoes a topological phase transition at 4.0% volume expansion.

The phase transition results in a Dirac semi-metal with band inversion.

LiMgBi is identified as a strong topological insulator with potential applications.

Abstract

Topological Insulators (TI) exhibit robust spin-locked dissipationless Fermion transport along the surface states. In the current study, we use \textit{first-principles} calculations to investigate a Topological Phase Transition (TPT) in a Half-Heusler (HH) compound LiMgBi driven by a Volume Expansive Pressure (VEP) which is attributed to the presence of, intrinsic voids, thermal perturbations and/or due to a phenomena known as cavity nuclei. We find that, the dynamically stable \textit{face-centred cubic} (FCC) structure of LiMgBi (which belongs to the F$\overline{4}$3m[216] space group), undergoes TPT beyond a critical VEP at 4.0\%. The continuous application of VEP from 0.0\% to 8.0\% results in a phase transition from a, band insulator to a Dirac semi-metal nature. Qualitatively, the Dirac cone formation and band inversion along the high symmetry point $\mathbf{\Gamma}$ in the Brillouin Zone (BZ) are analysed in terms of Electronic Band Structure (EBS) and Projected Local Density of States (LDOS). The TPT is further characterised by the $\mathbb{Z}_2$ invariant, ($\nu_0$, $\nu_1$ $\nu_2$ $\nu_3$) $\equiv$ (1, 0 0 0) along the (0001) surface which indicates quantitatively that, HH LiMgBi is a strong TI. We hence propose, HH LiMgBi (known for its piezoelectric, thermo-electric and semi-conducting applications) as a strong TI with potential multipurpose application in the field of electronics, spintronics and quantum computation.

Figures

Tables

Table 1. Table 1 : Theoretical and experimental values of lattice constant (a) and energy band gap.

	Work	LiMgBi
Lattice Constant (Å)	PBE	6.81
	Theory [11]	6.71
	Experiment [11, 40]	6.74
Energy Band Gap (eV)	PBE	0.35
	PBE+HSE	0.70
	Theory [11]	0.62
	Experiment	-

Table 2. Table 2 : Energy band gap from GGA and HSE calculations at 0.0%, 4.0% (critical VEP) and 6.5% VEP.

Functional	% VEP	Band Gap (eV)
GGA	0.0	0.35
	4.0	0.00
	6.5	0.001
HSE	0.0	0.70
	4.0	0.00
	6.5	0.10

Equations

$$(-1{)}^{{\nu }_i}=\prod _{i=1}^8\delta ({\Lambda }_i)$$

$$\nu_{0}=\Big{(}\mathbb{Z}_{2}\big{(}k_{i}=0\big{)}+\mathbb{Z}_{2}\big{(}k_{i}=0.5\big{)}\Big{)}mod2$$

$$\nu_{i}=\mathbb{Z}_{2}\big{(}k_{i}=0.5\big{)}$$

Full text

Introduction

Topological Insulating (TI) phase of matter is interesting due to the exotic quantum properties such as, spin-locked dissipationless transport of Fermions along the topologically protected surfaces. As a result of this, the spin helicity along the surface has several potential applications in the field of spintronics, quantum computation etc [1]. Till date, several three dimensional (3D) materials have been predicted theoretically and some are realised experimentally [2, 3, 4, 5, 6] to exhibit TI nature. A variety of 3D compounds (such as, Bi2Se3, $\beta$-As2Te3, AMgBi (A = Li, Na, K) and elemental tellurium (Te)) exhibit TI nature driven by strain [7, 8, 9]. These materials also have multi-purpose applications [10, 11, 12]. It is known that, specific electronic properties of material depend explicitly on their crystal structure and associated space group. This implies that, crystals of different structural class with previously known applications can be designed precisely to investigate for the Topological Phase Transition (TPT) by employing well known band engineering technique which is based on application of strain or pressure (compressive or expansive in nature).

Strain engineering of the electronic band structure is utilized to realize a strong TI [7, 8]. This is because, the application of strain / pressure enhances the strength of Spin Orbit Coupling (SOC) in some systems [8]. But, there are systems such as, 3D elemental Te (three fold screw symmetric helical chains of Te atoms) where, it is observed that, weak van der Waals interaction between the helical layered arrangement of Te atoms can be manipulated by application of shear strain (violates three fold screw symmetry of helices) and hydrostatic pressure leading to band inversion which is a characteristic property of a TI system [9]. This suggests that, without venturing into the computationally costly; SOC calculations, simple application of strain or pressure can yield the desired TPT in materials [9, 13]. One such system is, Anti-Ferroelectric TI (AFTI) material AMgBi (A = Li, Na, K) which exhibits strong TI nature characterized by the $\mathbb{Z}_{2}$ invariant $\nu$ = 1 in both Pnma and P63mc structures [13]. This is achieved under the influence of epitaxial strain and hydrostatic pressure [13]. Off these, the most promising candidate TI was found to be LiMgBi under epitaxial strain leading to AFTI of Type I (represents normal anti-polar and topological polar states) [13].

LiMgBi which belongs to Pnma and P63mc space groups and exhibits AFTI of Type I [13] also exists as a Half Heusler (HH) compound (in the form PQR where, P and Q are transition metals and R is a p-block element) with a variety of applications such as a semi-conducting, piezoelectric and thermo-electric material [14, 11, 12]. Theoretical search for TI nature in HH compounds has been pursued since a long time [15]. Concrete efforts have led to theoretical predictions of HH compounds as candidate TI’s [16, 17, 18, 19, 20, 21, 22, 23, 24]. Also, there are studies where, strain tuning led to realization of TI nature in materials [25, 26, 27, 23, 24]. We perform analogous calculations to investigate TPT in Face-Centred Cubic (FCC) HH LiMgBi described by the F$\overline{4}$3m[216] space group in Hermann-Mauguin representation. Experimentally it is observed that, LiMgBi in FCC phase exhibits trivial band insulating behaviour [11]. One study suggests that, when screening potential method is employed [28], a large number of candidate HH compounds can be classified as TI and Normal Insulator (NI). Based on this study it turns out that, thirty three HH compounds exhibit TI nature, but, LiMgBi exhibits NI nature at ambient conditions. Feng et.al. show that, strain tuning of band order in cubic systems can lead to TI nature [25].

With this background, we investigate the possibility of TPT in FCC LiMgBi by strain tuning the band order which is quite uncharted. We find that, under Volume Expansive Pressure (VEP), HH LiMgBi undergoes a TPT at a critical value of 4.0% and exhibits a strong TI nature characterized by the $\mathbb{Z}_{2}$ invariant as (1, 0 0 0). Thus, indicating that, HH LiMgBi is a candidate TI material for novel applications in the field of spintronics, quantum computation etc.

Methodology

We performed density functional theory based first-principles calculations to investigate TPT under the application of pressure. For this purpose, we obtained crystallographic information for LiMgBi in F$\overline{4}$3m[216] space group from MaterialsProject repository [29]. Optimization of lattice constant (a) was performed following the convergence test for total energy with respect to, wave function cut off and k-mesh. We used norm conserving pseudopotentials under Generalized Gradient Approximation (GGA) which considers 1s1, 3s2 and 6s26p3 orbitals of Li, Mg and Bi respectively. The pseudopotential method is based on Martins-Troullier with exchange correlation of Perdew-Burke-Ernzerhof (PBE) functional type [30]. The optimization was performed by, finding a global minima in terms of the total energy of the system and then narrowing down (using bisection method) to a local minima with a constrain that, the total pressure on the atoms is 0.00 kbar. The converged value of plane wave kinetic energy cut-off and charge density are 50 Ry and 200 Ry respectively. A uniform Monkhorst-Pack Grid [31] for k-points of 7 $\times$ 7 $\times$ 7 was used in the self-consistent calculations. For better prediction of energy band gap [32], we performed calculations with Heyd–Scuseria–Ernzerhof (HSE06) screened Coulomb hybrid functional [33] in Vienna ab initio simulation package (VASP) [34] with proper optimization. Phonon calculations were performed using Density Functional Perturbation Theory (DFPT) [35]. For phonon calculations, we used 10 $\times$ 10 $\times$ 10 q-mesh which was followed by plotting the dynamical matrices in the entire Brillouin Zone (BZ). Following the analysis of structural and dynamic stability, we performed electronic calculations pertaining to the Electronic Band Structure (EBS), Density of States (DOS) and Surface States (SS). Optimization and electronic properties were calculated using the QUANTUM ESPRESSO code [36] while the SS were calculated using WannierTools (WT) [37] package which employed the Wannier Function Centres (WFC) obtained from a post processing package Wannier90 (W90) [38]. The electronic properties led to qualitative results whereas for quantitative characterization of the TI nature, WT package was used for thorough $\mathbb{Z}_{2}$ analysis.

Results and Discussion

Structural Properties and Dynamic Stability

The optimized value of the lattice constant (a) under GGA approximation is 6.81 Å. This value is in descent agreement with the experimental (6.74 Å) and other theoretical (6.71 Å) values [11]. The optimized crystal structure is shown in Figure 1 which was found to be dynamically stable. The primitive cell vectors for HH LiMgBi are defined in terms of lattice constant ‘a’ as, v1 = (a$/$2)(-1,0,1), v2 = (a$/$2)(0,1,1), v3 = (a$/$2)(-1,1,0) with atoms Li, Mg and Bi occupying, 4b, 4c and 4a Wyckoff positions respectively [39]. This structure results in a direct Band Gap (BG) as seen from the EBS presented in Figure 2 which is quite small as compared to previous studies (Table 1). Since it is known that, the error in BG prediction is smaller in HSE calculations as compared to pure DFT calculations [32], we performed HSE calculations at 0.0%, 4.0% and 6.5% VEP. The BG from HSE calculations are in strong agreement (Table 1) with previous studies and better than the GGA calculations.

We increase lattice constant (a) in steps of 0.5% from 0.0% to 8.0% to investigate for TPT in HH LiMgBi. We refer to this increment as VEP which can be attributed to an effect created by the presence of an intrinsic void (exerting a volume isotropic pressure in the outward direction) in the crystal structure or due to extreme internal thermal perturbations. Also, it can be attributed to experimentally observed phenomena such as cavitation pressure [41] due to impurities known as cavity nuclei. Increment in unit cell volume (a.u.3) due to application of VEP (Figure 4) is interpreted in terms of volume strain as shown in Figure 4. Also, we analyse the difference in Total Energy ($\Delta$E = Ei - E0) of the system (in Ry units) with respect to % VEP (Figure 5) due to the volume strain and find that, due to increment in VEP, the system deviates from its equilibrium state leading to a quantum phase transition.

Phonon dispersion calculations were performed to check the lattice dynamical stability of the proposed crystal system LiMgBi belonging to the F$\overline{4}$3m[216] space group. This is important because, Phonon studies are fundamental in underpinning the vibrational dynamics in a phase transition and the practical feasibility to synthesize the material [42]. In the current study, Phonon Dispersion Curves (PDC) are obtained as shown in Figure 6 at 0% VEP. As there are three (3) atoms in the unit cell of LiMgBi, we obtain nine (9) Phonon branches in the dispersion relation (Figure 6). The PDC shows that, the system LiMgBi is, dynamically stable and has no imaginary modes of vibration. From Figure 6, we can see that, in the lower frequency regime there are 3 modes which correspond to the acoustic branches and in the higher frequency regime, there are six (3N - 3, here N = 3) optical modes [43]. The acoustic branch exhibits one (in plane) Longitudinal Acoustic (LA) and two (in plane) Transverse Acoustic (TA) modes of vibration which are degenerate along the directions, W to X and L to $\mathbf{\Gamma}$ in the BZ. The optical branch exhibits two (N - 1, here N = 3) Longitudinal Optical (LO) and four (2N - 2, here N = 3) Transverse Optical modes of vibration which dominate in the higher frequency regime [43]. The TO modes are degenerate along the directions $\mathbf{\Gamma}$ to K, K to X and X to $\mathbf{\Gamma}$ in the BZ. Similarly, the LO modes are degenerate at $\mathbf{\Gamma}$ in the BZ.

With proper understanding of the structural and dynamic stability of LiMgBi, we proceed for further analysis. We shall broadly divide our discussions here off in two perspectives, qualitative which revolves around the electronic properties and quantitative which revolves around the $\mathbb{Z}_{2}$ analysis.

Electronic Properties and Surface States

EBS analysis has proved to be a fruitful qualitative method to look for TI nature in materials [44, 45, 6, 9, 8, 7, 13]. We hence, analyse the EBS for band inversion by band engineering which is a characteristic signatures of TI nature. It is observed that, band inversion occurs in LiMgBi when subjected to VEP in steps of 0.5% at a critical value of 4.0%. From the EBS we observe a phase transition from a band insulator like behaviour to a Dirac semi-metal like behaviour. A direct BG (0.35 eV and 0.70 eV from PBE and HSE06 calculations respectively) exists at high symmetry point ($\mathbf{\Gamma}$) in the BZ at 0% VEP (Figure 7). This is compared to previous studies (Table 1) where the band gap is $\sim$ 0.62 eV [11]. BG from our HSE calculations match closely with the BG from previous studies [11]. For qualitative analysis we continue to discuss the electronic properties in terms of band gap under GGA and HSE approximation (tabulated in Table 2 for three different values of % VEP).

It is observed from the band structure that, the Valence Band (VB) maxima has a doubly degenerate state along the high symmetry point $\mathbf{\Gamma}$ whereas, the Conduction Band (CB) minima has single non-degenerate state along the high symmetry point $\mathbf{\Gamma}$. This feature is protected by the two-fold symmetry of the FCC structure. With the gradual increment of VEP, at 4.0% VEP, LiMgBi undergoes a phase transition from band insulator nature to a Dirac semi-metal nature (Figure 7) with BG 0.00 eV. We represent the eigen states along the high symmetry line $\mathbf{\Gamma}$ as $\mathbf{\Gamma}_{i}$ where ‘i’ is the band index without inclusion of spin degeneracy. The EBS for 0%, 4.0% and 6.5% are shown in Figure 7. At 0% VEP the band order is $\mathbf{\Gamma}_{6}$ (dominated by p orbital), $\mathbf{\Gamma}_{8}$ (dominated by p orbital) and $\mathbf{\Gamma}_{10}$ (dominated by s orbital). Here, $\mathbf{\Gamma}_{6}$ and $\mathbf{\Gamma}_{8}$ bands are degenerate as seen in Figure 7(a). With increment in VEP, at 4.0% a Dirac Cone is formed as shown in Figure 7(b). This feature exhibits a linear dissipationless transport of Fermions along the surface of LiMgBi. Also, the formation of a Dirac Cone leads to degeneracy between $\mathbf{\Gamma}_{8}$ and $\mathbf{\Gamma}_{10}$ bands as seen in Figure 7(b). The formation of Dirac cone also, lifts up the degeneracy between bands $\mathbf{\Gamma}_{6}$ and $\mathbf{\Gamma}_{8}$. By further increment in VEP, BG reopens with to 1.1 meV (under GGA) as shown in Figure 7(c) and 0.10 eV (under HSE). This reopening characterizes the band inversion in terms of the inversion in band order with, $\mathbf{\Gamma}_{8}$ (dominated by p orbital) being exchanged with $\mathbf{\Gamma}_{10}$ (dominated by s orbital) wherein, $\mathbf{\Gamma}_{10}$ feature is transferred to $\mathbf{\Gamma}_{6}$ (dominated by p orbital) band. Also, as the CB crosses Fermi level it indicates that, LiMgBi exhibits conduction of Fermions along the surface.

In order to justify the band inversion, we discuss the orbital contributions of Li, Mg and Bi. For this purpose, projected LDOS were calculated to compliment the EBS. We refrain our discussions to the regions in the vicinity of Fermi energy level (EF) which is important for interpretation of band inversion. At 0% VEP, the LDOS is as shown in Figure 8. VB in the vicinity of Fermi level in Li is dominated by a small contribution from unhybridized ‘s’ orbital, while the CB is dominated majorly by unhybridized ‘s’ orbital. Similarly, the VB in Mg has small contributions from the ‘s’ orbitals and the CB has major contribution from the ‘s’ orbital, both being unhybridized. In both the cases (Li and Mg), the ‘s’ orbital contribution from the CB crosses Fermi energy level (EF) as seen from Figure 8. As compared to this, Bi has a strong contribution solely due to the ‘p’ orbital in VB with a strong hybridization with the ‘s’ orbital in CB. Also, in the VB a mild contribution from the ‘s’ orbital is observed but, near Fermi energy level (EF) the ‘s’ orbital has major contribution as seen in Figure 8. Beyond a critical value of VEP at 4.5% the orbital contributions for the Dirac cone formation at the Dirac point on the Fermi energy level (EF) are interpreted from the LDOS as shown in Figure 9. The VB and CB of Li and Mg are dominated by the ‘s’ orbital with a increment in its orbital contribution at EF indicating the formation of Dirac cone (Figure 9). Although, the CB in both (Li and Mg) retains the ‘s’ orbital dominant feature. Similarly, in Bi, the major contribution for Dirac cone formation comes from the VB ‘p’ orbital crossing the Fermi energy level (EF) which has a weak hybridization with ‘s’ orbital, the strong hybridization persists and dominates the deeper regions of CB (Figure 9). At 6.5% of VEP when the band reopens, the orbital contributions change as compared to the 0% and 4.5% case above. Figure 10 shows the LDOS at 6.5% VEP. From the LDOS for Li and Mg, we observe that, the orbital contribution due to ‘s’ orbital increases in the VB region (Figure 10). Similarly from the LDOS for Bi, the ‘p’ orbital contribution reduces in the VB and increases in the interior regions of CB with a strong hybridization with the ‘s’ orbitals, whereas, the ‘s’ orbital contribution is almost negligible in VB but increases slightly in the interior of VB. The changes observed in the orbital contributions are contemplated as an indicator of the band order rearrangement observed in the EBS, which is a characteristic signature of band inversion. The major contribution for the band inversion is due to the Bi which has strong share in orbital rearrangements.

Further analysis for a deeper realization of TI behaviour was performed by calculating the SS. Figure 11 is the computational Angle Resolved Photo Emission Spectra (ARPES) [46, 47, 6] showing the SS [48] at 0% and 4.5% VEP respectively. Beyond the critical value of VEP (i.e., beyond 4.0% VEP) we observe dissipationless Fermi Surface along the edge states in the BZ. This is a qualitative evidence of TI nature characterized by, a single Dirac cone SS (similar to previous work [49]). This predicts that, experimental realization of FCC LiMgBi can be characterized by performing ARPES to examine the SS along the (0001) surface.

$\mathbb{Z}_{2}$ Analysis

The $\mathbb{Z}_{2}$ invariant is essential to characterize a material into a $\mathbb{Z}_{2}$ topological class. There are various methods to obtain the $\mathbb{Z}_{2}$ invariant [50, 51, 52, 53, 2] for systems with and without inversion symmetry. Generally, for systems with inversion symmetry, the product of parities of eigen values is calculated at eight (8) Time Reversal Invariant Momenta (TRIM) points for a 3D system (equation 1) and at four (4) TRIM points for a two dimensional (2D) system [2]. The $\mathbb{Z}_{2}$ invariant $\nu_{0}=1$ and $\nu_{0}=0$ corresponds to strong and weak TI respectively [2, 6]. Evolution of $\mathbb{Z}_{2}$ invariant along the application of pressure is shown in Figure 12. Since, the system is FCC, which lacks inversion symmetry, the $\mathbb{Z}_{2}$ analysis was performed by calculating the $\mathbb{Z}_{2}$ invariants \big{(}\mathbb{Z}_{2}\big{)}_{\pi} and \big{(}\mathbb{Z}_{2}\big{)}_{0}, for two different planes of momentum (k) in the BZ resulting in the $\mathbb{Z}_{2}$ invariant $\nu_{0}$ as, \nu_{0}=\Big{[}\big{(}\mathbb{Z}_{2}\big{)}_{\pi}-\big{(}\mathbb{Z}_{2}\big{)}_{0}\Big{]}\big{(}mod2\big{)}. These calculations were performed using the WannierTools (WT) package [37] which utilizes the information regarding the Wannier Function Centres (WFC) obtained by post processing package Wannier90 (W90) [38]. Maximally Localised Wannier Functions (MLWF) were obtained after PWSCF calculations using Marzari and Vanderbilt (MV) method in W90 package. The MV method minimizes the gauge dependent spread $\tilde{\Omega}$ with respect to the Bloch states U*(k)* which are obtained from the PWSCF calculations [36]. The WFC are then used for computation of Wilson Loops in WT. WT was used to obtain Tight Binding (TB) model for LiMgBi using the MLWF from W90. Using this TB model, we compute the $\mathbb{Z}_{2}$ invariant by calculating the Wilson loops around the Wannier Charge Centres (WCC). This method is one of the best methods for systems without inversion symmetry [37]. Computational APRES is obtained using WT along the (0001) crystal direction. WT is a reliable software used quite often for such studies [51, 54, 55, 56, 57].

[TABLE]

[TABLE]

[TABLE]

The underlying computation in WT is carried out along the six Time Reversal Invariant Planes (TRIP) kx = 0,$\pi$, ky = 0,$\pi$ and kz = 0,$\pi$ in the BZ. The $\mathbb{Z}_{2}$ indices ($\nu_{0}$, $\nu_{1}$ $\nu_{2}$ $\nu_{3}$) are computed using equation 2 and 3 [37]. Based on these calculations, LiMgBi under 0% VEP is characterized by $\mathbb{Z}_{2}$ invariant $\nu_{0}$ = 0 which implies that LiMgBi is a weak TI. But, beyond the critical value of VEP, at 4.5% the $\mathbb{Z}_{2}$ indices are obtained to be ($\nu_{0}$, $\nu_{1}$ $\nu_{2}$ $\nu_{3}$) $\equiv$ (1, 0 0 0) which implies that, LiMgBi is a strong TI. Thus we conclude quantitatively that, LiMgBi undergoes a TPT under the application of VEP.

Conclusion

In summary, we performed density functional theory based first-principles calculations and find that, HH compound LiMgBi which was previously studied for its semi-conducting, piezoelectric and thermo-electric properties undergoes a TPT driven by VEP. The FCC structure of LiMgBi was subjected to VEP isotropically by increment of the lattice constant (a). At equilibrium, the EBS has a direct BG of $\sim$ 0.35 eV (under GGA) and $\sim$ 0.70 eV (under HSE06). With gradual increment in VEP, we find a critical value at which, LiMgBi undergoes a phase transition from a trivial band insulator to a Dirac semi-metal material. With further increment in VEP, the band reopening is observed which is quite small; with a BG of the order of $\sim$ 1.1 meV (under GGA) and $\sim$ 0.10 eV (under HSE) indicating, band inversion due to exchange of orbital contribution in the VBM and CBM. This was complemented by thorough analysis of projected LDOS in order to understand the orbital contribution of the constituent elements. In order to indicate experimental characterization of the TI nature, SS were calculated using computational ARPES along the direction (0001) which exhibits the presence of a surface Dirac cone. These results were further quantified by studying the evolution of $\mathbb{Z}_{2}$ invariant along the application VEP. The strong TI nature of HH LiMgBi is characterized by the $\mathbb{Z}_{2}$ index, ($\nu_{0}$, $\nu_{1}$ $\nu_{2}$ $\nu_{3}$) $\equiv$ (1, 0 0 0). We thus conclude that, HH compound LiMgBi with FCC crystal structure undergoes a TPT driven by VEP. These calculations will open up diverse perspectives for multi-purpose applications of LiMgBi as a strong TI apart from the known semi-conducting, piezoelectric and thermo-electric applications.

Acknowledgement

Authors acknowledge, Department of Science and Technology (DST), Government of India, for providing financial assistance.

Figures