Tunable electronic structure and topological properties of LnPn (Ln=Ce, Pr, Sm, Gd, Yb; Pn=Sb, Bi)

Recently, the rare-earth monopnictide compounds LnPn have attracted considerable attention in condensed matter physics studies due to their possible topological properties. We have performed systematic first principles study of the electronic structure and band topology properties of LnPn (Ln=Ce, Pr, Sm, Gd, Yb; Pn=Sb, Bi). Assuming the f-electrons are well localized in these materials, both hybrid functional and modified Becke–Johnson calculations yield electronic structure in good agreement with experimental observations, while generalized gradient approximation calculations severely overestimate the band inversions. From Ce to Yb, a systematic reduction of band inversion with respect to the increasing Ln atomic number is observed, and Z2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\cal Z}_2$$\end{document} for CePn and YbPn are [1;000] and [0;000], respectively. In both hybrid functional and modified Becke–Johns calculations, a topologically non-trivial to trivial transition is expected around SmSb for the antimonides and around DyBi for the bismuthides. Such variation is related with lanthanide contraction, but is different from simple pressure effects. Rare-earth monopnictide compounds have attracted attention for diverse physical properties such as low temperature antiferromagnetism. Here, the electronic structure of rare-earth monopnictide compounds is studied using density-functional based first-principles methods to identify the emergence of topological properties.

T he rare-earth monopnictide compounds were discovered decades ago, including other members in the same family (e.g. CeN, CeP, etc). All the compounds crystallize in simple rocksalt face-centered cubic structure except YbBi, but they exhibit a diversity of physical properties. Most of these compounds are antiferromagnetic at low temperature [1][2][3][4] , except for the non-magnetic yittrium and lanthanum compounds, as well as the praseodymium compounds that are Van Vleck type paramagnetic 5 . Angular resolved photoemission spectroscopy (ARPES) and de Haas-van Alphen (dHvA) measurements have been performed [6][7][8][9][10][11][12][13][14] to probe the electronic structure of these compounds. Both dynamical mean-field theory (DMFT) and LDA + U calculations have suggested that the 4f-electrons are almost fully localized in Pn=Sb/Bi compounds [15][16][17][18] , especially for the bismuth compounds. For CePn, the Ce-4f orbitals slightly mixes with the Ce-5d orbitals (d-f mixing) and Pn-p orbitals (p-f mixing) 6,19 , and the former leads to Kondo-like behaviors in this low density carrier system, which almost vanishes for CeBi 16 .
In addition, LaPn system was proposed to be a topological system by first principles calculation 20 . Although the initial proposal argued that the LaN is a 3D Dirac-semimetal and all other LaPn are topological insulators, it was later argued that the band inversion was severely overestimated in the calculation, thus only LaBi is topologically non-trivial in the series 21 . Both ARPES and transport measurements have confirmed that LaSb is a topologically trivial electron-hole compensated semimetal 22 similar to YSb 23 ; while LaBi is topologically non-trivial with clear evidence of protected surface states [24][25][26][27][28] . It has also been proposed from first-principles method that topological transition from trivial to non-trivial can be realized in LaSb by applying hydraulic pressure 29 . More recently, ARPES measurements have been performed with CeSb and CeBi compounds to investigate their topological properties [30][31][32][33][34][35] , and transport evidence for possible Weyl fermion has been found in CeSb at ferromagnetic state 36 . Large magnetoresistance has also been reported in NdSb 2 . Therefore, the renewed interests in LnPn system focus on possible topological properties of these materials. However, a systematic first-principles study of rare-earth monopnictide compound in this aspect is still missing.
In this article, we present our latest first-principles study of LnPn (Ln=Ce, Pr, Sm, Gd, Yb; Pn=Sb, Bi). In particular, we compare the electronic structure obtained using Perdew-Burke-Ernzerhof (PBE), modified Becke-Johnson (mBJ), and Heyd-Scuseria-Ernzerhof hybrid (HSE06) functionals, as well as the band inversion features and the associated Z 2 indices. We find good agreement between the mBJ/HSE06 results and quantum oscillation measurements, while PBE results usually yields substantially larger β and γ band frequencies due to overestimation of band inversions. From Ce to Yb, a systematic reduction of band inversion related with lanthanide contraction is observed, and the topologically non-trivial to trivial transition takes place around SmSb for antimonides and DyBi for bismuthides, respectively. Finally, in contrast to lanthanide contraction, simple pressure effects will increase the band inversion gap.

