Signature of topological crystalline insulating behavior in new B2X2Zn (X=Ir, Rh, Co) compound from first-principles Computation

Recent attempts at topological materials have revealed a large class of materials that show gapless surface states protected by time-reversal symmetry and crystal symmetries. Among them, topological insulating states protected by crystal symmetries, rather than time-reversal symmetry are classified as topological crystalline insulators. We computationally predict the signature of new three-dimensional topological crystalline insulating compounds of space group 139(I/4mmm). After conducting a full volume optimization process by allowing to rearrange of atomic positions and lattice parameters, the first principles calculation with a generalized gradient approximation is utilized to identify multiple Dirac-type crossings around X and P symmetric points near Fermi energy. Importantly the band inversion at point P is recognized. Further, We investigate the compound for topological crystalline insulating behavior and identify metallic surface states on high-symmetry crystal surfaces with the projection to the plane (001). Additionally, we performed formation energy, elastic properties, and phonon modes calculations to verify the structural, mechanical, and dynamical stability of the compounds. Therefore, we suggest the compounds for further investigation and experimental realization.


I. INTRODUCTION
Following the discovery of a two-dimensional topological insulating phase in graphene with a significant effect of spin-orbit coupling (SOC) 1,2 rigorous efforts have been committed to exploring new topological materials. Among them, the identification of topological insulators (TIs) and topological semimetals (TSMs) with unusual physical properties and potential applications lead to enthralling discoveries both theoretically and experimentally. TIs presence of a bulk insulating gap and gapless conducting surface states that form a Dirac point 3,4 while topological semimetals (Dirac semimetals, Weyl semimetals, and nodal line semimetals) exhibit gapless bulk states and band crossings in momentum space which are inverted beyond the crossing point (or line) [5][6][7] protected by topology and symmetry.
The gapless surface states across the band gap in both two-and three-dimensional TIs characterized by a Z 2 topological invariant, are topologically protected by time-reversal symmetry 3,4,8,9 . The binary alloys Bi 2 Se 3 , Bi 2 T e 3 , and Sb 2 T e 3 that hold topological surface states are identified as strong TIs [10][11][12] . Topological crystalline insulators (TCIs) originated as a counterpart of topological insulators without the significant effect of SOC. The crystal symmetries in TCIs play a vital role in the topology of the bands showing gapless surface states across the insulating 9,13,15-17 . TCI phase was first theoretically predicted and experimentally realized by angle-resolved photoemission spectroscopy in SnTe, P b 1−x Sn x Se, and P b 1−x Sn x Te 13 . An integer topological invariant known as the mirror Chern number arising from the mirror symmetry describes the non-trivial topology of these materials 9 . Additionally, SnS, SeTe, Ca 2 As family, and antiperovskites have also been reported to be TCIs 14,[18][19][20][21] .
Also, PbSe, PbTe, and PbS are predicted to be TCIs with a suitable combination of applied pressure or strain 22 .
The surface states of TCIs against the magnetic field strength seem to be more powerful than that of the TIs due to the mirror symmetry preservation without preserving the time reversal, but it may be ample to break the mirror symmetry and to be a trivial insulator 15,23 . Therefore, there are plus and minus points with the crystalline protection of TCIs 13,14,18-21 , but it is desirable to find new topological materials 24-27 with pronounced transport properties and a wide range of controllability and functionalities for growing potential applications: quantum computing and spintronic devices.
In this work, a detailed study of electronic band structures and density of states (DOS) of new ternary B 2 X 2 Zn (X=Ir, Rh, Co) compound in the tetragonal crystal structure of symmorphic space group 139 [I4/mmm] is presented. To the best of our knowledge, these compounds have not yet been reported either experimentally, or theoretically. In this paper, we report calculated electronic band structures of B 2 X 2 Zn after completing the full volume optimization process by allowing to rearrange of atomic positions and lattice parameters to minimize the energy. We predict the formation of type I multiple Dirac crossings around X and P points on the Brillouinzone (BZ) near Fermi energy. With the careful investigation of electronic band structures and DOS of three compounds with and without the spin-orbit coupling (SOC), topological features are identified and discussed in detail. Further, we theoretically demonstrate that B 2 X 2 Zn has a signature of topological crystalline insulating features with mirror symmetry. We study the surface states on the crystal plane (001) by projecting XΓ to surface M Γ. Additionally, elastic constants and formation energy arXiv:2302.06056v1 [cond-mat.mtrl-sci] 13 Feb 2023 calculations are performed to verify the mechanical and structural stability of the new compounds. Further, the dynamical stability of compounds is assured by conducting phonon mode calculations.

