% \section{ELECTRIC POLARIZATION AND ORBITAL MAGNETIZATION}
% \label{sec:pol}

First-principles calculations of electric dipoles and orbital magnetic
moments of molecular systems are straightforward.  The electric dipole
is
%
\beq
{\bf d}=-e\sum_j \matel{\psi_j}{\rr}{\psi_j}
\label{eq:finite-p}
\eeq
%
and the orbital moment is
%
\beq
{\bf m}= -\frac{e}{2c}\sum_j\matel{\psi_j}{\rr\times\vv}{\psi_j} \,,
\label{eq:finite-m}
\eeq
%
where the sum is over occupied Hamiltonain eigenstates $\ket{\psi_j}$,
$\rr$ is the position operator, $\vv=(i/\hbar)[H,\rr]$ is the
velocity operator, and Gaussian units are used.
However, these formulas cannot easily be generalized to the case of
crystalline systems, because the Hamiltonian
eigenstates take the form of Bloch functions $\ket{\psi\bnk}$ that
extend over all space.  The problem is that matrix
elements such as $\matel{\psi\bnk}{\rr}{\psi\bnk}$ and
$\matel{\psi\bnk}{\rr\times\vv}{\psi\bnk}$ are ill-defined for such
extended states \cite{nenciu-rmp91}.

To deal with this problem, the so-called ``modern theory of polarization''
\cite{resta-f92,king-smith-prb93,vanderbilt-prb93,resta-rmp94}
was developed in the 1990's, and a corresponding ``modern theory of
magnetization'' in the 2000's
\cite{xiao-prl05,thonhauser-prl05,ceresoli-prb06,shi-prl07,souza-prb08}.
Useful reviews of these topics have appeared \cite{resta-jpcm00,
vanderbilt-cccms06,resta-pf07,resta-jpcm10}.

These theories can be formulated either in terms of Berry phases and
curvatures, or equivalently, by working in the Wannier representation.
The basic idea of the latter is to consider a large but finite sample
surrounded by vacuum and carry out a unitary transformation from the
set of delocalized Hamiltonian eigenstates $\psi_j$ to a set of Wannier-like
localized molecular orbitals $\phi_j$.  Then one can use
\equ{finite-p} or \equ{finite-m}, with the $\psi_j$ replaced
by the $\phi_j$, to evaluate the electric or orbital magnetic
dipole moment per unit volume in the thermodynamic limit.  In doing
so, care must be taken to consider whether any surface contributions
survive in this limit.

In this section, we briefly review the modern theories of electric
polarization and orbital magnetization and related topics.
The results given in this section are 
valid for any set of localized WFs; maximally localized
ones do not play any special role.  Nevertheless, the
close connection to the theory of polarization has been one of
the major factors behind the resurgence of interest in WFs.
Furthermore, we shall see that the use of
MLWFs can provide a very useful,
if heuristic, local decomposition of polar properties in a an
extended system.  For these reasons, it is appropriate to
review the subject here.

%===========
\subsection{Wannier functions, electric polarization, and localization}
\label{sec:pol-theo}
%===========

%------------
\subsubsection{Relation to Berry-phase theory of polarization}
\label{sec:bp}
%------------

Here we briefly review the connection between the Wannier
representation and the Berry-phase theory of polarization
\cite{king-smith-prb93,vanderbilt-prb93,resta-rmp94}.
Suppose that we have constructed via \equ{blochlike} a set of
Bloch-like functions $\ket{\tilde{\psi}\bnk}$ that are smooth
functions of $\kk$.  Inserting these in place of $\ket{\psi\bnk}$
on the right side of \equ{wanniertransform}, the WFs in
the home unit cell $\RR$=$\0$ are simply
%
\beq
\ket{\0 n}=\ibz d\kk\; \ket{\tilde{\psi}\bnk} \,.
\label{eq:homewf}
\eeq
%
To find their centers of charge, we note that
%
\beq
\rr\,\ket{\0 n}=\ibz d\kk\; (-i\nabla_\kk e^{i\kk\cdot\rr}) \,\ket{\tilde{u}\bnk}\,.
\label{eq:rhomewf}
\eeq
%
Performing an integration by parts and applying $\bra{\0 n}$ on the left,
the center of charge is given by
%
\beq
\rr_n = \matel{\0 n}{\rr}{\0 n}=
  \ibz d\kk\; \matel{\tilde{u}\bnk}{i\nabla_\kk}{\tilde{u}\bnk} \,,
\label{eq:wfc}
\eeq
%
which is a special case of \equ{rmatel}.  Then, in the home unit cell, in
addition to the ionic charges $+eZ_\tau$ located at positions $\rr_\tau$,
we can imagine electronic charges $-e$ located at positions
$\rr_n$.\footnote{In these formulas, the sum over $n$ includes a
   sum over spin.  Alternatively a factor of 2 can be inserted to
   account explicitly for spin.}
Taking the dipole moment of this imaginary cell and dividing by the
cell volume, we obtain, heuristically 
%
\beq
\PP= \frac{e}{V} \Big( \sum_\tau Z_\tau \rr_\tau - \sum_n \rr_n \Big)
\label{eq:pol-real}
\eeq
%
for the polarization.

