Design and predict tetragonal van der Waals layered quantum materials of MPd5I2 (M=Ga, In and 3d transition metals)

Design and discovery of novel 2D van der Waals (vdW) layered topological and magnetic materials containing 3d transition metals (TM) are essential for advancing both fundamental science and applied technologies. For example, understanding and manipulating itinerant 3d TM magnetism in intermetallics are crucial to understand many emergent states, including superconductivity^1,2,3,4,5. At the same time, rare-earth-free permanent magnets with 3d TM having large coercivity and saturation moment are sought after for renewable energy technologies^{6,7,8,9,10,11}. When combining these two features in the same material, it is rare to find families of magnetic 2D vdW materials with only a handful examples such as CrX₃ (X = I, Br, Cl)^{12,13,14,15,16,17}, VI₃¹⁸, Cr₂Si₂Te₆¹⁹, Cr₂Ge₂Te₆^20,21, Fe₃GeTe₂²², and MnBi_2nTe_3n+1(MBT)²³. The MBT systems are particularly interesting because most of them are intrinsic antiferromagnetic (AF) topological insulators (TI) that can be exfoliated to a few layers. The discovery of these magnetic 2D vdW materials have galvanized extensive research activities to study their unique properties.

The recent discovery^{24,25,26,27,28,29} of magnetic materials with the structural motif of -MPd₅– and -MPt₅– slabs in the tetragonal anti-CeCo₅In structure has attracted wide attention for incorporating 3d TM magnetism of Cr, Mn, and Fe into slabs of Pd and Pt with large spin-orbit coupling (SOC). Ferromagnetic (FM) CrPt₅P²⁸ and MnPd₅P²⁹ have shown large magneto-anisotropy energy (MAE) for possible rare-earth-free permanent magnet applications, however, the easy axis is in-plane. These 1-5-1 compounds with separation of three-atomic-layered -MPd₅– and -MPt₅– slabs by a single layer of P or As are not exfoliable due to the strong bonding of P or As to Pd or Pt on both sides. But it hints at a new way to design 2D vdW materials by inserting more anion layers to separate the well-structured atomic metal slabs. On the other hand, non-magnetic (NM) 2D vdW materials that can be exfoliated to a single layer are in high demand to realize a range of emergent quantum states from superconductivity³⁰ to fractional quantum anomalous Hall effect^31,32,33 by twisting the few-layer systems, as well as charge density wave (CDW) induced quantum spin Hall (QSH) effect³⁴ by gating. These new tetragonal vdW-layered materials predicted here can potentially provide a platform for the twisted few-layer systems with a tetragonal lattice³⁵ besides the twisted hexagonal and orthorhombic lattices.

Here we report the design and prediction of 2D vdW-layered I-separated -MPd₅– slabs in the body-centered tetragonal crystal structure of MPd₅I₂ with (M= Ga, In, and 3d TMs) by first-principles calculations. The phase stability calculations based on density functional theory^36,37 (DFT) show that the NM ones are on the ground state (GS) hull and thermodynamically stable, similar to the already existing AlPd₅I₂³⁸, and a few magnetic ones are close to the GS hull and thus metastable. These materials have either topologically non-trivial band structures or interesting magnetic properties. In the two band manifolds near the Fermi energy (E_F) separated by the top valence band from above and below, respectively, the NM compounds host a pair of bulk Dirac points (BDPs) just above and likely an additional strong TI state below the top valence band. While TiPd₅I₂ is a Dirac semimetal (DSM) with a relatively clean Fermi surface (FS), InPd₅I₂ has a half-filled top valence band with the surface Dirac point (SDP) from the strong TI state appearing at the E_F. This combined feature results in the surface states with both the BDP projection and SDP appearing at the Γ point in a narrow energy window of 0.1–0.2 eV around the E_F in TiPd₅I₂ and InPd₅I₂. The single layer (1 L) TiPd₅I₂ is a QSH insulator with an indirect global band gap of 30 meV or a narrow-gap semiconductor with a small direct gap of 100 meV depending on the DFT exchange-correlation (XC) functional. The 1 L InPd₅I₂ hosts two QSH states just above and below the half-filled top valence band with the band inversions involving pieces of flat bands. Among the magnetic ones, CrPd₅I₂ stands out for having the highest MAE of 2.88 meV/f.u. with the FM easy axis along the c-axis giving the desirable out-of-plane magnetic anisotropy needed for a rare-earth-free permanent magnet. The large MAE remains for the 1 L FM CrPd₅I₂. From the many band crossings around E_F in both bulk and 1 L CrPd₅I₂, there are Weyl nodal lines protected by the horizontal mirror plane as the FM magnetic moment is along the c-axis.