Results
We show the crystal structure and its first Brillouine zone of LnPn in Fig. 1a, b, respectively. In Fig. 1c, we compare the lattice constants from structural relaxation with those obtained in experiments. The calculated lattice constants are slightly larger than the experimental observation, since generalized gradient approximation (GGA) tends to underestimate binding. Nevertheless, all the calculated values are within 3% error bar of density functional theory (DFT) calculations. In both calculation and experiments, the lattice constants of either bismuth compounds or antimony compounds decrease as the lanthanide element moves toward the end of the periodic table, i.e. the lanthanide contraction. Overall, the calculations correctly reproduce the trend of lattice constant change as the lanthanide element changes, demonstrating the validity of the DFT calculation. It is worth noting that to the best of our knowledge, there is no report on the rocksalt structured YbBi compound in the literature, possibly due to the small Yb 3+ radius.
Electronic band structure. To investigate the electronic structure of LnPn compounds, we first take the non-magnetic band structure of CeBi to study some general features. Without spinorbit coupling, the electron states near the Fermi level are dominated by the pnictogen p-orbitals and lanthanide t 2g -orbitals (Fig. 2a). Normally, the energy level of lanthanide t 2g -orbitals is higher than that of the pnictogen p-orbitals. However, such order can be reversed after considering the spin-orbit coupling in some of these LnPn compounds at X (π, π, 0), leading to a band crossing below E F along Γ-X (Fig. 2b). With the spin-orbit coupling, the 3 degenerate p-orbitals further split into a doubly degenerate Γ 8 state and a Γ 6 state at Γ point; while the 3 degenerate t 2g -orbitals split into a doubly degenerate Γ 8 state and a Γ 7 state at Γ. In addition, the spin-orbit coupling will gap out the band crossing between Γ and X due to the band inversion. As a result, one can define a partial gap at the electron occupation in the whole Brillouin zone for these materials. Therefore, the Z 2 indices are well defined for these materials. With spin-orbit coupling, there are 3 bands crossing the Fermi level, among which two of them are hole-like near Γ and another electron-like one near X point. As the electronic states are close to E F only around Γ and X, and the band topology is related with the anti-crossing between Γ and X, we shall focus on the band dispersion along this direction. We have also performed calculations with itinerant f-electron in both non-magnetic state and antiferromagnetic state for CeBi.
By taking the f-orbitals into valence, they become very itinerant and dominate the density of states near E F . This is a known artifact of local-density approximation due to lack of proper electron-correlation effect when it is applied to open-shell f compounds. To obtain the correct band structure, we have applied LDA + U method with U f = 6.1 eV in the antiferromagnetic calculations (Fig. 2c). The resulting band structure is very close to the non-magnetic band structure we obtained with f-core calculations, except that the antiferromagnetic band structure is folded (since the antiferromagnetic Brillouin zone is half of the original one), and that an extra flat f-band appears at~3 eV below E F . Both the inversion gap and the anti-crossing feature are present, and the size of the inversion gap is similar.
In Fig. 3, we show the electronic band structures of LnSb compounds. For CeSb, the band inversion occurs for t 2g orbitals and p orbitals, resulting in a large inversion gap of ≈400 meV at X under PBE approximation. However, as PBE band energies are not quasiparticle energies, and that PBE band widths are General features of LnPn band structure. a, b Electronic band structure of CeBi in non-magnetic state a without and b with spin-orbit coupling. The orbital character weights are represented by the size of the circles. The blue circles are Bi-6p orbitals, red circles are Ce-5d t 2g orbitals, and green circles are Ce-5d e g orbitals. c Comparison between f-core non-magnetic state CeBi band structure in the folded Brillouin zone (red lines) and LDA + U antiferromagnetic (AFM) state CeBi band structure (blue lines). Notice that in AFM state, the Brillouin zone is folded, so the band structure is more complicated. Nevertheless, the inversion gap is still present and the anti-crossing feature can be identified along Γ-X (red-circled area) overestimated, the inverted gap sizes in PBE are also overestimated. Therefore, we have also performed calculations in mBJ and HSE06 hybrid functional as well. Due to the overwhelming amount of calculation load for HSE06, the band structures of HSE06 hybrid functionals were obtained by Wannier fitting. The inversion gap size is substantially reduced (142 meV in mBJ and 92 meV in HSE06), but the band anti-crossing feature remains. From CeSb to YbSb, the inversion gap size reduces with respect to the increasing atomic number of the lanthanide element, and eventually disappears for YbSb for all functionals. For the mBJ method, the transition occurs at SmSb, where the Sm t 2g orbitals and Sb p-orbitals are nearly degenerate (<1 meV). For the HSE06 functional, the inversion gap sizes are generally smaller than those in mBJ calculations, and therefore this transition occurs closer to PrSb. While the decreasing trend of band inversion gap with heavier rare-earth element is in line with the ARPES results, quantitative deviations are still present in LnSb compounds. For example, ARPES measurements have shown that both CeSb and PrSb are close to the boundary of topologically trivial to nontrivial transitions, while SmSb is clearly on the trivial side with sizable non-inverted gap 33,35 .
Similar variation of inversion gap with respect to the increasing atomic number is also present in the bismuth compounds (Fig. 4). Due to the increased spin-orbit coupling in bismuth, the inversion gap is substantially enlarged. Under PBE approximation, the band inversion was so strong that the band from porbitals would cross the Fermi level around X in CeBi and PrBi, creating an additional hole pocket around X. However, this hole pocket is absent in either mBJ or HSE06 calculations. Using the HSE06 functional, the gap size is 550 meV in CeBi, 495 meV in PrBi, 270 meV in SmBi, 140 meV in GdBi and eventually disappears in YbBi. Therefore, a transition from anti-crossing to non-crossing band also exists in the bismuth compounds, presumably between GdBi and YbBi. In fact, we have located the position of this transition closest to DyBi, where the direct gap between Dy-t 2g and Bi-6p orbitals is ≈8 meV (5 meV) in mBJ (HSE06) calculations.
In addition to the variation of inversion gap, the energy level of the Γ 6 state at Γ point exhibits some systematic changes. Under PBE approximation, the Γ 6 state is slightly below E F in CeSb and slightly above E F in PrSb. It rises as one move toward YbSb. Therefore, an additional hole pocket would appear and gradually increases as one increases the atomic number under PBE. However, for both mBJ and HSE06 calculations, the Γ 6 state is much lower than it is in PBE, and is always below E F . Nevertheless, the trend remains visible, and the Γ 6 state in GdSb is much closer to E F than it is in CeSb. Such change in Γ 6 state is also present in bismuth compounds. However, the Γ 6 state is more than 1 eV below E F in these compounds, and therefore is irrelevant. Both mBJ and HSE06 results agree very well with ARPES measurements for CeBi, PrBi, SmBi, and GdBi 34 .
Using the maximally projected Wannier function method, we have also calculated the angular dependence of dHvA frequencies  of these materials 37 (Table 1) and compared with available experimental observed values. It is apparent that PBE considerably overestimates all frequencies; while the mBJ and HSE06 results agree well with the experimental observations at least for α and γ bands. Noticing that these bands are indeed the electrontype pockets due to t 2g orbitals close to X, and are directly relevant to the band inversion features. Moreover, the frequencies for β bands obtained in mBJ and HSE06 calculations also agree with experimental values except for CeSb and PrSb. It is worth noting that the β bands are hole-type pockets due to p-orbitals close to Γ, and are very heavy in CeSb (renormalized ≈10 times for β 4 ). In short, the assumption of localized f-electron appear to be quite rational for these compounds. calculations. These results agree with our band structure analysis in previous section. In particular, we note that SmSb and DyBi are on the edge to the transition from topologically non-trivial to trivial in more accurate mBJ and HSE06 calculations. A crucial difference between topological trivial and non-trivial state lies at their boundary, where topologically protected surface state exists for the non-trivial state. Employing the surface Green's function method, we have obtained the band structure and surface state for the semi-infinite [001] surface of LnPn compounds. In Fig. 5a, b, we compare two typical cases CeBi and YbSb for non-trivial and trivial case, respectively. The topological surface states are prominent in CeBi (Fig. 5a), where 2 Diraccones due to these states are formed at X point. The two Dirac points are very close in energy, mainly because the valence difference between surface atoms and bulk atoms are ignored in our surface Green's function calculations. Nevertheless, these two points are not energetically degenerate, and they will further separate if the valence difference is taken into consideration (e.g. in slab model calculations 34 ). However, slab model calculation with mBJ is flawed due to technical reasons, and HSE06 is extremely expensive and beyond our current calculation capability. In contrast, only bulk states are visible in YbSb compound (Fig. 5b). We have also performed calculations of the CeBi [001] surface at different energy (Fig. 5c-e). At the zone center (Γ), only bulk state is visible at all energies. The surface states can be observed around the zone corner X. At E F , the surface states are mixed with bulk states, and cannot be easily distinguished. Below E F − 0.1 eV, the bulk states are void in the band inversion gap, and the surface states become prominent. At E F − 0.2 eV, the surface states dominates the zone corner.
Finally, we discuss the inverted gap and Z 2 variations with respect to the rare-earth atomic numbers. In LnPn, there are two effects due to replacing the rare-earth element. First, when X is replaced with heavier rare-earth element, the spin-orbit coupling percentage increase of Δ is proportional to 2|Z X − Z Ce |/Z Ce 38 . In fact, the t 2g splitting is 0.12 eV in CeSb, 0.15 eV in SmSb and 0.18 eV in YbSb, roughly following the estimation. Second, when Ln is replaced, it causes lanthanide contraction. However, it is not equivalent to chemical pressure effect. In order to demonstrate this, we show the band structure of DyBi along Γ-X at ambient pressure (6.31 Å), as well as with YbBi lattice constant (6.25 Å, or equivalent to 2.1 GPa) (Fig. 6). It is apparent that DyBi at ambient pressure is topologically trivial (also confirmed by Z 2 calculation), whereas it is non-trivial at 2.1 GPa (Z 2 is [1;000]). It is worth noting that similar pressure-induced topological transition has also been proposed for LaSb between 3 and~4 GPa 29 . The reduction of lattice constant increases the hopping between orbitals, and therefore band dispersion. In LnPn, the increase of Ln-t 2g band dispersion enhances the inversion gap (Fig. 6c) while the increase of Pn-p band dispersion reduces the gap (Fig. 6d). In fact, the largest nearest neighboring hopping between Bi-6p orbitals increases 77% (from 0.13 eV in CeBi to 0.23 eV in YbBi); while the largest nearest neighboring hopping between Ln-t 2g orbitals increases only 5% (from 0.43 eV in CeBi to 0.45 eV in YbBi). This in turn reflects the dispersion changes, which suggests that p-band dispersion is more sensitive to the lanthanide contraction than d-band. On the contrary, when a pure pressure is applied, the largest nearest neighboring hopping between Ln-t 2g orbitals increases much faster than that between Pn-p orbitals, while the spin-orbit coupling strength Δ remains relatively constant. For example, if we compress that lattice constant of DyBi by 5%, the largest nearest neighboring hopping between Dyt 2g orbitals increases 20% (from 0.44 to 0.53 eV); whereas the largest nearest neighboring hopping between Bi-6p orbitals only increases 10% (from 0.19 to 0.22 eV). In short, the inverted gap in LnPn systems increases with respect to increasing spin-orbit coupling strength Δ, increasing Ln-t 2g band dispersion, or decreasing Pn-p band dispersion.