This argument can be put on somewhat firmer ground by imagining a
large but finite crystallite cut from the insulator of interest,
surrounded by vacuum.  The crystallite is divided into an ``interior''
bulk-like region and a ``skin'' whose volume fraction vanishes
in the thermodynamic limit.  The dipole moment is computed from
\equ{finite-p}, but using LMOs $\phi_j$ in place of the
Hamiltonian eigenfunctions $\psi_j$ on the right-hand side.
Arguing that the contribution of the skin to $\bf d$ is negligible
in the thermodynamic limit and that the interior LMOs
become bulk WFs, one can construct an
argument that arrives again at \equ{pol-real}.

If these arguments still seem sketchy, \equ{pol-real} can be
rigorously justified by noting that its second term
%
\beq
\PP_{\rm el}=-\frac{e}{(2\pi)^3}\sum_n \int_{\rm BZ} d\kk\;
   \matel{\tilde{u}\bnk}{i\nabla_\kk}{\tilde{u}\bnk} \,,
\label{eq:Pel}
\eeq
%
is precisely the expression for the electronic
contribution to the polarization in the Berry-phase theory
\cite{king-smith-prb93,vanderbilt-prb93,resta-rmp94},
which was derived by considering the flow of
charge during an arbitrary adiabatic change of the crystalline
Hamiltonian.

The Berry-phase theory can be regarded
as providing a mapping of the distributed quantum-mechanical
electronic charge density onto a lattice of negative point charges
of charge $-e$, as illustrated in Fig.~\ref{fig:mapping}.
Then, the change of polarization
resulting from any physical change, such as the displacement of one
atomic sublattice or the application of an electric field, can
be related in a simple way to the displacements of the Wannier centers
$\rr_n$ occurring as a result of this change.

\begin{figure}
\begin{center}
\includegraphics[width=7cm]{figs/rho}
\caption{Illustration of mapping from physical crystal
  onto equivalent point-charge system with correct dipolar
  properties.  (a) True system composed of point ions (+) and
  charge cloud (contours).  (b) Mapped system in which charge
  cloud is replaced by quantized electronic charges ($-$).
  In the illustrated model there are two occupied bands, i.e.,
  two Wannier functions per cell.}
\label{fig:mapping}
\end{center}
\end{figure}


A well-known feature of the Berry-phase theory is that the polarization
is only well-defined modulo a quantum $e\RR/V$, where $\RR$ is a
real-space lattice vector.  Such an indeterminacy is immediately
obvious from \equ{pol-real}, since the choice of which WFs
are assigned to the home unit cell ($\RR$=$\0$) --
or, for that matter, which ions are assigned to it -- is arbitrary.
Shifting one of these objects by a lattice vector $\RR$ merely
changes $\PP$ by the quantum.
Correspondingly, it can be shown that an arbitrary change of
gauge can shift individual Wannier centers $\rr_n$ in arbitrary
ways, except that the sum $\sum_n\rr_n$ is guaranteed to remain
invariant (modulo a lattice vector).

%------------
\subsubsection{Insulators in finite electric field}
\label{sec:efields}
%------------

The theory of crystalline
insulators in finite electric field $\EE$ is a subtle one;
the electric potential $-\EE\cdot\rr$ does not obey the conditions
of Bloch's theorem, and moreover is unbounded from below,
so that there is no well-defined ground state.  In practice one
wants to solve for a long-lived resonance in which the charge
density and other properties of the insulator remain periodic,
corresponding to what is meant experimentally by an insulator
in a finite field.  This is done by searching for
local minima of the electric enthalpy per cell
%
\beq
F=E_{\rm KS}-V\EE\cdot\PP
\label{eq:enthalpy}
\eeq
%
with respect to both the electronic and the ionic degrees of freedom.
$E_{\rm KS}$ is the ordinary Kohn-Sham energy as it would
be calculated at zero field, and the second term is the coupling of
the field to the polarization
as given in \equ{pol-real} \cite{nunes-prl94} or via the
equivalent Berry-phase expression \cite{nunes-prb01,souza-prl02,
umari-prl02}.  This approach is now standard for
treating periodic insulators in finite electric fields in 
density-functional theory.

%------------
\subsubsection{Wannier spread and localization in insulators}
\label{sec:wanloc}
%------------

We touch briefly here on another interesting connection to the
theory of polarization.  Resta and coworkers have defined a measure
of localization \cite{resta-prl99,sgiarovello-prb01,resta-jpcm02,
resta-jcp06} that distinguishes insulators from metals quite
generally, and have shown that this localization measure reduces,
in the absence of two-particle interactions or disorder, to the
invariant part of the spread functional $\Omega_{\rm I}$ given
in \equ{OmegaI}.  Moreover, \textcite{souza-prb00b}
have shown that this same
quantity characterizes the root-mean-square quantum fluctuations
of the macroscopic polarization. 
Thus,
while the Wannier charge centers are related to the {\it mean} value of
$\PP$ under quantum fluctuations, their invariant quadratic spread
$\Omega_{\rm I}$ is related to the corresponding {\it variance} of $\PP$.