DFT^36,37 calculations have been performed with different XC functionals using a plane-wave basis set and projector augmented wave method³⁹, as implemented in the Vienna ab initio simulation package^40,41 (VASP). Besides PBEsol⁴², for vdW interaction, we have used vdW density functional (vdW-DF) of optB86b⁴³ and the most recent r2SCAN + rVV10⁴⁴. Band structures have been calculated with SOC and the results have also been checked with modified Becke-Johnson⁴⁵ (mBJ) and HSE06⁴⁶ exchange functionals. We have used a kinetic energy cutoff of 400 eV, Γ-centered Monkhorst-Pack⁴⁷ with a k-point density of 0.025 1/Å and a Gaussian smearing of 0.05 eV. The ionic positions and unit cell vectors are fully relaxed, with the remaining absolute force on each atom being less than 1 × 10^–2 eV/Å. For the 1 L structures, ionic relaxation is allowed in all directions, while the lattice vectors are only relaxed along the in-plane directions (x–y) with a 20 Å vacuum inserted along the out-of-plane (z) direction. In magnetic systems, calculations are initialized with a magnetic moment of 5 μ_B for transition metals and 0 μ_B for other elements in both FM and AF configurations. MAE calculations are performed by changing the global spin quantization axis from the z to x direction. Phase stability analysis is conducted using the convex hull algorithm^48,49, which is implemented in the Pymatgen package^50,51. The high throughput calculations on electronic structure and thermodynamics have been carried out in the workflow of high throughput electronic structure pakage (HTESP)⁵². To calculate Wilson loop and surface spectral functions, maximally localized Wannier functions (MLWF)^53,54 and the tight-binding model have been constructed to reproduce closely the band structure within E_F±1 eV by using Group III sp, TM sd, and I p orbitals. The surface spectral functions have been calculated with the surface Green’s function methods^55,56 as implemented in WannierTools⁵⁷. The phonon band dispersions are calculated with the small displacement method as implemented in Phonopy⁵⁸.

Structural motif and phase stability

Figure 1a summarizes the design of layered compounds based on the -MPd₅– slabs. Without layer separation, the -MPd₅– slabs with shared plane boundary is in the Cu₃Au-type structure in space group (SG) Pm-3m (221) with M being surrounded by 12 nearest-neighbor Pd in a locally close-packed face-centered cubic (FCC) structure. As shown by the successful synthesis of MPd₅P and MPd₅As compounds^27,29, an anionic layer can be inserted between the -MPd₅– slabs to make new structures in the SG P4/mmm (123). With P/As in the nominal valence of 3–, it is possible to have two I to replace P, and more importantly, it becomes vdW layered, as shown by the existence of AlPd₅I₂³⁸, the only reported compound in this structure so far. It replaces P with two I and shifts the -MPd₅– slabs in-plane for every other slab, which results in a body-centered tetragonal (tI16) structure in SG I4/mmm (139) with a large distance of 3.89 Å between the I atoms of neighboring slabs. A very recent experimental study⁵⁹ on single crystal AlPd₅I₂ has shown that, indeed it can be exfoliated to a single AlPd₅I₂ slab with interesting properties. Here we design and predict new vdW-layered compounds in this structure by replacing Al with other group III elements and also 3d TMs via high throughput DFT calculations.

a Crystal structures of MPd₅ motif in the three-atomic-layer slab with increasing distance between slabs as in MPd₃ of Pm-3m (221) in Cu₃Au-type, MPd₅P of P4/mmm (123) and MPd₅I₂ of I4/mmm (139). The atomic species are shown in different colors and labeled accordingly. b Projected density of states (PDOS) on atomic orbitals of non-magnetic InPd₅I₂, TiPd₅I₂, and ferromagnetic (FM) CrPd₅I₂. c Convex hull energy (E_h) of MPd₅I₂ and MPd₅I (M=Al, Ga, In, Ti, V, Cr, Mn, Fe, Co, and Ni) calculated in density functional theory (DFT) with PBEsol and optB86b exchange-correlation functionals. The ground state magnetic configurations with easy axis are drawn with the listed moment size on the magnetic ions (from V to Ni) and magneto-anisotropy energy (MAE). FM and antiferromagnetic (AF) are shown in red and blue shaded squares, respectively.