II. COMPUTATIONAL METHOD
The first principle density functional theory (DFT) calculation implemented on Quantum ESPRESSO (QE) simulation package 28,29 has been performed. The plane wave pseudo-potential method formulated with generalized gradient approximation (GGA) and the Perdew-Burke-Ernzerhof (PBE) scheme has been employed with the ultra-soft pseudo-potentials from PSlibrary, including fully relativistic ultra-soft pseudo-potentials for SOC [30][31][32][33][34][35] . The k-mesh of 20 x 20 x 20, kinetic energy cutoff for wave functions of 80Ry, and kinetic energy cutoff for charge density and potential of 480Ry are held to receive extreme convergence of energy and charge to enhance the accuracy of the simulation. Additionally, WEIN2K simulation package 36,37 with PBE pseudo-potentials and plane-wave basis set with GGA is used to compare and verify the QE results. Phonopy software package interfaced with QE under harmonic approximation has been used to calculate the phonon spectrum 39 and ElaStic software package cooperated with QE has been used to calculate the full second-order elastic stiffness tensor 38 . The crystal momentum in units of k = (π/a, π/a, π/c) is used throughout the discussion unless otherwise specified.
The B 2 X 2 Zn structure is shown in Fig. 1 by denoting Zn, B, and X in blue, green, and red colored spheres respectively. The arrangement of atoms for Zn1, X1, and B1 layers is labeled as top-down views. The first BZ of the structure shows in Fig. 1 (b) with the high symmetric points and Γ(0, 0, 0) point located at the center. The (001) surface BZ shows in labeling with Γ, M , and X including surface BZ Γ at the center.

B. Volume Optimization
Since the volume that shows the lowest total energy can be identified as the ground state for the stable structure, first we optimize the crystals to find the cell volume that has the lowest total energy. In our calculations, we performed the geometry optimizations of the unit cell by allowing atomic positions to rearrange and the lattice parameters to move around together. We follow the Quantum Espresso volume optimization scheme with extreme convergence for non-relativistic ultra-soft pseudopotentials. The optimized lattice constants a and c for B 2 X 2 Zn are shown in the second and third columns of Table. I. Similarly calculated lattice positions for the optimized system is Zn(0, 0, 0), B1(0.25, 0.75, 0.50), B2(0.75, 0.25, 0.50), with a little differentiation of the Ir, Rh, and Co positions as Ir1(x, x, 0), Ir2(1-x, 1-x, 0), Rh1(y, y, 0), Rh2(1-y, 1-y, 0), Co1(z, z, 0), and Co2(1-z, 1-z, 0) with x=0.3725808923, y=0.3715044822, and z=0.3712045162. Hence, the crystal structures at optimized bulk lattice parameters and atomic positions are used for further investigation of B 2 X 2 Zn properties during the project.

C. Formation Energy
Since the compounds that show negative formation energies at thermal equilibrium with respect to their elemental phases are known to be stable, we calculate the formation energy to study the stability of the structures. In general, the formation energy per atom for ternary B 2 X 2 Zn can be calculated as where N B , N X , and N Zn are the numbers of B, X (Ir, Rh, Co), and Zn atoms in the unit cell, respectively. Since the B 2 X 2 Zn unit cell has 1 atom of Zn and 2 atoms of B and X, N Zn = 1 and N X = N B = 2 are taken. E B2X2Zn is the calculated total free energy of the compound, and E B , E X , and E Zn are the calculated total free energies per atom of the elemental phases of B, X, and Zn, respectively. During total energy calculation of B, X, and Zn, we use optimized structures of trigonal crystal of space group 166 (R3m) for B, a cubic crystal of space group 225 (Fm3m) for X, and hexagonal crystal of space group 194 (P 6 3 /mmc) for Zn. By using Eq. 1 we calculate the formation energy for B 2 X 2 Zn as -0.3005 eV/atom, -0.3245 eV/atom, and -0.2279 eV/atom, for X=Ir, Rh, and Co respectively, the last column of Table. I. Since all three compounds indicate negative formation energies with respect to their elemental phases, we identify those structures are theoretically stable.