%------------
\subsubsection{Many-body generalizations}
\label{sec:pol-manybody}
%------------

In the same spirit as for the many-body WFs discussed at
the end of Sec.~\ref{sec:manybody}, it is possible to generalize the
formulation of electric polarization and electron localization to the
many-body context.  One again considers $N$ electrons in a supercell;
for the present discussion we work in 1D and let the supercell have
size $L$.  The many-body theory of electric polarization was
formulated in this context by \textcite{ortiz-prb94}, and later
reformulated by \textcite{resta-prl98}, who introduced a ``many-body
position operator'' $\hat{X}=\exp(i2\pi\hat{x}/L)$ defined in terms of
the ordinary position operator $\hat{x}=\sum_{i=1}^N \hat{x}_i$.
While $\matel{\Psi}{\hat{x}}{\Psi}$ is ill-defined in the extended
many-body ground state $\ket{\Psi}$, the matrix element
$\matel{\Psi}{\hat{X}}{\Psi}$ is well-defined and can be used to
obtain the electric polarization, up to the usual quantum.  These
considerations were extended to the localization functional, and the
relation between localization and polarization fluctuations, by
\textcite{souza-prb00b}.  The variation of the many-body localization
length near an insulator-to-metal transition in 1D and 2D model
systems was studied using quantum Monte Carlo methods by
\textcite{hine-jpcm07}. Finally, the concept of electric enthalpy
was generalized to the many-body case by \textcite{umari_prl05},
allowing to calculate for the first time dielectric properties with
quantum Monte Carlo, and applied to the case of the polarizabilities
\cite{umari_prl05} and
hyperpolarizabilities \cite{umari_jcp09} of periodic hydrogen chains.

%===========
\subsection{Local polar properties and dielectric response}
\label{sec:pol-diel}
%===========

Is Sec.\ref{sec:bp} we emphasized the equivalence of the $\kk$-space
Berry-phase expression for the electric polarization, \equ{Pel},
and the expression written in terms of the locations of the
Wannier centers $\rr_n$, \equ{pol-real}.  The latter
has the advantage of being a real-space expression, thereby
opening up opportunities for localized descriptions and
decompositions of polar properties and dielectric responses.
We emphasize again that MLWFs have no
privileged role in \equ{pol-real}; the expression remains correct
for any WFs that are sufficiently well localized that the
centers $\rr_n$ are well defined.  Nevertheless, one may
argue heuristically that MLWFs provide the most natural
local real-space description of dipolar properties in crystals
and other condensed phases.

%------------
\subsubsection{Polar properties and dynamical charges of crystals}
\label{sec:pbulk}
%------------

Many dielectric properties of crystalline solids are
most easily computed directly in the $\kk$-space Bloch
representation.  Even before it was understood how to compute the
polarization $\bf P$ via the Berry-phase theory of \equ{Pel}, it was
well known how to compute derivatives of $\bf P$ using linear-response
methods \cite{gironcoli-prl89,resta-f92,baroni-rmp01}.  Useful derivatives include the electric
susceptibility $\chi_{ij}=dP_i/d\EE_j$ and the Born (or dynamical)
effective charges $Z^*_{i,\tau j}=VdP_i/dR_{\tau j}$,
where $i$ and $j$ are Cartesian labels and $R_{\tau j}$ is the
displacement of sublattice $\tau$ in direction $j$.  With the
development of the Berry-phase theory,
it also became possible to compute
effective charges by finite differences.  Similarly, with the
electric-enthalpy approach of \equ{enthalpy} it became possible to
compute electric susceptibilities by finite differences as well
\cite{souza-prl02,umari-prl02}.

The Wannier representation provides an alternative method for
computing such dielectric quantities by finite differences.  One computes
the derivatives $dr_{n,i}/d\EE_j$ or $dr_{n,i}/dR_{\tau j}$ of the
Wannier centers by finite differences, then sums these to get the desired
$\chi_{ij}$ or $Z^*_{i,\tau j}$.  An example of such a calculation
for $Z^*$ in GaAs was presented already in Sec.~VII of
\textcite{marzari-prb97},
%
and an application of the Wannier approach of \textcite{nunes-prl94}
in the density-functional context was used to compute $\chi$ by finite
differences for Si and GaAs \cite{fernandez-prb98}.
%
Dynamical charges were computed for several TiO$_2$ phases by
\textcite{cangiani-prb04} and \textcite{posternak-prb06}, and, as mentioned in 
Sec.~\ref{sec:bondcrystals},
observed differences between polymorphs were correlated with changes
in the chemical nature of the WFs associated with
OTi$_3$ structural units.
%
Piezoelectric coefficients, which are derivatives of $\bf P$ with
respect to strain, have also been carried out in the Wannier
representation for ZnO and BeO by \textcite{noel-prb02}.