Full size image

As plotted in Fig. 1b for the projected density of states (PDOS) of NM InPd₅I₂, NM TiPd₅I₂, and FM CrPd₅I₂, most of the I 5p orbitals hybridize with the bottom of Pd 4d orbitals from –6 to –4 eV to form the p-d bonding states. The p-d anti-bonding states are pushed to above 1 eV as empty states to gain more cohesion for the ternary compounds with additional bonding hybridization between Pd 5s (not shown) and I 5p orbitals in the same low-energy range. In contrast to I 5p, the p orbitals of group III elements at M site, for example In, mostly hybridizes with Pd 4 d states at a higher energy range (–4 to –2 eV), reflecting their electron positive character. But the states near the E_F in InPd₅I₂ are dominated by Pd 4d and I 5p orbitals. Next, moving to 3d TM, Ti 3d orbitals hybridize extensively with Pd 4d orbitals giving the broader and lower Pd 4d-derived bands than those in InPd₅I₂ below the E_F. There is also a large empty anti-bonding DOS peak just above the E_F due to the 3d-4d hybridization. Then, for Cr with two more 3d electrons, these empty states get partially filled and induce a large exchange interaction, as shown by the two splitting DOS peaks at 0 and +2 eV. This interaction gives a sizable magnetic moment of 2.80 ({mu }_{B}) on Cr to prefer FM, and importantly, the easy axis is along the c-axis with a large MAE of 2.88 meV/f.u. The PDOS for other magnetic 3d TM MPd₅I₂ are similar in terms of band hybridizations. The NM TiPd₅I₂ is an interesting case with the right number of valence electrons that the anti-bonding states are almost completely empty, giving a minimum DOS at E_F to form a semimetal, whose topological band structure will be detailed later.

The GS convex hull energy (E_h) for all the MPd₅I₂ compounds studied are plotted in Fig. 1c for PBEsol⁴² and also with the vdW exchange functional of optB86b⁴³. To confirm that the introduction of I prefers two anion layers instead of one, besides the MPd₅I₂ in the I4/mmm structure, we have also calculated the hypothetical MPd₅I in the P4/mmm structure. As shown in Fig. 1c, E_h for MPd₅I are all above 0.10 eV/atom, much higher than MPd₅I₂, confirming the qualitative argument that the broken metallic interactions in separating -MPd₅– slabs need to be compensated by a strong ionic interaction with enough anionic valence. For the group III elements, Al, Ga, and In, the E_h of MPd₅I₂ are all zero for PBEsol and slightly below 0.01 eV/atom for optB86b, showing they are all thermodynamically stable, because AlPd₅I₂ has already been found in experiment³⁸, although optB86b gives a small positive E_h of 0.006 eV/atom. InPd₅I₂ has the smallest E_h, thus the most stable among the group III compounds. For the 3d TMs, first TiPd₅I₂ is quite stable on the GS hull even for optB86b and we found it is a DSM with a clean FS. Next for V and Cr, E_h becomes positive, but smaller than 0.10 eV/atom. Then E_h decreases for Mn and Fe at the middle of 3d TM series, and increases again for the late 3d Co and Ni, which are the least stable among the 3d TM MPd₅I₂. Overall, optB86b gives a higher E_h than PBEsol, showing a systematic shift between the different XC functionals. But the trends for the variation in E_h across the whole series for both MPd₅I₂ and MPd₅I are the same for different XC functionals, showing the results are well converged.

The magnetic properties across the 3d TM series for MPd5I2 are quite interesting and tabulated at the top of Fig. 1c. Except for Ti being NM, the magnetic moment size increases first, starting with V, reaching the maximum of 3.92 ({mu }_{B}) for Mn before decreases at the end of the series for Co and Ni. With the gradual filling of the 3d orbitals, V, Cr, and Mn prefer FM, while Fe, Co, and Ni prefer AF. Importantly, both V and Cr prefer an easy axis along the c-axis, with Cr having the largest MAE of 2.88 meV/f.u., much higher than the 0.30 meV/f.u. for VPd₅I₂ and the rest. In contrast, FM Mn prefers the in-plane easy-axis, although with the largest moment. Then, for AF, first Fe moments prefer in-plane and then Co and Ni prefer out-of-plane directions.