Discussion
In conclusion, we have studied the electronic structure and band topology of LnPn compounds using state-of-art densityfunctionals. While PBE overestimates the band inversion gap, both mBJ and HSE06 functionals yield similar results. From CePn to YbPn, the atomic number increases but the band inversion gap size reduces. Therefore, a topological non-trivial to trivial transition is expected between PrSb and SmSb in LnSb compounds; while for LnBi compounds, the transition is expected around DyBi. Such variation in the inversion gap is related with lanthanide contraction, but is different from simple pressure effect. In particular, the d-orbitals respond differently in contraction than in pressure.

Methods
Calculation details. All the reported results were obtained with density functional theory (DFT) calculations as implemented in Vienna Abinitio Simulation Package (VASP) 39,40 . To ensure convergence, we employed plane-wave basis up to 480 eV, and 12 × 12 × 12 Γ-centered K-mesh so that the total energy converges to 1 meV per cell. The Perdew, Burke, and Ernzerhoff parameterization (PBE) of generalized gradient approximation (GGA) to the exchange correlation functional was employed 41 , with which the lattice constants were optimized until the internal stress less than 0.1 kbar. All the electronic structure calculations were obtained with optimized lattice constants unless specified otherwise. As PBE is known to overestimate the band inversions in solids due to the fact that it overestimates band widths, we have also performed calculations with both modified Becke-Johnson (mBJ) potentials as well as more expensive hybrid functional HSE06 as comparison.
Although the rare-earth elements in our calculations have an open f-shell, previous studies have shown that in antimonides and bismuthides the f-electrons are almost fully localized [15][16][17][18] , therefore the f-electrons are regarded as core state in all calculations. According to photoemission results, most lanthanide elements in LnPn compounds are 3+, therefore we have assumed Ln 3+ throughout the calculation.
Topological invariant and surface state. As LnPn compounds are centrosymmetric, the Z 2 indices can be calculated by evaluating the band parities at 8 timereversal invariant momenta (TRIM) 42 . The band structures obtained with PBE, mBJ, and HSE06 methods were fitted to a tight-binding (TB) Hamiltonian using the maximally localized wannier function (MLWF) method 43 with lanthanide t 2g orbitals and pnictogen-p orbitals. The resulting Hamiltonian were then used to calculate the Fermi surfaces as well as surface states using surface Green's function 44 .

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.