Some of the most extensive applications of this kind have been to
ferroelectric perovskites, for which the dynamical charges have been
computed in density-functional and/or Hartree-Fock contexts for
BaTiO$_3$, KNbO$_3$, SrTiO$_3$, and PbTiO$_3$
\cite{marzari-acp98,baranek-prb01,evarestov-ssc03,usvyat-ijqc04}.  In
these materials, partially covalent bonding associated with
hybridization between O $2p$ and Ti $3d$ states plays a crucial role
in stabilizing the ferroelectric state and generating anomalous
dynamical effective charges \cite{posternak-prb94,zhong-prl94}.
Recall that the dynamical, or Born, effective charge $Z^*$
is defined as $Z^*_{i,\tau j}=VdP_i/dR_{\tau j}$
and carries units of charge.  Naively, one might expect values around
$+4e$ for Ti ions and $-2e$ for oxygen ions in BaTiO$_3$ based on
nominal oxidation states, but instead one finds ``anomalous'' values
that are much larger.  For example, \textcite{zhong-prl94} reported
values of $+7.2e$ for Ti displacements, and $-5.7e$ for O
displacements along the Ti--O--Ti chains.

\begin{figure}
\begin{center}
\includegraphics[width=2.5in]{figs/wannier-perov}
\caption{(Color online) Amplitude isosurface plots of the maximally-localized
O$[2p_z]$--Ti$[3d_{z^2}]$ Wannier functions in
BaTiO$_3$. O is at center, surrounded by a plaquette of four
Ba atoms (green); Ti atoms (yellow, almost hidden) are above and below.
(a) Centrosymmetric structure.  (b) Ferroelectric structure
in which central O is displace upward relative to neighboring
Ti atoms.}
\label{fig:bto}
\end{center}
\end{figure}

This behavior can be
understood \cite{posternak-prb94,zhong-prl94} as arising from
hybridization between neighboring O $2p$ and Ti $3d$ orbitals that
dominate the valence and conduction bands, respectively.  This
hybridization, and the manner in which it contributes to an anomalous
$Z^*$, can be visualized by inspecting the changes in the MLWFs induced
by the atomic displacements.  Fig.~\ref{fig:bto}(a)
shows an 
O$[2p_z]$--Ti$[3d_{z^2}]$ MLWF in centrosymmetric
BaTiO$_3$ \cite{marzari-acp98}.
The hybridization to Ti $3d_{z^2}$ states appears in the form of
the ``donuts'' surrounding the
neighboring Ti atoms.  When the O atom moves upward relative to the
geometric center of the two neighboring Ti atoms as
shown in Fig.~\ref{fig:bto}(b), as it does in
ferroelectrically distorted BaTiO$_3$, the hybridization strengthens
for the upper O--Ti bond and weakens for the lower one, endowing the
WF with more Ti $3d$ character on the top than on the bottom.  As a
result, the center of charge of the WF shifts upward, and since
electrons carry negative charge, this results in a negative anomalous
contribution to the $Z^*$ of the oxygen atom.  The figure illustrates
this process for $\sigma$-oriented oxygen WFs, but a similar effect
occurs for the $\pi$-oriented oxygen WFs, and the total anomalous
dynamical charge can be accounted for quantitatively on the basis of
the distortion-induced changes of each kind of WF in the crystal
\cite{marzari-acp98}.

The above illustrates the utility of the MLWFs in providing a
{\it local} description of dielectric and polar responses in crystals.
This strategy can be carried further in many ways.  For example,
it is possible to decompose the $Z^*$ value for a given atom in
a crystal into contributions coming from various different
neighboring WFs, as was done for GaAs in Sec.~VII of
\textcite{marzari-prb97} and for BaTiO$_3$ by \textcite{marzari-acp98}.
Some chemical intuition is already gained by carrying out a
band-by-band decomposition of the $Z^*$ contributions
\cite{ghosez-prb95,ghosez-jpcm00}, but the WF analysis allows
a further spatial decomposition into individual WF contributions
within a band.  A deeper analysis that also involves the
decomposition of the WFs into atomic orbitals
has been shown to provide further insight into the anomalous
$Z^*$ values in perovskites \cite{bhattacharjee-pccp10}.

\begin{figure}
\begin{center}
\includegraphics[width=3.2in]{figs/nakhmanson}
\caption{(Color online) MLWFs for $\beta$-PVDF polymer chain.
(a) MLWF charge centers, indicated by diamonds.
(b)-(d) MLWFs localized on C--C, C--F, and C--H bonds,
respectively.  From \textcite{nakhmanson-prb05}.}
\label{fig:pvdf}
\end{center}
\end{figure}

Some insightful studies of the polar properties of polymer systems
in terms of MLWFs have also been carried out.  Figure~\ref{fig:pvdf},
for example, shows the WF centers and characters for the $\beta$
conformation of polyvinylidene fluoride ($\beta$-PVDF)
\cite{nakhmanson-prb05}, one of the more promising ferroelectric 
polymer systems.  An inspection of WF centers has also been
invoked to explain the polar properties of so-call ``push-pull''
polymers by \textcite{kudin-jcp07b} and of H$_2$O ice by
\textcite{lu-prl08}.  Finally, we note an interesting recent study
in which changes in polarization induced by corrugations in BN sheets
were analyzed in terms of WFs \cite{naumov-prl09}.