To study phase stability and construct GS hull, all the existing binary and ternary compounds, together with the elemental ones in the ternary phase diagrams, have been computed, and their stability are calculated via different possible reaction paths. Although PBEsol gives AlPd₅I₂ on the GS hull, agreeing with experiment in Fig. 1c, Fig. 2a shows that the PBEsol-calculated volume per atom is underestimated when compared to the available experimental data. Then as shown in Fig. 2b, this underestimation of volume can be improved by using optB86b exchange functional. Also, to explicitly include vdW interaction for the 1-5-2 compounds, we have chosen optB86b vdW exchange functional for the phase stability plots in Fig. 2, as well as the band structure and magnetic property calculations. The GS hulls with PBEsol are similar and can be found in Supplementary Fig. 1.

a Volume per atom of the fully relaxed elemental, binary and ternary compounds in PBEsol are compared to available experimental data with the mean absolute percentage errors (MAPE) listed. b Same comparison for optB86b. c–l The calculated phase stability and structural energies of MPd₅I₂ (M=Al, Ga, In, Ti, V, Cr, Mn, Fe, Co, and Ni) in optB86b. The compounds on the ground state hull are labeled as green dots. The MPd₅I₂ are indicated with arrows and their respective hull energies are listed in parenthesis.

Full size image

As shown in Fig. 2, for the Pd-I binary, there is only one stable line compound of PdI₂. For M-I binaries, there are many stable line compounds for Ga, In, Ti, and V. For the rest, there is only one stable binary M-I compound, including CrI₃ for Cr-I. For M-Pd binaries, Al, Ga, In, Ti, V, and Mn have many stable line compounds. Interestingly, for the -MPd₅– motif in MPd₃ or the Cu₃Au-type, this binary line compound structure exists for In, Ti, V, Cr, Mn, and Fe with Pd. But for Ni and Co, they form a random alloy or solid solution with Pd. Because both CoPd₃ and NiPd₃ are only slightly above the GS hull at 0.06 eV/atom, they can be used as good representatives for the binary solid solution, and we include them in the phase stability calculation of MPd₅I₂. Among the calculated MPd₅I₂, the NM compounds are more stable than the magnetic ones. Given that AlPd₅I₂ with the E_h of 0.006 eV/atom in optB86b has already been synthesized, we predict the existence of GaPd₅I₂, InPd₅I₂, and TiPd₅I₂ because they are on the GS hull. For magnetic ones, Mn and Fe have the smallest E_h above the GS hull and then followed by V and Cr. Considering the approximations used in DFT calculations, we predict these four magnetic ternaries are metastable, also because the binary MPd₃ line compounds with the -MPd₅– motif are stable and found in experiments. Lastly, for Co and Ni, they have the largest E_h above GS hull even with PBEsol, and due to the solid solution of MPd₃, it is also possible to form a solid solution for the ternary compounds. We predict these two are possible ternaries but with solid solution tendency, which needs to be further studied in the future.

From an experimental viewpoint, these MPd₅I₂ compounds are much more challenging to synthesize than MPd₅P and MPd₅As because of the higher vapor pressure or lower sublimation temperature of I than P and As. We have also studied the dynamical stability of these predicted compounds by calculating the phonon band dispersion. For example, the phonon band dispersions of both bulk and 1 L structures for TiPd₅I₂, InPd₅I₂, and CrPd₅I₂ are plotted in Supplementary Fig. 2. The absence of imaginary phonon modes shows that these predicted compounds in bulk and 1 L structures are dynamically stable.

Topological features of non-magnetic 1-5-2 compounds

For the band structures of NM MPd₅I₂, we chose TiPd₅I₂ and InPd₅I₂ to present the topological band features of both the bulk and 1 L structures. Bulk band structures of other MPd₅I₂ can be found in Supplementary Fig. 3. Figure 3a plots the band structure of bulk TiPd₅I₂ without SOC with the body-centered tetragonal Brillouin zone (BZ) and high symmetry k-points shown in Fig. 3c. The highest valence band (N) and band below (N–2) according to simple filling are shown in red and blue, respectively, with N being the number of valence electrons. Above the highest valence band, there is a sizable gap in most of the BZ, except for around the Z point in the (Gamma)-Z, Z-S₁, and Y₁-Z directions. Along the (Gamma)-Z direction, band N and N–2 are degenerate, giving triple degeneracy (or six-fold including spin) at the crossing point between the highest valence and lowest conduction band in middle of the (Gamma)-Z direction. From the triple degeneracy point to the Z point, a doubly degenerated nodal line segment appears, as also shown in Fig. 3d by plotting the zero-gap k-points in the whole BZ. The crossings along the Z-S₁ and Y₁-Z directions are parts of the nodal line loops around the Z point on the (110) and (1–10) planes, which are protected by the diagonal mirror symmetries.

