# Effects of self-consistent extended Hubbard interactions and spin-orbit
couplings on energy bands of semiconductors and topological insulators

Wooil Yang  Young-Woo Son  Email<EMAIL_ADDRESS>Korea Institute for Advanced
Study, Seoul 02455, Korea

###### Abstract

A first-principles computational method with self-consistent on-site and
inter-site Hubbard functionals is able to treat local and non-local Coulomb
interactions on an equal footing. To apply the method to understand solids
with strong spin-orbit coupling (SOC), we have extended a psuedohybrid
functional approach developed by Agapito-Curtarolo-Buongiorno Nardelli to
implement self-consistent extended Hubbard energy functionals for noncollinear
spin states. With this, energy bands of semiconductors with various SOC
strengths such as Si, Ge, GaAs, GaSb, CdSe and PdO are obtained, agreeing with
results from fully relativistic $GW$ approximation (FR-GWA) as well as
experiments. We also compute energy gaps of HgTe, CuTlS2, and CuTlSe2 and
assign them to be topological insulators correctly, unlike characteristic
failures for judging topological properties from typical hybrid functionals.
We demonstrate feasibility of our method to handle large systems by computing
surface bands of topological insulators, Bi${}_{\text{2}}$Se${}_{\text{3}}$
and Bi2Te3 with varying thickness up to eight quintuple layers. Considering
its low computational cost comparable to conventional ab intio methods and
improved accuracy to FR-GWA, we expect that our method provides an opportunity
to study large scale correlated systems with the strong SOC efficiently and
reliably.

##  I Introduction

First-principles calculations based on density functional theory (DFT) have
been a highly successful in simulating various properties of materials
efficiently [1, 2, 3, 4, 5]. Concomitantly, there has been a significant
progress in using DFT to understand electronic properties of systems having
heavy elements, in which the relativistic effects such as strong spin-orbit
coupling (SOC) should be taken into account. However, DFT calculations based
on approximate exchange-correlation (XC) functionals such as local density
approximation (LDA) [2] or generalized gradient approximation (GGA) [6] have
self-interaction error (SIE) related with heavy elements, posing a challenge
to capture their electronic properties [7, 8].

Among various methods to overcome such issues, the DFT+$U$ including the local
Coulomb repulsion of $U$ is widely employed [9, 10, 11, 12]. This is based on
the Hubbard model since SIEs occur primarily in localized $d$ and $f$
orbitals. Particularly, recent studies demonstrate accurate description of
ground electronic states for several materials with strong local Coulomb
repulsions and SOCs such as transition metal oxides [13], pyrochlore iridates
[14], americium monochalcogenides [15] and osmium-based double perovskites
[16]. We also note that recent developments of fully relativistic $GW$
approximations (FR-GWA) [17, 18, 19, 20, 21, 22, 23, 24, 25, 26] obtained
spin-orbit split energy bands of semiconductors in good agreements with
experiments.

Despite these successes of the DFT+$U$ method, it turns out to be insufficient
especially in case that orbital hybridization or inter-site Hubbard
interaction of $V$ becomes crucial [27, 28, 29, 30]. Then, the DFT+$U$+$V$
method based on the extended Hubbard model including both $U$ and $V$ has been
considered for several materials with negligible SOCs [27, 28, 29, 30],
demonstrating its improved capabilities in computing electronic and phonon
properties of solids [27, 28, 29, 30, 31, 32, 33, 34, 35, 36]. We also note
that a previous work [37] using a localized orbital basis set well
demonstrated an importance of including both inter-site Hubbard interactions
and SOCs for the nodal-line semimetal ZrSiSe. Therefore, we expect that the
DFT with both extended Hubbard interactions and SOC could be another efficient
computational tool to calculate reliable physical properties of solids having
heavy atomic elements ab initio.

In this paper, we investigate effects of self-consistent extended Hubbard
interactions of $U$ and $V$ on energy gaps and spin-orbit splittings of
semiconductors and topological insulators (TIs) with various SOC strengths.
For semiconductors, we compute spin-orbit split energy bands of Si, Ge, GaAs,
GaSb, CdSe and PdO. To compute $U$, $V$ and SOC on an equal footing, we extend
our previous work [28, 31] on self-consistent $U$ and $V$ for collinear states
to include noncollinear states. Particularly, for self-consistent $V$, we
apply the pseudohybrid density functionals by Agapito-Curtarolo-Buongiorno
Nardelli (ACBN0) [38] to compute the inter-site Hubbard parameters [28, 29]
for noncollinear spin states with a plane-wave (PW) basis set [39, 40]. For
TIs, we compute inversion energy gaps of HgTe, CuTlS2 and CuTlSe2 and confirm
them to be topologically nontrivial. This is in contrast to some erroneous
assignments from hybrid functionals [41]. Finally, we compute evolution of
topological surface states of Bi${}_{\text{2}}$Se${}_{\text{3}}$ and
Bi${}_{\text{2}}$Te${}_{\text{3}}$ as increasing their slab thickness up to
eight quintuple layers (QLs), demonstrating our capability to handle large
systems. We show that all results obtained by our new method agree with FR-GWA
as well as experiments, thus providing an efficient and reliable ab initio
method for large scale systems having heavy atoms.

The paper is organized as follows: in Sec. II, we present formalism for
DFT+$U$+$V$ method with self-consistent $U$ and $V$ for noncollinear states.
The detailed computational parameters for the calculations are presented in
Sec. III. In Sec. IV, we show computational results for atomic and electronic
structures of representative semiconductors and TIs with various SOCs and for
topological surface states of Bi2Se3 and Bi2Te3 with varying their thickness.
We conclude in Sec. V.

##  II Formalism

We briefly review the extended Hubbard functional formalism [27, 28, 29] with
atomic SOCs built on top of noncollinear DFT using the norm-conserving
pseudopotential (PP) and PW basis set. To obtain the generalized occupation
matrix for Hubbard interactions by projecting PW wavefunctions on atomic
states, the projector characterized by the total angular momentum of $j$ would
be required. However, we practically adopt the $j$-averaged radial components
of atomic wavefunctions from PPs [42, 43, 44, 45] such as,

$V({\bf r},{\bf
r^{\prime}})=V_{\text{loc}}+\sum_{l}\left[V^{\text{SR}}_{l}({\bf r},{\bf
r^{\prime}})+{\bf L}\cdot{\bf S}V^{\text{SO}}_{l}({\bf r},{\bf
r^{\prime}})\right],$ (1)

where ${\bf L}$ and ${\bf S}$ are angular and spin momentum, respectively and
$l$ an azimuthal quantum number. Here, the non-local potential can be split
into scalar-relativisitc (SR) and spin-orbit (SO) components as sums and
differences of relativistic potentials of $V^{\text{Rel}}_{j}({\bf r},{\bf
r^{\prime}})$ with $j=l+\frac{1}{2}\sigma$ ($j>0$ and spin index $\sigma=\pm
1$) [42, 43, 44, 45] such that

$\displaystyle V^{\text{SR}}_{l}$ $\displaystyle=$
$\displaystyle\frac{(l+1)V^{\text{Rel}}_{l+1/2}+lV^{\text{Rel}}_{l-1/2}}{2l+1},$
$\displaystyle V^{\text{SO}}_{l}$ $\displaystyle=$ $\displaystyle
2\frac{V^{\text{Rel}}_{l+1/2}-V^{\text{Rel}}_{l-1/2}}{2l+1}.$ (2)

Although the $j$-averaged atomic projector from potentials in Eqs. (1) and (2)
is not spin-dependent, we can practically assign spinor representation for the
projector following the previous studies [29, 46]. Then, the mean field
Hubbard energy (${\mathcal{E}}_{MF}$) including on-site ($U_{I}$) and inter-
site ($V_{IJ}$) electron-electron interactions for $I$-th and $J$-th atoms are
split into collinear (${\mathcal{E}}^{\text{c}}_{U}$ and
${\mathcal{E}}^{\text{c}}_{V}$) and noncollinear
(${\mathcal{E}}^{\text{nc}}_{U}$ and ${\mathcal{E}}^{\text{nc}}_{V}$)
contributions such that,

$\displaystyle{\mathcal{E}}_{\text{MF}}$ $\displaystyle=$
$\displaystyle{\mathcal{E}}^{\text{c}}_{U}+{\mathcal{E}}^{\text{nc}}_{U}+{\mathcal{E}}^{\text{c}}_{V}+{\mathcal{E}}^{\text{nc}}_{V}$
(3) $\displaystyle=$
$\displaystyle\sum_{I}\left[E^{\text{c}}_{U_{I}}+E^{\text{nc}}_{U_{I}}\right]+\sum_{\\{I,J\\}}\left[E^{\text{c}}_{V_{IJ}}+E^{\text{nc}}_{V_{IJ}}\right],$

where $\\{I,J\\}$ denotes a pair of different atoms $I$ and $J$ within a given
cutoff distance. Each term in Eq. (3) can be expressed as