%------------
\subsubsection{Local dielectric response in layered systems}
\label{sec:layer}
%------------

In a similar way, the theoretical study of dielectric properties of
ultrathin films and superlattices can also be enriched by a local
analysis.  Two approaches have been introduced in the literature.
In one, the local $x$-$y$-averaged electric field $\bar{\E}_z(z)$
is calculated along the stacking direction $z$, and then the local
dielectric permittivity profile $\varepsilon(z)=\bar{\E}_z(z)/\bar{D}_z$ or
inverse permittivity profile $\varepsilon^{-1}(z)=\bar{D}_z/\bar{\E}_z(z)$
is plotted, where $\bar{D}_z$ is the $x$-$y$-averaged electric displacement
field (constant along $z$ in the absence of free charge) determined
via a Berry-phase calculation of $P_z$ or by inspection of $\bar{\E}_z$
in a vacuum region.  Such an approach has been applied to study
dielectric materials such as SiO$_2$ and HfO$_2$ interfaced to Si
\cite{giustino-prl03,giustino-prb05,shi-prb06,shi-apl07}
and perovskite films and superlattices \cite{stengel-n06,stengel-prb06}.

\begin{figure}
\begin{center}
\includegraphics[width=2.6in]{figs/wu-prl-2006}
\caption{(Color online) Dispersion of WF center positions along $z$ as a function
of $(k_x,k_y)$ for a superlattice composed of alternating layers of
SrTiO$_3$ (sublayers 1A and 1B) and BaTiO$_3$ (sublayers 2A
and 2B).  From \textcite{wu-prl06}.}
\label{fig:wcband}
\end{center}
\end{figure}

The second approach is to use a Wannier analysis to assign a
dipole moment to each layer.  This approach, based on the concept
of hybrid WFs discussed in Sec.~\ref{sec:hybrid}, was pioneered by
\textcite{giustino-prl03,giustino-prb05} and used by them to study
Si/SiO$_2$ interfaces and related systems.  Later applications
to perovskite oxide films and superlattices have been fairly
extensive.  The essential observation is that, when studying a
system that is layered or stacked along the $z$ direction, one
can still work with Bloch functions labeled by $k_x$ and $k_y$ while
carrying out a Wannier construction or analysis only along $z$.  Since
the extraction of Wannier centers in 1D is rather trivial, even in the
multiband case \cite{marzari-prb97,sgiarovello-prb01,bhattacharjee-prb05}, it is not
difficult to construct a ``Wannier center band structure''
$\bar{z}(k_x,k_y)$, and use the planar-averaged values to assign
dipole moments to layers.  This approach was demonstrated by
\textcite{wu-prl06}, as shown in Fig.~\ref{fig:wcband}, and has since
been used to study perovskite superlattices and artificial
nanostructures \cite{wu-prl08,murray-prb09}.  In a related
development, \textcite{stengel-prb06} introduced a Wannier-based
method for computing polarizations along $z$ and studying electric
fields along $z$ that work even in the case that the stacking includes
metallic layers, as long as the system is effectively insulating along
$z$ \cite{stengel-prb07}.  This allows for first-principles
calculations of the nonlinear dielectric response of ultrathin
capacitors and related structures under finite bias, providing an
insightful avenue to the study of finite-size and dead-layer effects
in such systems \cite{stengel-nm09,stengel-np09,stengel-prb09}.

%------------
\subsubsection{Condensed molecular phases and solvation}
\label{sec:mol}
%------------

Wannier-function methods have also played a prominent role in the
analysis of polar and dielectric properties of dipolar liquids, mainly
H$_2$O and other H-bonded liquids.  While the dipole moment of an
isolated H$_2$O molecule is obviously well defined, a corresponding
definition is not easy to arrive at in the liquid phase where
molecules are in close contact with each other.  An influential
development was the proposal made by
\textcite{silvestrelli-prl99, silvestrelli-jcp99}
that the dipole moment of a solvated molecule could be defined
in terms of positive charges on ionic cores and negative charges
located at the centers of the MLWFs.  Using this definition,
these authors found that the water molecule dipole is somewhat
broadly distributed, but has an average value of about 3.0\,D,
about 60\% higher than in the gas phase.
These features were shown to be in conflict with the behavior
of widely-used empirical models.

Of course, such a definition in terms of the dipole of the
molecular WF-center configuration remains heuristic
at some level.  For example, this local measure of the dipole
does not appear to be experimentally measurable
even in principle, and clearly the use of one of the alternative
measures of maximal localization discussed in Sec.~\ref{sec:alt-loc}
would give rise to slightly different values.  Nevertheless,
the approach has been widely adopted.  For example,
subsequent work elaborated \cite{sagui-jcp04,dyer-jcp06} and
extended this analysis to water in
supercritical conditions \cite{boero-prl00,boero-potps00},
confined geometries \cite{coudert-c06,dellago-cpc05},
and with solvated ions present \cite{scipioni-jcp09}, and compared
the results obtained with different exchange-correlation
functionals~\cite{todorova-jpcb06}.