a Band structure of TiPd₅I₂ without spin-orbit coupling (SOC) with the highest valence band (N) and band below (N-2) shown in red and blue, respectively. Green shade stands for I p_z orbital projection. b Band structure of TiPd₅I₂ with SOC highlighting the top two valence bands and the same orbital projection as in (a). c Bulk Brillouin zone (BZ) with high symmetry k-points and those on (001) surface BZ are labeled. The bulk Dirac points (BDP) along (varGamma)-Z are indicated by red dots. d Nodal lines in TiPd₅I₂ without SOC between the highest valence and lowest conduction bands. e Bands zoomed in along (varGamma)-Z around the BDP between the highest valence and lowest conduction band. The green (gray) shade is for I p_z (Pd d_xz and d_yz) orbital. f Bands zoomed in along the (Z)–({S}_{1}) direction with a small gap opening. The orange (cyan) shade is for Pd d_yz (Ti d_xz) orbital. g Wilson loop of Wannier charge centers (WCC) on the k_z = 0.5 plane between band N–2 and N. h Surface spectral function along (bar{X})–(bar{varGamma })–(bar{M}) on (001). i (001) 2D surface Fermi surface at Fermi energy (E_F) + 0.045 eV and j E_F-0.185 eV with spin-texture shown in orange arrows.

Full size image

With SOC, as plotted in Fig. 3b, the orbital degeneracy for the top two valence bands along the (Gamma)-Z direction is lifted, and the nodal loops are all gapped out, except for the crossing between the highest valence and lowest conduction band along the (Gamma)-Z, forming a BDP as protected by the fourfold rotational symmetry. Because of the time-reversal symmetry (TRS) and inversion symmetry, each band is still doubly degenerated. The BDP is zoomed in Fig. 3e along the (Gamma)-Z direction showing the zero gap and the switching between I p_z and Pd d_xz/d_yz orbitals with the 2-dimensional irreducible representations of (Gamma)₉ and (Gamma)₆. The BDPs are at the momentum energy of (0, 0, ±0.1499 Å^–1; E_F + 0.0236 eV), as also shown by the red dots in Fig. 3c. The band dispersion of TiPd₅I₂ is zoomed along the Z-S₁ direction in Fig. 3f to show the small SOC-induced gap.

Additionally, for the highest valence band N (red), it is also gapped from below by the next valence band N–2 (blue). The lower branch of the band N–2 along the (Gamma)-Z direction forms a band inversion region around the Z point with the top valence band N. Wilson loop calculation shows that this N–2 band manifold hosts a strong TI (STI) state with the Fu-Kane⁶⁰ topological index of (1;001). The non-trivial Z₂ number is shown by the Wilson loop of Wannier charge centers (WCC) on the k_z = 0.5 plane in Fig. 3g with the odd number of crossings by the dashed line with WCC. Calculations with the more recent r2SCAN + rVV10⁴⁴ XC functional give similar band features (see Supplementary Fig. 4) for both DSM and STI in TiPd₅I₂. For mBJ⁴⁵ functional, the BDP still remains and predicts a DSM, despite most of the valence bands being pushed lower and conduction bands higher in energy. But the extra band inversion at the Z point is lifted between bands N–2 and N. The existence of this extra band inversion for STI or not in TiPd₅I₂ will be the features that need to be verified in the experiment.

To demonstrate the non-trivial band structure of TiPd₅I₂ with both a BDP above and an STI below the highest valence band, we have calculated the (001) surface spectral functions using the Wannier functions. On (001) surface, the projections of the BDP at ±k_z onto the same (bar{Gamma }) point stand out nicely around the E_F because of no overlap with other bulk band projection for the clean FS. There are topological surface states (TSS) stemming from the BDP projection, as seen in Fig. 3h at E_F + 0.03 eV. Below that at E_F–0.2 eV is the SDP from the STI of the N–2 band manifold also shown clearly even though on top of the other bulk band projections. The spin-texture of the surface Dirac cone is plotted at E_F–0.185 eV in Fig. 3j, confirming the spin-momentum locking of the surface Dirac cone. The spin-momentum locking of the TSS stemming from the BDP projection is shown in Fig. 3i at E_F + 0.045 eV. Thus, bulk TiPd₅I₂ is a DSM with a clean FS and also possibly hosts an STI just below the highest valence band.