$\displaystyle E^{\text{c}}_{U_{I}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{ij,\sigma\sigma^{\prime}}(Ii\sigma,Ii\sigma|Ij\sigma^{\prime},Ij\sigma^{\prime})n^{II\sigma\sigma}_{ii}n^{II\sigma^{\prime}\sigma^{\prime}}_{jj}$
(4)
$\displaystyle-\frac{1}{2}\sum_{ij,\sigma}(Ii\sigma,Ij\sigma|Ij\sigma,Ii\sigma)n^{II\sigma\sigma}_{ii}n^{II\sigma\sigma}_{jj},$
$\displaystyle E^{\text{nc}}_{U_{I}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{ij,\sigma}(Ii\sigma,Ij{\sigma}|Ij\bar{\sigma},Ii\bar{\sigma})n^{II\sigma\bar{\sigma}}_{ii}n^{II\bar{\sigma}\sigma}_{jj},$
(5) $\displaystyle E^{\text{c}}_{V_{IJ}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{ij,\sigma\sigma^{\prime}}(Ii\sigma,Ii\sigma|Jj\sigma^{\prime},Jj\sigma^{\prime})$
(6)
$\displaystyle\times\left[n^{II\sigma\sigma}_{ii}n^{JJ\sigma^{\prime}\sigma^{\prime}}_{jj}-\delta_{\sigma,\sigma^{\prime}}n^{IJ\sigma\sigma^{\prime}}_{ij}n^{JI\sigma^{\prime}\sigma}_{ji}\right],$
$\displaystyle E^{\text{nc}}_{V_{IJ}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{ij,\sigma}(Ii\sigma,Ii\sigma|Jj\bar{\sigma}Jj\bar{\sigma})n^{IJ\sigma\bar{\sigma}}_{ij}n^{JI\bar{\sigma}\sigma}_{ji}.$
(7)

Here, the generalized occupation number $n^{IJ\sigma\sigma^{\prime}}_{ij}$ can
be written as

$\displaystyle n^{IJ\sigma\sigma^{\prime}}_{ij}$ $\displaystyle\equiv$
$\displaystyle n^{I,n,l,\sigma,J,n^{\prime},l^{\prime},\sigma^{\prime}}_{ij}$
(8) $\displaystyle=$
$\displaystyle\sum_{m\mathbf{k}}w_{\mathbf{k}}f_{m\mathbf{k}}\langle\psi^{\sigma}_{m\mathbf{k}}|\phi^{I,n,l}_{i\sigma}\rangle\langle\phi^{J,n^{\prime},l^{\prime}}_{j\sigma^{\prime}}|\psi^{\sigma^{\prime}}_{m\mathbf{k}}\rangle,$

where $f_{m\mathbf{k}}$ is the Fermi-Dirac function of the PW Bloch state
$|\psi^{\sigma}_{m\mathbf{k}}\rangle$ of the m-th band at a momentum
$\mathbf{k}$ and $\textit{w}_{m\mathbf{k}}$ the $\mathbf{k}$-grid weight.
$|\phi^{I,n,l}_{i\sigma}\rangle$ is the Löwdin orthonormalized atomic
wavefunction [47] as a projector where $n$ is the principal quantum numbers.
$\bar{\sigma}$ denotes an opposite spin orientation to $\sigma$. We will use
an abbreviated notation of $I$, which specific principal and azimuthal qunatum
numbers of the $I$-th atom hereafter, $e.g.$,
$|\phi^{I}_{i\sigma}\rangle\equiv|\phi^{I,n,l}_{i\sigma}\rangle$. We note
that, when $I=J$, Eq. (8) becomes the on-site occupation, i.e.,
$n^{I}_{i\sigma}\equiv n^{II\sigma\sigma}_{ii}$. Here, $i,\cdots,l$ represent
the magnetic quantum number (note that $l$ is not azimuthal quantum number
hereafter). The Coulomb integrals are written using the notation,

$\displaystyle(Ii\sigma,Ij\sigma|Jk\sigma^{\prime},Jl\sigma^{\prime})$ (9)
$\displaystyle\equiv$ $\displaystyle\langle
Ii\sigma,Jk\sigma^{\prime}|V_{ee}|Ij\sigma,Jl\sigma^{\prime}\rangle$
$\displaystyle=$ $\displaystyle\int d{\bf r}\int d{\bf
r^{\prime}}\frac{\rho_{ij\sigma}({\bf r}-{\bf
R}_{I})\rho_{kl\sigma^{\prime}}({\bf r^{\prime}}-{\bf R}_{J})}{|{\bf r}-{\bf
r^{\prime}}-{\bf R}_{I}+{\bf R}_{J}|},$

where $\phi_{i\sigma}({\bf r}-{\bf R}_{I})=\langle{\bf r}-{\bf
R}_{I}|\phi_{i\sigma}^{I}\rangle$, $|\phi^{I}_{i\sigma}\rangle$ denotes an
associated atomic orbital state of the $I$-th atom, ${\bf R}_{I}$ its
position, $\rho_{ij\sigma}({\bf r}-{\bf R}_{I})=\phi^{*}_{i\sigma}({\bf
r}-{\bf R}_{I})\phi_{j\sigma}({\bf r}-{\bf R}_{I})$ and
$\rho_{kl\sigma^{\prime}}({\bf r^{\prime}}-{\bf
R}_{J})=\phi^{*}_{k\sigma^{\prime}}({\bf r^{\prime}}-{\bf
R}_{J})\phi_{l\sigma^{\prime}}({\bf r^{\prime}}-{\bf R}_{J})$.

Following Dudarev’s approach [11], we define the averaged on-site and inter-
site interactions [27] such as

$\displaystyle U_{I}$ $\displaystyle=$
$\displaystyle\frac{\sum_{ij}(Ii\sigma,Ii\sigma|Ij\sigma^{\prime},Ij\sigma^{\prime})}{(2l_{I}+1)^{2}},$
(10) $\displaystyle J_{I}$ $\displaystyle=$ $\displaystyle\frac{\sum_{i\neq
j}(Ii\sigma,Ij\sigma|Ij\sigma,Ii\sigma)}{2l_{I}(2l_{I}+1)},$ (11)
$\displaystyle V_{IJ}$ $\displaystyle=$
$\displaystyle\frac{\sum_{ij}(Ii\sigma,Ii\sigma|Jj\sigma^{\prime},Jj\sigma^{\prime})}{(2l_{I}+1)(2l_{J}+1)},$
(12)

where $l_{I}(l_{J})$ is an azimuthal quantum number of interested orbitals in
the $I(J)$-th atom. Using Eqs. (10)$\sim$(12), Eqs (4)$\sim$(7) can be
rewritten as

$\displaystyle E^{\text{c}}_{U_{I}}$ $\displaystyle=$
$\displaystyle\frac{U_{I}-J_{I}}{2}\sum_{i\neq
j,\sigma}n^{I}_{i\sigma}n^{I}_{j\sigma}+\frac{U_{I}}{2}\sum_{ij,\sigma}n^{I}_{i\sigma}n^{I}_{j\bar{\sigma}},$
(13) $\displaystyle E^{\text{nc}}_{U_{I}}$ $\displaystyle=$
$\displaystyle-\frac{U_{I}}{2}\sum_{i,\sigma}n^{I\sigma\bar{\sigma}}_{ii}n^{I\bar{\sigma}\sigma}_{ii}-\frac{J_{I}}{2}\sum_{i\neq
j,\sigma}n^{I\sigma\bar{\sigma}}_{ii}n^{I\bar{\sigma}\sigma}_{jj},$ (14)
$\displaystyle E^{\text{c}}_{V_{IJ}}$ $\displaystyle=$
$\displaystyle\frac{V_{IJ}}{2}\sum_{ij,\sigma\sigma^{\prime}}\left[n^{I}_{i\sigma}n^{J}_{j\sigma^{\prime}}-\delta_{\sigma\sigma^{\prime}}n^{IJ\sigma\sigma^{\prime}}_{ij}n^{JI\sigma^{\prime}\sigma}_{ji}\right],$
(15) $\displaystyle E^{\text{nc}}_{V_{IJ}}$ $\displaystyle=$
$\displaystyle-\frac{V_{IJ}}{2}\sum_{ij,\sigma}n^{IJ\sigma\bar{\sigma}}_{ij}n^{JI\bar{\sigma}\sigma}_{ji}.$
(16)

Following Mosey et. al. [48, 49] and ACBN0 [38, 50, 28, 29] approaches, we now
rewrite the mean field energies in Eqs. (4)$\sim$(7) using the noncollinear
renormalized density matrices of $P^{IJ\sigma\sigma^{\prime}}_{ij}$ that can
be written as

$P^{IJ\sigma\sigma^{\prime}}_{ij}=\sum_{m\mathbf{k}}w_{\mathbf{k}}f_{m\mathbf{k}}N^{IJ}_{{m\mathbf{k}}}\langle\psi^{\sigma}_{m\mathbf{k}}|\phi^{I}_{i\sigma}\rangle\langle\phi^{J}_{j\sigma^{\prime}}|\psi^{\sigma^{\prime}}_{m\mathbf{k}}\rangle.$
(17)

Here renormalized occupation number of $N^{IJ}_{{m\mathbf{k}}}$ for different
atoms ($I\neq J$) is defined as

$N^{IJ}_{{m\mathbf{k}}}\equiv\sum_{\\{I\\}}{\textrm{tr}}_{o}({\textrm{tr}}_{s}(N^{I\sigma\sigma^{\prime}}_{i,m{\bf
k}}))+\sum_{\\{J\\}}{\textrm{tr}}_{o}({\textrm{tr}}_{s}(N^{J\sigma\sigma^{\prime}}_{i,m{\bf
k}})),$ (18)

where $N^{I\sigma\sigma^{\prime}}_{i,m{\bf
k}}\equiv\langle\psi^{\sigma}_{m\mathbf{k}}|\phi^{I}_{i\sigma}\rangle\langle\phi^{I}_{i\sigma^{\prime}}|\psi^{\sigma^{\prime}}_{m\mathbf{k}}\rangle$,
${\textrm{tr}}_{o(s)}$ is a trace over orbitals (spins), and
$\sum_{\\{I(J)\\}}$ denote summations of all $I(J)$ atoms. For a pair of the
same atoms or for the on-site interactions, Eq. (18) reduces to
$N^{II}_{{m\mathbf{k}}}=\sum_{\\{I\\}}{\textrm{tr}}_{o}({\textrm{tr}}_{s}(N^{I\sigma\sigma^{\prime}}_{i,m{\bf
k}}))$. In Eq. (18), we consider traces over spins since
$N^{IJ}_{{m\mathbf{k}}}$ plays as a role of screening to charges that should
not depend on spins [38, 29, 28]. If including off-diagonal spin orientations
in Eq. (18), we confirm erroneous changes of interactions with varying spin
orientations. With introduction of Eq. (17), we rewrite Eqs. (4)$\sim$(7) in
most general forms as followings,

$\displaystyle
E_{U_{I}}=\frac{1}{2}\sum_{ijkl,\sigma\sigma^{\prime}}{P}^{I\sigma}_{ij}{P}^{I\sigma^{\prime}}_{kl}(Ii\sigma,Ij\sigma|Ik\sigma^{\prime},Il\sigma^{\prime})$
(19) $\displaystyle
E_{J_{I}}=-\frac{1}{2}\sum_{ijkl,\sigma\sigma^{\prime}}{P}^{I\sigma\sigma^{\prime}}_{ij}{P}^{I\sigma^{\prime}\sigma}_{kl}(Ii\sigma,Il\sigma|Ik\sigma^{\prime},Ij\sigma^{\prime})$
(20) $\displaystyle
E_{V_{IJ}}=\frac{1}{4}\sum_{ijkl,\sigma\sigma^{\prime}}\left[{P}^{I\sigma}_{ij}{P}^{J\sigma^{\prime}}_{kl}-{P}^{IJ\sigma\sigma^{\prime}}_{il}{P}^{JI\sigma^{\prime}\sigma}_{kj}\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times(Ii\sigma,Ij\sigma|Jk\sigma^{\prime},Jl\sigma^{\prime}),$
(21)

where $P^{I\sigma}_{ij}\equiv P^{II\sigma\sigma}_{ij}$ and
$P^{I\sigma\sigma^{\prime}}_{ij}\equiv P^{II\sigma\sigma^{\prime}}_{ij}$.
Then, by equating corresponding terms between Eqs. (13)$\sim$(16) and Eqs.
(19)$\sim$(21), respectively, we can obtain the pseudohybrid-type extended
Hubbard functionals for the noncollinear states [38, 50, 29] within
DFT+$U$+$V$ formalism [27, 28],

$\displaystyle{\mathcal{U}}_{I}$ $\displaystyle=$
$\displaystyle\sum_{ijkl,\sigma\sigma^{\prime}}\frac{{P}^{I\sigma}_{ij}{P}^{I\sigma^{\prime}}_{kl}}{{\mathcal{N}}_{U_{I}}}(Ii\sigma,Ij\sigma|Ik\sigma^{\prime},Il\sigma^{\prime}),$
(22) $\displaystyle{\mathcal{J}}_{I}$ $\displaystyle=$
$\displaystyle\sum_{ijkl,\sigma\sigma^{\prime}}\frac{{P}^{I\sigma\sigma^{\prime}}_{ij}{P}^{I\sigma^{\prime}\sigma}_{kl}}{{\mathcal{N}}_{J_{I}}}(Ii\sigma,Il\sigma|Ik\sigma^{\prime},Ij\sigma^{\prime}),$
(23) $\displaystyle{\mathcal{V}}_{IJ}$ $\displaystyle=$
$\displaystyle\sum_{ijkl,\sigma\sigma^{\prime}}\frac{{P}^{I\sigma}_{ij}{P}^{J\sigma^{\prime}}_{kl}-{P}^{IJ\sigma\sigma^{\prime}}_{il}{P}^{JI\sigma^{\prime}\sigma}_{kj}}{{\mathcal{N}}_{V_{IJ}}}$
(24)
$\displaystyle\times(Ii\sigma,Ij\sigma|Jk\sigma^{\prime},Jl\sigma^{\prime}),$

where

$\displaystyle{\mathcal{N}}_{U_{I}}$ $\displaystyle=$
$\displaystyle\sum_{i\neq
j,\sigma}n^{I}_{i\sigma}n^{I}_{j\sigma}+\sum_{ij,\sigma}n^{I}_{i\sigma}n^{I}_{j\bar{\sigma}}-\sum_{i,\sigma}n^{I\bar{\sigma}\sigma}_{ii}n^{I\bar{\sigma}\sigma}_{ii},$
$\displaystyle{\mathcal{N}}_{J_{I}}$ $\displaystyle=$
$\displaystyle\sum_{i\neq
j,\sigma}\left[n^{I}_{i\sigma}n^{I}_{j\sigma}+n^{I\sigma\bar{\sigma}}_{ii}n^{I\bar{\sigma}\sigma}_{jj}\right],$
$\displaystyle{\mathcal{N}}_{V_{IJ}}$ $\displaystyle=$
$\displaystyle\sum_{ij,\sigma\sigma^{\prime}}\left[{n^{I\sigma}_{i}n^{J\sigma^{\prime}}_{j}}-{n^{IJ,\sigma\sigma^{\prime}}_{ij}n^{JI,\sigma^{\prime}\sigma}_{ji}}\right].$

Considering the double counting corrections within the fully localized limit
as discussed before [51, 10, 52, 53, 54, 27, 55], Eqs. (13)$\sim$(16) have
following double counting energies (${\mathcal{E}}_{\text{dc}}$) for on-site
($E^{\text{dc}}_{U_{I}}$) and inter-site ($E^{\text{dc}}_{V_{IJ}}$) parts,

$\displaystyle E^{\text{dc}}_{U_{I}}$ $\displaystyle=$
$\displaystyle{\frac{U_{I}}{2}}(N^{2}_{I}-N_{I})-\frac{J_{I}}{2}\left(\frac{N^{2}_{I}+|{\bf
m}_{I}|^{2}}{2}-N_{I}\right)$ (25) $\displaystyle E^{\text{dc}}_{V_{IJ}}$
$\displaystyle=$ $\displaystyle{\frac{V_{IJ}}{2}}N_{I}N_{J}$ (26)

where
$N_{I}={\textrm{tr}}_{o}({\textrm{tr}_{s}}(n^{I\sigma\sigma^{\prime}}_{ij}))$
and ${\bf
m}_{I}={\textrm{tr}}_{o}({\textrm{tr}_{s}}({\bm{\sigma}}n^{I\sigma\sigma^{\prime}}_{ij}))$.

After subtracting Eqs. (13)$\sim$(16) by Eqs. (25) and (26) and then replacing
on-site and inter-site interactions in Eqs. (10), (11) and (12) with Eqs (22),
(23) and (24) respectively, the resulting energy functionals with SOC and
without double countings can be written as

$\displaystyle{E}_{U_{I}}$ $\displaystyle=$
$\displaystyle\frac{{\mathcal{U}}_{I}-{\mathcal{J}}_{I}}{2}\left[\sum_{i,\sigma}n^{I}_{i\sigma}-\sum_{ij,\sigma\sigma^{\prime}}{n^{II\sigma\sigma^{\prime}}_{ij}n^{II\sigma^{\prime}\sigma}_{ji}}\right],$
(27) $\displaystyle E_{V_{IJ}}$ $\displaystyle=$
$\displaystyle-{\frac{{\mathcal{V}}_{IJ}}{2}}\sum_{ij,\sigma\sigma^{\prime}}{n^{IJ\sigma\sigma^{\prime}}_{ij}n^{JI\sigma^{\prime}\sigma}_{ji}}.$
(28)

Compared with earlier results without SOC [38, 50, 28, 29], the forms of the
functionals are quite similar and have additional summations over all spin
orientations owing to noncollinear states.

##  III Computational details

Our ab initio calculations have been performed with Quantum Espresso (QE)
package [39, 40]. We employed LDA parametrized by Perdew and Zunger (PZ) [56],
and GGA by Perdew-Burke-Ernzehof (PBE) [6] and its revision for solids
(PBEsol) [57]. We used the optimized norm-conserving Vanderbilt PPs [44] with
parameters from the Pseudo Dojo database [58]. For scalar relativistic (SR)
calculations, we include relativstic mass correction and Darwin term together
with $V^{\text{SR}}_{l}$ in Eq. (1). For fully relativistic (FR) ones, we
further add $V^{\text{SO}}_{l}$ in Eq. (1) on top of the aformentioned ones.
Unless otherwise mentioned, we will use FR(SR)-LDA and FR(SR)-GGA for FR (SR)
calculations with FR(SR)-PZ and FR(SR)-PBEsol, respectively. Then, unless
otherwise mentioned, we call our method as FRL+$U$+$V$ if using FR-PZ
approximation for DFT+$U$+$V$ calculation and FRG+$U$+$V$ if using FR-PBEsol.

We employed our modified in-house QE package to obtain self-consistent
extended Hubbard parameters with PW basis sets [28] and SOCs. We use Löwdin
orthogonalized atomic orbitals for projectors as were done in our previous
studies [28, 31, 32]. In calculating $U$ and $V$ in Eqs. (22)$\sim$(24), we
use PPs without the semicore in valences because the occupation numbers in
valences using the $j$-averaged atomic projectors slightly vary whether the
semicore states are included or not. Although the difference due to semicore
is not significant, to reflect charge densities from semicore states reliably,
the Hubbard parameters are first obtained with SR approximations, if we use
LDA PPs having semicore states, and then proceed to FR ones with Eqs. (27) and
(28).

The on-site U for all s orbitals was set to be zero. The inter-site
interactions of $V$ were included up to the second nearest neighbors for Si,
Ge, GaAs, GaSb, Bi2Se3, and Bi2Te3 while for CdSe, PdO, HgTe, CuTlS2, and
CuTlSe2, ${\mathcal{V}}_{IJ}$ is considered for the adjacent atoms only, all
of which are confirmed to be converged [28]. Then, we determined the
equilibrium lattice constants and atomic positions with converged on-site and
inter-site Hubbard interactions [31, 59]. We use $15\times 15\times 15$
$k$-point grid of the Brillouin zone (BZ) sampling for Si, Ge, GaAs, GaSb and
HgTe, $15\times 15\times 10$ for CdSe and PdO, $11\times 11\times 11$ mesh for
CuTlS2 and CuTlSe2, and $11\times 11\times 1$ mesh for slab models of Bi2Se3
and Bi2Te3, respectively. We set the kinetic energy cutoff to 120 Ry for Si,
Ge, HgTe, CuTlS2, CuTlSe2, Bi2Se3, and Bi2Te3, 350 Ry for GaAs and GaSb, 200
Ry for CdSe and 100 Ry for PdO, respectively.

##  IV Results

###  IV.1 Si, Ge, GaAs, and GaSb

Table 1: Calculated $U$ and $V$ (in eV) with FRG+$U$+$V$ for a set of semiconductors. $U_{p}$, the on-site Hubbard interaction for $p$ orbital and, $V_{ss}$, $V_{sp}$, and $V_{pp}$, the inter-site Hubbard parameters between $s$-$s$, $s$-$p$, and $p$-$p$ orbitals of the adjacent atoms. For $U_{p}$ for GaX (X$=$As, Sb), the left and right values are for Ga and X $p$-orbitals, repsectively while, for $V_{sp}$, they are between Ga $s(p)$-orbital and X $p(s)$-orbital.  | $U_{p}$ | $V_{ss}$ | $V_{sp}$ | $V_{pp}$
---|---|---|---|---
Si | 3.49 | 0.88 | 0.72 | 1.85
Ge | 3.32 | 1.04 | 0.69 | 1.75
GaAs | 0.35, 1.96 | 0.84 | 1.30, 0.81 | 1.67
GaSb | 0.43, 1.43 | 0.9 | 1.07, 0.79 | 1.43
Figure 1: FR-LDA and FRL+$U$+$V$ band structures of Si, Ge, GaAs, and GaSb,
plotted along high-symmetry lines of the Brillouin zone. The valence band
maximum is set to zero.

We first consider Si, Ge, and GaAs that are prototypical semiconductors to
test our method and to compare calculation results with those from other
methods as well as many experiments. We also compute atomic and electronic
structures of GaSb of which spin-orbit splitting in its valence bands is
comparable to its band gap.

Table 1 shows the calculated self-consistent $U$ and $V$ values. We have
checked that for all of them, which are either diamond or zincblende
semiconductors, the calculated Hubbard parameters from SR- and FR-calculations
do not show much differences ($<0.01$ eV) agreeing with values from our
previous studies without SOCs [31]. With FRG+$U$+$V$, we determined the
equilibrium lattice constants (in Å) of 5.44, 5.67, 5.64 and 6.08 for Si, Ge,
GaAs and GaSb, respectively, in excellent agreements with experiments [60].

Table 2: The energy gap, interband gap, and spin-orbit splitting in eV for Si, computed using FR-LDA, FRL+U+V, FR-GGA, and FRG+U+V, compared to $GW$ calculations and experiments. | $E_{g}$ | $E^{\Gamma}_{6c-8v}$ | $\Delta^{\Gamma}_{v}$ | $\Delta^{\Gamma}_{c}$ | $\Delta^{\text{L}}_{v}$ | $\Delta^{\text{L}}_{c}$
---|---|---|---|---|---|---
$GW$+SOC111Reference [61] | 1.27 | 3.28 | 0.05 | 0.04 | 0.03 | 0.02
FR-$GW$222Reference [26] | 1.22 | 3.22 | 0.05 | 0.04 | 0.03 | 0.01
FR-LDA | 0.45 | 2.46 | 0.05 | 0.03 | 0.03 | 0.01
FRL+$U$+$V$ | 1.18 | 3.38 | 0.06 | 0.04 | 0.04 | 0.02
FR-GGA | 0.44 | 2.47 | 0.05 | 0.04 | 0.03 | 0.01
FRG+$U$+$V$ | 1.20 | 3.40 | 0.06 | 0.04 | 0.04 | 0.02
Exp | 1.22333Reference [62] | 3.34444Reference [63] | 0.044555Reference [60] | 0.03$\sim$0.045 | 0.035 | -
Table 3: The energy gap, interband gap, and spin-orbit splitting in eV for Ge, computed using FR-LDA, FRL+U+V, FR-GGA, and FRG+U+V, compared to $GW$ calculations and experiments. | $E_{g}$ | $E^{\Gamma}_{7c-8v}$ | $\Delta^{\Gamma}_{v}$ | $\Delta^{\Gamma}_{c}$ | $\Delta^{\text{L}}_{v}$ | $\Delta^{\text{L}}_{c}$
---|---|---|---|---|---|---
$GW$+SOC666Reference [61] | 0.54 | 0.38 | 0.32 | 0.24 | 0.20 | 0.12
FR-$GW$777Reference [26] | 0.74 | 0.96 | 0.30 | 0.21 | 0.19 | 0.08
FR-LDA | 0.13 | 0.14 | 0.31 | 0.22 | 0.19 | 0.10
FRL+$U$+$V$ | 0.53 | 0.68 | 0.40 | 0.29 | 0.25 | 0.11
FR-GGA | 0.00 |  $-$0.19 | 0.30 | 0.21 | 0.19 | 0.09
FRG+$U$+$V$ | 0.33 | 0.33 | 0.39 | 0.28 | 0.24 | 0.11
Exp888Reference [60] | 0.79 | 0.90 | 0.297 | 0.200 | 0.228 | -
Table 4: The energy gap and spin-orbit splitting in eV for GaAs, computed using FR-LDA, FRL+U+V, FR-GGA, and FRG+U+V, compared to $GW$ calculations and experiments. | $E_{g}$ | $\Delta^{\Gamma}_{v}$ | $\Delta^{\Gamma}_{c}$ | $\Delta^{\text{L}}_{v}$ | $\Delta^{\text{L}}_{c}$
---|---|---|---|---|---
$GW$+SOC999Reference [61] | 1.31 | 0.35 | 0.20 | 0.22 | 0.09
FR-$GW$101010Reference [26] | 1.49 | 0.34 | 0.17 | 0.21 | 0.07
FR-LDA | 0.54 | 0.32 | 0.19 | 0.20 | 0.08
FRL+$U$+$V$ | 1.34 | 0.37 | 0.19 | 0.22 | 0.08
FR-GGA | 0.31 | 0.34 | 0.19 | 0.21 | 0.08
FRG+$U$+$V$ | 1.03 | 0.39 | 0.19 | 0.23 | 0.08
Exp | 1.57111111Reference [64] | 0.34121212Reference [60] | 0.17112 | 0.2212 | 0.0512
Table 5: The interband gap, and spin-orbit splitting in eV for GaSb, computed using FR-LDA, FRL+U+V, FR-GGA, and FRG+U+V, compared to $GW$ calculations and experiments. | $E^{\Gamma}_{6c-8v}$ | $E^{\text{L}}_{6c}-E^{\Gamma}_{8v}$ | $\Delta^{\Gamma}_{v}$ | $\Delta^{\Gamma}_{c}$ | $\Delta^{\text{L}}_{v}$ | $\Delta^{\text{L}}_{c}$
---|---|---|---|---|---|---
$GW$+SOC131313Reference [61] | 0.70 | 0.85 | 0.73 | 0.21 | 0.42 | 0.12
FR-$GW$141414Reference [26] | 0.82 | 0.78 | 0.73 | 0.20 | 0.42 | 0.09
FR-LDA | 0.12 | 0.24 | 0.74 | 0.22 | 0.42 | 0.12
FRL+$U$+$V$ | 0.78 | 0.80 | 0.82 | 0.25 | 0.46 | 0.12
FR-GGA |  $-$0.23 | 0.13 | 0.72 | 0.22 | 0.42 | 0.12
FRG+$U$+$V$ | 0.43 | 0.68 | 0.81 | 0.24 | 0.45 | 0.12
Exp151515Reference [60] | 0.822 | 0.907 | 0.756 | 0.213 | 0.430 | 0.13

In Fig. 1, we display the band structures of Si, Ge, GaAs and GaSb calculated
with FR-LDA and FRL+U+V. It is immediately noticeable that all conduction
bands shift up significantly, increasing their band gaps. As discussed before
[27, 28], the inter-site interactions dramatically improve the band gap values
in these covalent materials, which are comparable to results from other
sophisticated methods. We also note that the energetic positions of the lowest
valence band obtained using FR-LDA and FRL+U+V are almost same to each other.
Furthermore, for these semiconductors, the bandwidths of the upper valence
bands increase with including $V$. This is slightly different from $GW$
calculations that increase the bandwidth of the upper valence bands only in Ge
and GaSb [26].

In Tables 2, 3, 4 and 5, we compare values of the interband gaps, spin-orbit
splittings and energy gaps of these semiconductors with $GW$ approximations
and experiments in details. Following nomenclatures in Ref. [26], we use
FR-$GW$ for full relativistic spinor $GW$ calculations while $GW$+SOC for
treating SOCs as perturbations to $GW$ approximation without SOC [61, 26].
Specifically, we compare energy gap ($E_{g}$), spin-orbit splittings of
conduction (valence) bands at high symmetric $\Gamma(\text{L})$-point
($\Delta^{\Gamma(\text{L})}_{c(v)}$), and interband gaps between conduction
and valence bands labeled according to their irreducible representations [26].

For Si in Table 2, the computed values using extended Hubbard interactions
with SOCs show good agreements with those from FR-$GW$. Considering negligible
contributions from $U$ [28], we confirm again that the inter-site $V$ plays a
pivotal role in computing band gaps and interband gap of $E^{\Gamma}_{6c-8v}$
regardless of XC functionals. Unlike the gaps, spin orbit splittings at
$\Gamma$ and L-points are more or less similar with each other, implying
marginal contributions of SOCs as expected.

For Ge, GaAs and GaSb in Tables 3, 4 and 5, the computed band gaps as well as
interband gaps from FRL+$U$+$V$ show excellent agreements with values from
FR-$GW$ as well as experiments. We also note that spin orbit splittings at
high symmetric points are overestimated with extended Hubbard interactions. As
was discussed in our previous study [28], the current formalism depends on
pseudopotentials and XC functionals of DFT since the projectors are obtained
from PPs. It is noticeable that, regardless of XC functionals used for DFT
calculations, the contribution from the extended Hubbard interactions to
energy band gaps and interband gaps are more or less similar. However, FR-GGA
calculations significantly underestimate the energy gaps of these materials,
if compared with those from FR-LDA, so that resulting gaps from FRG+$U$+$V$
are smaller than values from FR-$GW$, FRL+$U$+$V$ and experimental values.

###  IV.2 CdSe and PdO

We show effects of extended Hubbard interactions on CdSe and PdO, all of which
have significant SOCs. CdSe has a wurtzite structure with a large SOC of 407
meV [65]. DFT with local or semilocal XC approximations erroneously computes
the ground electronic property of semiconducting PdO to be a metal [66] so
that extended Hubbard interactions to correct SIEs for states associated with
the $d$ orbitals in the presence of the Pd-O bondings seem to be important. In
Table 6, we summarize the self-consistent Hubbard parameters obtained using
with FRG+$U$+$V$. With these interactions for CdSe, we obtain the lattice
constants of $a=4.30$ and $c=6.99$ in Å, agreeing with experimental values of
4.30 Å and 7.01 Å [67] very well.

Table 6: Calculated $U$ and $V$ (in eV) with FRG+$U$+$V$. $U_{d}$ and $U_{p}$ represent the on-site Hubbard parameters for $d$-orbitals in Cd and Pd, and and $p$-orbitals in Se and O, respectively. $V_{dp}$ represents the inter-site Hubbard parameters between $d$\- and $p$-orbitals of the adjacent atoms. | $U_{d}$ | $U_{p}$ | $V_{dp}$
---|---|---|---
CdSe | 13.78 | 3.90 | 2.09
PdO | 5.54 | 2.69 | 3.25
Figure 2: FR-LDA and FRL+$U$+$V$ band structure of CdSe, plotted along high-
symmetry lines of the Brillouin zone. The zero of the energy is at the valence
band maximum.

In Fig. 2, we plot energy bands of CdSe and confirm that the conduction bands
shift up significantly with extended Hubbard interactions from those obtained
by FR-LDA. The band widths of valence bands are overall slightly reduced with
Hubbard interactions compared with those from FR-LDA. In Table 7, we summarize
the energy gap, spin-orbit splitting of the valence band maximum and the
crystal field splitting of the valence bands at $\Gamma$ point. The energy gap
from FRG+$U$+$V$ agrees with those from FR-$GW$ and experiment while the gap
from FRL+$U$+$V$ is overestimated. As discussed for Ge, GaAs, and GaSb in the
previous subsection, our estimations of the spin-orbit splittings are also
slightly overestimated if compared with experiments. On the other hand, when
the extended Hubbard correction is included, the magnitude of crystal field
splitting is a little bit smaller than experiments.

Table 7: The band gap, spin-orbit splitting, and crystal field (CF) splitting of valence band at $\Gamma$ point ($\Delta^{\Gamma}_{\text{CF},v}$) in eV for CdSe with FR-LDA, FRL+$U$+$V$, FR-GGA, FRG+$U$+$V$ compared to $GW$ calculations and experiment.  | $E_{g}$ | $\Delta^{\Gamma}_{v}$ | $\Delta^{\Gamma}_{\text{CF},v}$
---|---|---|---
FR-$GW$161616Reference [26] | 1.85 | 0.405 | 0.026
FR-LDA | 0.58 | 0.373 | 0.037
FRL+$U$+$V$ | 2.33 | 0.481 | 0.023
FR-GGA | 0.42 | 0.350 | 0.031
FRG+$U$+$V$ | 1.76 | 0.481 | 0.017
Exp.171717Reference [65] | 1.82 | 0.407 | 0.025
Table 8: Lattice constants and band gap of PdO computed with FR-GGA and FRG+$U$+$V$ by using FR-PPs compared with PBE and hybrid functional (HSE06) computed with SR-PPs and experiments. |  $a$ (Å) |  $c$ (Å) | Band gap (eV)
---|---|---|---
FR-GGA | 3.04 | 5.33 | 0
FRG-$U$+$V$ | 3.05 | 5.33 | 0.95
PBE181818Scalar relativistic calculations, Reference [66] | 3.14 | 5.47 | 0
HSE0618 | 3.05 | 5.47 | 0.71
Exp. | 3.03191919Reference [68] | 5.3319 | 0.8202020Reference [69], 2.13212121Reference [70], 2.6721

We summarize lattice constants and energy gap of PdO in Table 8 and plot
energy bands obtained from FR-GGA and FRG+$U$+$V$ in Fig. 3. Like other
materials discussed so far, the computed lattice constants with extended
Hubbard interactions show excellent agreements with an experiment [68].
Furthermore, with non-negligible inter-site interactions of 3.25 eV between Pd
$d$-orbitals and O $p$-orbitals in Table 6, we obtain energy gap of 0.95 eV
showing qualitative agreements with an experiment [70], in sharp contrast to
metallic ground states from FR-GGA and a scalar relativistic PBE calculation
[66]. We also note that bandwidths of valence bands in Fig. 3 significantly
increase from those using GGA.

Figure 3: FR-GGA and FRG+$U$+$V$ band structures of PdO, plotted along high-
symmetry lines of the Brillouin zone. The Fermi energy of metallic PdO from
FR-GGA and the valence band maximum from FRG+$U$+$V$ are set to zero.

###  IV.3 HgTe, CuTlS2, and CuTlSe2

For TIs, we first choose HgTe, CuTlS2 and CuTlSe2 to test our method because
they are interesting materials for false-positive or false-negative
assignments of topological properties depending on computational methods [41].
It has been shown that a screened hybrid functional of Heyd–Scuseria–Ernzerhof
(HSE06) [71] results in inaccurate or uncertain determinations on topological
characteristics of HgTe and CuTlS2 [41].

Following a previous study [41], we use the order of bands to distinguish
three-dimensional (3D) bulk TIs from normal insulators exploiting their band
inversion. The magnitude of band inversion is defined through
$\Delta_{i}=\varepsilon_{\textit{s}}-\varepsilon_{\textit{{p,d}}}$ being
negative at the time reversal invariant momentum. Here,
$\varepsilon_{\textit{s}}$ and $\varepsilon_{\textit{{p,d}}}$ represent the
s-like conduction and the (p,d)-like valence band energy, respectively.

Although DFT has been widely adopted to evaluate the topological class in
various materials [72, 73, 74, 75, 76, 77], it sometimes yields incorrect
results due to the band gap underestimation [41]. Alternative computational
methods such as HSE06 and the modified-Becke-Johnson (mBJ) potential are often
employed to obtain improved band gaps compared with DFT-GGA or DFT-LDA [78,
79, 80]. However, in cases of HgTe and CuTlS2, HSE06 produces qualitatively
different results compared to standard DFT and self-consistent $GW$
calculation (sc-$GW$) [41]. Similarly, the mBJ+$U$ approach, which combines
mBJ potential with Hubbard $U$, shows positive inversion energy in CuTlS2
[80]. Meanwhile, HSE06 and mBJ+$U$ give the same assignment of topological
property for CuTlSe2 if compared with other methods.

In Table 9, we compare $\Delta_{i}$ from various methods. Albeit smaller
magnitudes of $\Delta_{i}$ than those from sc-$GW$, our estimations of
inversion energies for all three TIs remain to be negative, thus providing
qualitatively correct assessments of their topological properties. This
highlights important roles of site-dependent variations of on-site and inter-
site Hubbard interactions in our method in contrast to the constant mixing of
exchanges in typical hybrid approaches as was already highlighted in a
previous study [28] for electronic structures of the reconstructed Si(111)
surface.

Table 9: Calculated FR-GGA and FRG+$U$+$V$ inversion energies ($\Delta_{i}$) compared with HSE06 and sc-$GW$ results for HgTe, CuTlS2 and CuTlSe2. Here we use PBE as the XC functional instead of PBEsol. |  $\Delta_{i}$ (eV)
---|---
| FR-GGA | sc-$GW$ 222222Reference [41] | HSE06 22 | FRG+$U$+$V$
HgTe | $-1.13$ | $-0.26$ | $\sim 0$ | $-0.16$
CuTlS2 | $-0.72$ | $-0.41$ | 0.06 | $-0.2$
CuTlSe2 | $-1.02$ | $-0.85$ | $-0.4$ | $-0.43$

###  IV.4 Bi2Se3 and Bi2Te3

In Figs 4 and 5, we present bulk energy bands of Bi2Se3 and Bi2Te3, both of
which are well-known 3D TIs with strong SOCs [81]. We used the experimental
lattice parameters and atomic coordinates following previous works [82, 26,
83]. Here, the strong SOC in bismuth atom causes large spin-orbit splitting in
bands around the $\Gamma$ point. This brings band inversions between the
valence and conduction bands across the energy gap, mixing valence and
conduction states near the $\Gamma$ point, and ultimately resulting in a
nontrivial Z2 topological index [81]. The size of energy gap after band
inversions is proportional to the degree of mixing between the valence and
conduction states. The original band gaps without SOCs are underestimated in
standard DFT so that the mixing results in ‘camelback’-shaped DFT bands around
their charge neutral energies [81, 84]. This makes the DFT energy gap of 3D
bulk Bi2Se3 and Bi2Te3 be indirect as shown in DFT bands in Fig. 4(b) and
5(b).

Figure 4: FR-LDA and FRL+$U$+$V$ band structures of bulk Bi2Se3 plotted along
(a) high-symmetry lines of the Brillouin zone, (b) enlarged ones near the
$\Gamma$ point and (c) surface energy bands in the vicinity of the $\Gamma$
point as a function of slab thickness. (d) Band gap of topological surface
bands at the $\Gamma$ point as a function of slab thickness. $GW$+SOC and
experimental results are reproduced from Refs. [84] and [85], respectively.
The valence band maximum is set to zero in (a) and (b) while in (c), the Fermi
energy is set to 0 eV.

For Bi2Se3, however, our FRL+$U$+$V$ calculation obtains parabolic valence and
conduction bands as well as the direct band gap as shown in Figs. 4(b). This
agrees with previous $GW$+SOC [84] and FR-$GW$ [26] studies very well. The
computed band gap is 0.25 eV, smaller than 0.30 eV from $GW$+SOC [84] and 0.38
eV from FR-$GW$ [26]. Our computed gaps are comparable to various experimental
estimations, e.g., 0.332 eV obtained from angle-resoled photoemission
spectroscopy (ARPES) [86], 0.3 eV from scanning tunneling spectroscopy (STS)
[87] and 0.2 eV from optical measurements [88]. We also estimate the effective
masses from FRL+$U$+$V$ energy bands. The effective masses for the holes and
electrons are 0.10 $m_{e}$ and 0.09 $m_{e}$, respectively. These values are in
line with 0.19 $m_{e}$ (hole) and 0.14 $m_{e}$ (electron) from FR-$GW$ and are
consistent with the 0.14 $m_{e}$ for both electrons and holes from magneto-
optic measurements [88].

Figure 5: FR-GGA and FRG+$U$+$V$ band structures of bulk Bi2Te3 along (a)
high-symmetry lines of the Brillouin zone and (b) near the $\Gamma$ point and
(c) surface bands in the vicinity of the $\Gamma$ point as a function of slab
thickness. (d) Surface band gap at the $\Gamma$ point as a function of slab
thickness. $GW$ calculation results are reproduced from Ref. [84]. The valence
band maximum is set to zero in (a) and (b) while in (c), the Fermi energy is
set to 0 eV.

Since DFT+$U$+$V$ calculation requires moderate computational cost compared to
standard DFT, we test our method to calculate variation of topological surface
bands (TSBs) with varying thickness of slab geometries. To construct supercell
geometries of (111)-directional slabs with increasing thickness from one to
eight QLs, atomic positions are taken from experimental data [83] like bulk
calculations. Figures 4(c) show the evolutions of TSBs in the vicinity of the
$\Gamma$ point as a function of thickness. It is well known that as decreasing
thickness of TIs, the TSBs on opposite surfaces hybrid each other and open a
gap [89, 90, 91, 85]. Such a gap of massive TSBs in ultrathin Bi2Se3 decreases
as its thickness increase and eventually disappear to form massless TSBs
beyond a critical thickness [89, 90, 91, 85].

We note that within FR-LDA, TBSs of Bi2Se3 become to be massless beyond three
or four QLs [89, 84]. However, with extended Hubbard interactions, the gaps of
massive TSBs are greatly enhanced as shown in Fig. 4(d). Such a trend in
variation of surface band gaps agrees with $GW$+SOC calculations [84] and with
an experiment [85]. Notably, for slabs thicker than three QL, our FRL+$U$+$V$
calculations showing finite gaps closely match the experimental values. From
our method, the critical thickness for gap closure in TSBs of Bi2Se3 is eight
QLs of which gap is below 10 meV. For the massless dispersion, we estimate the
group velocity of 3.6 eV$\cdot$Å, being close to the experimental values, 3.6
eV$\cdot$Å [92] and 3.2 eV$\cdot$Å [93].

For Bi2Te3, our FRL+$U$+$V$ energy bands in Figs. 5 (a) and (b) show a
reduction of the dip in the valence bands at the $\Gamma$ point compared to
the band structures from FR-GGA. Unlike Bi2Se3, the valence band top is more
or less flat owing to a remnant corrugation of bands from mixing. The
resulting indirect gaps from FR-GGA and FRG+$U$+$V$ are 0.11 and 0.25 eV,
respectively, with the latter being overestimated if compared with 0.165 eV
from ARPES measurement [94] and 0.17 eV from $GW$+SOC [84]. For Bi2Te3, we use
PBE as a XC functional instead of PBEsol.

For thin slabs of Bi2Te3, the evolution of surface band gaps are similarly
dependent of computational methods as was discussed above. One notable
difference is an oscillatory decrease of gap [89, 90, 91] shown in Fig. 5 (d)
such that the gap of TSBs in six QLs is slightly higher than that for slabs
with five QLs. From the linear dispersion near the $\Gamma$ point in Fig. 5
(c), the group velocities of the TSBs are estimated to be 3.0 eV$\cdot$Å along
the $\Gamma$-M and 3.1 eV$\cdot$Å along the $\Gamma$-K direction. These values
are close to 2.2 (2.6) eV$\cdot$Å along the $\Gamma$-M direction and 2.7
eV$\cdot$Å along the $\Gamma$-K direction from experimental measurements [94,
95].

##  V Conclusion

We have implemented the fully relativistic extended Hubbard functionals to
calculate physical properties of materials having the spin-orbit coupling and
the local and non-local Coulomb interactions simultaneously. For on-site and
inter-site interactions, we determine them self-consistently within the DFT
routine by extending the ACBN0 approach [38, 28, 29] to construct inter-site
Hubbard functionals for noncollinear states.

We applied our method to calculate atomic and electronic structures of several
characteristic semiconductors and topological insulators with various
strengths of SOCs. Overall, all our results for those materials show quite
good agreements with fully relativistic $GW$ level calculations as well as
experiments. Though our new method significantly improves results from simple
XC approximations in DFT, in some cases, there are noticeable dependency of
computed energy gaps on XC functionals and pseudopotentials as discussed for
the cases of III-V semiconductors.

Considering that current formalisms without SOCs already have provided
fruitful results on phonon dispersions as well as electron-phonon interactions
[59, 31, 32, 36], our new extension to noncollinear DFT+$U$+$V$ is expected to
be a reliable and accurate method to calculate lattice dynamical properties
and electron-phonon interactions of materials with heavy atomic elements.

###### Acknowledgements.

We thank Bradford A. Barker for providing us with fully relativistic norm
conserving LDA pseudopotentials of Si, Ge, Ga, P, As, Bi, Se and Sb atoms and
Myung Joon Han for discussions. Y.-W.S. was supported by KIAS individual Grant
(No. CG031509). W.Y. was supported by KIAS individual Grant (No. 6P090103).
Computations were supported by the CAC of KIAS.

## References

  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Ihm [1988] J. Ihm, Total energy calculations in solid state physics, Rep. Prog. Phys. 51, 105 (1988).
  * Jones [2015] R. O. Jones, Density functional theory: Its origins, rise to prominence, and future, Rev. Mod. Phys. 87, 897 (2015).
  * Lejaeghere  _et al._ [2016] K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. D. Corso, S. de Gironcoli, T. Deutsch, J. K. Dewhurst, I. D. Marco, C. Draxl, M. Dułak, O. Eriksson, J. A. Flores-Livas, K. F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E. K. U. Gross, A. Gulans, F. Gygi, D. R. Hamann, P. J. Hasnip, N. A. W. Holzwarth, D. Iuşan, D. B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçükbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordström, T. Ozaki, L. Paulatto, C. J. Pickard, W. Poelmans, M. I. J. Probert, K. Refson, M. Richter, G.-M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunström, A. Tkatchenko, M. Torrent, D. Vanderbilt, M. J. van Setten, V. V. Speybroeck, J. M. Wills, J. R. Yates, G.-X. Zhang, and S. Cottenier, Reproducibility in density functional theory calculations of solids, Science 351, aad3000 (2016).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Cococcioni and de Gironcoli [2005] M. Cococcioni and S. de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the $\mathrm{LDA}+\mathrm{U}$ method, Phys. Rev. B 71, 035105 (2005).
  * Cohen  _et al._ [2008] A. J. Cohen, P. Mori-Sánchez, and W. Yang, Insights into current limitations of density functional theory, Science 321, 792 (2008).
  * Anisimov  _et al._ [1991] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Band theory and Mott insulators: Hubbard U instead of Stoner I, Phys. Rev. B 44, 943 (1991).
  * Anisimov  _et al._ [1993] V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, Density-functional theory and NiO photoemission spectra, Phys. Rev. B 48, 16929 (1993).
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B 57, 1505 (1998).
  * Himmetoglu  _et al._ [2014] B. Himmetoglu, A. Floris, S. de Gironcoli, and M. Cococcioni, Hubbard-corrected DFT energy functionals: The $\text{LDA}+\text{U}$ description of correlated systems, Int. J. Quant. Chem. 114, 14 (2014).
  * Kim  _et al._ [2008] B. J. Kim, H. Jin, S. J. Moon, J.-Y. Kim, B.-G. Park, C. S. Leem, J. Yu, T. W. Noh, C. Kim, S.-J. Oh, J.-H. Park, V. Durairaj, G. Cao, and E. Rotenberg, Novel ${J}_{\mathrm{eff}}=1/2$ Mott State Induced by Relativistic Spin-Orbit Coupling in ${\mathrm{Sr}}_{2}{\mathrm{IrO}}_{4}$, Phys. Rev. Lett. 101, 076402 (2008).
  * Wan  _et al._ [2011] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Topological semimetal and fermi-arc surface states in the electronic structure of pyrochlore iridates, Phys. Rev. B 83, 205101 (2011).
  * Zhang  _et al._ [2012] X. Zhang, H. Zhang, J. Wang, C. Felser, and S.-C. Zhang, Actinide topological insulator materials with strong interaction, Science 335, 1464 (2012).
  * Erickson  _et al._ [2007] A. S. Erickson, S. Misra, G. J. Miller, R. R. Gupta, Z. Schlesinger, W. A. Harrison, J. M. Kim, and I. R. Fisher, Ferromagnetism in the Mott Insulator ${\mathrm{Ba}}_{2}{\mathrm{NaOsO}}_{6}$, Phys. Rev. Lett. 99, 016404 (2007).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Marini  _et al._ [2009] A. Marini, C. Hogan, M. Grüning, and D. Varsano, Yambo: An ab initio tool for excited state calculations, Comp. Phys. Commun. 180, 1392 (2009).
  * Friedrich  _et al._ [2010] C. Friedrich, S. Blügel, and A. Schindlmayr, Efficient implementation of the $GW$ approximation within the all-electron FLAPW method, Phys. Rev. B 81, 125102 (2010).
  * Sakuma  _et al._ [2011] R. Sakuma, C. Friedrich, T. Miyake, S. Blügel, and F. Aryasetiawan, $GW$ calculations including spin-orbit coupling: Application to Hg chalcogenides, Phys. Rev. B 84, 085144 (2011).
  * Kühn and Weigend [2015] M. Kühn and F. Weigend, One-Electron Energies from the Two-Component GW Method, J. Chem. Theory Comput. 11, 969 (2015).
  * Scherpelz  _et al._ [2016] P. Scherpelz, M. Govoni, I. Hamada, and G. Galli, Implementation and Validation of Fully Relativistic GW Calculations: Spin–Orbit Coupling in Molecules, Nanocrystals, and Solids, J. Chem. Theory and Comput. 12, 3523 (2016).
  * Holzer and Klopper [2019] C. Holzer and W. Klopper, Ionized, electron-attached, and excited states of molecular systems with spin–orbit coupling: Two-component GW and Bethe–Salpeter implementations, J. Chem. Phys. 150, 204116 (2019).
  * Pashov  _et al._ [2020] D. Pashov, S. Acharya, W. R. Lambrecht, J. Jackson, K. D. Belashchenko, A. Chantis, F. Jamet, and M. van Schilfgaarde, Questaal: A package of electronic structure methods based on the linear muffin-tin orbital technique, Comp. Phys. Commun. 249, 107065 (2020).
  * Marsili  _et al._ [2021] M. Marsili, A. Molina-Sánchez, M. Palummo, D. Sangalli, and A. Marini, Spinorial formulation of the $GW$-BSE equations and spin properties of excitons in two-dimensional transition metal dichalcogenides, Phys. Rev. B 103, 155152 (2021).
  * Barker  _et al._ [2022] B. A. Barker, J. Deslippe, J. Lischner, M. Jain, O. V. Yazyev, D. A. Strubbe, and S. G. Louie, Spinor $GW$/Bethe-Salpeter calculations in BerkeleyGW: Implementation, symmetries, benchmarking, and performance, Phys. Rev. B 106, 115127 (2022).
  * Campo and Cococcioni [2010] V. L. Campo and M. Cococcioni, Extended DFT+$U$+$V$ method with on-site and inter-site electronic interactions, J. Phys.: Condens. Matter 22, 055602 (2010).
  * Lee and Son [2020] S.-H. Lee and Y.-W. Son, First-principles approach with a pseudohybrid density functional for extended Hubbard interactions, Phys. Rev. Research 2, 043410 (2020).
  * Tancogne-Dejean and Rubio [2020] N. Tancogne-Dejean and A. Rubio, Parameter-free hybridlike functional based on an extended Hubbard model: $\mathrm{DFT}+\text{U}+\text{V}$, Phys. Rev. B 102, 155117 (2020).
  * Timrov  _et al._ [2021] I. Timrov, N. Marzari, and M. Cococcioni, Self-consistent Hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formulations, Phys. Rev. B 103, 045141 (2021).
  * Yang  _et al._ [2021] W. Yang, S.-H. Jhi, S.-H. Lee, and Y.-W. Son, Ab initio study of lattice dynamics of group IV semiconductors using pseudohybrid functionals for extended Hubbard interactions, Phys. Rev. B 104, 104313 (2021).
  * Yang  _et al._ [2022] W. Yang, B. G. Jang, Y.-W. Son, and S.-H. Jhi, Lattice dynamical properties of antiferromagnetic oxides calculated using self-consistent extended Hubbard functional method, J. Phys.: Condens. Matter 34, 295601 (2022).
  * Ricca  _et al._ [2020] C. Ricca, I. Timrov, M. Cococcioni, N. Marzari, and U. Aschauer, Self-consistent $\mathrm{DFT}+U+V$ study of oxygen vacancies in ${\mathrm{SrTiO}}_{3}$, Phys. Rev. Res. 2, 023313 (2020).
  * Timrov  _et al._ [2022] I. Timrov, F. Aquilante, M. Cococcioni, and N. Marzari, Accurate Electronic Properties and Intercalation Voltages of Olivine-Type Li-Ion Cathode Materials from Extended Hubbard Functionals, PRX Energy 1, 033003 (2022).
  * Mahajan  _et al._ [2021] R. Mahajan, I. Timrov, N. Marzari, and A. Kashyap, Importance of intersite Hubbard interactions in $\beta\text{$-$}{\mathrm{MnO}}_{2}$: A first-principles $\mathrm{DFT}+U+V$ study, Phys. Rev. Mater. 5, 104402 (2021).
  * Jang  _et al._ [2023] B. G. Jang, M. Kim, S.-H. Lee, W. Yang, S.-H. Jhi, and Y.-W. Son, Intersite coulomb interactions in charge-ordered systems, Phys. Rev. Lett. 130, 136401 (2023).
  * Gatti  _et al._ [2020] G. Gatti, A. Crepaldi, M. Puppin, N. Tancogne-Dejean, L. Xian, U. De Giovannini, S. Roth, S. Polishchuk, P. Bugnon, A. Magrez, H. Berger, F. Frassetto, L. Poletto, L. Moreschini, S. Moser, A. Bostwick, E. Rotenberg, A. Rubio, M. Chergui, and M. Grioni, Light-Induced Renormalization of the Dirac Quasiparticles in the Nodal-Line Semimetal ZrSiSe, Phys. Rev. Lett. 125, 076401 (2020).
  * Agapito  _et al._ [2015] L. A. Agapito, S. Curtarolo, and M. Buongiorno Nardelli, Reformulation of $\mathrm{DFT}+U$ as a Pseudohybrid Hubbard Density Functional for Accelerated Materials Discovery, Phys. Rev. X 5, 011006 (2015).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Vidal  _et al._ [2011] J. Vidal, X. Zhang, L. Yu, J.-W. Luo, and A. Zunger, False-positive and false-negative assignments of topological insulators in density functional theory and hybrids, Phys. Rev. B 84, 041109 (2011).
  * Hybertsen and Louie [1986] M. S. Hybertsen and S. G. Louie, Spin-orbit splitting in semiconductors and insulators from the ab initio pseudopotential, Phys. Rev. B 34, 2920 (1986).
  * Hemstreet  _et al._ [1993] L. A. Hemstreet, C. Y. Fong, and J. S. Nelson, First-principles calculations of spin-orbit splittings in solids using nonlocal separable pseudopotentials, Phys. Rev. B 47, 4238 (1993).
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).
  * Garrity  _et al._ [2014] K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Pseudopotentials for high-throughput dft calculations, Comput. Mater. Sci. 81, 446 (2014).
  * Binci and Marzari [2023] L. Binci and N. Marzari, Noncollinear $\mathrm{DFT}+U$ and Hubbard parameters with fully relativistic ultrasoft pseudopotentials, Phys. Rev. B 108, 115157 (2023).
  * Löwdin [1950] P. Löwdin, On the Non‐Orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals, J. Chem. Phys. 18, 365 (1950).
  * Mosey and Carter [2007] N. J. Mosey and E. A. Carter, Ab initio evaluation of Coulomb and exchange parameters for $\mathrm{DFT}+\mathrm{U}$ calculations, Phys. Rev. B 76, 155123 (2007).
  * Mosey  _et al._ [2008] N. J. Mosey, P. Liao, and E. A. Carter, Rotationally invariant ab initio evaluation of Coulomb and exchange parameters for DFT+U calculations, J. Chem. Phys. 129, 014103 (2008).
  * Tancogne-Dejean  _et al._ [2017] N. Tancogne-Dejean, M. J. T. Oliveira, and A. Rubio, Self-consistent $\mathrm{DFT}+U$ method for real-space time-dependent density functional theory calculations, Phys. Rev. B 96, 245133 (2017).
  * Liechtenstein  _et al._ [1995] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators, Phys. Rev. B 52, R5467 (1995).
  * Czyżyk and Sawatzky [1994] M. T. Czyżyk and G. A. Sawatzky, Local-density functional and on-site correlations: The electronic structure of ${\mathrm{La}}_{2}$${\mathrm{CuO}}_{4}$ and ${\mathrm{LaCuO}}_{3}$, Phys. Rev. B 49, 14211 (1994).
  * Solovyev  _et al._ [1994] I. V. Solovyev, P. H. Dederichs, and V. I. Anisimov, Corrected atomic limit in the local-density approximation and the electronic structure of d impurities in Rb, Phys. Rev. B 50, 16861 (1994).
  * Bultmark  _et al._ [2009] F. Bultmark, F. Cricchio, O. Grånäs, and L. Nordström, Multipole decomposition of $\text{LDA}+U$ energy and its application to actinide compounds, Phys. Rev. B 80, 035121 (2009).
  * Ryee and Han [2018] S. Ryee and M. J. Han, The effect of double counting, spin density, and Hund interaction in the different DFT+U functionals, Sci. Rep. 8, 9559 (2018).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * van Setten  _et al._ [2018] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, Comput. Phys. Commun. 226, 39 (2018).
  * Timrov  _et al._ [2020] I. Timrov, F. Aquilante, L. Binci, M. Cococcioni, and N. Marzari, Pulay forces in density-functional theory with extended Hubbard functionals: From nonorthogonalized to orthogonalized manifolds, Phys. Rev. B 102, 235159 (2020).
  * Madelung [2004] O. Madelung,  _Semiconductors: data handbook_ (Springer Science & Business Media, 2004).
  * Malone and Cohen [2013] B. D. Malone and M. L. Cohen, Quasiparticle semiconductor band structures including spin–orbit interactions, J. Phys. Condens. Matter 25, 105503 (2013).
  * Cardona and Thewalt [2005] M. Cardona and M. L. W. Thewalt, Isotope effects on the optical spectra of semiconductors, Rev. Mod. Phys. 77, 1173 (2005).
  * Masovic  _et al._ [1983] D. R. Masovic, F. R. Vukajlovic, and S. Zekovic, Local-pseudopotential calculation for optical properties and photoemission valence-band spectrum of silicon, J. Phys. C: Solid State Phys. 16, 6731 (1983).
  * Lautenschlager  _et al._ [1987] P. Lautenschlager, M. Garriga, S. Logothetidis, and M. Cardona, Interband critical points of GaAs and their temperature dependence, Phys. Rev. B 35, 9174 (1987).
  * Pedrotti and Reynolds [1962] F. L. Pedrotti and D. C. Reynolds, Spin-Orbit Splitting in CdS: Se Single Crystals, Phys. Rev. 127, 1584 (1962).
  * Kansara  _et al._ [2016] S. Kansara, D. Singh, S. K. Gupta, and Y. Sonvane, Density functional studies of structural, electronic and vibrational properties of palladium oxide, Solid State Commun. 245, 36 (2016).
  * Hotje  _et al._ [2003] U. Hotje, C. Rose, and M. Binnewies, Lattice constants and molar volume in the system ZnS, ZnSe, CdS, CdSe, Solid State Sci. 5, 1259 (2003).
  * Waser  _et al._ [1953] J. Waser, H. A. Levy, and S. W. Peterson, The structure of PdO, Acta Crystallogr. 6, 661 (1953).
  * Nilsson [1979] P. O. Nilsson, Optical properties of PdO in the range of 0.5-5.4 eV, J. Phys. C: Solid State Phys. 12, 1423 (1979).
  * Rey  _et al._ [1978] E. Rey, M. R. Kamal, R. B. Miles, and B. S. H. Royce, The semiconductivity and stability of palladium oxide, J. Mat. Sci. 13, 812 (1978).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened Coulomb potential, J. Chem. Phys. 118, 8207 (2003).
  * Chadov  _et al._ [2010] S. Chadov, X. Qi, J. Kübler, G. H. Fecher, C. Felser, and S. C. Zhang, Tunable multifunctional topological insulators in ternary heusler compounds, Nat. Mater. 9, 541 (2010).
  * Lin  _et al._ [2010] H. Lin, L. A. Wray, Y. Xia, S. Xu, S. Jia, R. J. Cava, A. Bansil, and M. Z. Hasan, Half-Heusler ternary compounds as new multifunctional experimental platforms for topological quantum phenomena, Nat. Mater. 9, 546 (2010).
  * Xiao  _et al._ [2010] D. Xiao, Y. Yao, W. Feng, J. Wen, W. Zhu, X.-Q. Chen, G. M. Stocks, and Z. Zhang, Half-heusler compounds as a new class of three-dimensional topological insulators, Phys. Rev. Lett. 105, 096404 (2010).
  * Sun  _et al._ [2010] Y. Sun, X.-Q. Chen, S. Yunoki, D. Li, and Y. Li, New family of three-dimensional topological insulators with antiperovskite structure, Phys. Rev. Lett. 105, 216406 (2010).
  * Zhang  _et al._ [2011] H.-J. Zhang, S. Chadov, L. Müchler, B. Yan, X.-L. Qi, J. Kübler, S.-C. Zhang, and C. Felser, Topological insulators in ternary compounds with a honeycomb lattice, Phys. Rev. Lett. 106, 156402 (2011).
  * Yang  _et al._ [2017] H. Yang, J. Yu, S. S. P. Parkin, C. Felser, C.-X. Liu, and B. Yan, Prediction of Triple Point Fermions in Simple Half-Heusler Topological Insulators, Phys. Rev. Lett. 119, 136401 (2017).
  * Chen  _et al._ [2011] S. Chen, X. G. Gong, C.-G. Duan, Z.-Q. Zhu, J.-H. Chu, A. Walsh, Y.-G. Yao, J. Ma, and S.-H. Wei, Band structure engineering of multinary chalcogenide topological insulators, Phys. Rev. B 83, 245202 (2011).
  * Kim  _et al._ [2017] J. Kim, J. Kim, Y.-S. Song, R. Wu, S.-H. Jhi, and N. Kioussis, Weyl node assisted conductivity switch in interfacial phase-change memory with van der Waals interfaces, Phys. Rev. B 96, 235304 (2017).
  * Ke  _et al._ [2023] D. Ke, L. Li, and Y. Zhang, Relativistic effects in complex compounds containing heavy elements: ternary Cu–Tl–X (X = S/Se/Te), Phys. Chem. Chem. Phys. 25, 7445 (2023).
  * Zhang  _et al._ [2009] H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Topological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surface, Nat. Phys. 5, 438 (2009).
  * Nakajima [1963] S. Nakajima, The crystal structure of Bi2Te3-xSex, J. Phys. .Chem. Solids 24, 479 (1963).
  * Wyckoff [1964] R. W. G. Wyckoff,  _Crystal structures_ , 2nd ed., Vol. 2 (John Wiley and Sons, 1964).
  * Yazyev  _et al._ [2012] O. V. Yazyev, E. Kioupakis, J. E. Moore, and S. G. Louie, Quasiparticle effects in the bulk and surface-state bands of Bi2Se3 and Bi2Te3 topological insulators, Phys. Rev. B 85, 161101 (2012).
  * Zhang  _et al._ [2010] Y. Zhang, K. He, C.-Z. Chang, C.-L. Song, L.-L. Wang, X. Chen, J.-F. Jia, Z. Fang, X. Dai, W.-Y. Shan, S.-Q. Shen, Q. Niu, X.-L. Qi, S.-C. Zhang, X.-C. Ma, and Q.-K. Xue, Crossover of the three-dimensional topological insulator Bi2Se3 to the two-dimensional limit, Nat. Phys. 6, 584 (2010).
  * Nechaev  _et al._ [2013] I. A. Nechaev, R. C. Hatch, M. Bianchi, D. Guan, C. Friedrich, I. Aguilera, J. L. Mi, B. B. Iversen, S. Blügel, P. Hofmann, and E. V. Chulkov, Evidence for a direct band gap in the topological insulator Bi2Se3 from theory and experiment, Phys. Rev. B 87, 121111 (2013).
  * Kim  _et al._ [2011] S. Kim, M. Ye, K. Kuroda, Y. Yamada, E. E. Krasovskii, E. V. Chulkov, K. Miyamoto, M. Nakatake, T. Okuda, Y. Ueda, K. Shimada, H. Namatame, M. Taniguchi, and A. Kimura, Surface Scattering via Bulk Continuum States in the 3D Topological Insulator ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$, Phys. Rev. Lett. 107, 056803 (2011).
  * Orlita  _et al._ [2015] M. Orlita, B. A. Piot, G. Martinez, N. K. S. Kumar, C. Faugeras, M. Potemski, C. Michel, E. M. Hankiewicz, T. Brauner, i. c. v. Drašar, S. Schreyeck, S. Grauer, K. Brunner, C. Gould, C. Brüne, and L. W. Molenkamp, Magneto-Optics of Massive Dirac Fermions in Bulk ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$, Phys. Rev. Lett. 114, 186401 (2015).
  * Linder  _et al._ [2009] J. Linder, T. Yokoyama, and A. Sudbø, Anomalous finite size effects on surface states in the topological insulator ${\text{Bi}}_{2}{\text{Se}}_{3}$, Phys. Rev. B 80, 205401 (2009).
  * Liu  _et al._ [2010] C.-X. Liu, H. Zhang, B. Yan, X.-L. Qi, T. Frauenheim, X. Dai, Z. Fang, and S.-C. Zhang, Oscillatory crossover from two-dimensional to three-dimensional topological insulators, Phys. Rev. B 81, 041307 (2010).
  * Lu  _et al._ [2010] H.-Z. Lu, W.-Y. Shan, W. Yao, Q. Niu, and S.-Q. Shen, Massive Dirac fermions and spin physics in an ultrathin film of topological insulator, Phys. Rev. B 81, 115407 (2010).
  * Kuroda  _et al._ [2010] K. Kuroda, M. Arita, K. Miyamoto, M. Ye, J. Jiang, A. Kimura, E. E. Krasovskii, E. V. Chulkov, H. Iwasawa, T. Okuda, K. Shimada, Y. Ueda, H. Namatame, and M. Taniguchi, Hexagonally Deformed Fermi Surface of the 3D Topological Insulator ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$, Phys. Rev. Lett. 105, 076802 (2010).
  * Zhu  _et al._ [2011] Z.-H. Zhu, G. Levy, B. Ludbrook, C. N. Veenstra, J. A. Rosen, R. Comin, D. Wong, P. Dosanjh, A. Ubaldini, P. Syers, N. P. Butch, J. Paglione, I. S. Elfimov, and A. Damascelli, Rashba Spin-Splitting Control at the Surface of the Topological Insulator ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$, Phys. Rev. Lett. 107, 186405 (2011).
  * Chen  _et al._ [2009] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, Experimental Realization of a Three-Dimensional Topological Insulator, Bi2Te3, Science 325, 178 (2009).
  * Li  _et al._ [2010] Y.-Y. Li, G. Wang, X.-G. Zhu, M.-H. Liu, C. Ye, X. Chen, Y.-Y. Wang, K. He, L.-L. Wang, X.-C. Ma, H.-J. Zhang, X. Dai, Z. Fang, X.-C. Xie, Y. Liu, X.-L. Qi, J.-F. Jia, S.-C. Zhang, and Q.-K. Xue, Intrinsic Topological Insulator Bi2Te3 Thin Films on Si and Their Thickness Limit, Adv. Mat. 22, 4002 (2010).