It should be noted that the decomposition into Wannier dipoles
is closer to the decomposition of the charge density into static
(Szigeti) charges than to a decomposition into dynamical (Born)
charges.  The first one corresponds to a spatial decomposition of
the total electronic charge density, while the second is connected
with the force that appears on an atom in response to an applied
electric field.  As a counterpoint to the WF-based definition,
therefore, \textcite{pasquarello-prb03} have argued the a
definition based on these forces provides a more fundamential
basis for inspecting molecular dipoles in liquids.  In particular,
they define a second-rank tensor as the derivative of the
torque on the molecule with respect to an applied electric field,
and finding that this is typically dominated by its antisymmetric part,
they identify the latter (rewritten as an axial vector) as the
molecular dipole.  Surprisingly, they find that the magnitude of
this dipole vector is typically only about 2.1\,D in liquid water,
much less than the value obtained from the WF analysis.

Clearly the WF-based and force-based approaches to defining
molecular dipoles provide complementary perspectives, and
a more complete reconciliation of these viewpoints is the subject
of ongoing work.

Finally, we note that there is an extensive literature in which
Car-Parrinello molecular-dynamics simulations are carried out for H$_2$O
and other liquids, as already mentioned in Sec.~\ref{sec:aimd}
and surveyed in several reviews
\cite{tuckerman-jpcb00,tse-aropc02,tuckerman-jpcm02,kirchner-prsopl07}.
Using such approaches, it is possible to compute the dynamical
dipole-dipole correlations of polar liquids and compare the results
with experimental infrared absorption spectra.  While it is possible
to extract the needed information from the time-time correlation
function of the total polarization ${\bf P}(t)$ of the entire
supercell as calculated using the Berry-phase approach, methods which
follow the time-evolution of local dipoles, as defined via WF-based
methods, provide additional insight and efficiency
\cite{bernasconi-prl98,pasquarello-prb03,sharma-prl05,
  sharma-prl07,mcgrath-mp07,chen-prb08} and can easily be extended to
other kinds of molecular systems
\cite{gaigeot-mp07,gaigeot-jpcb03,gaigeot-jctc05,pagliai-jcp08}.  The
applicability of this kind of approach has benefited greatly from the
development of methods for computing WFs and their centers ``on the
fly'' during Car-Parrinello molecular dynamics simulations
\cite{sharma-ijqc03,iftimie-jcp04,wu-prb09}.

%===========
\subsection{Magnetism and orbital currents}
\label{sec:magnetism}
%===========

%------------
\subsubsection{Magnetic insulators}
\label{sec:mag}
%------------

Just as an analysis in terms of WFs can help clarify
the chemical nature of the occupied states in an ordinary insulator,
they can also help describe the orbital and magnetic ordering in
a magnetic insulator.

In the magnetic case, the maximal localization proceeds in the same
way as outlined in Sec.~\ref{sec:theory}, with trivial extensions
needed to handle the magnetic degrees of freedom.  In the case of
the local (or gradient-corrected) spin-density approximation,
in which spin-up and spin-down electrons are treated independently,
one simply carries out the maximal localization procedure
independently for each manifold.  In the case of a spinor calculation
in the presence of spin-orbit interaction, one instead implements
the formalism of Sec.~\ref{sec:theory} treating all wavefuntions as
spinors.  For example, each matrix element on the right-hand side
of \equ{Mkb} is computed as an inner product between spinors, and
the dimension of the resulting matrix is the number of occupied
spin bands in the insulator.

Several examples of such an analysis have appeared in the literature.
For example, applications to
simple antiferromagnets such as MnO \cite{posternak-prb02},
novel insulating ferromagnets and antiferromagnets \cite{ku-prl02,ku-jssc03},
and complex magnetic ordering in rare-earth manganates
\cite{yamauchi-prb08,picozzi-jpcm08} have proven to be illuminating.

%------------
\subsubsection{Orbital magnetization and NMR}
\label{sec:om}
%------------

In a ferromagnetic (or ferrimagnetic) material, the total
magnetization has two components. One arises from electron spin and is
proportional to the excess population of spin-up over spin-down
electrons; a second corresponds to circulating orbital currents.  The
spin contribution is typically dominant over the orbital one (e.g., by
a factor of 10 or more in simple ferromagnets such as Fe, Ni and Co \cite{ceresoli-prb10a}),
but the orbital component is also of interest, especially in unusual
systems cases in which it can dominate, or in the context of
experimental probes, such as the anomalous Hall conductivity, that
depend on orbital effects.  Note that inclusion of the spin-orbit
interaction is essential for any description of orbital magnetic
effects.