When the vdW-layered TiPd₅I₂ is exfoliated down to 1 L, the band structure is plotted in Fig. 4a. A large band gap exists in most of the BZ, while a small gap appears along the (Gamma)-X direction, which is projected from the Z-S₁ direction from the bulk band structure. Overall there is an indirect global band gap of 30 meV between the valence band maximum at X and conduction band minimum at (Gamma) point. The Wilson loop calculation of the highest valence band manifold (red) in Fig. 4b indicates a QSH with an odd number of crossings of WCC. In contrast, for the N–2 band manifold (blue), the even number of crossings in the Wilson loop in Fig. 4c shows it is topologically trivial. The edge spectral functions are plotted in Fig. 4d, e for the different TiPd- and PdI-terminations, respectively. The topological edge states connect the gapped valence with conduction band projections to form a TRS-protected edge Dirac points (EDP) at the (bar{Gamma }) point. While the EDP on the TiPd-termination is inside the QSH gap, that on PdI-termination is merged into the valence band projection at E_F–0.06 eV. With E_F cutting through the QSH gap, the spin-momentum locked edge states are unavoidable from TRS topological protection. Calculations with r2SCAN + rVV10 XC functional give a similar band inversion feature at the (Gamma) point for QSH (see Supplementary Fig. 4). In contrast, HSE06⁴⁶ functional pushes the valence band lower and conduction bands higher in energy, and lift the band inversion and changes the indirect band gap to a direct gap of just 100 meV at the (Gamma) point (see Supplementary Fig. 4d). Such a small gap size can be potentially tuned to close and reopen by strain to induce the band inversion for a topological phase transition to realize a critical 2D Dirac point at the (Gamma) point and then a QSH.

a Band structure of 1 L TiPd₅I₂ with the highest valence band (N) and band below (N-2) shown in red and blue, respectively. b Wilson loop of Wannier charge center (WCC) between band N and N + 2. c Wilson loop between band N–2 and N. d Edge spectral functions along (bar{X})–(bar{varGamma })–(bar{X}) on TiPd- and e PdI-terminations.

Full size image

Next for InPd₅I₂ as an example from group III MPd₅I₂, its bulk band structure with SOC is plotted in Fig. 5a. Due to TRS and inversion symmetry, each band is doubly degenerated. Because of the odd number of electrons, the highest valence band (red in Fig. 5a) is only half-filled, indicated by E_F sitting right in the middle of the bandwidth. But it is still meaningful to discuss the topological features of the band manifolds above and below this half-filled top valence band with a finite FS. Similarly to TiPd₅I₂, between the highest valence and lowest conduction band of InPd₅I₂, there is only one crossing along the (Gamma)-Z direction as protected by the fourfold rotational symmetry for a pair of BDP. The BDP is zoomed in Fig. 5b shown by the projection on I p_z and Pd d_xz/d_yz orbitals with the two-dimensional irreducible representations of (Gamma)₉ and (Gamma)₆. The BDP is at the momentum energy of (0, 0, ±0.0406 Å^–1; E_F + 0.0955 eV). In contrast, between the highest and the next valence band (blue), there is no band crossing. For such gapped band manifolds, the Fu-Kane topological index has been calculated as (1;001) showing a STI state with the Wilson loop of WCC at the k_z = 0.5 plane being plotted in Fig. 5c. To confirm the STI with surface spectral function on (001), the surface Dirac cone is shown clearly in Fig. 5d with the SDP right at the E_F inside the projected bulk gap. The spin-texture of the surface states at the E_F is plotted in Fig. 5e, confirming the spin-momentum locking topological feature without the overlap with bulk band projection. In contrast, the BDP projection at E_F + 0.10 eV on (001) surface in Fig. 5d is buried inside the other bulk band projection and shows no TSS, unlike those in TiPd₅I₂ with a clean FS. Thus, group III MPd₅I₂ hosts a BDP above and a STI below the half-filled highest valence band. For InPd₅I₂, both the SDP and BDP projection appear at the (bar{Gamma }) point on (001), and they are also within an energy window of 0.1 eV, with the SDP being right at the E_F and the BDP just above the E_F.

