\def\draftversion{false}
%\def\draftversion{true}
\RequirePackage{ifthen}
\ifthenelse{\equal{\draftversion}{true}}{
  \documentclass[aps,prb,10pt,galley,amsmath,amssymb,showpacs,
                superscriptaddress]{revtex4}
}{
  \documentclass[aps,prb,10pt,twocolumn,amsmath,amssymb,showpacs,
                 superscriptaddress]{revtex4-1}
}


%\documentclass[11pt,fleqn]{article}
%\usepackage[pdftex]{graphicx}
%\usepackage{epstopdf}
\usepackage{times}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage[usenames]{color} % colors
\usepackage{bm} % bold math
\usepackage{subfigure}
\newcommand{\angstrom}{\text{\normalfont\AA}}

\def\Z2{$\mathbb{Z}_2$}
\def\bisb{(Bi$_{1-x}$Sb$_{x})_2$Se$_3$}
\def\biin{(Bi$_{1-x}$In$_{x})_2$Se$_3$}
\def\bise{Bi$_2$Se$_3$}
\def\bite{Bi$_2$Te$_3$}
\def\inse{In$_2$Se$_3$}
\def\sbse{Sb$_2$Se$_3$}
\def\I{\uppercase\expandafter{\romannumeral 1}}
\def\II{\uppercase\expandafter{\romannumeral 2}}
\def\III{{\uppercase\expandafter{\romannumeral 3}}}
\def\IV{{\uppercase\expandafter{\romannumeral 4}}}
\def\V{{\uppercase\expandafter{\romannumeral 5}}}

% ----------------------------------------
% Better to add more macro definitions
% ----------------------------------------

\def\Pt{\widetilde{P}}
\def\Qt{\widetilde{Q}}

\def\ket#1{\vert#1\rangle}
\def\bra#1{\langle#1\vert}
\def\ip#1#2{\langle#1\vert#2\rangle}
\def\me#1#2#3{\langle#1\vert#2\vert#3\rangle}
\def\ev#1{\langle#1\rangle}
\def\dag{^\dagger}

\def\Tra#1{\textrm{Tr}[#1]}

\def\kk{\mathbf{k}}
\def\rr{\mathbf{r}}
\def\GG{\mathbf{G}}
\def\RR{\mathbf{R}}

% my personal notation is to use 'b' for 'suBscript' and 'p'
% for 'suPerscript' as the first letter of a macro name:

\def\bnk{_{n\kk}}
\def\bmk{_{m\kk}}
\def\bnpk{_{n'\kk}}
\def\bocc{_\mathrm{occ}}
\def\bso{_\mathrm{so}}
\def\bR{_\mathrm{R}}

\def\ut{\tilde{u}}
\def\wt{\tilde{w}}
\def\pt{\tilde{\psi}}
\def\Ct{\widetilde{C}}
% ----------------------------------------




\begin{document}

\title{Spin-orbit spillage as a measure of band inversion in insulators}
\date{July 4, 2014}

\author{Jianpeng Liu}
\affiliation{ Department of Physics and Astronomy, Rutgers University,
 Piscataway, NJ 08854-8019, USA }

\author{David Vanderbilt}
\affiliation{ Department of Physics and Astronomy, Rutgers University,
 Piscataway, NJ 08854-8019, USA }

\date{\today}

\begin{abstract}
We propose a straightforward and effective approach for quantifying
the band inversion induced by spin-orbit coupling in band
insulators.  In this approach we define a quantity as
a function of wavevector in the Brillouin zone measuring
the spillage between the occupied states of a system with and
without SOC.  Plots of the spillage throughout the BZ provide
a ready indication of the number and location of band inversions
driven by SOC.  To illustrate the method, we apply this approach to
the 2D Kane-Mele model, a 2D Bi bilayer with applied Zeeman field,
and to first-principles calculations of some
3D materials including both trivial and \Z2\ topological
insulators. We argue that the distribution of spillage in the BZ is
closely related to the topological indices in these systems.
Our approach provides a fresh perspective for understanding
topological character in band theory, and should be helpful in
searching for new materials with non-trivial band topology.
\end{abstract}

\pacs{71.70.Ej, 03.65.Vf}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\marginparwidth 2.7in
%\marginparsep 0.5in
%\def\dvm#1{\marginpar{\small DV: #1}}
%\def\jlm#1{\marginpar{\small JPL: #1}}
%\def\scr{\scriptsize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%-------- MARGIN COMMENTS --------------
\def\scr{\scriptsize}
\ifthenelse{\equal{\draftversion}{true}}{
  \marginparwidth 2.7in
  \marginparsep 0.5in
  \newcounter{comm} % counter for commentaries
  % increase counter
  \def\commnext{\stepcounter{comm}}
  % commentary in text
  \def\commtext{{\bf\color{blue}[\arabic{comm}]}}
  % commentary in margin
  \def\commmar{{\bf\color{blue}[\arabic{comm}]}}
  % comment commands for all authors
  \def\dvm#1{\commnext\marginpar{\small DV\commmar: #1}\commtext}
  \def\jlm#1{\commnext\marginpar{\small JPL\commmar: #1}\commtext}
  \def\mlab#1{\marginpar{\small\bf #1}}
  \def\tnewpage{\marginpar{\small Temporary newpage}\newpage}
}{
  \def\dvm#1{}
  \def\jlm#1{}
  \def\mlab#1{}
  \def\tnewpage{}
}
%----------------------------------


%-------- COLORS --------------
\def\Red#1{\textcolor{red}{#1}}
\def\Blue#1{\textcolor{blue}{[#1]}}
\def\Magenta#1{\textcolor{magenta}{#1}}
%----------------------------------

%------------------------------------------------
\section{Introduction}
\label{sec:intro}
%------------------------------------------------

Spin-orbit coupling (SOC) is a relativistic effect originating
from the interaction
between the spin and orbital motions of electrons.
It has played a key
role in various aspects of condensed-matter physics, including the
electronic structure of solids and the
transport properties in mesoscopic systems.
\cite{winkler2003,satofumi-negf} It has been known since the
1950s that SOC can induce anisotropic spin splitting
in some \III-\V\
semiconductors with the zinc-blende structure,
known as the Dresselhaus splitting.\cite{winkler2003}
In 2D and quasi-2D
systems, the SOC resulting from the electric
field perpendicular to the 2D plane
gives rise to a Rashba splitting linear in $k$ with
interesting ``helical'' spin textures.\cite{winkler2003,satofumi-negf}
The SOC is also crucial in determining the transport behavior of
low-dimensional electronic systems. One famous example
is the weak antilocalization
in spin-orbit-coupled 2D electronic systems,
where the backscattering
amplitudes interfere destructively due to a geometric Berry
phase\cite{berry84} associated with the intrinsic SOC,
leading to a suppressed resistivity when an
external magnetic field is absent.\cite{hln}
SOC is also responsible for
spin precession in 1D and quasi-1D
systems,\cite{satofumi-negf} the spin Hall effect in paramagnetic
metals,\cite{Hirsch99} and numerous other effects.

The SOC has received renewed attention recently because of its
central role in the physics of topological insulators (TIs) and
related topological states.  Typically, the transition from a
topologically trivial to a non-trivial phase is accomplished by
a SOC-driven inversion of states of different symmetry at the
conduction-band minimum (CBM) and valence-band maximum (VBM).
For example, such a SOC-driven topological band inversion
between $\Gamma_6$-derived ($s$-like) and $\Gamma_8$-derived
($p$-like) states at the zone center is responsible for the
quantum spin Hall (QSH) state observed in HgTe/CdTe quantum
wells.\cite{bhz06,konig06science}
Similarly, the Kane-Mele model of 2D graphene-like systems
\cite{kane-prl05-a,kane-prl05-b} enters the QSH state
when two band inversions occur at the K and K$'$ points as the
SOC strength is increased at a constant staggered potential.
In 3D band insulators with time-reversal (TR) symmetry,
a SOC-induced band inversion can transform the system from a
trivial insulator into a strong TI displaying an odd number of
gapless Dirac cones in the surface states, as occurs for
\bise\ and \bite.\cite{kane-rmp10,zhang-rmp11,yao-review12,zhang-np09}


In the case of a 3D strong TI with inversion symmetry such as
\bise, the strong \Z2\ index can be
uniquely determined by the parities of the occupied bands
at the TR-invariant momenta (TRIM) in the Brillouin zone
(BZ).\cite{fu-prb07} If the highest occupied states and lowest
unoccupied states at one of the TRIM possess opposite parities
without SOC, and they are inverted by turning on SOC, then the
system transforms from a normal to a topological insulator.
For example, in \bise, two pairs of Kramers-degenerate occupied states
at the BZ center ($\Gamma$) are inverted by SOC, resulting in
the nontrivial \Z2\ index.  For TIs without inversion symmetry,
the band inversion may happen at arbitrary points in the BZ,
instead of at the TRIM.  We can identify such band inversion
points as the points where a band touching occurs between valence
and conduction bands as the SOC is adiabatically turned on; TR
symmetry implies that an inversion at $\mathbf k_0$ will always
be accompanied by one at $-\mathbf k_0$.
Even in the absence of inversion
symmetry, therefore, a band inversion driven by SOC is typically
a hallmark of the non-trivial topology in TIs with TR symmetry.

The SOC also plays a crucial role in giving rise to the Chern insulator
(CI) state, also known as the quantum anomalous Hall state,
which can occur in 2D insulators lacking time-reversal symmetry.
The possibility of a CI state was first introduced by
Haldane,\cite{Haldane-model} who constructed an explicit model
that demonstrates the effect.  Although the Haldane model is a model
of spinless Fermions on a honeycomb lattice, its key feature is the
presence of complex second-neighbor hoppings, which can be regarded
as arising from intrinsic atomic SOC through a second-order
perturbation process in a more realistic system of spinor
electrons.\cite{hongbin-prb13} An example is a Bi bilayer with an applied 
Zeeman field, as will be discussed below.

The concept of topological band inversion
has been much discussed in the topological-insulator literature,
but in the absence of symmetry it may be difficult to recognize
when a band inversion has actually occurred.  The usual approach is to
look at the symmetry or orbital character at a high-symmetry
point in the BZ where a band inversion is suspected, but this
only works if sufficient symmetry is present.  Some authors
have tried to deduce the presence of band-inversion behavior by
studying other properties of the system, such as by looking at
the qualitative shape of the bands near the symmetry point,
\cite{klintenberg-archive10} or even more indirectly, by studying
the variation of the band-energy differences with
strain in the absence of SOC.\cite{yang-nm12} However, the
reliability of such methods is questionable, as they do not
give a direct and quantitative evaluation of the SOC-induced
band inversion.

In this paper, we propose that the calculation of spin-orbit
spillage, which measures the degree of mismatch between the
occupied band projection operators with and without SOC, provides
a simple and effective measure of SOC-driven band inversion
in insulators.  We demonstrate that the mapping of this spin-orbit spillage
in $k$-space easily allows a direct identification of any region
in the BZ where band inversion has occurred, and that the maximum
spillage is a useful indicator of topological character.
We illustrate the method in the context of both tight-binding
models and realistic first-principles calculations.

The paper is organized as follows. In Sec.~\II\ the
formal definition of SOC-induced spillage is introduced, and
the correspondence between topological indices
and spillage is also discussed.
In Sec.~\III\ the formalism is applied to various
systems, including the 2D Kane-Mele model, a Bi bilayer
with tunable SOC and exchange field, and realistic materials
including \bise, \inse, and \sbse.
In Sec.~\IV\ we make a brief summary.

%------------------------------------------------
\section{Formalism}
\label{sec:formal}
%------------------------------------------------

%------------------------------------------------
\subsection{Definitions}
\label{sec:def}
%------------------------------------------------

Mathematically, the spillage between two projection operators $P$
and $\Pt$, both of rank $N$, is defined as
%
\begin{equation}
\gamma=N-\Tra{P\Pt}=\Tra{P\Qt}=\Tra{Q\Pt}
\label{equa:spillage}
\end{equation}
%
where $Q=1-P$ and $\Qt=1-\Pt$ denote the complementary projections.
Clearly the spillage vanishes if $P=\Pt$ at one extreme, and rises
to $N$ at the other extreme if there is no overlap at all
between the subspaces associated
with $P$ and $\Pt$.  Thus, the spillage provides a measure of
the degree of mismatch between the two subspaces.

Here we apply this concept to the band projection operators
%
\begin{equation}
P(\kk)=\sum_{n=1}^{n\bocc} \ket{\psi\bnk}\bra{\psi\bnk}
\label{equa:bandproj}
\end{equation}
%
associated with a given wavevector $\kk$ in the BZ
of a crystalline insulator with $N\!=\!n\bocc$ occupied bands.
We assume
an effective single-particle Hamiltonian such as that appearing
in density-functional theory (DFT).\cite{dft1,dft2}
Then the SOC-induced spillage $\gamma(\kk)$ is defined as
%
\begin{equation}
\gamma(\kk)=\Tra{P(\kk) \Qt(\kk)}
\label{equa:def1}
\end{equation}
%
where $P$ and $\Pt$ (and their complements) refer to the
system with and without SOC respectively.
More explicitly,
%
\begin{align}
\gamma(\mathbf{k})& =n\bocc-\Tra{P(\mathbf{k})
\Pt(\mathbf{k})}\notag\\
& =n\bocc-\sum_{m,n=1}^{n\bocc}\vert M_{mn}(\kk)\vert^2
\label{equa:def2}
\end{align}
%
where
%
\begin{equation}
M_{mn}(\kk)= \ip{\psi\bmk}{\pt\bnk}
\label{equa:Mmn}
\end{equation}
%
is the overlap between occupied Bloch functions with and without SOC at the
same wavevector $\kk$.
Equivalently, this can be written as $M_{mn}(\kk)= \ip{u\bmk}{\ut\bnk}$
if one prefers to work in terms of the cell-periodic $\ket{u\bnk}$
defined as $u\bnk(\rr)=e^{-i\kk\cdot\rr}\psi\bnk(\rr)$.

In the case of realistic DFT calculations in a plane-wave basis,
the overlap matrix elements
are easily evaluated as
%
\begin{equation}
M_{mn}(\kk)=
\sum_{\GG}
\ip{\psi\bmk}{\kk+\GG} \ip{\kk+\GG}{\pt\bnk} \,,
\label{equa:plane-wave}
\end{equation}
%
where $\ket{\kk+\GG}$
is the plane wave $e^{i(\kk+\GG)\cdot\rr}$ for reciprocal vector $\GG$
normalized to the unit cell.
The evaluation should also be straightforward in other first-principles
basis sets.
For simple lattice models the Hamiltonian is typically
written in an orthonormal tight-binding basis, so that the
wavefunctions are
%
\begin{equation}
\ket{\psi\bnk}=\sum_j C_{nj,\kk}\,\ket{\chi_{j\kk}}
\label{equa:psichi}
\end{equation}
%
where $\ket{\chi_{j\kk}}$ are the Bloch basis states
%
\begin{equation}
\chi_{j\kk}(\rr)=\sum_{\RR} e^{i\kk\cdot\RR} \,\varphi_j(\rr-\RR)
\label{equa:chi}
\end{equation}
%
and $\varphi_j(\rr-\RR)$ is the $j$'th tight-binding basis orbital
in unit cell $\RR$.  Then the spillage is trivially computed using
\begin{equation}
M_{mn}(\kk)=\sum_j C^*_{mj,\kk} \Ct_{nj,\kk} \,.
\label{equa:tboverlap}
\end{equation}

Since the use of Wannier interpolation methods
\cite{MLWF-2,yates-prb07,MLWF-rmp}
is becoming increasingly frequent, we also comment on this case
here.  In this approach, the occupied Bloch states are again written
as in Eq.~(\ref{equa:psichi}), but this time the Bloch basis states are
%
\begin{equation}
\chi_{j\kk}(\rr)=\sum_{\RR} e^{i\kk\cdot\RR} \,w_j(\rr-\RR)
\label{equa:chiw}
\end{equation}
%
where $w_j(\rr-\RR)$ is the $j$'th Wannier function in unit cell $\RR$.
Then the spillage is again computed using Eqs.~(\ref{equa:def2}) and
(\ref{equa:tboverlap}).
This will be accurate as long as the WFs for the systems with and
without SOC are chosen to be the same, or as similar as possible.
As we shall see in the following section, the results from the
Wannier basis match those of the direct plane-wave calculation
very closely for the cases studied here.

%------------------------------------------------
\subsection{Relation to topological character}
\label{sec:topo}
%------------------------------------------------

Here we argue that the presence of non-trivial topological
indices will be reflected in certain features of the spillage
distribution in the BZ.

We first consider the relatively simple
case in which the SOC-driven band inversion is associated with
the crossing of highest valence and lowest conduction states
belonging to two different irreducible representations (irreps)
at a high-symmetry point $\kk\!=\!\Lambda_0$ in the BZ.
Since the states belonging to different irreps have no overlap with
each other,
the spillage at $\Lambda_0$ must be greater than or equal
to the irrep dimension.  In TR-invariant \bise, for example,
the four states around the Fermi level at $\Gamma$ consist of two
Kramers doublets of opposite parity.  In this case the dimension of the
irreps is two, so we expect a peak in $\gamma(\kk)$ centered at
$\Gamma$ whose height is $\gamma_{\rm max}\!\ge\!2$.
As we shall show in Sec.~\ref{sec:materials}, this is exactly what we find
in \bise.

Next, we argue that a correspondence between topological
character and spillage should also remain valid for more
general cases without special lattice symmetry.  Let us first
consider the case of CIs (i.e., with broken TR symmetry).
We assume the Bloch functions $\psi\bnk$ are those of a normal
system with Chern number $C\!=\!0$, while $\pt\bnk$ are topologically
nontrivial with a nonzero Chern number $\Ct$.  We argue that this
implies the existence of at least one point in the BZ where the
spillage is $\ge1$.  If we assume the contrary, i.e., $\gamma(\kk)<1$
everywhere in the BZ, then the determinant of the overlap matrix
of Eq.~(\ref{equa:Mmn}) between $\psi\bnk$ and $\pt\bnk$ obeys
$\det(M_\kk)>0$ everywhere, since
a singular $M$ would imply $\gamma\ge1$.  Because the system
$\ket{\psi\bnk}$ is topologically normal, we know it is possible
to choose a smooth and periodic gauge for it, and we assume without
loss of generality that this has been done.  But if $M_\kk$ is
nowhere singular, the $\ket{\psi\bnk}$ can be used as ``trial
functions'' to construct a smooth and periodic gauge for the
$\ket{\pt\bnk}$, as follows.  At each
$\kk$, carry out a singular value decomposition to express
$M=V^\dagger\Sigma W$ ($V$ and $W$ are unitary and $\Sigma$
is real positive diagonal), and then use the unitary matrix $V^\dagger W$
to transform the original $\pt\bnk$ to a new set $\pt'\bnk$.  Then
$M=V^\dagger\Sigma V$, i.e., it is Hermitian and positive definite.
Intuitively, this means that a smooth and periodic
gauge has been chosen for the
states $\pt'\bnk$ to make them ``maximally aligned'' with the
states $\psi\bnk$.  But a smooth and periodic gauge is inconsistent
with a nonzero Chern number, completing the proof by contradiction.
Thus, if $\gamma<1$ everywhere in the BZ, then $M_{mn,\kk}$ is
nonsingular everywhere, and the system $\ket{\pt\bnk}$ is normal.
Conversely, a topological system must have $\gamma(\kk)\ge1$
somewhere in the BZ, which provides both a signal for the topological
phase and an indication of where in the BZ the band inversion has
occurred.

For the TR-invariant \Z2\ TIs, similar arguments can be put forward
that work even in the absence of inversion symmetry. If the system
of $\ket{\psi\bnk}$ is in the \Z2-even\ phase, one can always make
a smooth gauge choice over the entire BZ that respects TR symmetry.
In the \Z2-odd\ case, however, such a gauge choice does not
exist.\cite{fu-prb06,soluyanov-prb12} Therefore, $\det(M_\kk)$ must
vanish somewhere in BZ, or else the smooth gauge could be transferred
to the $\ket{\pt\bnk}$, resulting in a contradiction.  Due to the
TR symmetry, $\det(M_\kk)=\det(M_{-\kk})$, so one would generically expect
$\gamma(\kk)\!\geq\!1$ at two points ($\kk_0$ and $-\kk_0$) in the
BZ.  For the case of inversion-symmetric TIs, $\kk_0$ and $-\kk_0$
merge at one of the TRIM, the two spillages add up, and one expects
$\gamma\ge2$ at one of the TRIM.

In the following section, we numerically test and confirm
the above arguments by applying the formalism to systems in different
topological phases.


%==============================================================
\section{Applications}
%==============================================================

%-------------------------------------------------
\subsection{Application to the Kane-Mele model}
\label{sec:kane-mele}
%-------------------------------------------------

The Kane-Mele model is a four-band TB model on a graphene lattice,
including nearest-neighbor (NN) spin-independent
hoppings and both NN and next-NN spin-dependent hoppings:
%
\begin{align}
H\!=\!& \sum_{\langle ij\rangle} t c^{\dagger}_{i} c_{j}+
\sum_{\langle\!\langle ij\rangle\!\rangle}i
\lambda\bso\nu_{ij}c^{\dagger}_{i}s_z c_{j}\notag\\
& +\sum_{\langle ij\rangle}i\lambda\bR c^{\dagger}_{i}
(\mathbf{s}\times\hat{\mathbf{d}}_{ij})_{z}c_{j}
+\sum_{i}{\epsilon(-1)^{i}c^{\dagger}_{i}c_{i}} \,.
\label{equa:kane-mele}
\end{align}
%
Here spin is implicit,
$t$ is the NN spin-independent hopping amplitude, $\lambda\bso$ is
the strength of the next-NN non-spin-flip SOC,
$\lambda\bR$ is the NN Rashba-like SOC amplitude, and $\epsilon$ is the
magnitude of on-site energy, with signs $\pm 1$ for A and B
sublattices respectively.  Also, $\nu_{ij}\!=\!\pm 1$ with the
sign depending on the chirality of the
next-NN bond from site $i$ to $j$, and $\hat{\mathbf{{d}}}_{ij}$ is
the unit vector pointing from site $i$
to its NN $j$.
In this model, $\lambda\bso$
competes with $\lambda\bR$ and $\epsilon$, in the sense that
$\lambda\bso$ tends to drive the system to the QSH phase
while $\lambda\bR$ and $\epsilon$ tend to retain the trivial band topology.

For simplicity, we first drop the Rashba coupling, so that
spin is a good quantum number.  The system
is in the QSH phase when $3\sqrt{3}\lambda\bso\!>\!\epsilon$,
and in the normal phase otherwise.  Without the Rashba term,
the Kane-Mele model can be considered as a superposition of two
copies of the Haldane model with opposite Chern numbers.\cite{Haldane-model}
If one calculates the 2D Chern numbers for spin-up and spin-down electrons
separately, one would find that the two Chern
numbers are $\pm1$ in the QSH phase.  While the Haldane-model system
goes from a normal insulator to a CI via a band inversion at
either the K or K$'$ point, the Kane-Mele model transitions to the
QSH state via simultaneous band inversions at both K and K$'$, 
but for opposite spins at these two points.

%
\begin{figure}
\centering
\subfigure{
\includegraphics[width=7.0cm,bb=18 251 559 552]{KM_spinspillage.eps}}
\subfigure{
\includegraphics[width=7.0cm,bb=44 256 536 544]{KM_spillage.eps}}
\caption{Spin-orbit spillage of the Kane-Mele model in the QSH phase,
with $t\!=\!1$, $\lambda\bso\!=\!0.1t$, and $\epsilon\!=\!0.1t$.
(a) Spin-resolved spillage without Rashba coupling; solid (green) and
dashed (red) lines denote spin-up and spin-down spillage.
Inset shows
$\Gamma$-M-K-M-K$^\prime$-$\Gamma$ path used here (blue) and
K-$\Gamma$-M-K path used in Fig.~\ref{fig:bi}
(magenta).
(b) Total spillage without (solid line) and with
(dashed line) Rashba coupling.}
\label{fig:km}
\end{figure}
%

The SOC-induced spillage without the Rashba term is
shown in Fig.~\ref{fig:km}(a).  In this case the spins act
independently, so the spin-up and spin-down spillages
%
$\gamma_\sigma(\mathbf{k})= n\bocc/2-\sum_{m,n=1}^{n\bocc/2}
\vert M_{n\sigma,m\sigma}(\mathbf{k})\vert^2$
(where $\sigma\!=\!\{\uparrow,\downarrow\}$)
%
are shown separately.  Clearly the spin-up band inversion at K is
responsible for $\gamma_{\uparrow}\!=\!1$, and conversely at
K$'$.  The total spillage
%
$\gamma(\mathbf{k})=\gamma_{\uparrow}(\mathbf{k})+
\gamma_{\downarrow}(\mathbf{k})$
%
is shown by the solid line in Fig.~\ref{fig:km}(b).  The symmetry between
the behavior at K and K$'$ has been restored by summing over spins.
Note that the peak values are $\gamma=\!1$ exactly; the fact that they
do not exceed one is an artifact of the simplicity of the model.

When the Rashba coupling is included, as shown by the dashed line
in Fig.~\ref{fig:km}(b), spin is no longer a good quantum number,
so that a spin decomposition is not well-defined.
As expected, adding the Rashba term does not significantly change the
results;\cite{comment_rashba}
one still finds that the spillage reaches unity at
K and K$'$ as before, providing an indication of the spin-Hall phase.

%-------------------------------------------------
\subsection{Application to Chern insulators}
\label{sec:ci}
%-------------------------------------------------

We now consider the case of broken TR symmetry, so that the
\Z2\ index is no longer well-defined, but the possibility of
CI phases appears.  As discussed in Sec.~\ref{sec:intro}, 
SOC is important here as well.  Here we study a buckled honeycomb
Bi bilayer with a Zeeman field applied normal to the plane, 
which can be regarded as having been cut from a
3D Bi crystal on a $(111)$ plane.
The Bi $(111)$ bilayer has been proposed as 
a candidate for QSH insulator.\cite{murakami-prl06}
If a Zeeman field is further applied,
it is possible to obtain CI phases with Chern numbers
$C\!=\!1$ or $C\!=\!-2$.\cite{hongbin-prb12,hongbin-prb13}
To describe this system we use a TB model based on Bi 6$s$ and
6$p$ orbitals,
where the first-neighbor $ss$,
$sp$, $pp\sigma$, and $pp\pi$ hoppings, as well as
the second-neighbor $pp\sigma$ hoppings, are included.
The hopping parameters
are taken from a TB model for 3D bulk Bi.\cite{bi-tb}
In order to obtain non-zero Chern numbers,
an on-site $p$-shell SOC ($\lambda_{\rm SOC}$)
and a Zeeman field ($H_z$) are further applied.
It turns out that if $H_z$ is fixed at 0.8\,eV,
then the phases with $C\!=\!-2$ and $+1$ are realized at
$\lambda_{\rm SOC}\!=\!2.4$\,eV and 0.6\,eV respectively.
If the SOC is completely turned off, $C\!=\!0$.

\begin{figure}
\centering
\subfigure{\includegraphics[height=4.0cm,bb=50 256 539 540]{bi_spillage.eps}}
\subfigure{\includegraphics[height=4.2cm,bb=1 233 625 566]{Bi_spillage_2Dcart.eps}
}
\caption{(a) Spin-orbit spillage of the Bi bilayer
for $C\!=\!1$ (dashed blue) and
$C\!=\!-2$ (solid red) phases, plotted along the K-$\Gamma$-M-K path
(magenta path in inset of Fig.~\ref{fig:km}(a)).
(b) Spillage for $C\!=\!1$ phase plotted in the 2D BZ
($k_x$ and $k_y$ in units of \angstrom$^{-1}$).
(c) Same for $C\!=\!-2$ phase.}
\label{fig:bi}
\end{figure}
%

The spillage for the Bi bilayer is shown along a high-symmetry
$k$-path in Fig.~\ref{fig:bi}(a), and as a distribution in the 2D BZ in
Figs.~\ref{fig:bi}(b-c), for the two parameter sets giving the
$C\!=\!1$ and $C\!=\!-2$ phases.
In both cases the spillage distribution is concentrated at $\Gamma$,
indicating a band inversion there, although it is much more
sharply peaked in the $C\!=\!1$ case.
Clearly the spillages provide a signature of the presence of
a Chern-insulator phase, including the location of the band inversion and
the magnitude (but not the sign) of the Chern number.  Here again
the peak values of the spillage are exactly equal to the magnitude
of the Chern number.  For more realistic systems with more bands 
included, the spillage can be expected to exceed these values slightly,
but a clear correlation between the peak values of spillage and the
Chern number is still expected.


%-------------------------------------------------
\subsection{Application to 3D topological insulators}
\label{sec:materials}
%-------------------------------------------------

%
\begin{figure}
\centering
\includegraphics[width=6cm]{lattice.eps}
\caption{(a) Lattice structure of \bise.
(b) The BZ of \bise; the spillage and bandstructures shown in
Fig.~\ref{fig:materials}(a) and Fig.~\ref{fig:proj} are plotted along the 
black path. }
\label{fig:lattice}
\end{figure}
%
In this subsection we apply our formalism to realistic first-principles
calculations of \bise, \inse\ and \sbse. \bise\ is a
well-known strong TI,\cite{zhang-np09} where the SOC-induced band inversion
takes place at $\Gamma$. We also consider \inse\ and \sbse\ in the same
crystal structure (known as $\beta$ phase for \inse\ and not realized
experimentally for \sbse), which are theoretically predicted
(and experimentally confirmed for \inse) to be trivial insulators.
\cite{zhang-np09,oh-prl12,armitage1,jpl-prb13}
Here it is interesting to note that even though Sb and
In have very similar atomic SOC strength, the substitution of In atoms
tends to drive \bise into a trivial-insulator phase
much faster than does Sb substitution, due to the existence of
In $5s$ orbitals.\cite{jpl-prb13}

As shown in Fig.~\ref{fig:lattice},
the considered structure is rhombohedral, with
two cations and three Se atoms in the primitive unit cell.
The five 2D monolayers are stacked in an
\hbox{$A$-$B$-$C$-$A$-...} sequence
along the (111) direction to form quintuple layers (QLs).  Experimentally
the in-plane hexagonal parameters are $a\!=\!4.138$ and 4.05\,\AA,
and the QL size is $c\!=\!9.547$ and 9.803\,\AA, for \bise\ and
\inse\ respectively.   In our calculations, we take the experimental 
lattice parameters of \bise\ and \inse, but relax their internal atomic
coordinates.  As for \sbse, because its rhombohedral structure 
is not adopted in nature, both the lattice parameters and atomic positions
are relaxed.
The ground state of rhombohedral \sbse\ is predicted to be a
trivial insulator with $a\!=\!4.11$\,\AA\ and $c\!=\!10.43$\,\AA.

We use the \textsc{Quantum ESPRESSO} package\cite{QE-2009} to carry out
first-principles calculations on these systems both with and without SOC.
The PBE generalized-gradient
approximation (GGA) is taken to treat the exchange-correlation functional,
\cite{pbe-1,pbe-2} and norm-conserving
pseudopotentials are generated from OPIUM package.\cite{opium-web,opium-paper}
The energy cutoff is taken as 65\,Ry for \inse\ and 55\,Ry for
\bise\ and \sbse, with an
8$\times$8$\times$8 Monkhorst-Pack $\bf k$ mesh.\cite{monkhorst-pack}
The wavefunctions defined in the plane-wave basis are extracted from these
calculations and Eq.~(~\ref{equa:plane-wave}) is applied to evaluate
the spillage.


As mentioned in Sec.~\ref{sec:intro}, the spillage can also be calculated
in the Wannier basis. Starting from the first-principles calculations,
we use the \textsc{Wannier90} package\cite{wannier90} to
construct Wannier functions (WFs) and a
corresponding realistic TB model\cite{comment_wannier90} for
each of the three materials.
The basis WFs are constructed by projecting 30 atomic $p$ trial
orbitals onto the Bloch subspace of $p$-like bands to generate a
30-band spinor model for \bise\ and \sbse, 
whereas four additional In $5s$ projectors and bands are
included in the model for \inse.
In order that they will retain their atomic-like identity as much
as possible, the projected WFs are not optimized to minimize the
spread functional.\cite{MLWF-rmp} We find that the
WFs generated by this projection method are almost the same for the 
systems with and without SOC,
so that the matrix elements $M_{mn}(\kk)$ defined in Sec.~\ref{sec:def}
can be evaluated with good accuracy.
%\tnewpage
\begin{figure}
\centering
\subfigure{
\includegraphics[height=4.7cm,bb=34 236 551 556]
{spillage_3.eps}}
\subfigure{\includegraphics[width=7.0cm,bb=86 246 515 545]{bi2se3_spillage2D.eps}
}
\caption{(a) Spin-orbit spillage of rhombohedral 
\bise, \sbse\ and \inse\ as indicated by blue, green and red lines
respectively.  Solid (dashed) lines show the spillage computed
from direct plane-wave (Wannier-based) calculations.
(b) Spillage of \bise\ in the
$(k_x,k_y)$ plane at $k_z\!=\!0$
(units of \angstrom$^{-1}$).}
\label{fig:materials}
\end{figure}


The spillage from the direct plane-wave calculations are shown as
the solid lines in Fig.~\ref{fig:materials}(a).
For \bise, the spillage $\gamma(\mathbf{k})$ has a peak value of 2.12
at $\Gamma$, which is slightly larger than 2, indicating that two
Kramers degenerate bands at $\Gamma$ have been inverted by SOC.
On the other hand, the effect of SOC in \inse\ and \sbse\ seems to be
negligible everywhere in the BZ, which is consistent with the fact
that they are both trivial insulators.

\begin{figure}
\centering
\includegraphics[width=8.5cm,bb=33 252 579 555]{proj_band.eps}
\caption{
(a) Wannier-interpolated bandstructure of \sbse.
(b) Same for \bise. Color coding indicates weight of
Sb $5p$ or Bi $6p$ orbitals. }
\label{fig:proj}
\end{figure}

The calculations carried out in the Wannier basis are shown by
the dashed lines in Fig.~\ref{fig:materials}(a).  The spillage
is typically slightly larger for the direct plane-wave
calculations, since the fact that the WFs have a slightly
different plane-wave representation with and without SOC is not
taken into account in the Wannier-based calculations.
Still, the qualitative features are the same,
showing that the Wannier-based approach can successfully provide the
same kind of information about the nature and location of the
topological band inversion. In
Fig.~\ref{fig:materials}(b) we also show the spillage
of \bise\ in the ($k_x, k_y$) plane at $k_z\!=\!0$, as calculated
in the Wannier basis,
which again indicates a highly localized band inversion near
$\Gamma$ and is fully consistent with the expected picture of
the band inversion in \bise.

To see the band inversion from another perspective, we plot
in Fig.~\ref{fig:proj}
the bulk bandstructures of \sbse\ and \bise\ projected onto
Sb $5p$ and Bi $6p$ orbitals respectively.  It is clear that
for \sbse, the Sb $5p$ orbitals are almost all concentrated in the
conduction bands, whereas in \bise\ there is a localized region
around $\Gamma$ where the corresponding Bi $6p$ orbitals
contribute mostly to the top valence band.  This is precisely
the region of the band inversion
corresponding to the peak at $\Gamma$ in Fig.~\ref{fig:materials}.

%==============================================================
\section{Summary}
%==============================================================

To summarize, we have introduced the SOC-induced spillage $\gamma(\kk)$
as a useful quantitative tool for evaluating the degree of band inversion
driven by SOC and mapping it as a function of $\kk$ in the BZ.
We have applied this method to the 2D Kane-Mele model and a
tight-binding model of a Bi bilayer with applied Zeeman field,
as well as to realistic materials including
both trivial and topological insulators. A clear correspondence
between non-trivial topological indices and non-trivial spillage
distributions is evident. In the Kane-Mele model,
one gets two peaks of spillage at K and K$'$ with the 
peak value of 1, which indicates that a single band is inverted
at these two points corresponding to an odd 2D \Z2\ index. 
In the Bi bilayer with applied Zeeman field,
a peak of spillage shows up at
$\Gamma$, with the peak value corresponding to the absolute value
of the Chern number. In \bise, the spillage is slightly
greater than 2 at one of the TRIM, namely $\Gamma$,
implying that two bands are inverted by SOC there and signaling
the presence of a nontrivial strong \Z2\ index.

As mentioned above, other methods exist for the direct computation
of topological Chern and \Z2\ indices, with or without inversion
symmetry,\cite{fukui-jpsj05,fu-prb07,soluyanov-prb11}
and we still recommend these if a direct and definitive determination
of the topological indices is needed.  However, the present spillage-based
approach has the advantage of providing a BZ map of the strength,
position, and degree of localization of the band inversion responsible
for the topological character, thus giving valuable physical intuition
about the origin of the topological properties of the material in question.
In addition, compared with
direct methods for topological index calculation,  the spillage
calculation only requires the evaluation of overlaps between two
wavefunctions at the same $\kk$ point, which is easy to implement
and numerically very efficient. Therefore, it is our hope that the
calculation of SOC spillage will prove to be a widely useful tool
that can be applied both for high-throughput screening for topological
materials and for obtaining a deeper understanding of the critical
features of the bandstructures in known topological materials.

\acknowledgments

This work was supported by NSF Grant DMR-10-05838.
We appreciate valuable discussions with Hongbin Zhang and Huaqing Huang.

\bibliography{spillage}

\end{document}