Naively one might imagine computing the orbital magnetization
${\bf M}_{\rm orb}$ as the thermodynamic limit of \equ{finite-m}
per unit volume for a large crystallite.  However, as we discussed
at the beginning of Sec.~\ref{sec:pol}, Bloch matrix elements
of the position operator $\rr$ and the circulation operator
$\rr\times\vv$ are ill-defined.  Therefore, such an approach is
not suitable.  Unlike for the case of electric polarization,
however, there is a simple and fairly accurate approximation
that has long been used to compute ${\bf M}_{\rm orb}$: one
divides space into muffin-tin spheres and interstitial regions,
computes the orbital circulation inside each sphere as a spatial
integral of $\rr\times{\bf J}$, and sums these contributions.
Since most magnetic moments are fairly local, such an approach
is generally expected to be reasonably accurate.

Nevertheless, it is clearly of interest to have available an
exact expression for ${\bf M}_{\rm orb}$ that can be used to
test the approximate muffin-tin
approach and to treat cases in which itinerant
contributions are not small.  The solution to this problem has
been developed recently, leading to a closed-form expression for
${\bf M}_{\rm orb}$ as a bulk Brillouin-zone integeral.
Derivations of this formula via a semiclassical
\cite{xiao-prl05} or long-wave quantum \cite{shi-prl07} approach
are possible, but here we emphasize the derivation carried out
in the Wannier representation
\cite{thonhauser-prl05,ceresoli-prb06,souza-prb08}.
For this purpose, we restrict our interest to insulating
ferromagnets.  For the case of electric polarization, the
solution to the problem of $\rr$ matrix elements was sketched in
Sec.~\ref{sec:bp}, and a heuristic derivation of \equ{pol-real}
was given in the paragraph following that equation.  A similar
analysis was given in the above references for the case of
orbital magnetization, as follows.

Briefly, one again considers a
large but finite crystallite cut from the insulator of interest,
divides it into ``interior'' and ``skin'' regions, and transforms
from extended eigenstates to LMOs $\phi_j$.
For simplicity we consider the case of a two-dimensional
insulator with a single occupied band.  The interior gives a
rather intuitive ``local circulation'' (LC) contribution to the
orbital magnetization of the form
%
\beq
M_{\rm LC} = -\frac{e}{2A_0c}\matel{\0}{\rr\times\vv}{\0}
\label{eq:MLC}
\eeq
%
where $A_0$ is the unit cell area, since in the interior the
LMOs $\phi_j$ are really just bulk WFs.
This time, however, the skin contribution does {\it not} vanish.
The problem is that $\matel{\phi_j}{\vv}{\phi_j}$ is nonzero
for LMOs in the skin region, and the pattern of these velocity
vectors is such as to describe a current circulating around the
boundary of the sample, giving a second ``itinerant circulation''
contribution that can, after some manipulations, be written
in terms of bulk WFs as
%
\beq
M_{\rm IC} = -\frac{e}{2A_0c\hbar}\sum_{\RR}
[R_x\matel{\RR}{y}{\0}-R_y\matel{\RR}{x}{\0}] \,.
\label{eq:MIC}
\eeq
%
When these contributions are converted back to the Bloch
representation and added together, one finally obtains
%
\beq
{\bf M}_{\rm orb}=\frac{e}{2\hbar c} \, {\rm Im}
\int\frac{d^2k}{(2\pi)^2}\,\matel{\partial_\kk u_\kk}
{\times\,(H_\kk+E_\kk)}{\partial_\kk u_\kk},
\label{eq:Mkspace}
\eeq
%
which is the desired $\kk$-space bulk expression for the orbital
magnetization~\cite{thonhauser-prl05}.\footnote{In the case of metals
  \equ{Mkspace} must be modified by adding a $-2\mu$ term inside the
  parenthesis, with $\mu$ the chemical potential \cite{xiao-prl05,
    ceresoli-prb06}.  Furthermore, the integration is now restricted
  to the {\it occupied} portions of the Brillouin zone.}  The
corresponding argument for multiple occupied bands in three dimensions
follows similar lines \cite{ceresoli-prb06,souza-prb08}, and the
resulting formula has recently been implemented in the context of
pseudopotential plane-wave calculations \cite{ceresoli-prb10a}.
Interestingly, it was found that the interstitial contribution
--~defined as the difference between the total orbital magnetization,
\equ{Mkspace}, and the muffin-tin result~-- is not always
negligible. In bcc Fe, for example, it amounts to more than 30\% of
the spontaneous orbital magnetization, and its inclusion improves the
agreement with gyromagnetic measurements.

The ability to compute the orbital magnetization is also of use
in obtaining the magnetic shielding of nuclei. This is responsible for
the chemical shift effect observed in nuclear magnetic resonance (NMR)
experiments. A first principles theory for magnetic shielding in
solids was established by examining the perturbative response to a
periodic magnetic field in the long wavelength limit
\cite{mauri-prl96,pickard-prb01}. An alternative perturbative approach
used a WF representation of the electronic structure
together with a periodic position operator
\cite{sebastiani-jpca01,sebastiani-cpc02,sebastiani-mplb03}. However,
magnetic shieldings can also be computed using a ``converse'' approach
in which one uses Eq.~(\ref{eq:Mkspace}) to compute the orbital
magnetization induced by a fictitious point magnetic dipole on the
nucleus of interest~\cite{thonhauser-jcp09,ceresoli-prb10b}.
The advantage of such approach is that it does not require 
linear-response theory, and so it is amenable to large-scale
calculations or complex exchange-correlation functionals (e.g. including
Hubbard U corrections, or Hartree-Fock exchange), albeit at the cost
of typically one self-consistent iteration for every 
nucleus considered. Such converse approach has then been extended also to the
calculation of the EPR $g$-tensor by \textcite{ceresoli-prb10a}.