a Band structure of bulk InPd₅I₂ with the highest valence band (N) and band below (N-2) shown in red and blue, respectively. b Bands zoomed in along (varGamma)–Z around the Dirac point (DP) between the highest valence and lowest conduction band. The green (gray) shade is for I p_z (Pd d_xz and d_yz) orbital. c Wilson loop of Wannier charge center (WCC) on the k_z = 0.5 plane between band N–2 and N. d Surface spectral function along (bar{X})–(bar{varGamma })–(bar{M}) on (001). e (001) 2D surface Fermi surface at Fermi energy (E_F) with spin-texture shown in orange arrows. f Band structure of 1 L InPd₅I₂ with the highest valence band (N) and band below (N-2) shown in red and blue, respectively. g Wilson loop of WCC between band N and N + 2. h Wilson loop between band N–2 and N. i Edge spectral functions along (bar{X})–(bar{varGamma })–(bar{X}) on InPd- and j PdI-terminations.

Full size image

The band structure of 1 L InPd₅I₂ is plotted in Fig. 5f. Again, due to the odd number of electrons, the top valence band is half-filled with E_F sitting in the middle of the bandwidth. But the valence band is continuously gapped from both below and above with band inversions, so topological properties can be calculated. The Wilson loop calculations of the band manifolds in Fig. 5g, h show the odd number of crossings of WCC, confirming it hosts two QSH states. The edge spectral functions are plotted in Fig. 5i, j for two different terminations, which are rather similar. The TRS-protected EDP at E_F + 0.35 eV is for the upper QSH, and the EDP at E_F is for the lower QSH. The topological edge states of the upper QSH around the (bar{Gamma }) point stand out, which is accessible, when the E_F can usually be tuned by doping and gating. Calculations with r2SCAN + rVV10, mBJ, and HSE06 XC functionals all give the similar two band inversions and topological features (see Supplementary Fig. 5), which are much less affected than TiPd₅I₂, because InPd₅I₂ bands are more metallic from the half-filled top valence band than TiPd₅I₂. So 1 L InPd₅I₂ hosts two QSH states in a tetragonal structure despite being a metal. Additionally, the band inversions for the two QSH states around the E_F in Fig. 5f involve pieces of flat bands from the orbital-decorated square lattice as recently proposed⁵⁹. Together with the QSH insulator with a small indirect band gap or a narrow-gap semiconductor in 1 L TiPd₅I₂, the few-layer tetragonal systems of these exfoliable 1-5-2 compounds will be an interesting playground for emergent quantum states in future studies.

Magneto-anisotropy of CrPd₅I₂

Among the magnetic 3d TM MPd₅I₂ compounds, CrPd₅I₂ has the largest MAE and also the easy-axis is along the c-axis. First, without SOC, the spin DOS of the bulk and 1 L CrPd₅I₂ are plotted in Fig. 6a, b, respectively. While the spin down (majority) forms a pseudo gap near the E_F, the spin up (minority) has a local DOS maximum at E_F. Via the hybridization between Cr-3d and Pd-4d orbitals to form bonding and anti-bonding states just above E_F, the exchange splitting interaction gives a sizable magnetic momentum of 2.8 ({mu }_{B}) on the Cr atom. The 1 L spin DOS in Fig. 6b is similar and has narrower and sharper peaks than the bulk in Fig. 6a due to the smaller band dispersion from the reduced interlayer interactions. To analyze the origin of the large MAE in CrPd₅I₂, we have calculated the k-point-resolved MAE over the entire BZ by fixing the magnetic charge density but rotating the magnetic axis from [001] to [100] with SOC. As shown in Fig. 6c, d for the (Delta)MAE = ± 0.03 meV/f.u., respectively, the positive MAE contribution (favoring the c-axis) in Fig. 6c is mostly around the (Gamma) and X points. In contrast, the negative MAE contribution (favoring in-plane) in Fig. 6d is mostly from the Z, S point and also half way between the (Gamma) and X points. Going to the 1 L CrPd₅I₂, the whole bandwidth is reduced, but most of the hybridization peaks remain in the same positions, which shows that the 1 L can retain the chemical stability. The magnetic moment does not change much, the MAE is still as high as 2.47 meV/f.u.

a, b Up (minority) and down (majority) density of state (DOS) for bulk and 1 L CrPd₅I₂ without spin-orbit coupling (SOC). c, d k-point resolved magneto-anisotropy energy (MAE) isosurface at ±0.03 meV/f.u. for bulk CrPd₅I₂ with switching the magnetic axis from [001] to [100] with SOC. e, f Band structure with SOC for bulk and 1 L CrPd₅I₂ with the highest valence band (N) shown in red. g, h Weyl nodal lines of bulk and 1 L CrPd₅I₂ between band N and N + 1.

