Tunable Electronic Structure and Topological Properties of $LnPn$ ($Ln$=Ce, Pr, Gd, Sm, Yb; $Pn$=Sb, Bi)

We have performed systematic first principles study of the electronic structure and band topology properties of $LnPn$ compounds ($Ln$=Ce, Pr, Gd, Sm, 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 $\mathcal{Z}_2$ for Ce$Pn$ and Yb$Pn$ are [1;000] and [0;000], respectively. In both hybrid functional and modified Becke-Johns calculations, a topologically nontrivial 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 effect.


INTRODUCTION
The 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 nonmagnetic 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 nontrivial 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 nontrivial 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 nontrivial 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 nontrivial 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 and 1b, 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% errorbar 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 towards 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 rock-salt 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 nonmagnetic band structure of CeBi to study some general features. Without spin-orbit 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 are 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 Γ-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. Since the electronic states are close to E F only around Γ and X, and the band topology is related with the anti-crossing between Γ-X, we shall focus on the band dispersion along this direction.
We have also performed calculations with itinerant f -electron in both nonmagnetic 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- 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 moves towards 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 are 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 [38] (TAB. I) 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 electron-type 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 worthy 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.  . 5a), where 2 Dirac-cones 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 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.
Firstly, when X is replaced with heavier rare-earth element, the spin-orbit coupling is increased. Since ∆ ∝ Z 2 X where Z X is the atomic number of X, the percentage increase of ∆ is proportional to 2|Z X − Z Ce |/Z Ce [39]. 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. Secondly, 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 nontrivial at 2.1 GPa (Z 2 is [1;000]). It is worthy noting that similar pressure-induced topological transition has also been proposed for LaSb between 3 ∼ 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 Dy-t 2g orbitals increases 20% (from 0.44 eV to 0.53 eV); whereas the largest nearest neighboring hopping between Bi-6p orbitals only in-7 creases 10% (from 0.19 eV 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 density-functionals. 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 nontrivial 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.

Calculation details
All the reported results were obtained with density functional theory (DFT) calculations as implemented in Vienna Abinitio Simulation Package (VASP) [40,41]. 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 exhange correlation functional was employed [42], 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. Since PBE is known to overesti- 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 lan-thanide elements in LnPn compounds are 3+, therefore we have assumed Ln 3+ throughout the calculation.

Topological invariant and surface state
Since LnPn compounds are centrosymmetric, the Z 2 indices can be calculated by evaluating the band parities at 8 time-reversal invariant momenta (TRIM) [43]. 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 [44] with lanthanide-t2g 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 [45].

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. Ref. [11], c Ref. [12], d Ref. [13], e Ref. [14], and f unpublished SdH results.