%------------
\subsubsection{Berry connection and curvature}
\label{sec:berrystuff}
%------------

Some of the concepts touched on in the previous section can be
expressed in terms of the $\kk$-space Berry connection
%
\beq
\label{eq:berry-connection}
{\bf A}_{n\kk}=\bra{u_{n\kk}}i{\nabla}_\kk \ket{u_{n\kk}}
\eeq
%
and Berry curvature
%
\beq
\label{eq:curv-pot}
\bm{\mathcal{F}}_{n\kk}={\nabla}_\kk\times {\bf A}_{n\kk}
\eeq
%
of band $n$.  In particular, the contribution of this band to the
electric polarization of Eq.~(\ref{eq:Pel}), and to the second term
in the orbital magnetization expression of Eq.~(\ref{eq:Mkspace}), are
proportional to the Brillouin-zone integrals of
${\bf A}_{n\kk}$ and $E_{n\kk}{\bm{\mathcal{F}}}_{n\kk}$,
respectively.  These quantities will also play a role in the
next subsection and in the discussion of the anomalous Hall
conductivity and related issues in Sec.~\ref{sec:curv-ahc}.

%------------
\subsubsection{Topological insulators and orbital magnetoelectric
response}
\label{sec:top}
%------------

There has recently been a blossoming of interest in so-called
topological insulators, i.e., insulators that cannot be adiabatically
connected to ordinary insulators without a gap closure.
\textcite{hasan-rmp10,hasan-arcmp11} provide excellent reviews
of the background, current status of this field, and provide
references into the literature.

One can distinguish two kinds of topological insulators.  First,
insulators having broken time-reversal ($T$) symmetry (e.g.,
insulating ferromagnets and ferrimagnets) can be classified by
an integer ``Chern invariant'' that is proportional to the
Brillouin-zone integral of the Berry curvature ${\bm{\mathcal{F}}}_{n\kk}$
summed over occupied bands $n$. Ordinary insulators are characterized
by a zero value of the invariant. An insulator with a non-zero
value would behave like an integer quantum Hall system, but without
the need for an external magnetic field; such systems are usually
denoted as ``quantum anomalous Hall'' (QAH) insulators.  While no
examples are known to occur in nature,
tight-binding models exhibiting such a
behavior are not hard to construct \cite{haldane-prl88}.  It can
be shown that a Wannier representation is not possible for a QAH
insulator, and \textcite{thonhauser-prb06} have explored the
way in which the usual Wannier construction
breaks down for model systems.

% \dvm{Revised paragraph to clarify meaning of ${\cal Z}_2$.}
%
Second, dependng on how their Bloch functions wrap the Brillouin
zone, nonmagnetic ($T$-invariant) insulators can be sorted into
two classes denoted as ``${\cal Z}_2$-even'' and ``${\cal
Z}_2$-odd'' (after the name ${\cal Z}_2$ of the group \{0,\,1\}
under addition modulo 2).  Most (i.e., ``normal'') insulators are
${\cal Z}_2$-even, but strong spin-orbit effects can lead to the
${\cal Z}_2$-odd state, for which the surface-state
dispersions are topologically required to display characteristic
features that are amenable to experimental verification.
Several materials realizations of ${\cal Z}_2$-odd insulators
have now been confirmed both theoretically and experimentally
\cite{hasan-rmp10,hasan-arcmp11}.

% \dvm{Revised paragraph to clarify meaning of axion.}
%
In a related development, the orbital magnetoelectric coefficient
$\alpha_{ij}=\partial M_{{\rm orb},j}/\partial \E_i$ was found to
contain an isotropic contribution having a topological character
(the ``axion'' contribution, corresponding to an $\EE\cdot\bm{B}$
term in the effective Lagrangian).  This term 
can be written as a Brillouin-zone integral of the Chern-Simons
3-form, defined in terms of multiband generalizations of the Berry
connection ${\bf A}_\kk$ and curvature ${\bm{\mathcal{F}}}_\kk$
introduced in the previous subsection \cite{qi-prb08,essin-prl09}.
The Chern-Simons magnetoelectric coupling has been evaluated from
first-principles with the help of WFs for both
topological and ordinary insulators~\cite{coh-prb11}.

A careful generalization of
Eq.~(\ref{eq:Mkspace}) to the case in which a finite electric field is
present has been carried out by \textcite{malashevich-njp10} in the
Wannier representation using arguments similar to those in
Secs.~\ref{sec:bp} and \ref{sec:om}, and used to derive a complete
expression for the orbital magnetoelectric response,
of which the topological Chern-Simons term is only one contribution
\cite{malashevich-njp10,essin-prb10}.