Full size image

The band structures of FM bulk and 1 L CrPd₅I₂ with SOC, are plotted in Fig. 6e, f, respectively. The band double degeneracies are all lifted with the top valence band shown in red. There are many bands crossing the E_F, indicating a more complicated FS than the NM TiPd₅I₂ and InPd₅I₂. These many crossings form twofold degenerated Weyl nodal lines as plotted in Fig. 6g, h. For the FM bulk CrPd₅I₂, besides the main Weyl nodal loops on the k_z = ±0.5 plane, there are also loops around the X points. For the FM 1 L CrPd₅I₂, there are three Weyl nodal loops, one around the X point and two around the M point. These Weyl nodal lines are within E_F ± 0.2 eV.

The high MAE in CrPd₅I₂ reflects the unique structural motif of the -MPd₅– slab, where each moment-bearing 3d TM atom is surrounded by Pd with much larger SOC strength. The distance among the 3d TM atoms is much larger than that in elemental solids. The magnetic coupling among the 3d TM atoms are through Pd with a larger SOC and itself is near the Stoner magnetic instability. Such a combination gives a range of magnetic configurations in MPd₅I₂. With the gradual filling of the 3d orbitals. V, Cr and Mn prefer FM, while Fe, Co and Ni prefer AF. Importantly, both V and Cr prefer easy axis along the c-axis with Cr having the largest MAE of 2.88 meV/f.u. In contrast, FM Mn prefers the in-plane easy axis, although with a larger moment. Then, for AF, first Fe prefers in-plane, and then Co and Ni prefer out-of-plane. Among the FM ones, with the largest MAE and easy-axis being out-of-plane, CrPd₅I₂ can give a large coercivity field, which is attractive for developing rare-earth-free permanent magnets.

Using high throughput density functional theory calculations, we have explored the phase stability, topological and magnetic properties of MPd₅I₂ compounds (M=Ga, In and 3d TM), a family of -MPd₅– slabs separated by two layers of I to design and predict new vdW-layered quantum materials with tetragonal structure. After confirming the existing AlPd₅I₂ is on the ground state (GS) hull, we find non-magnetic (NM) compounds with M=Ga, In, and Ti are also on the GS hull and thermodynamically stable. For the magnetic ones with 3d TM, we find V, Cr, Mn, and Fe are not far above the GS hull and metastable, given the existence of the binary structures with -MPd₅– slab in the cubic MPd₃ structure. For Co and Ni, the hull energy is the largest and also these MPd₃ form random alloys. Using TiPd₅I₂ and InPd₅I₂ as examples, we show that the NM MPd₅I₂ host a pair of bulk Dirac points for the band manifolds just above the highest valence band and also likely an extra strong topological insulator (TI) state from below, respectively. While TiPd₅I₂ is a Dirac semimetal with a mostly clean Fermi surface, InPd₅I₂ has a half-filled top valence band with the surface Dirac point from the strong TI appearing at the Fermi energy (E_F). This combination gives the (001) surface hosting both a surface Dirac point and a bulk Dirac projection just 0.1–0.2 eV separation at the (bar{Gamma }) point. From different exchange-correlation functionals, the 1 L TiPd₅I₂ is either a quantum spin Hall (QSH) insulator with an indirect global band gap of 30 meV or a narrow-gap semiconductor with a direct gap of 100 meV, which potentially can be tuned for a topological phase transition. In contrast, the 1 L InPd₅I₂ is always a metal from half band-filling and hosts two QSH states with pieces of the flat bands near E_F. For the magnetic MPd₅I₂ with 3d TMs, the preferred magnetic ground state changes with the gradual filling of the 3d orbitals. V, Cr, and Mn prefer a ferromagnetic (FM) ground state with less or at half-filling, while Fe, Co, and Ni prefer an antiferromagnetic configuration with more than half-filling of the 3d orbitals. Interestingly both VPd₅I₂ and CrPd₅I₂ have their FM moment easy axis along the out-of-plane c-axis, a desirable feature to develop rare-earth-free permanent magnets. CrPd₅I₂ has a large MAE of 2.88 meV/f.u. Thus, our calculations predict that a group of MPd₅I₂ are synthesizable tetragonal vdW-layered quantum materials with either non-trivial topological features or highly tunable magnetic properties.