D. Elastic Properties
First principles density functional calculation implemented on QE is used to explore the mechanical properties of the structures. Mechanical properties together with crystal stability and stiffness are easily investigated by calculating the elastic stiffness matrix C ij or flexibility matrix S ij (= [C ij ] −1 ). The bulk modulus, Young's modulus, shear modulus, and Poisson's ratio of polycrystals are calculated by Voigt-Reuss approximation methods 38 . Average poly-crystalline modules (Hill's average) are obtained by using the upper and lower limit of the actual effective modulus corresponding to the Voigt bound and Reuss bound. The Hill's average is said to be mostly agreed with the experimental result 38 . Since the structure is a tetragonal crystal, there are six independent non-zero elastic constants, namely C 11 , C 12 , C 13 , C 33 , C 44 , and C 66 . The calculated elastic constants and effective bulk, shear, and Young's modules for B 2 X 2 Zn are presented in Table. II.
The elastic constants calculated for the tetragonal crystal satisfy the following Born mechanical stability criteria as discussed in 40,41 : where C ij represent six independent non-zero elastic constants. The elastic tensor of the second order is calculated by using the expansion of the elastic energy in terms of the applied strain. The C ij results are obtained for large deformations with high-order polynomial fit by identifying the plateau regions which provides good reasonable results 38 . Calculated elastic constants for all three B 2 X 2 Zn compounds satisfy the Born mechanical stability criteria implemented in Eq. 2.

E. Phonon Frequencies
A collective excitation of a set of atoms that decomposed into different modes plays an essential role in material science. Therefore, we perform the first principles of phonon calculations with force constants which are said to be an important calculation for studying dynamical behaviors and thermal properties of the materials. Studying topological phonons and their properties is another frontier research field that we are not discussing here. Our focus here is to check the dynamical stability of the samples by observing non-imaginary phonons frequencies. The vibrational band structures of B 2 X 2 Zn do not have regions with the imaginary frequency which imply that the phonon stability of the samples. It should be noted that the phonon calculation is the most widely used in the stability analysis, therefore those compounds are likely to be realized experimentally.  Due to the negative formation energies with the fulfillment of the mechanical stability scheme and nonimaginary phonon spectra of all three B 2 X 2 Zn compounds, we conclude that all three compounds are mechanically, structurally, and dynamically stable. Electronic band structure calculation of B 2 Ir 2 Zn within GGA with and without the inclusion of the SOC along high-symmetry k-path Γ − X − Z − P − Γ are plotted by setting the Fermi level at 0 eV on energy scale as shown in Fig. 3. The top panel (a) shows the band structure of B 2 Ir 2 Zn compound with the inclusion of the SOC effect and the bottom panel (c) shows the calculated electronic band structure without the inclusion of the SOC effect in irreducible representation. Irreducible representation (symmorphic crystal symmetries) of band structure shows band symmetries by using different colored solid lines. There are few bands near the Fermi level. Interesting band features near the Fermi level are noted around X and P points. Since irreducible representation allows us to access each eigenvalue along the chosen kpath, we can identify connecting lines of bands and the symmetries by looking for the same colored bands for the same symmetry 42,43 .
The two linear upright crossings in the Γ-X-Z plane are identified as Dirac-like crossings. Since the SOC (some of the degenerate atomic levels split without magnetic field) have appeared as promising candidates for exotic band behaviors of Dirac materials, we perform SOC calculation to identify the topological features at the crossings. As shown in the top panel (a) in Fig. 3, it is clear that those crossings are gapped out with the inclusion of SOC. The Dirac point located at (0.136 eV) with the coordination of k=(0.499, -0.499, 0) is gaped out into two-fold degeneracy. Blue color represents the Γ 2 and Green color represents the Γ 4 in irreducible representations with space group C 2v . The Dirac point located at (0.110 eV) is gaped out into two-fold degeneracy. Blue color represents the Γ 2 and black color represents the Γ 1 in irreducible representations with space group c s . Therefore, both crossings are identified as type I twofold degenerate Dirac points. The lines at crossing points display the same and opposite slopes around ±7 eVÅ. These Dirac points are protected by the absence of SOC with the predicted electron velocity of around 1.0 × 10 16 A /s calculated by (1/ )dE/dk, which is similar to the experimentally measured velocity of Cd 3 As 2 44 . The two crossings located at the Z-P-Γ plane (left at -0.128eV and right at -0.061eV) are gapped out again into two-fold degeneracy with band inversion. At the left crossing, black and blue bands represent the Γ 1 and Γ 2 in an orderly space group C s . At the right crossing, black and blue colors represent the Γ 1 and Γ 2 , respectively with space group C s . Therefore, both crossings at around -0.128eV and -0.061ev are identified as type I two-fold degenerate Dirac points. The lines at crossing points display the same and opposite slopes around ±5 eVÅ. These Dirac points are protected by the absence of SOC but gap out by the presence of SOC. It is also observed that there is a band inversion during the gap out. All the band crossing near Fermi energy are identified mainly as Ir−d orbitals with barely hybridized with all−p orbitals by using the fat band orientations.
The results of total and partial DOS of B 2 Ir 2 Zn provide valuable information about the origin of bands with contributions from each atom and each orbital. Total and atom-projected DOS calculations without SOC effect for B 2 Ir 2 Zn display in Fig. 3(c). Black solid lines represent the total DOS from all atoms of the B 2 Ir 2 Zn compound. Red, blue, and green represent the atom-projected total DOS for Ir, B, and Zn atoms respectively. The total DOS at the Fermi level is around 0.25 states per eV per unit cell and is dominated by Ir atom DOS. Further, total DOS at the Fermi level is dominated by Ir−d orbitals denoted by red dotted lines. This is agreed with the fat band orientations we discussed in the bands' diagram in Fig. 3.
The calculated Fermi surfaces display in Fig. 3 (d) by indicating cone-shaped Fermi pockets at the corner of the Fermi surface within the BZ. This work is extended to the classification of band structures in a different direction including crystal point group symmetries. Recent studies of topological materials reveal a large class of materials with gapless surface states. Among them, topological insulating states protected by crystal symmetries, rather than time-reversal symmetry are introduced as topological crystalline insu-lators. This motivates us to investigate topological crystalline insulators' behaviors of the compound. As we show in Fig. 1 (c), the projection of N-X-Γ in bulk BZ represents from X-M -Γ in (001) surface BZ. In Fig. 1, (a) again represents the (001) surface BZ. Calculated surface states on k-Path labeled in green display on the right top to bottom by using red solid lines for different termination layers as in (b-i), (c-i), and (d-i) of Fig. 4. To study the surface states compared to bulk bands, we plot the bulk band on top of the surface bands by using blue solid lines. All the calculations were performed for relativistic potentials. We found that all B, Ir, and Zn termination layers show the same features of surface states although they show different band densities. Additionally, we investigate clear surface states on the X-N plane since there is no crossing observed in bulk bands. The surface states on the X-Γ plane are barely recognized due to small gap openings of bulk bands. As you can see we are unable to project the interested exact bulk band k-path in Fig. 3 to any surface directly. That was the reason to look at the (001) surface plane with the projection of N-X-Γ. We further studied the B 2 Ir 2 Zn by investigating the dependence of the band gap on the lattice constants. We investigate the band gap as a function of the percentage of lattice parameters. The behavior of the two crossings identified as Dirac points on Γ-X-Z plane in Fig 3 is investigated as shown in the top panel of Fig 5. It displays that there is no manifest effect to change of lattice parameters even by compressing or stretching. We were barely able to identify the gap closing of the X-Γ plane around 5% of lattice compressing, but no effect was identified on the X-Z plane. The behavior of the two crossings identified as Dirac points with band inversion Z-P-Γ path in It displays that there is a considerable effect on band gap with respect to lattice parameters by stretching. As the lattice parameters increase, the band gap of B 2 Ir 2 Zn decreases to zero and then re-opens. This gap closing signals a topological crystalline insulator at that ambient pressure which is expected to close 13% of stretching the lattice. The energy gap as a function of the percentage change of lattice parameters is displayed in Fig 5 (i) by denoting purple solid color for Z-P and brown solid color for P-Γ paths. Continuously tunable band gaps may lead to wide-ranging applications in thermoelectrics, infrared detection, and tunable electronics.   We perform the electronic band structure and DOS calculations for B 2 Rh 2 Zn and B 2 Co 2 Zn by using the same structure as B 2 Ir 2 Zn. All the calculations are done by using the optimized lattice parameters from Table . I and including the SOC effect as shown in (a) and (b) of Fig. 6. We conclude that all three compounds display the same band characteristics as discussed above for B 2 Ir 2 Zn. The energy locations of four crossing points are different than the B 2 Ir 2 Zn. It is also identified that band structure is more compressed toward the Fermi energy in Co and Rh than Ir. The DOS features around the Fermi level are investigated in three samples with and without SOC effect as shown in (c) of Fig. 6. We did not identify the manifest SOC effect in DOS by comparing solid and dotted color lines represented with and without relativistic potentials. Fig. 6 shows that Co has more DOS near Fermi energy than Rh and Rh has more DOS than Ir.

IV. CONCLUSION
In summary, newly designed B 2 X 2 Zn (X = Ir, Rh, Co) compounds show interesting topological properties and the hallmark of the topological crystalline insulator. Dirac band behaviors near the Fermi level with the effect of spin-orbit coupling can be used to predict theoretical and experimental significance on material properties. Linear band crossings near the Fermi level are discovered. Type I Dirac crossings are identified by investigating the SOC effect. We perform the surface states calculation on high symmetric crystal surfaces by projecting to the (001) plane. The gaping behavior is quantitatively investigated by applying stress and strain to the lattice volume. It shows that gap is closing by stretching the lattice and also re-opening for further stretching which is identified as the topological crystalline insulating behavior. Further, calculated elastic constants, formation energy, and phonon spectra predict the mechanical, structural, and dynamical stabilities of the compounds. The predicted electronic structures of B 2 X 2 Zn compounds, their important topological properties, and stability criteria will be useful for investigating further studies. Since these materials show high mobilities and a wide range of functionalities, it creates an extremely valuable platform for exploring advanced technology.

V. ACKNOWLEDGMENT
NH and KH acknowledge the Extreme Science and Engineering Discovery Environment (XSEDE), supported by grant number TG-PHY190050. KH and JH acknowledge the financial support from Undergraduate Prestigious Fellowships from Seton Hall University.