\section*{Abstract}
{\bf
For a long time,  band theory of solids has been focused on the energy spectrum, or Hamiltonian eigenvalues. Recently, it was realized that the collection of eigenvectors also contains important physical information. The local geometry of eigenspaces determines the electric polarization, while their global twisting gives rise to the metallic surface states in topological insulators. These phenomena are central topics of the present notes.  The shape of eigenspaces is also responsible for many intriguing physical analogies, which have their roots in the theory of vector bundles. We give an informal introduction to the geometry and topology of vector bundles and describe various physical models from this mathematical perspective.
}

\vspace{10pt}
\setcounter{tocdepth}{2}
\noindent\rule{\textwidth}{1pt}
\tableofcontents\thispagestyle{fancy}
\noindent\rule{\textwidth}{1pt}
\vspace{10pt}

\section*{Preface}
\subsection*{Why topological insulators?}
Topological insulators are unique materials that do not conduct electricity in the bulk, but support metallic states on their boundary \cite{Hasan2010}. These states have unusual dispersion relation and cannot be removed from the surface or realized as bulk states of a stand-alone crystal. There is no order parameter responsible for these properties, and topological insulators fall outside the classification of phases based on symmetries. Instead, they are characterized by topological invariants, which take integral values and remain unchanged during any variations of the Hamiltonian that do not close the bulk gap.

Topological properties of an insulator can be captured even by the simplest quantum model of a crystal, the Bloch theory of non-interacting electrons. However, the nature of topological invariant is strikingly different from that of the other observables available in the framework of band theory. For example, consider a two-dimensional crystal with a single occupied band. Each individual eigenvector $|\psi_k\rangle$ is determined up to a phase factor, which is known to be physically irrelevant. But is it possible to choose the phases of \emph{all} eigenstates in a smooth fashion over the whole momentum space?  Clearly, one can do it in the case of $k$-independent Hamiltonian that describes a ``crystal'' constructed from decoupled atoms or molecules. It turns out that for some Hamiltonians, the eigenspaces determined by $|\psi_k\rangle$ fit together in a twisted way that inevitably leads to the appearance of singularities in the phase field. These are the Hamiltonians of topological insulators, and the number of singularities is given by the topological invariant.

While the properties of topological insulators are essentially quantum mechanical phenomena of condensed matter, the underlying principles are universal and have been successfully applied to various physical systems. Examples include topological acoustics \cite{Ma2019} and photonics \cite{Ozawa2019}, topoelectric circuits \cite{Lee2018} and topological interpretation of geophysical waves \cite{Delplace2020}. The field of topological matter is the subject of active current research, both in the solid state physics and in other contexts. The theory of topological insulators provides a decent starting point for anyone who wants to learn about the subject.

\subsection*{Why vector bundles?}
It is rather impractical to compute the topological invariant of an insulator by counting singularities in the phases of $|\psi_k\rangle$. Instead, one uses the Hamiltonian eigenstates to define a certain effective gauge field  in the momentum space. The gauge transformations correspond to the choice of the phases of eigenstates. The field strength, known as Berry curvature, depends only on the geometry of eigenspaces $|\psi_k\rangle$ and is gauge-invariant. The integral of the Berry curvature over the whole momentum space, divided by $2\pi$, equals the topological invariant. So, the global question about singularities in the phase field can be answered by integrating the local quantity which is independent of the phase choice.

This situation shares some similarities with quantum description of an electron moving in the magnetic field, if we allow for monopole configurations. Imagine a closed surface with magnetic monopoles sitting inside and consider the wave function of an electron, which can move freely around the monopoles.  Recall that position of each monopole is a starting point of the nodal line, or line of zeroes of the electron's wave function. If we restrict the wave function to the surface, there will be a singularity at each point where the surface is pierced by a nodal line. Here, the topological invariant is given by the total magnetic charge contained inside the surface. It determines the number of singularities in the wave function and can be computed by integrating magnetic field strength over the surface.

Moreover, an analogous construction arises in a purely mathematical context. There is a problem of classification of all closed two-dimensional surfaces. The answer is given by the genus of the surface, or the number of ``handles''. For example, the sphere has genus zero, while for the torus the genus is one. The genus is a topological invariant in the sense that it does not change under smooth deformations of the surface. Now consider a tangent vector field on the surface. On the torus, there exists a continuous nowhere-vanishing field. But on the sphere, any tangent field will contain singularities. The number of singularities of a tangent vector field on a surface is determined by its genus. It turns out that this number is given by the integral of the Gaussian curvature over the whole surface, divided by $2\pi$.

What exactly do topological insulators, magnetic monopoles, and two-dimensional surfaces have in common? The answer lies in the mathematical theory of vector bundles. These objects first appeared in physics in the context of gauge fields describing fundamental interactions. Later it was realized that vector bundles play an important role in many other physical situations. In fact, adjectives in ``geometric phase'' and ``topological insulator'' refer to the two ways one can characterize the shape of a vector bundle. Geometry is focused on the local quantities that can be measured in the vicinity of a given point. In physical terms, geometry of a vector bundle manifests itself as an effective gauge field. In particular, this allows one to interpret the Berry curvature of Bloch eigenstates as a ``magnetic field in momentum space''. Topology, on the other hand, is concerned with the global perspective and deals with properties that characterize vector bundle as a whole. In physics, such properties give rise to precisely quantized observables, which stay invariant under the wide range of perturbations.

\subsection*{About these notes}
The goal of these notes is to present the basic theory of topological insulators with emphasis on the mathematical structures that lie behind the physical models. We start with an informal introduction to geometry and topology of vector bundles in Sections \ref{Chap:ConnVB}--\ref{Chap:TopVB}. Our discussion is far from being mathematically rigorous: we will often give definitions by examples and use intuitive arguments instead of proofs. The idea is to learn just enough mathematical language to be able to spot certain patterns in physical models and to build analogies while understanding their origin and limitations.

In the second part of the notes, Sections \ref{Chap:ModPol}--\ref{Chap:Semi}, we consider a number of geometric and topological phenomena in condensed matter physics. These include geometric phase formula for electric polarization, topological invariant of adiabatic charge pump, and bulk-boundary correspondence for Chern insulators. We will also  discuss basic classification of topological matter with symmetries and the role of vector bundles in the theory of  topological semimetals. We will not focus on the foundations of physical models or their limits of applicability --- for example, we simply take tight-binding approximation for granted. Instead, we will consider a ``minimal working example'' of each phenomenon and  highlight its geometric aspects.

The notes are based on the material of one-semester course taught by the author in the Moscow State University. The intended audience includes beginning graduate students and advanced undergraduates. Prerequisites are kept to the bare minimum:  the reader must be familiar with vector calculus in mathematics and Dirac formalism in quantum mechanics. The text includes exercises, many of which are straightforward checks, but some can be more challenging. It is important to solve or at least attempt all of them before moving to the next topic. The \defn{boldface} font is used to indicate the first mention of a term, accompanied by its definition.

While most of the material is not new, there are some novelties in presentation and interpretation, such as:
\begin{itemize}
\item Vector bundle picture for acceleration of particle on a ring in Sec.~\ref{Sec:ParticleRing}.
\item Construction of spectrum for dimerized charge pump from the level (anti)crossing spectra of individual dimers in Sec.~\ref{Sec:EndStatesDimer}.
\item Interpretation of Bloch Hamiltonian as a section of a non-trivial real vector bundle used to explain symmetry-enforced band crossings in Sec.~\ref{Sec:TopoSymCons}.
\end{itemize}

\subsection*{Sources and further reading}
The first part of the notes is loosely based on the textbook \cite{Frankel2011}, which provides a physicist-oriented introduction to the machinery of differential geometry and basic algebraic topology. This is a good place to start learning about manifolds and differential forms, disguised as surfaces and certain well-behaved integrands in our notes. Another friendly account on these topics is given in Ref.~\cite{Baez1994}. The standard reference \cite{Nakahara2003} covers additional material, but is more formal and rather advanced. For the discussion of vector bundles aimed at mathematicians, see Ref.~\cite{Hatcher2003}.

The physical part of the notes shares the general story line with textbooks \cite{Asboth2016, Vanderbilt2018}. Ref.~\cite{Asboth2016} provides a concise example-driven presentation of the theory of topological insulators. In textbook \cite{Vanderbilt2018}, discussion revolves around the concept of Berry phase and includes wide range of topics, combining analytical and computational perspectives. Our notes aim to complement these excellent sources with mathematical point of view. The reader may also consult the introductory article \cite{Fruchart2013}, which uses the language of vector bundles. The paper \cite{Hasan2010}, besides the pedagogical exposition, provides an overview of early theoretical developments and experimental results on topological insulators and superconductors. Further physical insights can be found in the online course \cite{TCM} featuring video mini-lectures by pioneers of the field. A thorough research-level presentation of the subject is given in textbook \cite{Bernevig2013}. A review of Weyl semimetals can be found in Ref.~\cite{Armitage2018}.

There are two major omissions in our discussion. First, we will not consider the integer quantum Hall effect, which can be thought of as a common ancestor of all topological materials. We refer interested reader to the comprehensive lecture notes \cite{Tong2016}. Second, we will not discuss connections between topological insulators and high-energy physics. In particular, we will not use the method of finding the edge modes based on Dirac equation. The reason is that it requires the linearization of Bloch Hamiltonian at some point of momentum space, which obscures relation with global properties of the vector bundle. On the other hand, this construction can be found in almost any introductory text. For the detailed derivation, see Chapter 7 of Ref.~\cite{Asboth2016}. The high-energy perspective on topological insulators is developed in lecture notes \cite{Witten2015}, introductory article \cite{Cayssol2013}, and textbook \cite{Shen2017}.

\subsection*{Acknowledgments}
I am grateful to O.G.\,Kharlanov, K.V.\,Antipin, and A.A.\,Markov for many fruitful discussions. I would like to thank all students who attended the lectures for their tricky questions and keen interest in the subject. I thank the organizers of Topological Quantum Matter conference at Nordita, Stockholm in August 2019 for the hospitality. The travel support by the Foundation for the Advancement of Theoretical Physics and Mathematics ``BASIS'', grant No.\,19-36-013, is gratefully acknowledged. I would like to thank J.\,van Wezel and A.\,Bouhon for insightful discussions during the conference. The work was partially supported by RFBR grant No.\,19-02-00828 A.

\section{Connection on a vector bundle}\label{Chap:ConnVB}
Vector field is one of the fundamental mathematical objects used throughout all areas of physics. But sometimes the properties of the space where the field takes values can be as important as the specific configuration of the field itself. Such spaces are known as vector bundles. In this section, we will consider several examples of vector bundles and introduce some terminology.

One also needs to be able to take directional derivatives of vector fields. Generalization of such derivative to the realm of vector bundles gives rise to the concept of connection and has surprisingly far-reaching physical consequences, as we will see later.

\subsection{From vector fields to vector bundles}
\subsubsection{Idea of a vector bundle}\label{Sec:IdeaVB}
Consider a planar vector field $\vec v$ over the plane $\mathbb{R}^2$, shown in Fig.~\ref{Fig:VBExamples} on the left.  There is a vector space $V_p$ at each point $p\in \mathbb{R}^2$. The value of the field $\vec v$ at the point $p$ is given by the vector $\vec v(p)\in V_p$.  The collection of all vector spaces
\begin{equation}
V = \{V_p\mid p\in \mathbb{R}^2\}
\end{equation}
is our first example of a \defn{vector bundle}. In a sense, this is the space where all possible planar vector fields live. More generally, a vector bundle over a surface can be thought of as a collection of vector spaces attached to the points of the surface. The surface is called \defn{base space} of the bundle, and vector space at some point is known as a \defn{fiber}. All fibers are assumed to have the same dimension and to ``vary smoothly'' along the base space. Such a description is too vague to be useful for mathematicians, but it will suffice for our needs. In this example, the base space is the plane $\mathbb{R}^2$ and the fiber at $p$ is the vector space $V_p$. A \defn{section} $\vec v$ of a vector bundle is a smooth choice of an element $\vec v(p)$ in each fiber $V_p$. For the bundle $V$, section is just another term for vector field.

\begin{figure}\begin{center}
\psfrag{V}{$V$}
\psfrag{Vp}{$V_p$}
\psfrag{vp}{$\vec v(p)$}
\psfrag{p}{$p$}
\psfrag{R2}{$\mathbb{R}^2$}
\psfrag{TS1}{$TS^1$}
\psfrag{S1}{$S^1$}
\psfrag{VS1}{\hspace{-0.3cm}\textcolor{half-gray}{$V|_{S^1}$}}
\psfrag{R3}{$\mathbb{R}^3$}
\psfrag{S2}{$S^2$}
\psfrag{TpS2}{$T_pS^2$}
\psfrag{Wp}{\textcolor{half-gray}{$W_p$}}
\includegraphics{img/1bundles}
\caption{Examples of vector bundles. \label{Fig:VBExamples}}
\end{center}
\end{figure}

Starting from the bundle $V$, we can obtain new bundles by cutting subsets of the base space or vector subspaces of the fiber. Suppose we have a circle $S^1\subset \mathbb{R}^2$ on the plane. Then we can \defn{restrict} the bundle $V$ to the circle by considering only those fibers that are attached to the points of $S^1$:
\begin{equation}
V\rvert_{S^1} = \{V_p\mid p\in S^1\}.
\end{equation}
One can further select a vector subspace in each fiber, thus specifying a \defn{subbundle}. For example, choose in each fiber $V_p$ of $V\rvert_{S^1}$ a line that is tangent to the circle. In this way one obtains the \defn{tangent bundle} to the circle $TS^1 \subset V\rvert_{S^1}$, as shown in the middle panel of Fig.~\ref{Fig:VBExamples}. Its fiber at point $p$ is denoted $T_pS^1$.

In a similar way, we define the tangent bundle to the sphere $TS^2$ (Fig.~\ref{Fig:VBExamples}, right). Start from the vector bundle
\begin{equation}
W= \{W_p\mid p\in R^3\}
\end{equation}
of three-dimensional vector spaces over a three-dimensional space $\mathbb{R}^3$. Then restrict it to the two-dimensional sphere $S^2\subset \mathbb{R}^3$ and further select planes that are tangent to the sphere. Alternatively, one can define these planes in terms of normal vectors $\vec n$ to the sphere. At point $p$, define the plane $T_pS^2$ as an orthogonal complement to $\vec n(p)$ in $W_p$.

\subsubsection{How to differentiate a section?}
Our next goal is to take a directional derivative of a section of a vector bundle at some point $p$ of the base space. We wish the differentiation to be performed ``inside the vector bundle'': the derivative of a section at $p$ must belong to the fiber $V_p$ at this point.

\begin{figure}[ht]
\begin{center}
\psfrag{v}{$\vec v$}
\psfrag{R2}{$\mathbb{R}^2$}
\psfrag{S2}{$S^2$}
\psfrag{x}{$x$}
\psfrag{y}{$y$}
\psfrag{ex}{$\vec e_x$}
\psfrag{ey}{$\vec e_y$}
\psfrag{ph}{\textcolor{half-gray}{$\varphi$}}
\psfrag{th}{\textcolor{half-gray}{$\theta$}}
\psfrag{eph}{$\vec e_\varphi$}
\psfrag{eth}{$\vec e_\theta$}
\includegraphics{img/2Bases}
\caption{Coordinate basis sections for $V$ and $TS^2$.\label{Fig:Bases}}
\end{center}
\end{figure}

Consider again the bundle $V$ over the plane $\mathbb{R}^2$. We introduce Cartesian coordinates $(x, y)$ on the plane, which allows us to define basis sections $\{\vec e_x, \vec e_y\}$ that give a basis $\{\vec e_x(p), \vec e_y(p)\}$ in each fiber $V_p$, as shown in Fig.\ref{Fig:Bases} on the left. Thus the vector $\vec v(p)$ can be decomposed as
\begin{equation}
\vec v(p) = \vec e_i(p) v_i(p),
\end{equation}
where $i=x, y$ and summation over repeated indices is assumed. From now on, we will omit specification of the point $p$ in expressions like this, for clarity. One can easily define a derivative, say, in $x$ direction:
\begin{equation}
\partial_x \vec v =\partial_x (\vec e_i v_i) =\vec e_i \partial_x v_i.
\end{equation}
Note that we assumed that $\partial_x \vec e_i=0$, or that the basis vectors are constant. This need not be true for other basis choices, but one can always express another basis in terms of the Cartesian one. In components, we have
\begin{equation}
v(p) = \begin{pmatrix}
v_x (x,y) \\ v_y (x,y)\\
\end{pmatrix}.
\end{equation}
That is, the section $\vec v$ can be represented as a $\mathbb{R}^2$-valued function, with differentiation defined component-wise.

Now consider a section of the tangent bundle to the sphere $TS^2$ (Fig.~\ref{Fig:Bases}, right). Suppose we wish to take the derivative along a circle of latitude, $\partial_\varphi \vec v$. As above, we introduce a basis $\{\vec e_\theta, \vec e_\varphi\}$ in the tangent plane $V_p$ at each point $p\in S^2$ and decompose the vector in its terms. However, the formula for the derivative
\begin{equation}
\partial_\varphi \vec v =\partial_\varphi (\vec e_\alpha v_\alpha) =\vec e_\alpha \partial_\varphi v_\alpha, \quad \alpha = \theta, \varphi
\end{equation}
fails to work in this case. To see this, note that locally the sphere  resembles a plane, so we can apply our intuitive understanding of the derivative on a flat surface.

Near the equator the pair $\{\vec e_\theta, \vec e_\varphi\}$ looks like a Cartesian basis, and the condition $\partial_\varphi \vec e_\alpha=0$ seems to hold. But near each of the poles, this basis corresponds to the polar coordinate system, so its vectors cannot be declared to be constant. And it is not the unfortunate basis choice that causes the problem. There is a global constraint that forces any tangent vector field on a sphere to have singularities. Thus it is impossible to define a global basis, let alone to find a constant one. Note also that these singularities make a section of $TS^2$ something very different from a function $S^2\to \mathbb{R}^2$, which need not vanish anywhere.

\subsection{Covariant derivative}
\subsubsection{Projection from ambient space}
To find a basis-independent way to differentiate sections of $TS^2$, recall that this bundle can be thought of a part of the bundle $W$ of three-dimensional vector spaces over $\mathbb{R}^3$. Thus one can consider a vector $\vec v(p)\in T_pS^2$ as an element of $W_p$. As such, it can be decomposed as
\begin{equation}
\vec v= \vec e^{\,a}_i v^{\,a}_i, \quad i=x, y, z,
\end{equation}
where $\{\vec e_i^{\,a}\}$ is a Cartesian basis for $W_p$, and superscript $^a$ indicates that we work in this ``ambient space''. As in the case of the flat bundle $V$, we can assume that these Cartesian bases are constant, and calculate the following derivative:
\begin{equation}
\partial_\varphi^{\,a} \vec v = \vec e_i^{\,a} (\partial_\varphi v_i^{\,a}).
\end{equation}
However, the result of this differentiation at a point $p$ need not belong to the plane $T_pS^2$. For example, consider a velocity vector of a traveler who moves around the Earth along the equator. The derivative $\partial_\varphi^{\,a}$ of this field is perpendicular to the fibers of $TS^2$. This problem can be solved by using projection. Define the following operator:
\begin{equation}
\nabla_\varphi \vec v = \Proj (\partial^{\,a}_\varphi \vec v),
\end{equation}
where $\Proj$ is the orthogonal projection\footnote{We assume that the fiber $W_p$ has the standard inner product, which also induces an inner product in $T_pS^2$.} to the plane $T_pS^2$ at the point $p$ where we evaluate the derivative. We will call $\nabla_\varphi$ a \defn{covariant derivative} in $\varphi$ direction. In a similar way, one defines the covariant derivative in $\theta$ direction. Thus, we have constructed a differential operator on $TS^2$ that is immune to the problems described above. Several remarks are in order:
\begin{enumerate}
\item Applied to a section $\vec v$ at point $p$, the operator $\nabla_\varphi$ measures its rate of change inside the fiber $T_pS^2$.
\item By construction, the value of $\nabla_\varphi \vec v$ does not depend on the basis choice for $T_pS^2$.
\item We use ambient space to define $\nabla_\varphi$, so this operator is an external structure for the bundle $TS^2$.
\end{enumerate}

Since projection is a linear operator, $\nabla_\varphi$ inherits properties of the derivative $\partial^{\,a}_\varphi$:
\begin{align}
\text{Linearity: }&\nabla_\varphi (\vec v + \lambda \vec w) =
\nabla_\varphi \vec v + \lambda \nabla_\varphi \vec w, \\
\text{Leibniz rule: }&\nabla_\varphi (f \vec v ) = (\partial_\varphi f) \vec v + f\nabla_\varphi \vec v,
\end{align}
where $\lambda$ is a constant and $f$ is a function.

\subsubsection{Connection coefficients}
Suppose that a basis in fibers of $TS^2$ is chosen. For definiteness, we choose the standard basis $\{\vec e_\theta, \vec e_\varphi\}$. Then we can find the expression for $\nabla_\varphi \vec v$ in terms of components of $\vec v$, as follows.

First, decompose the vector $\vec v$ in terms of basis vectors and use the Leibniz rule:
\begin{equation}
\nabla_\varphi \vec v = \nabla_\varphi(\vec e_\alpha v_\alpha) = \vec e_\alpha \partial_\varphi v_\alpha + (\nabla_\varphi \vec e_\alpha) v_\alpha, \qquad \alpha=\theta, \varphi.
\end{equation}
Since covariant derivative of each basis vector belongs to $V_p$, it can be decomposed in terms of this basis, too:
\begin{equation}
 \nabla_\varphi \vec e_\alpha =
 \vec e_\beta (\nabla_\varphi \vec e_\alpha)_\beta \equiv
 \vec e_\beta \omega_{\varphi \alpha}^\beta,
 \end{equation}
 where $\omega_{\varphi \alpha}^\beta$ are  functions called \defn{connection coefficients}. Finally, we obtain the following expression for the covariant derivative:
\begin{equation}\label{Eq:CDreal}
\nabla_\varphi \vec v =
\vec e_\beta (\partial_\varphi v_\beta + \omega_{\varphi \alpha}^\beta v_\alpha ).
\end{equation}

To find the connection coefficients, we first need to express basis vectors as a three-dimensional vectors. Let $\vec r =(r \sin \theta \cos \varphi, \,r \sin \theta \sin \varphi, \,r \cos \theta)^T$ be a radius vector for a point on the sphere. Then define $\{\vec e_\theta, \vec e_\varphi\}$ as its normalized derivatives $\partial_\alpha \vec r$:
\begin{equation}\label{Eq:Basis3D}
\vec e_\theta = \begin{pmatrix}
\cos\theta\cos\varphi\\
\cos\theta\sin\varphi \\
-\sin\theta \\
\end{pmatrix}
\qquad
\vec e_\varphi =\begin{pmatrix}
-\sin\varphi \\
\cos \varphi \\
0 \\
\end{pmatrix}
\end{equation}
According to the definition given above, to find $\nabla_\varphi \vec e_\varphi$ we project the derivative $\partial_\varphi^{\,a}\vec e_\varphi$ to the tangent plane using the unit normal $\vec n = \vec r/r$:
\begin{equation}
\nabla_\varphi \vec e_\varphi = \partial_\varphi^{\,a}\vec e_\varphi
-\vec n \langle \vec n, \partial_\varphi^{\,a}\vec e_\varphi \ke
\end{equation}

\exr {Calculate $\nabla_\varphi \vec e_\varphi$ and $\nabla_\varphi \vec e_\theta$. Express these vectors  in terms of the basis $\{\vec e_\theta, \vec e_\varphi\}$. }

By results of the exercise, the only non-zero connection coefficients with respect to the basis $\{\vec e_\theta, \vec e_\varphi\}$ are
\begin{equation}
\omega_{\varphi\theta}^\varphi = \cos \theta \qquad
\omega_{\varphi\varphi}^\theta = -\cos\theta.
\end{equation}
With these functions at hand, one need not calculate projections to find the covariant derivative of a section $\vec v$.

\subsubsection{Complex plane notation}
One can further simplify the expression for $\nabla_\varphi\vec v$ by considering a fiber $T_pS^2$ with an orthonormal basis as a complex vector space.

Let $\vec v$ be an element of a two-dimensional real vector space $V$. Once a basis is chosen, $\vec v$ can be represented as a column of components $(v_1, v_2)^T$ with $ v_i\in \mathbb{R}$. Now consider a complex one-dimensional space, or a \defn{complex line} $W$. A choice of basis vector $\vec 1$  allows one to identify any vector $\vec w$ with a complex number: $\vec w = w \vec 1$, where $w\in\mathbb{C}$. We wish to describe the real plane as a complex line. While their (real) dimensions coincide, the latter carries extra structure: multiplication by $i$. In real terms, this is given by a  linear operator $I$ that squares to minus identity and can be described as a $\frac{\pi}{2}$ rotation. We only need to specify the sense of rotation, which is fixed by an orientation of the real plane (a definition of the term ``clockwise''). One says that the choice of $I$ endows $V$ with a \defn{complex structure}.

\begin{figure}[h]
\begin{center}
\psfrag{eph}{\textcolor{half-gray}{$\vec e_\theta$}}
\psfrag{eth}{\textcolor{half-gray}{$\vec e_\varphi$}}
\psfrag{phi}{$\varphi$}
\psfrag{npep}{\hspace{-0.2cm}$\nabla_\varphi \vec e_\theta$}
\psfrag{npet}{$\nabla_\varphi \vec e_\varphi$}
\psfrag{npv1}{\hspace{-0.1cm}$\nabla_\varphi \vec 1$}
\psfrag{arr}{$\Rightarrow$}
\psfrag{v1}{\textcolor{half-gray}{$\vec 1$}}
\includegraphics{img/3BasRotComp}
\caption{Covariant derivatives of orthonormal basis vectors in real and complex notation. \label{Fig:3BasRotComp}}
\end{center}
\end{figure}

Consider a collection of planes attached to points of some path with an orthonormal basis chosen in each plane. The only degree of freedom these vectors have inside the fiber is rotation. Thus the covariant derivatives $\nabla_\varphi \vec e_\alpha$ of the basis vectors  are perpendicular to them and have equal magnitude, as shown in Fig.~\ref{Fig:3BasRotComp}. Let us denote this magnitude by $\omega_\varphi$. Then
\begin{equation}
\nabla_\varphi \vec e_\theta =\omega_\varphi \vec e_\varphi \qquad
\nabla_\varphi \vec e_\varphi =- \omega_\varphi \vec e_\theta.
\end{equation}
In complex notation with $\vec e_\theta =\vec 1$ and $\vec e_\varphi = i\vec 1$ this becomes
\begin{equation}
\nabla_\varphi \vec 1 = \omega_\varphi i\vec 1 = i\omega_\varphi \vec 1 \qquad
\nabla_\varphi (i \vec 1) = -\omega_\varphi  \vec 1.
\end{equation}
It follows that $\nabla_\varphi (i\vec 1) =i\nabla_\varphi \vec 1$, so that $\nabla_\varphi$ is a complex linear differential operator.

The expression (\ref{Eq:CDreal}) for covariant derivative takes the form
\begin{equation}\label{Eq:CD_1}
\nabla_\varphi \vec v = (\partial_\varphi v)  \vec 1 + v i\omega_\varphi \vec 1 = (\partial_\varphi +i \omega_\varphi) \vec v,
\end{equation}
where last equality defines the action of $\partial_\varphi$ on $\vec v$. Thus, when one restrict oneself to use orthonormal bases in real planes, only one connection coefficient is needed to determine the covariant derivative. Informally, $\omega_\varphi$ is the angular velocity of rotation of the basis vector inside the fiber.

Note that the last expression should be used with care since it does not contain the basis vector $\vec 1$ explicitly, and the dependence of both $\omega_\varphi$ and $\partial_\varphi$ on the basis choice is obscure.
\subsubsection{Transformation laws for connection coefficients}\label{Sec:ConnCoefTransform}
\begin{figure}[h]
\begin{center}
\psfrag{ph}{$\varphi$}
\psfrag{v1}{$\vec 1$}
\psfrag{v1pr}{\textcolor{half-gray}{$\vec 1'$}}
\psfrag{bph}{$\beta(\varphi)$}
\psfrag{gph}{\textcolor{half-gray}{$\varphi$}}
\psfrag{th}{\textcolor{half-gray}{$\theta$}}
\psfrag{ga}{$\gamma$}
\includegraphics{img/4Transform}
\caption{Different choice of basis section and change of coordinates for the base space.\label{Fig:Transform}}
\end{center}
\end{figure}

Let us examine how connection coefficients transform under change of basis. Consider two bases related by a position-dependent phase rotation: $\vec 1' = e^{i\beta(\varphi)}\vec 1$, as shown in the left panel of Fig.~\ref{Fig:Transform}. Then
\begin{equation}
\nabla_\varphi \vec 1 = i \omega_\varphi \vec 1
\qquad \nabla_\varphi \vec 1' = i \omega'_\varphi \vec 1'
\end{equation}
On the other hand,
\begin{equation}
\nabla_\varphi (e^{i\beta} \vec 1) = i(\partial_\varphi \beta) e^{i\beta}\vec 1 +e^{i\beta} \nabla_\varphi \vec 1  = i(\omega_\varphi +\partial_\varphi	\beta)e^{i\beta}\vec 1
\end{equation}
by the Leibniz rule. It follows that transformation law for connection coefficients has the form
\begin{equation}\label{Eq:ConnBasis}
\omega'_\varphi = \omega_\varphi+\partial_\varphi\beta.
\end{equation}

Another transformation law is associated with the change of coordinates on the surface. Let us find connection coefficients along some curve, which is defined parametrically as two functions $\theta(\gamma), \varphi(\gamma)$. From the definition of the covariant derivative one finds that
\begin{equation}\label{Eq:ConnCurve}
\omega_\gamma = \frac{d\varphi}{d\gamma}\omega_\varphi +\frac{d\theta}{d\gamma}\omega_\theta.
\end{equation}
This allows one to calculate the covariant derivative in any direction given by a curve parameterized by $\gamma$.

More generally, one can use the functions $\theta(x_1, x_2), \varphi(x_1, x_2)$ to define new coordinates on the sphere. Indeed, for a fixed value of the second coordinate $x_2=x_2^0$,  the pair of functions  $ \theta(x_1, x^0_2)$ and $\varphi(x_1, x^0_2)$ define a coordinate curve parameterized by $x_1$. The connection coefficients along new coordinate curves are given by
\begin{equation}\label{Eq:ConnRep}
\omega_i = (\partial_i \alpha) \omega_\alpha,
\end{equation}
 where $i=1,2$ and $\alpha =\theta, \varphi$.

\subsection{Parallel transport}
Recall that our first attempt to differentiate a section of $TS^2$ failed because we could not find an appropriate constant basis. Now that we have a differential operator, we can find a section that is locally constant along some curve parameterized by $\gamma$. It is given by a solution to the equation
\begin{equation}\label{Eq:DefPT}
\nabla_\gamma \vec v_{PT} = 0.
\end{equation}
If we fix some vector $\vec v_0\in T_pS^2$ at the starting point $p$ of the curve as a boundary condition, the solution is unique. The section $\vec v_{PT}$ is called the \defn{parallel transport} of the vector $\vec v_0$ along the curve.

\subsubsection{Physical interpretation}
As we saw above, the differentiation is based on the idea of a constant basis. One can describe such a basis in physical terms by  by employing a conservation law.

Consider a gyroscope with rotation axis that can move freely. If it has angular momentum $\vec J$, one can transport it an any direction and define a constant vector field in this way, since $\partial_t \vec J =0$. This can be thought of as a definition of a parallel transport for the bundle $W$ of three-dimensional spaces over $\mathbb{R}^3$.

\begin{figure}[h]
\begin{center}
\psfrag{J}{$\vec J(p)$}
\psfrag{gV}{\textcolor{half-gray}{$T_pS^2$}}
\includegraphics{img/5GyroFP}
\caption{Restricted gyroscope and parallel transport in $TS^2$.}
\end{center}
\end{figure}

Now we place the gyroscope on the sphere. Suppose that  internal motion of the gyroscope is restricted in such a way that its axis is always lies in the tangent plane to the sphere, $\vec J(p)\in T_pS^2$ (here we neglect the height of the gimbal). Then if we transport the gyroscope along the surface, the plane will rotate, and generally  $\partial_t \vec J$ will be non-zero and will not belong to $T_pS^2$. But the angular momentum will be conserved inside the plane:
\begin{equation}
\Proj (\partial_t \vec J)=0,
\end{equation}
 where $\Proj$ is the projection to $T_pS^2$. In this way one can define a parallel transport for the tangent bundle to the sphere.

We are using projection again, so it should not be surprising that this parallel transport corresponds to the covariant derivative defined above. Indeed, if we treat the time $t$ as a parameter for a curve on the sphere, we have
\begin{equation}
\Proj (\partial_t \vec J)=0 \quad \Leftrightarrow \quad
\nabla_t \vec J =0.
\end{equation}

Parallel transport in $TS^2$ defined by constrained gyroscope has another famous physical realization. This is rotation of the swing plane of the Foucault pendulum. Indeed, the normal vector to the swing plane is always tangent to the surface of the Earth. The pendulum oscillates freely, and direction of its (oscillating)  angular momentum is locally conserved inside the tangent plane.

\subsubsection{Parallel transport angle}\label{Sec:PTangle}
To apply the developed machinery to a real-life situation, let us find the daily rotation angle of the Foucault pendulum which is located at latitude $\varphi$.

We choose a complex basis vector $\vec 1 = \vec e_\theta$ in each tangent plane along the circle of latitude. The complex connection coefficients are $\omega_\varphi = \cos \theta$ and $\omega_\theta =0$. We expect that the parallel transport preserves length of the transported vector. Indeed, the change of length of the vector happens inside the fiber, so it would lead to a non-zero covariant derivative. Thus, we express the vector field as  $\vec v_{PT} = e^{i\alpha(\varphi)} \vec 1$ with initial condition $\alpha(0)=0$. From the parallel transport equation,
\begin{equation}
 \nabla_\varphi \vec v_{PT} = 0 \quad \Rightarrow \quad
 (\partial_\varphi + i \omega_\varphi)e^{i\alpha(\varphi)}=0,
 \end{equation}
one finds that
\begin{equation}
\partial_\varphi \alpha = - \omega_\varphi \quad \Rightarrow \quad
\alpha(\varphi) = \int_0^{\varphi} (-\omega_{\varphi'}) d\varphi' = -\varphi \cos \theta.
\end{equation}
Let us reflect on this integral for a moment. Covariant derivative vanishes when applied to parallel transported vectors. Thus, connection coefficient defined by $\nabla_\varphi \vec 1= i \omega_\varphi \vec 1$ measures the angular velocity of rotation of $\vec 1$ with respect to the parallel transport. Then its negative $-\omega_\varphi$ gives the angular velocity of rotation of $\vec v_{PT}$ with respect to $\vec 1$. Here, ``velocity'' is measured per unit of coordinate. So, the integration of angular velocity over the curve gives the total angle of rotation.

Finally, for a full circle of latitude,
\begin{equation}
 \Delta \alpha = - 2\pi \cos \theta,
 \end{equation}
which agrees with the observed daily rotation of the Foucault pendulum swing plane. This result shows that in general, the parallel transported vector need not return to itself after traveling along a closed curve.

\exr {Consider the motion of the Foucault pendulum on three circles of latitude: equator and two small circles near the poles. Find corresponding daily rotation angles and explain their magnitudes and signs.}\label{Exr:Foucault3c}

We generalize this situation as follows. Let $\mathcal{C}$ be any closed oriented contour on a sphere, parameterized by $\gamma$ (which increases in positive direction). Assume that some smooth basis $\vec 1$ is chosen at each point of the contour, and let $\omega_\gamma$ be the corresponding connection coefficient. We define the \defn{parallel transport angle} as
\begin{equation}
\Delta \alpha (\mathcal{C}) = \int_\mathcal{C} (-\omega_\gamma)d\gamma.
\end{equation}
We know that $\omega_\gamma$ depends on the parametrization of the path and basis choice in the fibers. Let us examine $\Delta \alpha(\mathcal{C})$ from this point of view. First, the value of $\Delta\alpha$ does not depend on the choice of parametrization $\gamma$. Indeed, let $t$ be another parameter along the curve $\mathcal{C}$. Then
\begin{equation}
 \omega_t = \frac{d\gamma}{dt} \omega_\gamma \quad
 \Rightarrow \quad
 \int_\mathcal{C} (-\omega_t) dt =
 \int_\mathcal{C} \frac{d\gamma}{dt} (-\omega_\gamma) dt=
 \int_\mathcal{C} (-\omega_\gamma)d\gamma,
 \end{equation}
where we again assume that $t$ increases in positive direction.

Moreover, the quantity
\begin{equation*}
\Delta \alpha \bmod 2\pi
\end{equation*}
is independent of the choice of basis $\vec 1$ along the curve $\mathcal{C}$.  First, let $\vec 1'$ be a basis obtained from $\vec 1$ by a smooth variation, $\vec 1' = e^{i\beta(\gamma)}\vec 1$. Then nothing changes:
\begin{equation}
\Delta \alpha' = \int (-\omega'_\gamma)d\gamma =
\int (-\omega_\gamma - \partial_\gamma \beta)d\gamma = \Delta \alpha,
\end{equation}
since $\int (\partial_\gamma \beta) d\gamma =0$ for any smooth function $\beta(\gamma)$ over a closed path. However, the basis $\vec 1'$ can make a full revolution with respect to $\vec 1$, which will be described by a $2\pi$ discontinuity in $\beta$ and will lead to a $2\pi$ shift of $\Delta\alpha$.

One can interpret this invariance geometrically by thinking of vectors $\vec v_{PT}$ as being engraved on the planes along the path $\mathcal{C}$. This vector field itself depends only on the geometry of embedding of the planes into the ambient space. It is clearly unaffected by reparametrization of the path. The angle between initial and final vectors is determined up to $2\pi$, which is reflected in the ambiguity of $\Delta \alpha$ discussed above.

\subsection{Connection on a vector bundle}
\subsubsection{Covariant derivative, parallel transport, and connection}
The notions of parallel transport and covariant derivative are closely related to each other. Indeed, Eq.~\eqref{Eq:DefPT} defines the parallel transport from the covariant derivative. In the reverse direction, given a parallel transport one can \emph{define} a covariant derivative $\widetilde\nabla_\gamma$ as
\begin{equation}\label{Eq:CDfromPT}
\widetilde\nabla_\gamma \vec v = \vec 1_{PT} \partial_\gamma v,
\end{equation}
where $\vec 1_{PT}$ is basis section obtained by the parallel transport in $\gamma$ direction. It follows that this derivative vanishes at parallel transported vectors, $\widetilde \nabla_\gamma \vec 1_{PT}=0$, and we arrive again at Eq.~\eqref{Eq:DefPT}.

Now let us re-examine the expression for the covariant derivative
\begin{equation}\label{Eq:CD}
\nabla_\gamma \vec v = (\partial_\gamma + i\omega_\gamma) \vec v.
\end{equation}
Both the derivative $\partial_\gamma \vec v = (\partial_\gamma v)\vec 1$ and connection coefficient depend on the basis choice, while their sum does not. To understand this, note that  $\nabla_\gamma$ measures the rate of change of $\vec v$ with respect to the parallel transport. For simplicity, suppose that the vectors $\vec v$ have unit length. Then the last equation simply describes addition of relative angular velocities. To find the angular velocity of $\vec v$ with respect to the parallel transport, we need to sum the velocity of $\vec v$ with respect to the basis $\vec 1$ and the velocity of $\vec 1$ with respect to the parallel transport.

\defn{Connection} on a vector bundle is an external structure that allows one to take derivatives of sections. This structure can be equivalently thought of as a covariant derivative or as a parallel transport operator. Connection for $TS^2$ introduced above is far from unique. Imagine, for example, declaring another basis in $W$ to be constant, instead of the Cartesian one. However, this connection is distinguished in that it is defined via projection from another bundle with connection. We will call such structures \defn{projected connections}.

\subsubsection{Complex line bundle with connection}\label{Sec:CLBundles}
In the next sections, we will need the \defn{plane bundle} with a real plane as a fiber. The construction is similar to that of $TS^2$ given in Sec.~\ref{Sec:IdeaVB}, but without tangency condition.

Consider a smooth, non-vanishing three-dimensional vector filed $\vec m$ defined on a surface $\mathcal{B}\subset \mathbb{R}^3$. One can think of it as a section of the bundle $W\rvert_\mathcal{B}$. At each point $p\in \mathcal{B}$ the vector $\vec m(p)$ determines the subspace $M_p\subset W_p$ such that $\vec m(p) \perp M_p$. This defines a bundle $M$ of real planes over $\mathcal{B}$.  In a particular case when the vectors $\vec m$ are surface normals, we obtain the tangent bundle $T\mathcal{B}$ to the surface $\mathcal{B}$.

\begin{figure}[ht]
\begin{center}
\psfrag{R3}{$\mathbb{R}^3$}
\psfrag{M}{$M$}
\psfrag{Surf}{$\mathcal{B}$}
\psfrag{vm}{$\vec m$}
\includegraphics{img/6PlaneBundle}
\caption{A plane bundle $M$ over a two-dimensional surface $\mathcal{B}$.}
\end{center}
\end{figure}

A crucial difference between the bundle $M$ and the tangent bundle is that sections of $M$ are not related to the base space. Recall that we constructed the basis sections $\{\vec e_\theta, \vec e_\varphi\}$ for $TS^2$ from velocities $\partial_\alpha \vec r$ of a point moving on the sphere along the coordinate curves. In contrast, the fibers and sections of the bundle $M$ are independent of the geometry and coordinates of the base space $\mathcal{B}$.

Each plane $M_p$ is oriented by its normal vector $\vec m$ and the right-hand rule (i.e., the orientation of three-dimensional space $W_p$). Thus it can be identified with complex one-dimensional space, making $M$ into \defn{complex line bundle} over $\mathcal{B}$. Note that the switching the direction of normals, $\vec m\to -\vec m$ leaves real planes invariant, but changes their orientation. In complex terms, this amounts to replacement $i \to -i$, or complex conjugation. As we will see, it is important to distinguish between a complex line bundle and its conjugate version.

As a final ingredient, we define the projected connection on $M$. This gives us the way to perform the parallel transport of any vector along any curve on $\mathcal{B}$. Once a basis section $\vec 1$ is chosen, a connection coefficient $\omega_i$ can be calculated along each coordinate curve $x_i$ in $\mathcal{B}$. This allows one to find the covariant derivative of a section $\vec v$ of $M$. The whole construction can be  easily generalized to base spaces of other dimensions.

\subsubsection{Connections in general}

In mathematics, one often uses properties of some concrete object as an axiomatic definition of its abstract version. In this way, a group of transformations becomes an abstract group, and an embedded surface inspires the definition of an abstract manifold. Vector bundles and connections can also be defined abstractly. For example, one can describe a fiber of a tangent bundle in terms of velocity vectors of points moving through a given point of the base. Connection is defined as a differential operator with certain properties, such as Leibniz rule. But in the present notes, we will not discuss these abstract versions. In the cases that are most important to us, vector bundles arise naturally as subbundles of other bundles and are equipped with projected connection.

In our definition of the covariant derivative, we used the orthogonal projection, and thus the inner product. The smooth choice of the inner product in each fiber of a tangent bundle is known as \defn{metric}. Once the tangent bundle is endowed with metric, the base space becomes a  \defn{Riemannian manifold}. The machinery of connections originated in the context of Riemannian geometry, which focuses on the metric. For example, one can express the connection coefficients in terms of the metric tensor. However, we will not develop this point of view, since our main interest lies in the geometry of non-tangent bundles.

Finally we note that connection describes a local, or geometric property of a vector bundle. Indeed, the covariant derivative is a local differential operator, and parallel transport is defined between nearby points along some curve. Despite this local nature, connections can be used to determine the global properties of vector bundles. Later we will discuss this phenomenon in detail.

\section{Electromagnetic field and curvature of connection}\label{Chap:EMF}
Vector bundles provide the natural language for description of gauge fields. In this section, we will use a complex line bundle defined in Sec.~\ref{Sec:CLBundles} to visualize the geometry of electromagnetic field. We will also introduce the concept of curvature, an important local characterization of connection on a vector bundle.

\subsection{Particle in electromagnetic field}
\subsubsection{From classical to quantum description}
Consider a point-like classical particle of an elementary electric charge $q=|e|$ in electromagnetic field. It  feels the Lorentz force:
\begin{equation}
\vec  F= q(\vec{\emf E} + [\vec v \times \vec {\emf B}]).
\end{equation}
The corresponding Hamiltonian reads
\begin{equation}
H=\frac{1}{2m} \sum_j (p_j-q\emf A_j)^2 +q\phi
\end{equation}
where $\phi$ and $\emf A_j$ are electric and magnetic potentials of the field and $j = x, y, z$.

If the particle is microscopic, quantum mechanics prescribes to turn momentum into the differential operator  $\hat p_j = -i\hbar \partial_j$, so that the quantum Hamiltonian becomes
\begin{equation}
\hat H = \sum_j\frac{1}{2m}(-i\hbar \partial_j -q \emf A_j)^2 +q\phi
\end{equation}
Then a wave function $\psi$ describing our particle obeys the Schr\"odinger equation
\begin{equation}
i\hbar \partial_t \psi = \hat H \psi,
\end{equation}
which can be rewritten as
\begin{equation}
\biggl(\partial_t + i\frac{q}{\hbar} \phi\biggr) \psi =
\sum_j \frac{(-i\hbar)^2}{2m} \biggl(\partial_j-i\frac{q}{\hbar} \emf A_j\biggr)^2 \psi.
\end{equation}

\subsubsection{Vector potential and connection coefficients}
Comparing the last equation with Eq.~\eqref{Eq:CD}, we note that the combinations in parentheses can be interpreted as coordinate expressions of a covariant derivative. Then, ``multiplying'' by a basis section $\vec 1$ from the right, we obtain:
\begin{equation}
\nabla_t\vec \psi =\sum_j \frac{1}{2m} (-i\hbar \nabla_j)^2 \vec \psi.
\end{equation}
This suggests the following geometric interpretation. We have an abstract complex line bundle over the four-dimensional spacetime. There is a connection on this bundle that represents the electromagnetic field and enables covariant differentiation. State of the particle is given by a section $\vec \psi$ of this bundle. Once the basis section $\vec 1$ is chosen, the section $\vec \psi$ can be described as a complex-valued wave function $\psi(x, y, z, t)$ via $\vec \psi = \psi \vec 1$. The connection coefficients with respect to $\vec 1$ are gauge potentials of the electromagnetic field:
\begin{equation}
\omega_t = \frac{q}{\hbar}\phi, \qquad \omega_j = -\frac{q}{\hbar}\emf A_j.
\end{equation}
Finally, local rotations of basis correspond to gauge transformations: for another basis section $\vec 1' = e^{i\beta}\vec 1$ one has
\begin{eqnarray}
\vec \psi = \psi e^{-i\beta} e^{i\beta} \vec 1 & \qquad \Rightarrow \qquad &\psi' = \psi e^{-i\beta}, \\
\omega'_j = \omega +\partial_j \beta & \qquad \Rightarrow \qquad &\emf A'_j = \emf A_j - \frac{\hbar}{q}\partial_j\beta.
\end{eqnarray}
Thus, a gauge transformation manifestly does not affect the state of the particle and the configuration of the electromagnetic field. It merely changes the way of description of the same physical and geometrical situation.

However, beside the conceptual appeal, it is not clear whether this geometric picture can be useful for a physicist. In a world without Dirac monopoles we are probably living in, one can simply choose a global section $\vec 1$ and forget about vector bundle altogether. But it turns out that the structures similar to gauge fields arise in other physical systems, and those fields do allow monopole configurations, which will be important for us later. Before discussing monopoles, we need to understand the geometry of gauge fields in more detail. To do so, we will combine geometric picture of the complex line bundle described in Sec.~\ref{Sec:CLBundles} with the physics of electromagnetic field, in the following way. We start with the bundle $M$ over a three-dimensional space and assume that its section  $\vec \psi$ represents a wave function of a particle. Then one can ask: what kind of magnetic field is described by the projected connection on $M$? Below, we will explore this question, and establish the relation between the shape of the bundle and the configuration of the corresponding electromagnetic field.

\subsection{Magnetic flux from parallel transport angle}
\subsubsection{Flux ambiguity}\label{Sec:FluxAmb}
Let $M$ be a complex line  bundle over $\mathbb{R}^3$ with coordinates $x_j = x, y, z$. Assume that some basis section $\vec 1$ is specified. The connection coefficients $\omega_j$ correspond to the components of magnetic vector potential $\emf A_j$. Let $\mathcal{C}\subset R^3$ be an oriented closed contour defined parametrically by three functions $x_j(\gamma)$. Then connection coefficient with respect to the parameter $\gamma$ is given by the sum $\omega_\gamma = (\partial_\gamma x_j) \omega_j$.

Let us find the parallel transport angle for the contour $\mathcal{C}$:
\begin{equation}
\Delta \alpha(\mathcal{C}) = \int_\mathcal{C} (-\omega_\gamma) d\gamma = \int_\mathcal{C} (-\omega_j) dx_j = \frac{q}{\hbar}\int_\mathcal{C} \vec{\emf A}\cdot d\vec l = 2\pi \frac{\Phi(\mathcal{B})}{\Phi_0},
\end{equation}
where $\Phi(\mathcal{B})$ is the magnetic flux threading some surface $\mathcal{B}$ bounded by the contour and $\Phi_0 = \frac{h}{e}$ is the magnetic flux quantum. Thus, we can express the flux through a surface $\mathcal{B}$ in terms of the angle of the parallel transport  along its boundary:
\begin{equation}\label{Eq:FluxDa}
\Phi(\mathcal{B}) = \frac{\Delta \alpha(\mathcal{C})}{2\pi}\Phi_0.
\end{equation}
Recall that $\Delta\alpha$ depends on the choice of a basis section $\vec 1$. It follows that Eq.~(\ref{Eq:FluxDa}) determines the flux through the surface $\mathcal{B}$ only up to a flux quantum. In order to deal with this ambiguity, we need to determine the flux at each point of $\mathcal{B}$ from the local geometry of the vector  bundle.

\subsubsection{A note on orientation}
To fix the signs of line and surface integrals we need to define an \defn{orientation} of spaces in a coherent way. One says that a contour is oriented, if there is a preferred direction. This can be fixed by specifying a tangent vector at some point. For a two-dimensional surface, orientation means a preferred sense of rotation. This is given by an ordered pair of tangent vectors. To orient a three-dimensional space, one chooses handedness, or an ordered triple of vectors.
We will deal only with \defn{orientable} spaces, meaning that they admit a consistent global choice of orientation in all tangent spaces.

\begin{figure}[h]
\begin{center}
\psfrag{e1}{$\vec e_1$}
\psfrag{e2}{\hspace{-0.1cm}$\vec e_2$}
\psfrag{e3}{\hspace{-0.1cm}$\vec e_3$}
\includegraphics{img/7Orient}
\caption{Orientation of spaces given by ordered sets of tangent vectors.}
\end{center}
\end{figure}

Once a given space is oriented, there is a standard  way to orient its subspaces. For example, let us orient $\mathbb{R}^3$ by the right-hand rule, and let the triple $\{\vec e_1, \vec e_2, \vec e_3\}$ be positively oriented. To orient a surface $\mathcal{B}\subset \mathbb{R}^3$, we need to specify a normal $\vec n$ (if the surface is the boundary of some region, it is common to choose the outer normal). Then a pair $\{\vec e_1, \vec e_2\}$ of tangent vectors is positively-oriented if the triple $\{\vec n,  \vec e_1, \vec e_2\}$ is right-handed in $\mathbb{R}^3$.

In a similar way, one can orient any contour $\mathcal{C}$ on the surface. Suppose the contour is a boundary of some region $\Sigma\subset \mathcal{B}$, so that $\mathcal{C}=\partial \Sigma$. Consider a pair of vectors tangent to the surface, $\{\vec e_n, \vec e_1\}$, where $\vec e_1$ is tangent to the contour, and $\vec e_n$ is an outer normal of the region $\Sigma$. Then $\vec e_1$ defines positive direction in $\mathcal{C}$, if the pair is positively-oriented.

We will also need the notion of \defn{positively-oriented coordinates}. To any coordinate curve one can associate a velocity vector of a point that moves in positive direction. Thus, the orientation defined in terms of ordered sets of tangent vectors similarly applies to the ordered sets of coordinates.

\subsubsection{Flux through a coordinate plaquette}
For some inspiration from Cartesian coordinates, consider a small oriented plaquette  $\boxdot$ lying on the $xy$ plane, with a contour $\partial \,\boxdot$ as a boundary:
\begin{equation}
\partial \,\boxdot = (x_0,y_0) \rightarrow (x_0+\Delta x, y_0) \rightarrow
(x_0+\Delta x, y_0+\Delta y)
\rightarrow(x_0, y_0 +\Delta y) \rightarrow  (x_0,y_0)
\end{equation}
We wish to compute a magnetic flux through it as $\Phi(\boxdot)=\int_{\partial \,\boxdot }  \vec{\emf A}\cdot d\vec l$. For the second edge, in linear approximation, one has:
\begin{equation}
\int_{y_0}^{y_0+\Delta y} \emf A_y(x_0+\Delta x, y) dy \approx \emf A_y(x_0+\Delta x, y_0)\Delta y \approx  \emf A_y(x_0, y_0) \Delta y +\tfrac{\partial \emf A_y}{\partial x}\Big\rvert_{(x_0,y_0)} \Delta x \Delta y
\end{equation}
Adding together analogous expressions for the other edges of an infinitesimal contour, we obtain a familiar result:
\begin{equation}\label{Eq:FluxBoxdot}
\Phi(\boxdot)= (\partial_x \emf A_y - \partial_y \emf A_x)dxdy = \emf F_{xy}dxdy = \emf B_z dS
\end{equation}
where $\emf F_{xy}$ is a component of a field strength tensor, or equivalently the $z$ component of the magnetic field strength $\vec {\emf B}$.

\subsection{Curvature of connection}
\subsubsection{Curvature components}
Consider an oriented surface $\mathcal{B}\subset \mathbb{R}^3$ with positively-oriented coordinates $(x_1, x_2)$. Choose a basis section $\vec 1$ for the bundle $M\rvert_\mathcal{B}$.  By the same token as above, we calculate the parallel transport angle for the boundary of coordinate plaquette  $\boxdot$,
\begin{equation}\label{Eq:CurvatureDef}
\Delta\alpha(\partial \,\boxdot)= \bigl[\partial_1(-\omega_2) - \partial_2 (-\omega_1)\bigr] dx_1 dx_2\equiv f_{12} dx_1 dx_2
\end{equation}
where $f_{12}$ stands for the $12$-component of the \defn{curvature of connection}.
Importantly, this calculation assumes that connection coefficients are smooth functions of coordinates, that is, the basis section $\vec 1$ over $\boxdot$ is smooth. Note also that the sign of $\Delta \alpha$ depends both on the orientation of the base space and the orientation of the fiber (the latter being given here by the complex structure).

If the base space has dimension more than two, there will be a curvature component for each pair of coordinate indices. Since the field strength is invariant under gauge transformations, one can expect that the curvature is independent of the basis choice.

\exr {Consider another basis section, $\vec 1' = \vec 1 e^{i \beta(x_1, x_2)}$, where $\beta$ is some smooth function. Show that $f_{12}' = f_{12}$.}

\exr{One can define curvature coefficient in a manifestly gauge-invariant way through the commutator of covariant derivatives. Prove that
\begin{equation}
[\nabla_1, \nabla_2] \vec v = \frac{1}{i}f_{12} \vec v.
\end{equation}
\label{Ex:CurvComm} }

\subsubsection{Curvature and parallel transport angle for $TS^2$}
As an example, we calculate the curvature of the projected connection on the tangent bundle over the sphere $TS^2$. Recall that with the basis choice $\vec 1 = \vec e_\theta$ the only non-zero connection coefficient is $\omega_\varphi = \cos \theta$, so that
\begin{equation}
\Delta \alpha(\partial \,\boxdot)=(\partial_\theta (-\omega_\varphi ) - \partial_\varphi (-\omega_\theta)) d\theta d\varphi  = \sin\theta d\theta d\varphi.
\end{equation}
Thus, the $\theta\varphi$-component of the curvature reads
\begin{equation}
f_{\theta\varphi} = \sin \theta.
\end{equation}
Observe that for elementary contour, the parallel transport angle is given by an elementary solid angle: $\Delta\alpha(\partial \,\boxdot) = d\Omega$.

\begin{figure}[h]
\begin{center}
\psfrag{gdO}{\textcolor{half-gray}{$d\Omega$}}
\psfrag{dSig}{$\partial \Sigma$}
\psfrag{da}{$\Delta\alpha$}
\psfrag{gdt}{\textcolor{half-gray}{$d\theta$}}
\psfrag{gdp}{\textcolor{half-gray}{$d\varphi$}}
\psfrag{vpt}{$\vec v_{PT}$}
\includegraphics{img/8OmegaDalpha}
\caption{Parallel transport angle for an elementary coordinate contour on a sphere.}
\end{center}
\end{figure}

Consider now a region $\Sigma\subset S^2$ with a boundary $\partial \Sigma$ parameterized by $\gamma$. Assume that the section $\vec 1 $ is smooth over $\Sigma$. Let us break $\Sigma$ into small plaquettes $\boxdot$ and add together all contour integrals over their boundaries  $\partial \,\boxdot$. Since the contributions from the internal links cancel out, we have for the parallel transport angle
\begin{equation}\label{Eq:SumPlaq}
\int_{\partial \Sigma} (-\omega_\gamma) d\gamma
=\sum_{\boxdot}\int_{\partial \,\boxdot} (-\omega_j)dx_j =
\sum_\boxdot \Delta\alpha(\partial \,\boxdot) = \int_\Sigma d\Omega.
\end{equation}
It follows that the parallel transport angle for the boundary $\partial\Sigma$ is given by the solid angle spanned by the region $\Sigma$:
\begin{equation}
\Delta\alpha(\partial \Sigma) = \Omega(\Sigma).
\end{equation}
This is a purely geometric quantity, which is manifestly independent of the basis choice. If the region $\Sigma$ is negatively-oriented, so are the constituent plaquettes, and the parallel transport angle changes its sign.

Now we need to reconcile this result with the $2\pi$ ambiguity of $\Delta\alpha$. It turns out that there is no ambiguity once we demand that the basis section is smooth over the whole region $\Sigma$. Indeed, let $\vec 1$ be such a section. Then the basis changes along $\partial \Sigma$ that lead to $2\pi$ jumps inevitably create singularities of $\vec 1$ in the interior of $\Sigma$, and thus are not allowed.

\exr { \begin{enumerate}[a.]
\item Use the last equation to find the daily rotation angle of the Foucault pendulum for the three cases from Ex.~\ref{Exr:Foucault3c}. Explain the discrepancy in the results.
\item Introduce a new section $\vec 1_N$, obtained by parallel transporting some vector from the north pole along all meridians. Compute connection coefficients (to avoid projection, use transformation laws). Find $\Delta \alpha$ for three circles as contour integrals.
\end{enumerate}}\label{Exr:1N}

\subsubsection{Stokes' theorem}
The approach described above is not restricted to $TS^2$ and can be applied to any complex line bundle $M$. For a two-dimensional region $\Sigma$ of a base space $\mathcal{B}$, the parallel transport angle over the boundary contour $\partial \Sigma$ is related to the curvature of the bundle over $\Sigma$ as
\begin{equation}\label{Eq:DaStokes}
\Delta \alpha(\partial \Sigma) = \int_{\partial \Sigma} (-\omega_\gamma) d\gamma  =
\int_\Sigma f_{12} dx_1 dx_2,
\end{equation}
and there is no ambiguity in the first integral if the basis section $\vec 1$ is smooth over $\Sigma$.

In terms of electromagnetism and vector calculus, this is a version of the \defn{Stokes' theorem} that relates circulation of the vector potential with the integral of flux piercing the surface:
\begin{equation}
\int_{\partial\Sigma} \vec{\emf A}\cdot d\vec l
 = \int_\Sigma [\vec \nabla \times \vec{\emf A}]\cdot d\vec S.
 \end{equation}
Note that the last equation works only for a three-dimensional vector field $\vec{\emf A}$ over $\mathbb{R}^3$. In contrast, Eq.~\eqref{Eq:DaStokes} can be applied to bundles with base spaces of higher dimension.
It has even more simple and general formulation
\begin{equation}
\int_{\partial\Sigma} \omega = \int_\Sigma d\omega
\end{equation}
in the language of differential forms. Without going into any further detail, we note that the differential operator $d$ includes gradient, divergence, and curl as special cases.

\subsection{Magnetic flux from geometry of vector bundle}\label{Sec:FluxGeom}
In this section, we will generalize the results for $TS^2$ to a general vector bundle $M$ with projected connection. Then we will find expressions for magnetic flux and field strength described by this bundle.
\subsubsection{Curvature and Jacobian}\label{Sec:CurvJacobian}
Consider a plane bundle $M$ over surface $\mathcal{B}$ defined by a vector field $\vec m$ of unit normals, $|\vec m|=1$. Let $(x_1, x_2)$ be coordinates on the surface near a point $p\in \mathcal{B}$. The unit vector $\vec m(p)$ associates $p$ with a point $m(p)$ on a unit sphere, as shown in Fig.~\ref{Fig:mCurvature}. In this way, the field $\vec m$ gives rise to a map $m: \mathcal{B}\to S^2$ from the surface to the sphere. In coordinates, this map is described as two functions $\theta(x_1, x_2)$ and $\varphi(x_1, x_2)$. In a ``good'' situation (to be specified in a moment), one can  consider them as a definition of the new coordinates on the sphere near $m(p)$, as discussed in Sec.~\ref{Sec:ConnCoefTransform}.  Let us express curvature of $TS^2$ in these coordinates.

\begin{figure}[h]
\begin{center}
\psfrag{m}{$m$}
\psfrag{p}{$p$}
\psfrag{vm}{$\vec m$}
\psfrag{x1}{$x_1$}
\psfrag{x2}{$x_2$}
\psfrag{mp}{$m(p)$}
\psfrag{Surf}{$\mathcal{B}$}
\psfrag{S2}{$S^2$}
\includegraphics{img/9mCurvature}
\caption{The image of a coordinate contour under the map to the sphere defined by a vector field.\label{Fig:mCurvature}}
\end{center}
\end{figure}

Using Eq.~(\ref{Eq:ConnRep}), we have:
\begin{equation}
f_{12} = \partial_1(-\omega_2) - \partial_2 (-\omega_1) = \partial_1(-(\partial_2 \alpha)\omega_\alpha) - [1\leftrightarrow 2],
\end{equation}
where $\alpha = \theta, \varphi$. By the chain rule,
\begin{equation}
\partial_1 \omega_\alpha = (\partial_\beta \omega_\alpha)(\partial_1\beta).
\end{equation}
\exr{Show that the curvature component with respect to the coordinates $(x_1, x_2)$ on the sphere is
\begin{equation}\label{Eq:CurvatureJ}
f_{12} = J f_{\theta \varphi},
\end{equation}
where $J$ stands for the Jacobian determinant
\begin{equation}
J=\det \begin{pmatrix}
\partial_1 \theta & \partial_2\theta\\
\partial_1\varphi & \partial_2\varphi
\end{pmatrix}.
\end{equation}\label{Ex:f12}
}

The matrix in the last equation is known as the \defn{differential} $Dm$ of the map $m$, a linear transformation that acts from the tangent space at $p$ to the tangent space at $m(p)$.  The differential maps velocity vector of a point that moves through $p$ to the velocity vector of its image.  Thus, if the images of velocity vectors associated with coordinate curves are linearly independent, the images of these curves can define coordinates near $m(p)$. This happens when $Dm$ is non-degenerate, that is, $J$ is non-zero. It is the ``good'' situation we referred to above.

Now we need to return from the sphere to the surface $\mathcal{B}$. For each point $p\in \mathcal{B}$ one can identify $M_p$ and $T_{m(p)}S^2$ as subspaces of $\mathbb{R}^3$. This allows us to transfer section of $TS^2$ over some neighborhood of $m(p)$ to section of $M$ near $p$. The curvature of $M$ over $\mathcal{B}$ is determined by the parallel transport angle $\Delta \alpha (\partial \,\boxdot)$ for an infinitesimal contour $\partial \,\boxdot$ near the point $p$. Now observe that the planes over this contour, the basis section, and parametrization are by construction equivalent to those of the contour $m(\partial \,\boxdot)$ on the sphere. Note also that in the degenerate case $J=0$, the coordinate contour is mapped to a line or even a point, so the angle $\Delta \alpha (\partial \,\boxdot)=0$, and the curvature vanishes. We conclude that $f_{12}$ given by Eq.~\eqref{Eq:CurvatureJ} is the desired curvature of the projected connection on the bundle $M$.

\subsubsection{Flux as total solid angle}\label{Sec:FluxOmega}
Now, to find the flux described by the bundle over $\Sigma$, or the parallel transport angle, we integrate  over $\Sigma$ the curvature given by Eq.~\eqref{Eq:CurvatureJ}.  Since the curvature of the tangent bundle $TS^2$ is $f_{\theta\varphi}=\sin\theta$, we have
\begin{equation}\label{Eq:IntCurvJ}
\Delta \alpha (\partial \Sigma )=
\int_\Sigma f_{12} dx_1dx_2 =
\int_\Sigma J \sin\theta dx_2dx_2=
\int_{m(\Sigma)} d\Omega.
\end{equation}
In the last equality we used change of variables formula for multi-dimensional integrals. The result is, in a sense, the \defn{total solid angle} $\Omega(m(\Sigma))$ spanned by the image $m(\Sigma)$ on the sphere. Note that integration takes the sign of $J$ into account.

\begin{figure}[h]
\begin{center}
\psfrag{m}{$m$}
\psfrag{vm}{$\vec m$}
\psfrag{S1}{$S^1$}
\psfrag{gL}{\textcolor{half-gray}{$L$}}
\psfrag{gga}{\large\textcolor{half-gray}{$\gamma$}}
\psfrag{th}{$\theta$}
\psfrag{gmL}{\textcolor{half-gray}{$m(L)$}}
\psfrag{b}{$\beta$}
\includegraphics[width=0.6\textwidth]{img/10LowDimJ}
\caption{Vector field on a line segment defines a map to the circle.\label{Fig:LowDimJ} }
\end{center}
\end{figure}

To clarify the meaning of this result, let us consider a similar situation in lower dimension. Let $\vec m$ be a planar unit vector field on a line segment $L$, as shown in Fig.~\ref{Fig:LowDimJ}.  The field defines a map $m: L \to S^1$ to the circle. Assume that $L$ and $S^1$ are oriented, and $\gamma$ and $\beta$ are their respective positively-oriented coordinates. In coordinates, the map $m$ is given by a function $\beta(\gamma)$. We are interested in the total angle $\theta$ spanned by the image $m(L)$:
\begin{equation}
\theta = \int_{m(L)} d\beta = \int_L \frac{d\beta}{d\gamma}d\gamma.
\end{equation}
Essentially, we are integrating the angular velocity and obtain the resulting angle of rotation. Despite the folding of the image, different layers have opposite signs of the derivative, and their contributions cancel each other out. Note that the sign of the derivative also tells us if the map preserves the orientation at a given point. Finally,  the value of $\theta$ is determined by the image of the boundary $m(\partial L)$, that is, endpoints of the line segment. If the image $m(L)$ covers the circle more that once, this value is shifted by an integer multiple of $2\pi$.

These observations are readily generalized to the two-dimensional situation. Here, the Jacobian plays the role of the derivative above. Once the coordinates $(x_1, x_2)$ and $(\theta, \varphi)$ are positively-oriented, the sign of $J(p)$ indicates whether the map $m$ preserves orientation at $p$ (one can deduce this from the geometric meaning of the differential). If the image $m(\Sigma)$ develops folding, the layers with different signs of $J$ do not contribute to the integral. The parallel transport angle $\Delta \alpha(\partial\Sigma)$ is determined by the solid angle spanned by the image of the boundary $m(\partial\Sigma)$. This value can be shifted by $4\pi$ if the image covers the full sphere.

Thus the magnetic flux described by the bundle over $\Sigma$ is proportional to the total solid angle spanned by the image of $\Sigma$ under $m$:
\begin{equation}
\Phi(\Sigma) = \Phi_0\frac{\Omega(m(\Sigma))}{2\pi},
\end{equation}
which is a basis-independent and unambiguous expression.

\begin{figure}[h]
\begin{center}
\psfrag{m}{$m$}
\psfrag{vm}{$\vec m$}
\psfrag{Sig}{$\Sigma$}
\psfrag{S2}{$S^2$}
\psfrag{mSig}{$m(\Sigma)$}
\psfrag{Jp}{$J>0$}
\psfrag{Jn}{$J<0$}
\psfrag{Jz}{$J=0$}
\includegraphics[width=0.7\textwidth]{img/11CurvExamples}
\caption{Examples of vector fields $\vec m$ over a surface $\Sigma$ that describe positive, negative, and zero curvature.\label{Fig:CurvExamples}}
\end{center}
\end{figure}
Figure \ref{Fig:CurvExamples} shows three vector fields over a surface $\Sigma$ and corresponding images $m(\Sigma)$ on a sphere (we assume that no folding occurs here). Arrow on the boundary contour $\partial\Sigma$ indicates the orientation of the surface. On the left, the map $m$ preserves orientation, and the curvature of the bundle over $\Sigma$ is positive. In the middle, the image of the boundary is the same, but the whole image is negatively-oriented on the sphere, hence the curvature is negative. On the right, we have a constant map with image collapsed into a single point. The map is degenerate, $J=0$, and curvature vanishes.

The first two cases also illustrate that the total solid angle $\Omega(m(\Sigma))$ is determined by the image of the boundary $m(\partial \Sigma)$, up to $4\pi$. Suppose that in the first case the solid angle spanned by the image is $\Omega_1$. Then the solid angle for the second case is $\Omega_2 = -(4\pi -\Omega_1)$, where the first minus sign is due to the the orientation reversal. Thus, $\Omega_1 = \Omega_2+4\pi$, so the two images ``differ by a full sphere''.

\subsubsection{Field strength and Gaussian curvature}
Recall that the definition of the curvature component in Eq.~\eqref{Eq:CurvatureDef} was inspired by the electromagnetic field tensor in Eq.~\eqref{Eq:FluxBoxdot}. The notion of field strength is slightly different, since it is related to the area element. Magnetic flux through a surface $\Sigma$ can be expressed as
\begin{equation}
\Phi(\Sigma) = \int_\Sigma \emf F_{12} dx_1 dx_2  = \int_\Sigma \vec {\emf B} \cdot d\vec S.
\end{equation}
In the Cartesian coordinates for $\mathbb{R}^3$, the area element $dS_{ij}$ is given simply by $dx_idx_j$, so that
\begin{equation}
\emf F_{xy} = \emf B_z \quad \emf F_{xz} =-\emf B_y \quad \emf F_{yz}=\emf B_x
\end{equation}
where $\emf B_y$ has a negative sign since if one orients $xz$ plane by specifying the normal  $\vec n =\vec e_y$, the pair $(\vec e_x, \vec e_z)$ has the negative orientation. Geometrically, these components  are proportional to the three curvature components $\emf F_{ij} = \frac{\Phi_0}{2\pi}f_{ij}$. In a similar fashion, the components of electric field $\emf E_j$ can be thought of as curvatures over the elementary space-time contours $(x_j, t)$.

In general, magnetic field strength can be found by inverting the formula for a flux through a surface. Let $\mathcal{C}$ be a boundary of a small surface element  with normal $\vec n$ near point $p$. Then the normal component of $\vec {\emf B}$ is
\begin{equation}
\emf B_n(p) = \lim_{\mathcal{C}\rightarrow p} \frac{\Phi(\mathcal{C})}{S(\mathcal{C})},
\end{equation}
where $\Phi(\mathcal{C})$ is the magnetic flux through the surface element and $S(\mathcal{C})$ is its area.

Note that since magnetic flux is determined by the solid angle swept by $\vec m$, it depends on how ``curved'' the field $\vec m$ is. The relation to geometric curvature becomes more direct for the tangent bundle $M=T\mathcal{B}$ to a surface $\mathcal{B}\subset \mathbb{R}^3$. Here $\vec m=\vec n$ defines the \defn{normal map} $n: \mathcal{B}\to S^2$ from the surface to the unit sphere.  Defining a quantity similar to the field strength $\emf B_n$, one obtains the \defn{Gaussian curvature} of the surface:
\begin{equation}
\kappa(p)=\lim_{\mathcal{C}\rightarrow p}	\frac{\Omega(n(\mathcal{C}))}{S(\mathcal{C})},
\end{equation}
where $\Omega(n(\mathcal{C}))$ is the solid angle enclosed by the image $n(\mathcal{C})$ on the sphere and $S(\mathcal{C})$ is the area of the region enclosed by $\mathcal{C}$ on the surface. It turns out that the Gaussian curvature  can also be expressed as
\begin{equation}\label{Eq:kRadii}
\kappa(p)=\frac{1}{R_1 (p) R_2(p)},
\end{equation}
where $R_j$ are \defn{principal radii of curvature}, which are defined as follows. Consider a plane that contains a surface normal at $p$. The intersection of this plane with the surface gives some plane curve. Then $R_1$ and $R_2$ are minimal and maximal values of radius of curvature at $p$ of such curve.
For example, a sphere of radius $r$ has  $R_1=R_2 = r$, and $\kappa=\frac{1}{r^2}$. Then the parallel transport angle for some contour $\mathcal{C} = \partial \Sigma$ is, not surprisingly,
\begin{equation}
\Delta\alpha(\mathcal{C}) = \int_\Sigma \kappa dS = \int_\Sigma \frac{dS}{r^2} = \Omega(\mathcal{C}).
\end{equation}

\subsubsection{Curvature and parallel transport on a cone}

\begin{figure}[h]
\begin{center}
\psfrag{C1}{$\mathcal{C}_1$}
\psfrag{gC2}{\textcolor{half-gray}{$\mathcal{C}_2$}}
\psfrag{nC1}{$n(\mathcal{C}_1)$}
\psfrag{gnC2}{\hspace{-0.1cm}\textcolor{half-gray}{$n(\mathcal{C}_2)$}}
\psfrag{n}{$n$}
\includegraphics{img/12Cone}
\caption{Images of contours on a surface of a cone with a rounded tip under the normal map.  \label{Fig:Cone} }
\end{center}
\end{figure}

Now let $\mathcal{B}$ be a conical surface shown in Fig.~\ref{Fig:Cone}.  The tip is rounded, so the field $\vec n$ is smooth everywhere on $\mathcal{B}$. We are interested in the parallel transport angles for contours that lie away form the tip, that is, belong to the lateral surface.  Under the normal map $n$, any point of the lateral surface is mapped to a point in a certain circle of constant latitude (the cone would touch the sphere along this circle). Let $\Omega_0$ be the solid angle enclosed by the circle. Then for any contour on the lateral surface of the cone one has $\Delta\alpha = N\Omega_0$, where $N\in\mathbb{Z}$ counts how many times  the contour encircles the tip. In the figure, $N=1$ for $\mathcal{C}_1$ and $N=0$ for $\mathcal{C}_2$.

The Gaussian curvature $\kappa$ of lateral surface vanishes, since taking the limit $\mathcal{C}\rightarrow p$ requires $N=0$. One can also deduce that $\kappa=0$ from Eq.~\ref{Eq:kRadii}: since the generator of a cone is a straight line, the maximal radius of curvature is $R_2\to\infty$. Whenever a contour lies in the zero-curvature region, the resulting $\Delta\alpha$ is independent of the shape of the contour.

\begin{figure}[h]
\begin{center}
\psfrag{g1}{\raisebox{-0.1cm}{\textcolor{half-gray}{$\vec 1$}}}
\psfrag{gph}{\raisebox{-0.1cm}{\textcolor{half-gray}{$\widetilde\varphi$}}}
\psfrag{gpi}{\textcolor{half-gray}{$\pi$}}
\psfrag{g0}{\raisebox{-0.1cm}{\textcolor{half-gray}{$0$}}}
\psfrag{ph}{$\varphi$}
\psfrag{pi}{$\pi$}
\psfrag{0}{$0$}
\psfrag{be}{\raisebox{-0.1cm}{$\beta$}}
\psfrag{vpt}{$\vec v_{PT}$}
\includegraphics{img/FlatCone}
\caption{Parallel transport on a flattened cone.\label{Fig:FlatCone}}
\end{center}
\end{figure}

Vanishing of the curvature of the cone almost everywhere gives an interesting way to find a parallel transport vector field $\vec v_{PT}$. This time, consider a cone with point-like tip, which touches a sphere along a circle with azimuthal angle $\theta$. The parallel transport on a cone along this circle is the same as in $TS^2$, since the tangent spaces coincide. Now we can cut the cone along its generating line and lay it flat on the plane. Note that this deformation preserves lengths of paths that do not cross the cut: the surface bends without stretching. There is no distortion of the coordinate system (in contrast with maps of the Earth). This deformation is an \defn{isometry}, and preserves lengths of tangent vectors and angles between them.

Let us see what happens after such flattening to the vectors $\vec v_{PT}$ that were parallel transported along the circumference of the cone, as shown in Fig.~\ref{Fig:FlatCone}. Recall from Sec.~\ref{Sec:PTangle} that this vector field is expressed as $\vec v_{PT} = e^{i \alpha}\vec 1$, where $\alpha(\varphi) = - \varphi \cos\theta$, and the basis section is $\vec 1 = \vec e_{\theta}$. On the plane, the shape of the cone is characterized by the opening angle $\beta$. Form the elementary geometry one finds that
\begin{equation}
\cos \theta = 1 - \frac{\beta}{2\pi}.
\end{equation}
Thus the parallel transport equation reads
\begin{equation}
\alpha = - \varphi \biggl( 1- \frac{\beta}{2\pi}\biggr) = - \widetilde \varphi,
\end{equation}
where $\widetilde \varphi$ is the polar angle on the plane. We conclude that after flattening, the parallel transported vectors on the cone become parallel on the plane, in the usual sense. Folding the cone back makes $\vec v_{PT}$ discontinuous at the cut. The parallel transport angle equals the opening angle $\beta$, which therefore measures the amount of curvature concentrated in the tip of the cone.

\section{Geometry of quantum states}\label{Chap:GeomQuant}
In this section, we will discuss several quantum effects whose description includes covariant differentiation or parallel transport on a relevant vector bundle.

\subsection{Quantum particle on a ring}\label{Sec:ParticleRing}
\subsubsection{Spectrum and flux}\label{Sec:ParticleSpec}
Consider a charged particle that moves freely on a ring $S^1$ surrounding a solenoid with magnetic flux $\Phi$. Assume that all of the magnetic field strength is concentrated inside the coil. Then the (angular) momentum operator reads
\begin{equation}
\hat p
= \frac{\hbar}{i} \biggl(\partial_\varphi  - i\frac{q}{\hbar} \emf A_\varphi\biggr).
\end{equation}
Vector potential along the ring must satisfy $\int \emf A_\varphi d\varphi = \Phi$, and we choose it to be constant, $\emf A_\varphi = \frac{\Phi}{2\pi}$.  Let the particle have elementary charge $q=e=|e|$, then
\begin{equation}
\hat p = \frac{\hbar}{i}\partial_\varphi -e \emf A_\varphi = \frac{\hbar}{i}\partial_\varphi -\hbar \frac{\Phi	}{\Phi_0},
\end{equation}
Eigenstate wave function $\psi_n$ and corresponding momentum eigenvalue $p_n$ are
\begin{equation}
\psi_n = e^{in\varphi} \quad \Rightarrow \quad
\hat p \psi_n = \hbar\biggl( n-\frac{\Phi}{\Phi_0}\biggr)\psi_n.
\end{equation}
and the spectrum of the Hamiltonian $\hat H = \frac{\hat p^2}{2m}$ is
\begin{equation}
\varepsilon_n = \frac{p_n^2}{2m}=\frac{\hbar^2}{2m} \biggl( n-\frac{\Phi}{\Phi_0}\biggr)^2.
\end{equation}
The spectrum as a function of $p$ consists of the set of points on the parabola determined by the condition $n\in \mathbb{Z}$, as shown in Fig.~\ref{Fig:Spectrum}.

\begin{figure}[h]
\begin{center}
\psfrag{S1}{$S^1$}
\psfrag{e}{$e$}
\psfrag{ph}{$\varphi$}
\psfrag{P}{$\Phi$}
\psfrag{ep}{$\varepsilon$}
\psfrag{PP0}{$\frac{\Phi}{\Phi_0}$}
\psfrag{phb}{$\frac{p}{\hbar}$}
\psfrag{Pez}{$\Phi=0$}
\psfrag{Pg}{$\Phi\uparrow$}
\psfrag{PeP0}{$\Phi=\Phi_0$}
\psfrag{1}{$1$}
\psfrag{0}{$0$}
\psfrag{-1}{$-1$}
\includegraphics{img/13Spectrum}
\caption{Particle on a ring pierced by magnetic flux and its spectrum for different values of the flux.  \label{Fig:Spectrum} }
\end{center}
\end{figure}

When $\Phi=0$, the momentum of the particle is $n\hbar $, where $n$ is a ``winding number'' of the wave function $\psi_n$, and the wave function determines the physical state of the particle unambiguously. The situation becomes different in the presence of magnetic flux. The variation of the field strength, or \defn{flux insertion},  changes the momentum and energy eigenvalues smoothly while the eigenfunction remains the same. Thus, the presence of magnetic field alters the relationship between the wave function and the state of the particle. As a result, the concept of the winding number as a measure of the momentum is flawed.

Note that one can make a ``large gauge transformation''
\begin{equation}\label{Eq:LargeGT}
\psi'_n = e^{ik\varphi} \psi_n  \qquad \emf A' = \emf A+k\frac{\Phi_0}{2\pi} \qquad \hat p'\psi' = \hat p \psi,
\end{equation}
which simultaneously changes the winding number of the wave function and shifts the flux by an integer number of flux quanta, while keeping the momentum eigenvalue invariant. Strictly speaking, this is not a gauge transformation since it changes the field strength inside the solenoid. However, in a situation when the flux is unknown, one can think of it as a mapping between various ways to describe a given state of the particle.

\subsubsection{Flux insertion as a bundle deformation}\label{Sec:FluxDef}
Let us see what the flux insertion looks like in the geometric picture. We consider a plane bundle $M$ defined by the vector field $\vec m$ over the plane $\mathbb{R}^2$ that contains the ring $S^1$ and intersects a solenoid in a disc $D$. The case of $\Phi=0$ is described by a constant field, say, $\vec m = \vec e_z$. We need to find a bundle configuration that would represent a non-zero flux. That is, the curvature must be  non-zero inside $D$ and vanish outside. Note that the tangent bundle to the cone with rounded tip has exactly this form. We just need to move the field of normals down to the plane. The result is shown in Fig.~\ref{Fig:BundleDef} on the left.

\begin{figure}[h]
\begin{center}
\psfrag{S1}{$S^1$}
\psfrag{gD}{\textcolor{half-gray}{$D$}}
\psfrag{vm}{$\vec m$}
\psfrag{psi0}{$\vec \psi_0 = \vec 1$ }
\psfrag{P}{$\Phi$}
\psfrag{ep}{$\varepsilon$}
\psfrag{P0}{$\Phi_0$}
\psfrag{2P0}{$2\Phi_0$}
\psfrag{p0hb}{\hspace{0.3cm}$\frac{p}{\hbar}$}
\psfrag{-2}{$-2$}
\psfrag{0}{$0$}
\psfrag{-1}{$-1$}
\psfrag{t}{$t$}
\includegraphics[width=0.8\textwidth]{img/14BundleDef}
\caption{Left: deformation of the plane bundle that describes flux insertion. The field $\vec m$ has cylindrical symmetry about vertical axis. Right:  acceleration of the particle. Fibers of the plane bundle along the circle $S^1$ are merged into bands for clarity. Arrows show the section $\vec 1$ that describes the state of the particle.
  \label{Fig:BundleDef} }
\end{center}
\end{figure}

In the process of flux insertion, increasing amount of the curvature contained inside $D$ leads to the rotation of planes outside the solenoid, by continuity of $\vec m$. In this way, the field strength affects the geometry of the bundle in zero-field region. Moreover, this picture incorporates the Faraday's law of induction. Consider the time-dependent bundle defined by $\vec m(r, \varphi, t)$ as a bundle over space-time. Then variation of the curvature in the plane inside the solenoid (magnetic field strength) leads to the non-zero curvature component $f_{\varphi t}$, which describes the strength of the circular electric field.

\subsubsection{Wave function as a section}
Now that we know the evolution of the bundle, we can find the wave function of the eigenstate as its section. First, we need to choose such a basis section for $M$ that connection coefficient corresponds to the constant potential $\emf A_\varphi = \frac{\Phi}{2\pi}$. Consider parts of the bundle $M$ that lie over concentric circles inside the ring $S^1$. For each circle, we can identify the part of the bundle with restriction of $TS^2$ to some circle of latitude. By this identification, we transport the section $\vec 1_N$ of $TS^2$ defined in the Exercise \ref{Exr:1N} to the bundle inside the ring $S^1$. The resulting section $\vec 1$ of $M$ fulfills our requirement. Indeed, by results of the Exercise, connection coefficient $\omega_\varphi$ does not depend on $\varphi$. The section over the region inside $S^1$ is smooth, and the parallel transport angle computed as contour integral $\int_{S^1}(-\omega_\varphi) d\varphi$ equals the solid angle $\Omega(m(S^1))$. Thus, over the ring $S^1$, the connection coefficient with respect to $\vec 1$ is
\begin{equation}
\omega_\varphi = -\frac{\Omega(m(S^1))}{2\pi}.
\end{equation}

Geometrically, the momentum operator is given by the covariant derivative:
\begin{equation}
\hat p = \frac{\hbar}{i}\nabla_\varphi,
\end{equation}
and we are looking for the eigenstate section $\vec \psi_n$ such that
\begin{equation}
\nabla_\varphi \vec \psi_n = i\frac{p_n}{\hbar} \vec \psi_n.
\end{equation}
For $n=0$, we have:
\begin{equation}
\nabla_\varphi \vec \psi_0 = i \biggl(-\frac{\Phi}{\Phi_0}\biggr) \vec \psi_0.
\end{equation}
Now observe that the basis section $\vec 1$ satisfies this equation:
\begin{equation}
\nabla_\varphi \vec 1= i \omega_\varphi \vec 1  = i \biggl( -\frac{\Omega(m(S^1))}{2\pi}\biggr) \vec 1 = i \biggl(-\frac{\Phi}{\Phi_0}\biggr) \vec 1.
\end{equation}
We conclude that for general $n$, the eigenstate section is given by
\begin{equation}
\vec \psi_n = e^{in\varphi} \vec 1.
\end{equation}

The section $\vec \psi_0 = \vec 1$ is shown in the Fig.~\ref{Fig:BundleDef} on the right. When there is no magnetic field, $\Phi=0$, the bundle is flat, and all fibers lie in the horizontal  plane. The particle is at rest, with corresponding section given by a constant vector field. The flux insertion leads to the rotation of the fibers, which inevitably changes the basis section. In the process, the eigenvalue of $\hat p$ changes smoothly. For the values of the flux that are not multiples of the flux quantum $\Phi_0$, there is no eigenstate section with zero momentum. Indeed, it is impossible to define a covariantly constant vector field along a contour on a cone that encloses the tip (cf. right panel of Fig.~\ref{Fig:FlatCone}).

In the final state with $\Phi= 2\Phi_0$, the fibers along $S^1$ again lie in the horizontal plane (with opposite orientation, since $\vec m$ changed sign). However, the basis section winds now twice. In this way, the geometric picture allows one to change the ``covariant winding number'' of a vector field by a \emph{continuous} deformation. Note that during the process, the wave function remains constant, $\psi(\varphi, t) = 1$. Bundle deformation affects the basis section, which changes the state of the particle. In classical terms, this represents the acceleration of the particle under the action of the circular electric field.

\exr {Find a bundle deformation that starts form the same initial state $\vec m = \vec e_z$, but leads to the momentum value $p = \hbar$. \label{Ex:BundleDef}}

\subsection{Aharonov-Bohm effect}
Recall that  the covariant derivative is closely related to the parallel transport. If we treat wave function as a section of a vector bundle, the parallel transport looks like a phase rotation of a state vector characterizing a point-like particle during the motion along some path. At first sight, this is not compatible with quantum physics, since quantum particles are delocalized and do not have a definite trajectory. However, there is a way to formulate quantum mechanics that includes this picture of the parallel transport: it is the Feynman path integral approach. Here we use this framework to discuss effects of the gauge potential on the interference pattern of electrons in the double-slit experiment.
\subsubsection{Path integral formulation}
Suppose that electron is emitted from the source at the point $a$ at the moment $t=0$,  as shown in Fig.~\ref{Fig:AB}. Let us find the probability of detecting electron at the point $b$ of the screen at $t=T$. It is given by the modulus squared of the wave function, $|\psi(b,T)|^2$. This can be obtained from the initial state $|\psi(0)\rangle = |a\rangle$ by applying unitary evolution operator:
\begin{equation}
\psi(b,T) = \langle b|\psi(T)\rangle = \langle b| \hat U(T) |a\rangle.
\end{equation}
The last expression is called the \defn{amplitude} of going form $a$ to $b$ in time $T$. According to Feynman, it is proportional to the sum
\begin{equation}\label{Eq:PathintSum}
\sum_{\mathcal{P}_i} \exp \biggl(\frac{i}{\hbar} S(\mathcal{P}_i)\biggr)
\end{equation}
over all possible paths $\mathcal{P}_i$ connecting points $a$ and $b$. Here $S(\mathcal{P}_i)$ is the classical action along the path $\mathcal{P}_i$. Each path can be thought of as a map from the time interval to the space:
\begin{equation}
\mathcal{P}: [0, T] \to \mathbb{R}^3, \quad \mathcal{P}(0)=a,\quad \mathcal{P}(T)=b.
\end{equation}
It turns out that the ``interference'' of these phases leads to the formation of the fringes in the probability distribution $|\psi(b,T)|^2$ on the screen. Even if the source emits only one electron at a time, there are certain points on the screen that will never be hit by an electron.

\begin{figure}[h]
\begin{center}
\psfrag{P}{$\Phi$}
\psfrag{P1}{\textcolor{half-gray}{$\mathcal{P}_1$}}
\psfrag{P2}{\textcolor{half-gray}{$\mathcal{P}_2$}}
\psfrag{a}{$a$}
\psfrag{b}{$b$}
\includegraphics{img/15AB}
\caption{Left: classical and non-classical trajectories of an electron in double-slit diffraction experiment. Right: shift of the interference pattern in the Aharonov-Bohm effect.
  \label{Fig:AB} }
\end{center}
\end{figure}

\subsubsection{Parallel transport}
Now let us insert between the slits a solenoid with magnetic flux. In presence of magnetic field, one adds the term $L^B= q \dot x_j \emf A_j$ to the Lagrangian, where $j=x,y,z$ and $\emf A_j$ are the vector potential components. Corresponding change in the action reads
\begin{equation}
S^B = q\int \dot x_j \emf A_j dt = q \int_\mathcal{P} \vec{\emf A} \cdot d\vec x
\end{equation}
According to the relation between the vector potential components and connection coefficients,
\begin{equation}
\exp \biggl(q\frac{i}{\hbar} \int_\mathcal{P} \vec{\emf A} \cdot d\vec x  \biggr) =
\exp \biggl(i \int_\mathcal{P} (-\omega_j) dx_j  \biggr) = e^{i\alpha}
\end{equation}
For a given path and a unit vector $\vec 1(0)$ at the starting point, the last expression describes a complex coordinate of the vector $e^{i\alpha}\vec 1$ that is parallel transported along the path.
Here connection coefficients are determined with respect to the basis section $\vec 1$.

The additional ``magnetic'' phase difference between any two paths that are going through the different slits is
\begin{equation}
\Delta \phi_{12}^B  = \int_{\mathcal{P}_2} (-\omega_j) dx_j-\int_{\mathcal{P}_1} (-\omega_j) dx_j = \\= \int_{\mathcal{P}_2-\mathcal{P}_1} (-\omega_j) dx_j = \int_\mathcal{C} (-\omega_j) dx_j
\end{equation}
where $\mathcal{C}$ is a closed contour formed by the difference of the two paths. It follows that
\begin{equation}
\Delta \phi_{12}^B = 2\pi \frac{\Phi}{\Phi_0}
\end{equation}
where $\Phi$ is the magnetic flux through the solenoid. By varying the flux one can continuously shift the diffraction pattern, so that insertion of one flux quantum shifts the pattern by one period. This phenomenon is known as the \defn{Aharonov-Bohm effect}.

\subsection{Geometric phase}
So far, all vector bundles we have looked at were formed by real vector spaces, perhaps described in terms of complex numbers and interpreted quantum mechanically. Now we are going to consider inherently quantum mechanical complex vector bundles. They arise as collections of eigenspaces of a quantum system over a \defn{parameter space}. Each point of this space corresponds to a set of values of some parameters that control  the Hamiltonian of the system. Then variation of parameters can be thought of as traversing a path in the parameter space. Adiabatic evolution of an  eigenstate over a cyclic path will result in the appearance of a phase, similar to the daily rotation angle of the Foucault pendulum. Moreover, both phenomena stem form the restriction of the evolution to some subspace, and we will discuss this analogy in detail.

\subsubsection{Vector bundle over parameter space}
Denote $R$ the space of control parameters of a Hamiltonian $\hat H$ acting on a finite-dimensional space of states $\mathcal H$. The parameter space can assume many shapes: for example, a single real parameter gives $R=\mathbb{R}$, a direction in three-dimensional space is described by a point on  a sphere $R=S^2$ and two periodic variables give $R$ the form of a torus $T^2$.

We attach a copy $\mathcal H_r$ of $\mathcal H$ to each point $r\in R$, forming the \defn{Hilbert space bundle}.  Let $|n_r\rangle\in\mathcal H_r$ be an $n$-th Hamiltonian eigenstate separated by an energy gap from other levels for all combinations of control parameters. At each $r$, the eigenvector determines the eigenspace
\begin{equation}
V^n_r =\{a |n_r\rangle, \, a\in\mathbb{C}\}
\end{equation}
 as a complex one-dimensional subspace of $\mathcal H_r$. Taken together, these subspaces form a complex line bundle over $R$, called an $n$-th \defn{eigenspace bundle} $V^n$.

\begin{figure}[h]
\begin{center}
\psfrag{R}{$R$}
\psfrag{gH}{\textcolor{half-gray}{$\mathcal H_r$}}
\psfrag{r}{$r$}
\psfrag{nr}{\hspace{-0.1cm}$|n_r\rangle$}
\psfrag{Vnr}{$V^n_r$}
\includegraphics{img/16GP}
\caption{Eigenspace bundle over parameter space.
  \label{Fig:GP} }
\end{center}
\end{figure}

\subsubsection{Adiabatic evolution and parallel transport}\label{Sec:BerryPT}
Assume that at the moment $t=0$ the system is prepared in the eigenstate $|n_{r(0)}\rangle$ of the Hamiltonian $\hat H(r(0))$. Let $r(t)$ describe an adiabatically slow variation of control parameters, such that the system remains in the eigenspace $V^n_{r(t)}$ during the process. Note that $|n_r\rangle$ as a complex vector is a specific element of the one-dimensional space $V^n_r$, and $e^{i\alpha}|n_r\rangle$ for some $\alpha\in\mathbb{R}$ describes the same physical state. We wish to consider, after Berry \cite{Berry1984}, the adiabatic evolution of the state in this general form.

Let $e^{i\alpha(t)} |n_{r(t)}\rangle$ be an eigenstate of $\hat H(r(t))$. From the Schr\"odinger equation
\begin{equation}
i\hbar \partial_t (e^{i\alpha} |n\rangle) =\hat H e^{i\alpha} |n\rangle
\end{equation}
one finds
\begin{equation}
\partial_t \alpha = i\langle n|\partial_t n\rangle - \frac{\varepsilon_n}{\hbar}
\end{equation}
The last term describes the ordinary time evolution, or dynamical phase, which can be eliminated by considering $\hat H_n = \hat H - \varepsilon_n\hat {\mathbb{I}}$ or equivalently by setting $\varepsilon_n=0$.

The adiabatic evolution defines a way to transport a state vector $|n(0)\rangle$ through the set of eigenspaces  over a path $r(t)\subset R$. Geometrically, we can use this as a definition of a parallel transport in the bundle $V^n$, and so we define the \defn{Berry connection}. Let us find the form of the covariant derivative associated with this parallel transport. To do this, we compare the present setting with that of the Foucault pendulum.

\begin{center}
\begin{tabular}{| C{4cm} | R{4cm} | L{5cm} | }
\hline
Geometry & Foucault pendulum & Quantum state \\ \hline
Base space & surface of the Earth $S^2$& parameter space $R$  \\
``Ambient bundle'' & $W\rvert_{S^2}$ (see Sec.~\ref{Sec:IdeaVB})  & Hilbert space bundle $\{\mathcal H_r\}$\\
Subspace bundle & tangent planes $TS^2$ & eigenspaces $V^n$ \\
Basis section & $\vec 1$ & $|n\rangle$\\
Parallel transport &$\vec v_{PT} = e^{i\alpha}\vec 1$
 & $e^{i\alpha}|n\rangle$ \\
PT equation & $\partial_\varphi \alpha = -\omega_\varphi$ & $\partial_t \alpha = i\langle n|\partial_t n\rangle$ \\
\hline
\end{tabular}
\end{center}

 It follows that the connection coefficient is $\omega_t = -i\langle n|\partial_t n\rangle$. Negative of this,
\begin{equation}
A_t=i\langle n|\partial_t n\rangle,
\end{equation}
is called \defn{Berry potential}, similarly to a component of the electromagnetic vector potential. The  corresponding parallel transport angle for a closed contour $\mathcal{C}$
\begin{equation}
\Delta \alpha(\mathcal{C}) = \int_\mathcal{C} A_\gamma d\gamma
\end{equation}
is known as the \defn{Berry phase}.

We proceed with the analogy. The definition of $\omega_\gamma$ translates as

\begin{center}
\begin{tabular}{| C{4cm} | R{4cm} | L{5cm} | }
\hline
Connection coefficient & $\nabla_\gamma \vec 1 = i\omega_\gamma \vec 1$  & $\nabla_t |n\rangle = \langle n|\partial_t n\rangle |n\rangle$ \\
\hline
\end{tabular}
\end{center}
The last expression can be rewritten as $|n\rangle\langle n|\partial_t n\rangle$. Introducing projection operator $\Proj^n = |n\rangle\langle n|$ on the $V^n_r$ subspace at each $r$, we find
\begin{center}
\begin{tabular}{| C{4cm} | R{4cm} | L{5cm} | }
\hline
Covariant derivative &
$\nabla_\gamma \vec 1 = \Proj (\partial ^{\,a}_\gamma \vec 1)$
 &
 $\nabla_t |n\rangle = \Proj^n |\partial_t n\rangle$\\
 Derivative in the ambient space & $\partial_\gamma^{\,a} \vec v = \vec e_i^{\,a} \partial_\gamma v_i^{\,a}$ &  $|\partial_t n\rangle = \sum_\alpha |\alpha \rangle \partial_t n_\alpha $ \\
Constant basis & $\{\vec e^{\,a}_x, \vec e^{\,a}_y, \vec e^{\,a}_z\}$ in $W$ & $\{|\alpha\rangle\}$ in $\mathcal H_r$\\ \hline
\end{tabular}
\end{center}
We conclude that the evolution of the quantum state, subject to the combination of the Schr\"odinger equation and the adiabatic theorem, is mathematically equivalent to the rotation of the Foucault pendulum plane (up to the quantum dynamical phase). In both cases the parallel transport is defined by projection and depends only on geometry of certain subspaces.

There is, however, an important difference. In the case of the Foucault pendulum, the conservation law of the angular momentum gives us the parallel transport in the ambient space $W$, the corresponding constant basis, and the derivative $\partial^{\,a}$. Then we obtain the parallel transport in the subspaces by demanding that the projection of this derivative vanishes. But in the case of the Hilbert space bundle $\{\mathcal H_r\}$, there is no analogue of the ``conservation law''. Thus there is no preferred choice of the constant basis. Later we will see how different basis choices in $\{\mathcal H_r\}$ lead to different geometric phases for the same eigenspace bundle $V^n$.

\subsubsection{Eigenspace bundle for two-level quantum system}\label{Sec:Two-level}
Consider a quantum system with two-dimensional Hilbert space $\mathcal H$, which can be identified with $\mathbb{C}^2$ once a basis $\{|\spinup\rangle, |\spindown\rangle\}$ is chosen. A general form of the Hamiltonian matrix is
\begin{equation}
H=h_0 \mathbb{I} + \sum_\alpha h_\alpha \sigma_\alpha = h_0 \mathbb{I} + H_{\vec h},
\end{equation}
where $\sigma_\alpha$ are Pauli matrices, which serve as basis elements in the space of traceless Hermitian $2\times 2$ matrices:
\begin{equation}
\sigma_x = \begin{pmatrix}
0&1\\1&0
\end{pmatrix}
\quad
\sigma_y = \begin{pmatrix}
0& -i \\ i&0
\end{pmatrix}
\quad
\sigma_z = \begin{pmatrix}
1&0 \\0&-1
\end{pmatrix}.
\end{equation}

The eigenvalues  of $H$ are
\begin{equation}\label{Eq:TwoLevelEn}
\varepsilon_\pm = h_0 \pm |\vec h|.
\end{equation}
We are mostly interested in eigenstates, so we set $h_0=0$. Then an eigenstate is determined by the direction of $\vec h = (h_x, h_y, h_z)^T$, which is well-defined everywhere except the origin, and we have the space of parameters $R = \mathbb{R}^3 \setminus \{0\}$.

Following the general procedure described above, we attach a copy of $\mathcal H$ with a specified basis $\{|\spinup\rangle, |\spindown\rangle\}$ to each point  $\vec h\in R$ and choose a high-energy state eigenvector $|\psi_{\vec h+}\rangle\in \mathcal H$, so that $H_{\vec h}\psi_{\vec h+}=\psi_{\vec h+}$. Then we have a subspace $V^+_{\vec h}=\{c\psi_{\vec h+} | \, c\in \mathbb{C}\}\subset \mathbb{C}^2$ attached to each point $\vec h$  of the parameter space $R$ . These eigenspaces form a complex line bundle over $R$, which we will call the \defn{monopole bundle} and  denote $D_3^+$. Sometimes it is convenient to consider the restriction of this bundle to the unit sphere $S^2\subset R$:
\begin{equation}
D^+= D^+_3\rvert_{S^2}
\end{equation}
 As a base space for $D^+$, this sphere is known as the \defn{Bloch sphere} in the context of spin-$\frac{1}{2}$ particle in the magnetic field, or as the \defn{Poincar\'e sphere} when it is used  to describe states of circularly polarized light.

Now we need to find the eigenstates $\psi_{\vec h+}$ on the sphere. To this end, define the unitary matrix $U(\theta, \vec w)$ as
\begin{equation}
U(\theta, \vec w) = \cos\frac{\theta}{2} \mathbb{I} +
\sin\frac{\theta}{2}\sum_\alpha\frac{\sigma_\alpha}{i}w_\alpha,
\end{equation}
where $\vec w$ is a three-dimensional vector. Under conjugation by $U(\theta, \vec w)$, the Hamiltonian $H_{\vec h}$ is rotated in the parameter space $R$:
\begin{equation}
U(\theta, \vec w) H_{\vec h} U(\theta, \vec w)^{-1} = H_{\mathcal{R}(\theta, \vec w)\vec h},
\end{equation}
where $\mathcal{R}(\theta, \vec w)$ is $\theta$ rotation about $\vec w$ axis.  If $\psi_{\vec h}$ is an eigenstate of $H_{\vec h}$, then the same-energy eigenstate of the rotated Hamiltonian reads
\begin{equation}
\psi_{\mathcal{R}(\theta, \vec w)\vec h} = U(\theta, \vec w)\psi_{\vec h}.
\end{equation}

Let us apply this unitary rotation to the eigenstate $\psi=(1,0)^T$ of $H_{\vec h}=\sigma_z$, so $\vec h=(0,0,1)^T$, and choose $\vec w = (-\sin\varphi, \cos\varphi,0)^T$. Then the $\theta$ rotation of $\vec h$ about $\vec w$ axis  gives the point with coordinates $(\theta, \varphi)$ on the sphere. The corresponding eigenstate is
\begin{equation}\label{Eq:PsiPlus}
\psi_+(\theta, \varphi) = \begin{pmatrix}
\cos\frac{\theta}{2}\\
e^{i\varphi}\sin\frac{\theta}{2}\\
\end{pmatrix}
\end{equation}
Note that this expression is conceptually similar to the column of components of $\vec e_\theta$ as a three-dimensional vector, Eq.~(\ref{Eq:Basis3D}).

Thus, we have defined the eigenspace bundle for a two-level quantum system and specified a basis section $|\psi_+\rangle$.

\subsubsection{Geometry of eigenspace bundles}
Now we consider geometry of the monopole bundle with the Berry connection. Connection coefficients are
\begin{equation}
\omega_\varphi = -i\langle \psi_+|\partial_\varphi \psi_+\rangle = \sin^2\frac{\theta}{2}
\end{equation}
and $\omega_\theta=0$. Since the eigenstate is unchanged in the radial direction, $\omega_h=0$.

From the connection coefficients on the parameter space $R$ one can find the \defn{Berry curvature} component for each pair of coordinates of $R$. For the bundle $D^+_3$ over $\mathbb{R}^3\setminus \{0\}$, we have
\begin{equation}
f_{\theta \phi}=-\frac{1}{2}\sin \theta, \qquad
f_{\theta h} = f_{\varphi h} =0.
\end{equation}
Thus the Berry phase for a boundary contour $\partial\Sigma$ of  a surface  $\Sigma\subset R$  is
\begin{equation}\label{Eq:DalphaBerry}
\Delta\alpha(\partial\Sigma) = -\frac{1}{2}\Omega(\Sigma)
\end{equation}
where $\Omega$ is a solid angle enclosed by $\Sigma$ in the space of rays $\vec h$.

\exr{Employ the unitary transformation $U(\theta, \vec w)$ to find the negative-energy eigenstates $\psi_-$. Calculate connection coefficients and curvature components for the  corresponding bundle $D^-_3$.}\label{Ex:PsiMinus}

Let us find the curvature of $D^+_3$ in terms of Cartesian coordinates $h_\alpha = h_x, h_y, h_z$. The solid angle spanned by an infinitesimal coordinate rectangle $dS_{\alpha \beta}$ with radius-vector $\vec h$ is given by
\begin{equation}
d\Omega_{\alpha\beta} = \cos \nu \frac{dS_{\alpha\beta}}{h^2},
\end{equation}
where $h=|\vec h|$, and $\nu$ is the angle between the normal to the rectangle and $\vec h$. Since the element $dS_{\alpha \beta}$ is orthogonal to the third coordinate $h_\gamma$, the cosine is given by $\frac{h_{\gamma}}{h}$. Using Levi-Civita anti-symmetric symbol to take the orientation of the rectangle into account, we have
\begin{equation}
f_{\alpha \beta} = -\frac{1}{2} \frac{\varepsilon_{\alpha\beta\gamma} h_\gamma}{h^3}.
\end{equation}

Note that the situation here is essentially the same as for the curvature of real plane bundles discussed in Sec.~\ref{Sec:FluxGeom}. We have a standard vector bundle with a known curvature, and describe any other bundle by a map to the base space of the standard one. Thus we can  find the Berry curvature of a general vector bundle formed by eigenspaces of a two-level Hamiltonian over some parameter space $R$, as follows.  Consider a surface $\mathcal{B} \subset R$ with coordinates $(x_1, x_2)$ and let $\vec h(x_1, x_2)$ be a vector field that defines the Hamiltonian at each point of $\mathcal{B}$. We wish to find the curvature of the eigenspace bundle $V^+$ over $\mathcal{B}$. As we did in Sec.~\ref{Sec:CurvJacobian}, we consider the vector field $\vec h$  as a map from $\mathcal{B}$ to $\mathbb{R}^3\setminus\{0\}$, and interpret it as a change of coordinates. The result of Exercise~\ref{Ex:f12} can be written as
\begin{equation}
f_{12}=(\partial_1 \alpha)(\partial_2 \beta) f_{\alpha\beta},
\end{equation}
with $f_{\alpha\beta}$ understood now as curvature components of $D^+_3$.  Coordinates $\alpha$ are now denoted as $h_\alpha$, so the curvature component is
\begin{equation}
f_{12} = -\frac{1}{2} \frac{\varepsilon_{\alpha\beta\gamma} h_\gamma}{h^3} (\partial_1 h_\alpha) (\partial_2 h_\beta) = -\frac{1}{2h^3} \vec h\cdot [\partial_1 \vec h \times \partial_2 \vec h].
\end{equation}

\subsubsection{Dirac monopole and emergent electrodynamics}
Finally, let us interpret the monopole bundle $D^+_3$ over $\mathbb{R}^3\setminus \{0\}$ with Berry connection as a description of the magnetic field. Consider a region $\Sigma \subset S^2$ of a sphere of radius $r$. From Eq.~\eqref{Eq:DalphaBerry}, the magnetic flux through $\Sigma$ is given by
\begin{equation}
\Phi(\Sigma) = -\Phi_0 \frac{\Omega(\Sigma)}{4\pi}.
\end{equation}
Thus the total flux through the sphere is $\Phi(S^2) = - \Phi_0$, which means that the sphere contains a source of the magnetic field. We know that the only non-zero curvature component is $f_{\theta \varphi}$, which gives the radial component of the magnetic field strength. In the limit of region $\Sigma$ shrinking to a point $p$ on the sphere,
\begin{equation}
\emf B_r = \lim_{\Sigma \to p} \frac{\Phi(\Sigma)}{S(\Sigma)} = - \frac{1}{r^2} \frac{\Phi_0}{4\pi}  =
-\frac{1}{r^2} \frac{\hbar }{2e}.
\end{equation}
This is the field configuration of the \defn{Dirac monopole} with magnetic charge $-\frac{\hbar }{2e}$, placed at the origin.

In some situations, this electromagnetic interpretation of the Berry connection is taken literally by electrons themselves. Consider a \defn{skyrmion}, a vortex-like magnetic texture of magnetization $\vec m(x,y)$ on a plane, shown in Fig.~\ref{Fig:SkED} on the left. We assume that the magnetization has constant modulus $|\vec m|$. Such structure is said to have a topological charge, since the image of the plane under the map $m: \mathbb{R}^2\to S^2 $ defined by magnetization vector $\vec m$ covers the whole sphere. But for now we are interested in the local geometric effects of the spatial inhomogeneity of the magnetization.

\begin{figure}[h]
\begin{center}
\psfrag{Sk}{Sk}
\psfrag{gP}{\textcolor{half-gray}{$\mathcal{P}$}}
\psfrag{gmP}{\textcolor{half-gray}{$m(\mathcal{P})$}}
\psfrag{vm}{$\vec m$}
\psfrag{10}{$\begin{psmallmatrix}
1\\0
\end{psmallmatrix}$}
\includegraphics{img/17SkED}
\caption{Electron moving through a skyrmion feels the emergent magnetic field.
  \label{Fig:SkED} }
\end{center}
\end{figure}

Suppose that the skyrmion exists in a magnetic metal, and conductance electrons are free to flow through the structure. We further assume  that the motion is semi-classical and adiabatic, in the following sense. First, electron is localized in a region with approximately constant  direction of $\vec m$. Second, electron's spin is always aligned with the magnetization. Under these assumptions, any trajectory traced out by the electron in the plane is mapped to the space of spin states, or the Bloch sphere. So the evolution of the spin state is equivalent to the Berry phase problem of the spin in a slowly rotating field.

Now let us look at the motion of electron form the path integral perspective. Recall that the sum of phases along all possible path appears as a result of decomposition of the the unitary evolution operator.
Here, such evolution splits into a real-space part and a spin-space part. Thus, besides the phase \eqref{Eq:PathintSum}, each trajectory will be weighted by the Berry phase. But regardless of the origin of the phase shift, it leads to the real Lorentz force felt by electrons. In a similar way, time-dependent magnetic structures can create effective electric field. These effects are manifestations of \defn{emergent electrodynamics} \cite{EverschorSitte2014}.

\exr{Estimate the strength of the emergent magnetic field, given that the typical diameter of a skyrmion is $100~\text{nm}$.}

\section{Topology of vector bundles}\label{Chap:TopVB}
Let $M$ be a complex line bundle with a closed surface as a base space. One can associate to $M$ an integer $c(M)$, which will stay invariant under any smooth deformations of $M$. It turns out that the curvature of a connection gives a way to calculate $c(M)$, despite that locally, deformations will change the values of curvature. In physics, the integer invariant $c(M)$ is the source of topological quantization, which restricts possible values of magnetic charge of the hypothetical Dirac monopoles and provides the basic classification of topological phases in condensed matter.

\subsection{Trivial and non-trivial bundles}
\subsubsection{What is topology about?}
\defn{Topology} is a branch of mathematics studying global properties of objects that do not change under smooth deformations. In other words, such properties must be independent of local geometric details, like angles and distances. For example, one can define a circle as a one-dimensional space, such that a sufficiently long journey in any direction will bring you to the starting point. This definition is based on a global property and  does not require any local measurements. Such ``circle'' can assume many shapes, and even be represented as an interval with  identified endpoints.

This topological point of view might seem to contain too much freedom to have any physical application. Consider, however, an electron in a periodic potential of a one-dimensional crystal. States of the electron are labeled by crystal momentum $k$. Periodicity of the potential makes the states at $k$ and $k+G$ physically equivalent, where $G$ is a reciprocal lattice vector. Therefore the Brillouin zone formed by the physically distinct values of $k$ is the interval $[-G/2, G/2]$ with endpoints identified, that is, a circle.

Another essentially topological statement is the intermediate value theorem and its consequences. In particular, any continuous function that takes positive value at one end of an interval and a negative value at the other must have at least one root on the interval (or an odd number of them). As we will see in a moment, this is closely related to the topology of real line bundles.

For complex line bundles we studied above, examples of local measurements are covariant derivative and curvature, subjects of differential geometry. And what are the global properties? First, we need to specify what we mean by a smooth deformation of a vector bundle. Of course, such deformations must preserve vector bundle structure. For a bundle defined by a vector field as described in Sec.~\ref{Sec:CLBundles}, one can imagine smooth deformations of both the base space and the vector field. Soon we will be able to compute a topological number which stays invariant under such deformations.

\subsubsection{Topology of line bundles}
Let us first consider a real line bundle $V$ over the circle $S^1$, constructed as follows. We attach a real line to each point of the circle in such a way that the lines form the surface of the cylinder, as shown in Fig.~\ref{Fig:RTopoBun}. If we think of the circle as an interval with identification, we also need to identify the fibers over endpoints, which is indicated by gray arrows.

\begin{figure}[h]
\begin{center}
\psfrag{V}{$V$}
\psfrag{M}{$M$}
\psfrag{p}{$p$}
\psfrag{Vp}{$V_p$}
\psfrag{S1}{$S^1$}
\psfrag{ve}{$\vec e$}
\includegraphics{img/18RTopo}
\caption{Cylinder and M\"obius band as examples of  trivial and non-trivial real line bundles over the circle.
  \label{Fig:RTopoBun} }
\end{center}
\end{figure}

The choice of a basis vector in a fiber allows us to represent it as a real line:
\begin{equation}
V_p \xrightarrow[\text{basis } \vec e(p)]{\text{choice of}} \mathbb{R}.
\end{equation}
Globally, this yields an identification of the bundle with a product space:
\begin{equation}
V \xrightarrow[\text{section } \vec e]{\text{choice of}} S^1 \times \mathbb{R},
\end{equation}
once a non-vanishing basis section $\vec e$ is chosen.  Note also that in this case any other section $\vec v = v \vec e$ is described by a real function
\begin{equation}
v: S^1 \to \mathbb{R}.
\end{equation}
Bundles that can be identified with a product space by a choice of global non-vanishing section are called \defn{trivial}.

Now suppose that we reverse one of the fibers over the endpoints in the planar picture, and obtain a new bundle $M$. In three dimensions, this gives the bundle a shape of the M\"obius band (of infinite width). This twisted identification forces any smooth section to pass through zero, and makes the bundle \defn{non-trivial}. Such  bundle cannot be represented as a product space, and its sections cannot be described as functions on the base space.

Note that \emph{locally} bundles $V$ and $M$ look similar. By choosing a non-vanishing section along an interval $I\subset S^1$, one can identify the part of either bundle over $I$ with the product $I\times \mathbb{R}$. The difference between bundles $V$ and $M$ can only be detected if we look at the whole circle $S^1$. Also, one cannot transform one bundle into another by a smooth deformation. These observations indicate that (non-)triviality is a global topological property of a vector bundle.

\begin{figure}[h]
\begin{center}
\psfrag{TT2}{$TT^2$}
\psfrag{TS2}{$TS^2$}
\psfrag{Tp}{$T_pT^2$}
\psfrag{vs}{$\vec 1_N$}
\psfrag{gvv}{\textcolor{half-gray}{$\vec v$}}
\includegraphics{img/18CTopo}
\caption{Tangent bundles to the torus and the sphere as examples of  trivial and non-trivial complex line bundles.
  \label{Fig:CTopoBun} }
\end{center}
\end{figure}

In a similar way, a complex line bundle over base $\mathcal{B}$ can be identified with $\mathcal{B}\times \mathbb{C}$ if there exists a smooth non-vanishing section $\vec 1$ over all of $\mathcal{B}$.  If no such section exists, the bundle is non-trivial. More generally, a bundle with $n$-dimensional complex vector space as a fiber is said to be trivial if it can identified with a product $\mathcal{B}\times \mathbb{C}^n$. To make such an identification, one needs $n$ global linearly independent sections.

Consider the tangent bundle to the torus as a complex line bundle (Fig.~\ref{Fig:CTopoBun}). Since there is a smooth non-vanishing field of tangent vectors that goes along one of generating circles, the bundle is trivial. The tangent bundle to the sphere $TS^2$ has different nature. Its section $\vec 1 = \vec e_\theta$ has two singularities at the poles. One can try to construct a smooth section $\vec 1_N$ by transporting some vector from the north pole along the meridians. The resulting section will be smooth at the north pole, but there will be a singularity at the south pole, which winds twice around it, as shown in Fig.~\ref{Fig:CTopoBun}. We will see that any section of $TS^2$ must have singularities, and the bundle is topologically non-trivial.

Thus, the global twisting of the bundle space leads to the appearance of singularities in its sections. In all our examples, both real and complex, the singularities were zero-dimensional, or point-like. Generically, this happens when the dimensionality of the base space coincides with that of the fiber (as a real space). We wish to use these point-like singularities to characterize the non-trivial topology of a vector bundle. Note that the base space must also be closed, since otherwise the singularities can be pushed away through the boundary, and the bundle is always trivial. We are interested in topology of complex line bundles, so the base spaces will be closed two-dimensional surfaces.

\subsubsection{Index of a singularity}
We start with a formal definition of the winding number of a singularity. Suppose that a section $\vec s$ of a complex line bundle  has a \defn{singularity} at $p$. This means that $\vec s(p)=0$ and $\vec s$ is non-zero in some neighborhood of $p$.  Let $\vec v$ be some smooth non-vanishing section near $p$. Choose a positively-oriented contour $\mathcal{C}$ near $p$ such that both sections are non-zero on it.
Then along the contour two sections are related by $\vec s = r e^{i\beta}\vec v$, where $r$ and $\beta$ are  smooth functions of the coordinate $\gamma$ on $\mathcal{C}$. We define the \defn{index of a singularity} at $p$ as
\begin{equation}
\ind [\vec s(p)] = \frac{1}{2\pi}\int_\mathcal{C} (\partial_\gamma \beta)d\gamma.
\end{equation}
For example, the singularity of $\vec 1_N$ shown on the right in Fig.~\ref{Fig:CTopoBun}, has index $+2$.

The index is an integer independent of the choice of the smooth section $\vec v$ (compare with the low-dimensional example in Sec.~\ref{Sec:FluxOmega}). Note that the sign of the index depends on the orientation of the surface, which determines that of the contour. There is also an implicit dependence on the orientation of the fiber, encoded in the complex structure.

\subsection{Chern number}
\subsubsection{Singularities and integral of curvature}
Consider a complex line bundle $M$ over a closed two-dimensional surface $\mathcal{B}$. Let $\vec 1$ be a basis section with singularities at points $p_i\in \mathcal{B}$. Let us try to find the integral of curvature $f_{12}$ over $\mathcal{B}$ using the Stokes' theorem. However, the theorem relates the integral of curvature with the integral of connection coefficient along the boundary, and also requires a smooth basis section. In our case $\mathcal{B}$ has no boundary and the section $\vec 1$ is not smooth. We can mend the situation by cutting a hole around each singularity.

\begin{figure}[h]
\begin{center}
\psfrag{M}{$M$}
\psfrag{Sig}{$\Sigma$}
\psfrag{pi}{$p_i$}
\psfrag{B}{$\mathcal{B}$}
\psfrag{Di}{$D_i$}
\psfrag{dDi}{$\partial D_i$}
\includegraphics{img/19Chern}
\caption{Removing disks at the points with singularities from the surface $\mathcal{B}$.
  \label{Fig:Chern} }
\end{center}
\end{figure}

Let $\{D_i\}$ be a set of small discs near points $\{p_i\}$. Consider a region $\Sigma\subset \mathcal{B}$ that is obtained by removing these discs from $\mathcal{B}$:
\begin{equation}
\Sigma = \mathcal{B} \setminus \cup_i D_i.
\end{equation}
The section $\vec 1$ is smooth over $\Sigma$. The  boundary of $\Sigma$ is given by
\begin{equation}
\partial \Sigma = \cup_i (-\partial D_i),
\end{equation}
since the boundary of each hole has opposite orientation to that of the boundary of the disc. Now we apply the Stokes' theorem to the region $\Sigma$:
\begin{equation}\label{Eq:StokesDiscs}
\frac{1}{2\pi}\int_\Sigma f_{12} dx_1dx_2 = \frac{1}{2\pi} \sum_i \int_{-\partial D_i} (-\omega_{\gamma_i}) d\gamma_i,
\end{equation}
where $\gamma_i$ parameterizes $\partial D_i$.

Consider one of the integrals from the sum. Let $\vec 1'$ be a smooth section over the disc $D_i$. The singular section can be expressed as $\vec 1 = e^{i\beta(\gamma_i)}\vec 1'$ along $\partial D_i$. Using the transformation rule of connection coefficient (\ref{Eq:ConnBasis}), one finds:
\begin{equation}
\int\limits_{\partial D_i} \omega_{\gamma_i} d\gamma_i =
\int\limits_{\partial D_i} \omega'_{\gamma_i} d\gamma_i
+\int\limits_{\partial D_i} (\partial_{\gamma_i} \beta) d\gamma_i = -\int\limits_{D_i} f_{12} dx_1 dx_2 + 2\pi \ind[\vec 1(p_i)].
 \end{equation}

Finally, we take the limit $D_i\rightarrow p_i$. Then integrals of curvature over discs vanish $\int_D f \to 0$, while the region $\Sigma$ approaches the whole surface $\mathcal{B}$. Thus the limit of Eq.~(\ref{Eq:StokesDiscs}) is
\begin{equation}\label{Eq:ChernStokes}
\frac{1}{2\pi}\int_\mathcal{B} f_{12}dx_2dx_2 = \sum_i \ind [\vec 1(p_i)].
\end{equation}
This is a remarkable result. Left-hand side is the integral of a local geometric quantity that does not depend on the choice of a basis section $\vec 1$, while on the right we have a sum of integers. It follows that
\begin{enumerate}
\item For any section of $M$ with point-like singularities, the sum of their indices is the same.
\item The integral of curvature over $\mathcal{B}$ is an integer multiple of $2\pi$.
\end{enumerate}
Moreover, consider a smooth deformation of the bundle $M$. Locally the values of curvature will change, and the integral will vary smoothly. But any smooth integer-valued function is constant. It follows that the integer is a topological invariant of a complex line bundle. It is called the \defn{Chern number} of the bundle, and is denoted $c(M)$.

Another important consequence of Eq.~\eqref{Eq:ChernStokes} is that the integral of  curvature on the left does not depend on the choice of connection. Recall that connection is an external structure on a vector bundle. It is not unique, and different connections may have different curvatures. In contrast, the sum of indices on the right is an intrinsic property of the vector bundle. Thus, the integral of the curvature of a connection has the same value for any connection on $M$.

\subsubsection{Chern numbers of familiar bundles}
\newpar{Tangent bundle $TS^2$.}
The curvature of the projected connection is $f_{\theta \phi}=\sin \theta$. The integral of the curvature gives
\begin{equation}
 \frac{1}{2\pi} \int_{S^2} d\Omega = 2.
\end{equation}
The section $\vec 1=\vec e_\theta$ has two singularities on the poles, each with index $+1$:
\begin{equation}
\ind[\vec 1(S)] + \ind[\vec 1 (N)] = 2.
\end{equation}
The section $\vec 1_N$ has only one singularity on the south pole, but
\begin{equation}
\ind [\vec 1_N (S)] =2.
\end{equation}
We conclude that the Chern number is $c(TS^2)=2$. Note that the different orientation of the sphere would change signs of the curvature, indices, and Chern number $c(TS^2)$. On the other hand, we could define $i\vec e_\theta$ to be $-\vec e_\phi$ instead of $\vec e_\phi$, resulting in the conjugation of the complex coordinate and reversing the sign of the Chern number.

\newpar{Monopole bundle and topological quantization.}
For the eigenstate bundle $D^+$ over the Bloch sphere described in Sec.~\ref{Sec:Two-level}, one has
\begin{equation}
c(D^+) = \frac{1}{2\pi}\int_{S^2} -\frac{1}{2}d\Omega = -1.
\end{equation}
To obtain this in terms of singularities, consider the section given by Eq.~(\ref{Eq:PsiPlus}). At the north pole, $\theta=0$, the section is smooth. Near the south pole it has the form
\begin{equation}
\lim_{\theta\to \pi} \psi_+ (\theta, \varphi) = \begin{pmatrix}
0 \\ e^{i\varphi}\\
\end{pmatrix}.
\end{equation}
For a positively-oriented contour near the south pole, $\varphi$ increases in negative direction. It follows that the singularity has index $-1$, as expected.

Recall that in terms of electromagnetic field, $D^+_3$ over $\mathbb{R}^3\setminus \{0\}$ describes the magnetic field configuration of the Dirac monopole placed at the origin. For any closed surface $\mathcal{B}\subset \mathbb{R}^3$ that encloses the origin, the magnetic flux through $\mathcal{B}$ equals $-\Phi_0$, and we have $c(D^+_3\rvert_\mathcal{B})=-1$. Thus, a wave function $\vec \psi$ of an electron in the vicinity of the monopole has a singularity on $\mathcal{B}$.  By smoothness of $\vec \psi$, these singularities form a \defn{nodal line} of zeroes of wave function. The same is true for the basis section $\vec 1$. This leads to the line of singularities of the vector potential, known as the \defn{Dirac string}. Topology also puts constraints on the possible values of the  monopole charge. In terms of magnetic field, Eq.~\eqref{Eq:ChernStokes} tells that the magnetic flux through any closed surface is quantized in the units of $\Phi_0$. In other words, any magnetic monopole must have magnetic charge $n\frac{\hbar}{2e}$, where $n\in\mathbb{Z}$. This is known as the \defn{Dirac quantization condition}.

Note that if a complex line bundle has an odd Chern number, it is impossible to model it with real planes. By continuity of $\vec m$, the total solid angle $\Omega(m(\mathcal{B}))$ must be a multiple of $4\pi$, and the same is true for the integral of curvature over $\mathcal{B}$. Thus the Chern number of a complex line bundle constructed from real planes, as described in Sec.~\ref{Sec:CLBundles}, must be even.

\newpar{Gaussian curvature of the torus.}
Now consider the tangent bundle to the torus $T^2$. We already know that this bundle is trivial, so its Chern number is zero. Thus the integral of curvature of any connection must vanish:
\begin{equation}
 \int_{T^2} f_{12}dx_2dx_2=0.
 \end{equation}
In particular, we have
 \begin{equation}
 \int_{T^2} \kappa dS = 0
 \end{equation}
for the Gaussian curvature. Moreover, this will be true for the tangent bundle to any closed surface obtained from the torus by a smooth deformation. Thus, whenever a surface has a single hole, the integral of the curvature will vanish. Such surfaces are said to be of \defn{genus} $g=1$ (for the sphere $S^2$, $g=0$). All closed orientable two-dimensional surfaces are classified by the genus, up to a smooth deformation. From the results for the sphere and the torus, we have:
\begin{equation}
\frac{1}{2\pi}\int_{\mathcal{B}_0} \kappa dS = 2, \qquad
\frac{1}{2\pi}\int_{\mathcal{B}_1} \kappa dS =0,
\end{equation}
where $\mathcal{B}_g$ denotes a genus $g$ surface. Later we will find the Chern number of the tangent bundle $T\mathcal{B}_g$ for all possible $g$.

\subsection{Pullback construction and topology}
Recall that both real plane bundle $M$ and eigenspace bundle for two-level quantum system $V$ can be defined by a vector field. For the former it is the field of plane normals $\vec m$, while for the latter it is the field $\vec h$ that encodes the Hamiltonian by a  point in the space of Pauli matrices. In both cases, the vector field on the base space $\mathcal{B}$ gives rise to the map from $\mathcal{B}$ to the sphere $S^2$. Here we consider a general way to define vector bundles by a map to the base space of some other bundle.
\subsubsection{Bundles induced by a map}

Let $m: \mathcal{B} \to \mathcal{S}$ be a smooth map between surfaces, and let $E$ be a vector bundle over $\mathcal{S}$. We wish to define from this data a vector bundle over the surface $\mathcal{B}$. This can be done as follows: at each point $p\in \mathcal{B}$, define the fiber of the new bundle as
\begin{equation}
(m^*E)_p= E_{m(p)}.
\end{equation}
In this way, a vector space is associated with every point of $\mathcal{B}$. The bundle $m^* E$ is called \defn{pullback} bundle, or \defn{induced} bundle. In a similar way one defines a \defn{pullback section} $m^*\vec 1$ of $m^* E$:
\begin{equation}
(m^*\vec 1) (p) = \vec 1 (m(p)),
\end{equation}
where $\vec 1$ is some section of $E$. Fibers of $m^*E$ at different points are assumed to be independent\footnote{See Ref.~\cite{Hatcher2003} for the precise description of the pullback construction.}, so that there are sections of the pullback bundle that are not pullback sections.

For a simple example,  consider the pullback of $TS^2$ under the constant map $m_0: S^2\to S^2$, which sends the whole sphere to a single point of itself. The result is a trivial bundle $m_0^*TS^2$ over $S^2$.

\subsubsection{Curvature of pullback bundle and degree of a map}
Recall that in Sec.~\ref{Sec:CurvJacobian} we computed the curvature of the real plane bundle $M$. The bundle was described by a vector field of the plane normals $\vec m$, or equivalently by the map $m: \mathcal{B} \to S^2$. To find the curvature, we used the identification of the vector spaces $M_p$ and $T_pS^2$ as subspaces of $\mathbb{R}^3$. The pullback bundle is a more abstract version of the same construction, with the map $m$ playing now a central role.  The whole argument can be rephrased in terms of pullbacks, leading to the same expression for the curvature:
\begin{equation}
(m^*f)_{12} = J f_{\alpha \beta},
\end{equation}
where $m^*f$ stands for the curvature of the pullback bundle, $J$ is Jacobian determinant, and $f_{\alpha \beta}$ is the curvature of $E$ over $\mathcal{S}$ in coordinates $(x_\alpha, x_\beta)$.

We are interested in topology of induced bundles, and from now on we assume that surfaces $\mathcal{B}$ and
$\mathcal{S}$ are closed.  Similarly to  Eq.~(\ref{Eq:IntCurvJ}), we obtain for  the Chern number of $m^*E$:
\begin{equation}
c(m^*E) = \frac{1}{2\pi}\int_{m(\mathcal{B})} f_{\alpha\beta} dx_\alpha dx_\beta = c(E)\degree(m),
\end{equation}
where \defn{degree of a map} $\degree(m)$ is an integer that shows how many times the image of $\mathcal{B}$ under $m$ covers the sphere $\mathcal{S}$. The integrality of the degree follows from the fact that continuous maps preserve boundaries. Since $\partial \mathcal{B}=\varnothing$, the image $m(\mathcal{B})$ also does not have a boundary and covers all of $\mathcal{S}$, perhaps, several times. Let us consider some examples.

\newpar{Pullbacks of $TS^2$ under reflection and inversion.} Let $\sigma_x: S^2 \to S^2$ be the reflection in the $yz$ plane (the sphere is centered at the origin). This map preserves areas of regions on the sphere but reverses orientation of contours. Thus the curvature of the tangent bundle $TS^2$ changes its sign under pullback by $\sigma_x$, and we have\footnote{The bundle $\sigma_x^*TS^2$ may be useful in the context of Ex.~\ref{Ex:BundleDef}.}
\begin{equation}
c(\sigma_x^* TS^2) = -2.
\end{equation}
In a similar fashion, the inversion map $\mathcal{I}$ that sends $\vec r \to -\vec r$, gives
\begin{equation}
c(\mathcal{I}^*TS^2)=-2.
\end{equation}

\newpar{Pullback of $D^+$ over the Bloch sphere under inversion.} By the same token as above, the Chern number is
\begin{equation}
c(\mathcal{I}^*D^+) = 1.
\end{equation}
This bundle has a simple physical interpretation. We attach a vector space defined by $|\psi_{\vec h+}\rangle$ to the point $-\vec h$. Note that
\begin{equation}
\hat H_{-\vec h} |\psi_{\vec h +}\rangle= -\hat H_{\vec h} |\psi_{\vec h +}\rangle = -|\psi_{\vec h +}\rangle,
\end{equation}
that is, the positive-eigenvalue state of $\hat H_{\vec h}$ at a point $\vec h$ coincides with the negative-eigenvalue state of $\hat H_{-\vec h}$ at the opposite point of the Bloch sphere. It follows that
\begin{equation}
\mathcal{I}^* D^+=D^-,
\end{equation}
where $D^-$ is the bundle of low energy eigenspaces (compare with results of  Ex.~\ref{Ex:PsiMinus}.).

 At each point $\vec h$ of the Bloch sphere, the two eigenspaces, taken together, form a copy of the Hilbert space:
\begin{equation}
V^+_{\vec h}\oplus V^-_{\vec h} = \mathcal H.
\end{equation}
Globally, we have the direct sum of vector bundles:
\begin{equation}
D^+\oplus D^- = S^2\times \mathcal H.
\end{equation}
Thus, two eigenstate bundles form together a trivial bundle with fiber $\mathcal H$, which can be identified with $\mathbb{C}^2$ once the basis $\{|\spinup\rangle, |\spindown\rangle\}$ is chosen. Note that for Chern numbers we have
\begin{equation}\label{Eq:DChernSum}
c(D^+) + c(D^-) =0,
\end{equation}
which seems to agree with the previous formula, if we assume that zero on the right reflects triviality of the corresponding bundle. In the theory of vector bundles, one can associate a Chern number (in fact, several Chern numbers of different ``order'') with any complex bundle. For a trivial bundle, all Chern numbers vanish. Under taking a direct sum of vector bundles, these numbers behave according to the Whitney sum formula \cite{Nakahara2003}. In particular, for a decomposition of a trivial bundle $M$ as a direct sum of line bundles $L_i$, this formula implies that\footnote{Beware that in mathematical literature on the subject one usually considers Chern \emph{classes}, the objects to be integrated over the base space in order to find the Chern numbers. Our Chern number corresponds to the integral of the first Chern class denoted $c_1(V)$, while $c(V)$ stands for the total Chern class, or the sum of all Chern classes of a given bundle.
}
\begin{equation}\label{Eq:ChernSum}
0= c(M) = c\biggl(\bigoplus_i L_i\biggr) = \sum_i c(L_i),
\end{equation}
and the previous equation is a particular case.

\begin{figure}[h]
\begin{center}
\psfrag{T2}{$T^2$}
\psfrag{C1}{$\mathcal{C}_1$}
\psfrag{gC2}{\textcolor{half-gray}{$\mathcal{C}_2$}}
\psfrag{nC1}{\hspace{-0.3cm}$n(\mathcal{C}_1)$}
\psfrag{gnC2}{\textcolor{half-gray}{$n(\mathcal{C}_2)$}}
\psfrag{n}{$n$}
\includegraphics{img/20TorusCurv}
\caption{Tangent bundle to the torus has curvatures of opposite signs on the outer and inner surfaces.
  \label{Fig:TorusCurv} }
\end{center}
\end{figure}

\newpar{The tangent bundle to the torus.} This bundle can be thought of as a pullback of $TS^2$ under the normal map $n: T^2 \to S^2$. We already know that the integral of curvature vanishes. Let us look at the curvature at some points of the torus. Choose two positively-oriented contours, as shown in Fig.~\ref{Fig:TorusCurv}. The contour $\mathcal{C}_1$ on the outer surface of the torus is mapped to a positively-oriented contour on the sphere, while the orientation of $\mathcal{C}_2$  is reversed. Thus the curvature has opposite signs on the outer surface and the surface around the hole. Heuristically, the image of the outer surface covers the sphere once, while the image of the inner part is ``turned inside out'', so the total degree is zero.

Note that there are two circles that divide regions with opposite signs of curvature. At these circles, the Jacobian vanishes. This happens because the map is degenerate, since the circles are mapped to the poles of the sphere.

\subsubsection{Indices of singularities in pullback section}
Now let us see how the index of a singularity of a section of $E$ is transformed under the pullback.

Suppose that a section $\vec s$ of $E$ has a singularity at $q\in \mathcal{S}$ with index $\ind[\vec s(q)]$. Let  $p\in \mathcal{B}$ be some \defn{pre-image} of $q$ under $m$, that is, a point that satisfies $m(p)=q$.  We also demand that the map $m$ is invertible near $p$, so that Jacobian $J(p)\neq 0$. Then at $p$ we will have a singularity of the pullback section with index $\ind[m^*\vec s(p)]$.

Choose a contour $\mathcal{C}$ near $p$ and denote its image near $q$ by $m(\mathcal{C})$. We choose the orientation of $\mathcal{C}$ in such a way that $m(\mathcal{C})$ is positively-oriented in $\mathcal{S}$.  Then, if $J(p)>0$, the map preserves orientation of the surface at $p$, and the two indices will be the same. In the case of $J(p)<0$, the contour $\mathcal{C}$ has negative orientation with respect to the surface, and the index of the singularity changes its sign. Thus
\begin{equation}
\ind [m^*\vec s (p)] = \sign[J(p)] \ind [\vec s(q)].
\end{equation}

\exr {To show this formally, define positively-oriented coordinates $\varphi$ for $\mathcal{C}$ and $\gamma$ for $m(\mathcal{C})$. The map $m$ at the contour is given by the function $\gamma(\varphi)$. Introduce a smooth section $\vec v$ of $E$ near $q$, so that $\vec s = e^{i\beta(\gamma)}\vec v$ on $m(\mathcal{C})$. Then pull both sections back to $\mathcal{C}$ and find the index of $m^*\vec s$ at $p$.}

Now we can relate the Chern number of the pullback bundle $c(m^*E)$ with $c(E)$, as follows. Choose a point $q\in \mathcal{B}$ and consider the set of all its pre-images $\{p_i\}$, which we denote as   $m^{-1}(q)$. Suppose that the determinant $J\neq 0$ in some neighborhood of each pre-image $p_i$. It is a safe assumption, since generically a smooth real function of two variables takes zero value along some one-dimensional curves. If some of the pre-images lie on the curve where  $J=0$, one can move the point $q$ or deform the map $m$ until $J\neq 0$ for all pre-images. Let  $\vec 1 $ be a section of $E$ that has a single singularity at $q$ with index $\ind [\vec 1(q)] = c(E)$.

Then the pullback section $m^*\vec 1$ will have singularities near $p_i$, and there will be no other singularities, since pullback preserves smoothness of a section by construction. Thus, for the sum of indices of the singularities, we have:
\begin{equation}\label{Eq:ChernPBSum}
c(m^*E) = c(E) \sum_{p_i\,\in \,m^{-1}(q)} \sign [J(p_i)].
\end{equation}
Besides the Chern number of the pullback bundle, this formula provides a useful expression for the degree of a map as a discrete sum.

\exr {Plot the singularities of $\sigma_x^* \vec 1_N$ and $\mathcal{I}^* \vec 1_N$, where $\vec 1_N$ is a section of $TS^2$ with a single singularity at the south pole. Convince yourself that in both cases, singularities of pullback sections have index given by the Chern number.}

\begin{figure}[h]
\begin{center}
\psfrag{TS2}{$TS^2$}
\psfrag{TBg}{$T\mathcal{B}_g$}
\psfrag{p2}{$+2$}
\psfrag{p}{$p$}
\psfrag{m2}{$-2$}
\psfrag{gg}{\textcolor{half-gray}{$g$}}
\psfrag{n}{$n$}
\includegraphics{img/21GB-alt}
\caption{Counting of singularities of pullback section of the tangent bundle to the genus $g$ surface.
  \label{Fig:GB-alt} }
\end{center}
\end{figure}

Finally, we consider the tangent bundle to the genus  $g$ surface $\mathcal{B}_g$. This bundle can be described as a pullback
\begin{equation}
T\mathcal{B}_g = n^*TS^2
\end{equation}
under the normal map $n: \mathcal{B}_g \to S^2$. On the sphere, we move the singularity of $\vec 1_N$ to the point $p$, as shown in Fig.~\ref{Fig:GB-alt}. In the figure, the point $p$ has $g+1$ pre-images on the surface $\mathcal{B}_g$. From discussion of the curvature of the torus we know that the curvature inside the holes is negative, so the corresponding $g$ singularities  of the pullback section must have index $-2$. We also have one pre-image of $p$ in the region with positive curvature. It follows that the integral of Gaussian curvature of $\mathcal{B}_g$ is given by
\begin{equation}
\frac{1}{2\pi}\int_{\mathcal{B}_g} \kappa dS = 2 (1-g),
\end{equation}
which is the celebrated \defn{Gauss-Bonnet theorem}. Note that the result does not depend on the particular shape of $\mathcal{B}_g$ used in the calculation because of the topological stability of the Chern number.

\section{Modern theory of electric polarization}\label{Chap:ModPol}
Starting from this section, we will focus on the applications of the mathematical ideas developed above to the condensed matter physics. We will use the tight-binding approximation, a simple quantum model of a crystal that qualitatively describes many physical properties, including those related to the geometry and topology. We begin with a review of the tight-binding formalism, and then discuss the relation between a certain geometric phase and electric polarization of a crystal\footnote{Our discussion of the geometric phase formula for polarization follows Ref.~\cite{Sergeev2018}. For the original derivation based on the quantum mechanical calculation of the adiabatic current, please see textbook \cite{Vanderbilt2018} and lecture notes \cite{Resta2020}, both written by the authors of the theory. }.

\subsection{Tight-binding models and Bloch theory}
In the tight-binding approximation, one uses a discrete set of orthogonal orbitals as a basis for the Hilbert space of states of the crystal, so the operators are represented as matrices. Then the spectrum of the single-particle problem is populated according to the Pauli exclusion principle.
\subsubsection{Dimer as a two-level system}\label{Sec:Dimer}
First we consider a finite system, a molecule with two orbitals, $a$ and $b$. We denote the corresponding electronic states as $|a\rangle$ and $|b\rangle$. The Hamiltonian takes the form
\begin{equation}
\hat H = U_a |a\rangle\langle a| + U_b |b\rangle\langle b|+
t|b\rangle\langle a| + \overline t |a\rangle\langle b|,
\end{equation}
where $U_\alpha$ with $\alpha=a, b$ are on-site potentials, and $t$ is hopping amplitude from $a$ to $b$ site. Hermiticity requires that the hopping amplitude in the opposite direction is complex conjugate of $t$.

For simplicity, we assume that $U_a= -U_b=\Delta$ and $t\in \mathbb{R}$. Then in basis $\{|a\rangle, |b\rangle\}$ the Hamiltonian matrix can be expressed in terms of Pauli matrices:
\begin{equation}
H =   \sigma_x t+\sigma_z  \Delta.
\end{equation}

\subsubsection{Diatomic chain in real space}

\begin{figure}[h]
\begin{center}
\psfrag{tin}{$t_{in}$}
\psfrag{gm}{\textcolor{half-gray}{$m$}}
\psfrag{am}{$|\alpha_m\rangle$}
\psfrag{tex}{$t_{ex}$}
\psfrag{ta}{$\tau_\alpha$}
\psfrag{Ua}{$U_\alpha$}
\includegraphics{img/22TB}
\caption{Tight-binding model of a diatomic chain. $t_{in}$ and $t_{ex}$ are internal and external hopping amplitudes, the index $\alpha = a, b$ enumerates orbital inside each unit cell, $U_\alpha$ stands for on-site potential, $m$ is the number of a unit cell, $\tau_\alpha$ is the distance between the origin of the unit cell and the orbital $\alpha$. Tight-binding orbitals are denoted $|\alpha_m\rangle$.
  \label{Fig:TB} }
\end{center}
\end{figure}

Now we construct an infinite periodic crystal from dimers as unit cells, as shown in Fig.~\ref{Fig:TB}. Besides the internal hopping $t_{in}$ inside each unit cell, we add external hopping $t_{ex}$ between neighboring atoms of two adjacent cells. Thus
\begin{equation}
\hat H  = \sum_{\alpha m} U_\alpha |\alpha_m\rangle\langle \alpha_m| +
\sum_m \bigl[t_{in} |b_m\rangle\langle a_m|+ t_{ex} |a_{m+1}\rangle\langle b_m| + h.c.\bigr],
 \end{equation}
where $m$ is unit cell index an $\alpha=a,b$ enumerates orbitals. There are $N$ unit cells, and we impose periodic boundary conditions, $|\alpha_m\rangle = |\alpha_{m+N}\rangle$, so the system is invariant under translation by the lattice constant.

In the basis $\{|\alpha_m\rangle\}$ the Hamiltonian is represented by $2N\times 2N$ matrix that describes $N$ coupled dimers. Generally, it has the form:
\begin{equation}\label{Eq:GenRSH}
\hat H = \sum_{\substack{\alpha \beta \\ mn}}
|\alpha_{m+n}\rangle H^{\alpha\beta}_n \langle \beta_m |.
\end{equation}
Here $H^{\alpha \beta}_n$ describes the amplitude $|\beta_m\rangle\to |\alpha_{m+n}\rangle$. Note that this amplitude does not depend on $m$ because of periodicity.

For our chain, the matrix $H_n$ reads
\begin{equation}\label{Eq:ChainRSH}
H_n=\begin{pmatrix}
U_a \delta_{n,0} & \overline{ t_{in}} \delta_{n,0} + t_{ex} \delta_{n,1}\\
t_{in} \delta_{n,0} + \overline {t_{ex}} \delta_{n,-1} & U_b \delta_{n,0}\\
\end{pmatrix},
\end{equation}
where Kronecker delta $\delta_{ij}$ allows us to express $2N\times 2N$ matrix succinctly as a single $2\times 2$ block.

\subsubsection{Bloch waves and momentum space}
The problem of diagonalizing the Hamiltonian can be greatly simplified by introducing a new basis that respects the discrete translational symmetry of the crystal lattice. It is the \defn{Bloch wave basis}:
\begin{equation}
|\alpha_k\rangle = \frac{1}{\sqrt{N}} \sum_m e^{imk} |\alpha_m\rangle
\end{equation}
where we set lattice constant to unity. The inverse Fourier transform gives the expression of real-space orbitals in terms of Bloch waves:
\begin{equation}
|\alpha_m\rangle = \frac{1}{\sqrt{N}} \sum_m e^{-imk} |\alpha_k\rangle
\end{equation}
Periodic boundary conditions demand that $e^{ikN}=1$, so the \defn{crystal momentum} $k$ takes discrete values in the interval $[0, 2\pi)$, called \defn{Brillouin zone}, with increment $\Delta k = \frac{2\pi}{N}$.

The states $|\alpha_k\rangle$ form an orthonormal basis, which follows from $\langle \alpha_m|\beta_n\rangle = \delta_{\alpha\beta} \delta_{mn}$ and the delta function identity:
\begin{equation}
\frac{1}{N}\sum_{m} e^{im(k_i-k_j)} = \delta_{ij}.
\end{equation}

We substitute the expression for  $|\alpha_m\rangle $ as an inverse transform of $|\alpha_k\rangle$ into Eq.~(\ref{Eq:GenRSH}) and obtain
\begin{equation}
 \hat H=\sum_{\alpha \beta k} |\alpha_k\rangle
 \biggl[ \sum_n e^{-ink} H^{\alpha\beta}_n \biggr] \langle \beta_k| \quad\equiv\quad
 \sum_{\alpha \beta k}|\alpha_k\rangle
 H^{\alpha\beta}_k \langle \beta_k| = \sum_k \hat H_k,
 \end{equation}
where $\hat H_k$ is the \defn{Bloch Hamiltonian}. Thus in the Bloch wave basis the Hamiltonian of the diatomic chain describes $N$ independent two-level systems parameterized by the crystal momentum $k$. The matrix $H_k$ is given by Fourier transform of Eq.~(\ref{Eq:ChainRSH}) and reads
\begin{equation}
H_k = \begin{pmatrix}
U_a & \overline{ t_{in}} + t_{ex} e^{ik}\\
t_{in} + \overline{ t_{ex}} e^{-ik} & U_b\\
\end{pmatrix}.
\end{equation}

Again we assume that hopping amplitudes are real and there is no overall shift of on-site potentials, so the Hamiltonian matrix in the $\{|a_k\rangle, |b_k\rangle\}$ basis is
\begin{equation}
H_k = \sigma_x (t_{in}+ t_{ex}\cos k) +
\sigma_y t_{ex}\sin k +\sigma_z \Delta = \sum_i h_i \sigma_i.
\end{equation}
The eigenstates are solutions of $\hat H_k|\psi_k\rangle = \varepsilon_k |\psi_k\rangle$ and can be thought of as the states of spin-$\frac{1}{2}$ particle in the magnetic field $\vec h = (h_x, h_y, h_z)$. According to Eq.~\eqref{Eq:TwoLevelEn}, the energies of the two bands are
\begin{equation}\label{Eq:BlochBands}
\varepsilon_{k {\pm}} = \pm \sqrt{(t_{in} + t_{ex}\cos k)^2 + (t_{ex}\sin k)^2 + \Delta^2}.
\end{equation}
We will be mostly concerned with low-energy, or valence band eigenstates.

\subsubsection{Another momentum-space basis } \label{Sec:Bases}
So far Bloch theory has not included real space positions of orbitals. Denote $\tau_\alpha$ the position of the orbital of type $\alpha$ inside a unit cell. Then position operator acts on $|\alpha_m\rangle$ as follows:
\begin{equation}
\hat x |\alpha_m\rangle = (m+\tau_\alpha)|\alpha_m\rangle
\end{equation}
Later we will need the basis in momentum space that is aware of positions of the orbitals:
\begin{equation}
|\widetilde{\alpha_k}\rangle = e^{ik\tau_\alpha} |\alpha_k\rangle.
\end{equation}
The Bloch Hamiltonian eigenstates can be expressed in terms of either basis:
\begin{equation}
|\psi_k\rangle = \sum_\alpha \psi_{\alpha k} |\alpha_k\rangle = \sum_\alpha u_{\alpha k} |\widetilde{\alpha_k}\rangle.
\end{equation}
Note that while states $|\alpha_k\rangle$ obey
\begin{equation}
|\alpha_{k+2\pi}\rangle = |\alpha_k\rangle,
\end{equation}
the new basis vectors $|\widetilde{\alpha_k}\rangle$ are not periodic in momentum space because of the additional phase factor $e^{ik\tau_\alpha}$. For our chain, $\tau_a=0$, $\tau_b=\frac{1}{2}$, and so $|\widetilde{b_{k+2\pi}}\rangle = -|\widetilde{b_k}\rangle$. Since eigenstates $|\psi_k\rangle$ are $k$-periodic, the components $u_{\alpha k}$ are not.

\subsection{Difficulties with polarization}
Electric polarization is commonly associated with the spatial distribution of charge density $\rho(x)$. We will see that, for periodic systems, such a definition is problematic both on classical and quantum levels.

\subsubsection{Classical polarization and experiment}\label{Sec:ClassPol}
Consider a finite system, such as a molecule of length $a$. Its polarization can be readily defined as
\begin{equation}\label{Eq:Pdip}
P_{dip} = \frac{1}{a}\int_0^a x \rho(x) dx.
\end{equation}
One can try to define the polarization of a crystalline solid in a similar way, with the elementary unit cell playing the role of the molecule. However, we will see that even the perfect knowledge of the cell-periodic charge density is not sufficient to determine the value of polarization.

First, consider how electric polarization is measured in experiment. The standard method is to take a crystal in a non-polar state, and put electrodes on its surfaces. Then the crystal is squeezed, which leads to the appearance of the bulk polarization and surface charges. These charges are used as a measure for the bulk polarization. They can be found from the current that flows between the electrodes in the external circuit during the process of deformation.

Now we need to predict the results of the experiment form the microscopic, cell-periodic data. Suppose we know all microscopic parameters of a unit cell as functions of time. How can we find the surface charges? The problem is that microscopic theory cannot include the surface phenomena. Let us see what happens with the average dipole moment $P_{dip}$ of a unit cell when charge density is time-dependent. Taking into account continuity equation $\partial_t \rho = -\partial_x j$, one finds
\begin{equation}
\partial_t P_{dip} = \frac{1}{a}\int_0^a x (\partial_t \rho) dx = -\frac{1}{a} \int_0^a x (\partial_x j) dx = -\frac{1}{a} \bigl(xj\bigr)\big\rvert_0^a + \frac{1}{a}\int_0^a j dx.
\end{equation}
Here $j$ is the current density (or simply the current, since the model is one-dimensional). The last term describes the average current through the unit cell. Observe that this is exactly what determines the charge accumulated on the surfaces, and the bulk polarization. Thus it is natural to associate the rate of change of polarization with
\begin{equation}\label{Eq:ClassPol}
\dot P = \frac{1}{a} \int j dx = j(a) + \dot P_{dip}.
\end{equation}
This means that changes in polarization come both from the current through the unit cell boundary $j(a)$ and from the changes of the dipole moment of a unit cell. As we will see shortly, the current $j(a)$ cannot be determined from the cell-periodic charge density $\rho(x, t)$.

We take the last equation as a defining property of the electric polarization $P$. This formula does not fix the absolute value of $P$ for a given state $f$ of a crystal. However, it can be recovered as a  change of polarization from a non-polar initial state $i$ as
\begin{equation}
\Delta P = \int_{t_i}^{t_f} \dot P dt,
\end{equation}
in line with experimental methods.

\subsubsection{Quantum systems}\label{Sec:QuantPol}
Our main goal is to obtain the expression for the electric polarization of a crystal from its Bloch Hamiltonian. We start with attempts to relate current and charge density in simple tight-binding models.

First, consider the model of a diatomic molecule described in Sec.~\ref{Sec:Dimer}:
\begin{equation}
H = \Delta(p) \sigma_z+ t(p) \sigma_x,
\end{equation}
where on-site potentials and hopping amplitudes are now functions of the control parameter $p\in[0, \pi]$. Let Hamiltonian parameters vary as shown in Fig.~\ref{Fig:DimerPump}.

\begin{figure}[h]
\psfrag{t}{$t_0$}
\psfrag{gtp}{$\textcolor{half-gray}{t(p)}$}
\psfrag{p}{$p$}
\psfrag{pi}{$\pi$}
\psfrag{pi2}{$\frac{\pi}{2}$}
\psfrag{mD}{$-\Delta_0$}
\psfrag{Dp}{$\Delta(p)$}
\psfrag{D}{$\Delta_0$}
\begin{center}
\begin{tabular}{cc}
\raisebox{-1.5cm}{\includegraphics[scale=1]{img/23DimerPump}} &
\begin{tabular}{| l|c|c| c|}\hline
$p$ & $H$ & $|\psi\rangle$ & $e|\psi|^2$ \\ \hline
$0$ & $-\Delta_0\sigma_z$ & $ |a\rangle$ &
\includegraphics[scale=0.8]{img/24D1} \\
$\frac{\pi}{2}$ & $t_0\sigma_x$ & $ \frac{1}{\sqrt 2}(|a\rangle - |b\rangle)$&
\includegraphics[scale=0.8]{img/25D2}\\
$\pi$ & $\Delta_0\sigma_z$ & $ |b\rangle$&
\includegraphics[scale=0.8]{img/26D3}\\\hline
\end{tabular}
\end{tabular}
\end{center}
\caption{Charge pumping in a single dimer. Left: variation of Hamiltonian parameters  $t$ and $\Delta$ as functions of control parameter $p$. Right: Hamiltonian matrices, ground state wave functions, and electric charge densities for three values of $p$. Color of circles indicate on-site potentials. \label{Fig:DimerPump}}
\end{figure}

In the initial state, $p=0$, there is no hopping, and the ground state is $|a\rangle$ with energy $-\Delta_0$. Then we turn the hopping on and gradually reverse on-site potentials. At the point $p=\frac{\pi}{2}$ the Hamiltonian matrix is proportional to $\sigma_x$,  and the ground state is given by  $\frac{1}{\sqrt 2}(|a\rangle - |b\rangle)$. Finally, we arrive in another state with decoupled orbitals, with ground state $|b\rangle$. If the process is slow enough, the adiabatic theorem asserts that the system remains in the ground state at each stage. Note that the charge density $e|\psi|^2$ gets shifted from $a$ site to $b$ site. This is the simplest example of an adiabatic charge pump.

Now consider a one-dimensional crystal made of such dimers and set $t_{in}(p) = t(p)$ and $t_{ex}=0$ (left panel of Fig.~\ref{Fig:Currents}). The current flows inside cells from left to right. But if the roles of internal and external hopping amplitudes are interchanged, $t_{in} =0$ and $t_{ex}(p) = t(p)$, the charge flows between cells from right to left (Fig.~\ref{Fig:Currents}, right). For both protocols, the evolution of charge density $e|\psi|^2$ is exactly the same. It is the periodic nature of the crystal that makes it possible to go in opposite directions while moving from $a$ to $b$ sublattice in both cases.

\begin{figure}[h]
\begin{center}
\psfrag{j}{$j$}
\psfrag{ep}{$e|\psi|^2$}
\psfrag{p}{$p$}
\psfrag{pi}{$\pi$}
\psfrag{0}{$0$}
\includegraphics[scale=1]{img/27Currents}
\caption{Two charge-pumping protocols with identical  charge density evolution and opposite directions of current. Color of circles indicate on-site potentials.
  \label{Fig:Currents} }
\end{center}
\end{figure}

This example shows that one cannot determine electric polarization from the charge density alone. Classically, the difference stems from the term that describes the current through the unit cell boundary $j(a)$, which vanishes in the first case and is non-zero in the second case. Note that the Hamiltonian $H_k$ and thus eigenstates $|\psi_k\rangle$ do contain information about the current, but it is lost when we take modulus squared of the state vector. It should not be surprising, since in general the probability current is sensitive to the phase of the wave function.

In this extreme limit of decoupled dimers, we were able to track positions of electrons by inspection. In a more realistic situation, both hopping amplitudes are non-zero, and electrons are delocalized. We still can find the cell-periodic charge density and $P_{dip}$, but finding the current through the boundary $j(a)$ becomes a non-trivial problem. Fortunately, one can extract this information from the eigenstates  $|\psi_k\rangle$ by using the inverse Fourier transform.

\subsection{Wannier functions and geometry}
\subsubsection{Definition and basic properties}
We define $n$-th \defn{Wannier function} as inverse Fourier image of Bloch eigenstates:
\begin{equation}
|w^n\rangle = \frac{1}{\sqrt N} \sum_k e^{-ikn} |\psi_k\rangle
\end{equation}
These functions are not Hamiltonian eigenstates, but they are orthonormal and span the same Hilbert space as Bloch functions. For simplicity, we consider a single occupied band with eigenstates $|\psi_k\rangle$, which is the case for our diatomic chain at half-filling.

Since Bloch states are periodic in the real space, one can expect that Wannier functions are localized and thus well-suited to be acted on by position operator. The expectation value
\begin{equation}
x^n= \langle w^n |\hat x |w^n\rangle
\end{equation}
is called $n$-th \defn{Wannier center}.

The real space components of Wannier function $|w^n\rangle$ depend only on the difference $m-n$,
\begin{equation}\label{Eq:wpsi}
w^n_{\alpha m} = \frac{1}{N}\sum_k e^{ik(m-n)} \psi_{\alpha k}
\end{equation}
so that Wannier functions with different $n$ are related by lattice translation.

Let us find the cell-periodic charge density. For the diatomic chain at the half-filling, each unit cell contributes one electron, and charge density must satisfy
\begin{equation}
\sum_{\alpha m}\rho_{\alpha m} = Ne.
\end{equation}
From the orthonormality conditions, we have
\begin{equation}
\sum_n \langle w^n|w^n\rangle = \sum_k \langle \psi_k|\psi_k\rangle = N.
\end{equation}
It follows that the charge density can be expressed in terms of either set of functions:
\begin{equation}\label{Eq:ChargeDens}
\rho_{\alpha m} = \frac{e}{N} \sum_k |\psi_{\alpha k}|^2 = e \sum_n |w^n_{\alpha m}|^2.
\end{equation}

\exr {Consider two Bloch Hamiltonians for the diatomic chain, defined in terms of parameters $(\Delta, t_{in}, t_{ex})$ as
\begin{equation}
H^0_k: (0, 1, 0) \qquad H^1_k: (0, 0, 1).
\end{equation}
Find corresponding eigenstates $|\psi_k\rangle$ and Wannier functions $|w^n\rangle$ as their inverse Fourier transforms.  Check that two expressions for charge density given by Eq.~\eqref{Eq:ChargeDens} agree in both cases. Compute positions of Wannier centers and note that they coincide with centers of charge of the dimer molecules.}

\subsubsection{Zak phase}
Let us calculate the coordinate of zeroth Wannier center $x^0$. Acting with position operator on the Wannier function, we have
\begin{equation}
\hat x|w^0\rangle = \frac{1}{N}\sum_{m\alpha k} \psi_{\alpha k} e^{imk} \hat x |\alpha_m\rangle=\\
= \frac{1}{N}\sum_{m\alpha k} u_{\alpha k} \frac{1}{i}	\frac{\partial}{\partial k} \bigl(e^{i(m+\tau_\alpha)k}\bigr) |\alpha_m\rangle,
\end{equation}
where we used equations from Sec.~\ref{Sec:Bases}.

If $N$ is large, one can consider momentum space  sum as an integral over the Brillouin zone: $\frac{1}{N}\sum_k\rightarrow \frac{1}{2\pi}\int_{BZ}dk$. Integration by parts yields:
\begin{equation}
\hat x |w^0\rangle =
\sum_{\alpha m}|\alpha_m\rangle \frac{1}{2\pi i}\biggr[
u_{\alpha k } e^{i(m+\tau_\alpha)k }\big\rvert_0^{2\pi} - \int_{BZ} \bigl(\partial_k u_{\alpha k}\bigr) e^{i(m+\tau_\alpha)k} dk \biggl].
\end{equation}
We then perform summation over $m$, which turns $e^{imk}|\alpha_m\rangle$ into the Bloch wave $\sqrt{N}|\alpha_k\rangle$, and we have
\begin{equation}
\hat x|w^0\rangle = \frac{\sqrt{N}}{2\pi i} \sum_\alpha \psi_{\alpha k }\big\rvert_0^{2\pi} |\alpha_k\rangle +
\frac{\sqrt{N}}{2\pi} \sum_\alpha \int_{BZ} i \bigl(\partial_k u_{\alpha k}\bigr) |\widetilde{\alpha_k}\rangle dk.
\end{equation}
Observe that the first term vanishes because of periodicity of $\psi_{\alpha k}$. Finally, we multiply on the left by
\begin{equation}
\langle w^0| = \frac{1}{\sqrt N}\sum_{\beta k'} \overline{u_{\beta k'}}\langle \widetilde{\beta_{k'}}|
\end{equation}
and obtain
\begin{equation}\label{Eq:ZakPhase}
\langle w_0|\hat x |w_0\rangle = \frac{1}{2\pi} \int_{BZ} i\sum_\alpha \overline{u_{\alpha k}}\,\, \partial_k u_{\alpha k} dk \equiv \frac{\gamma}{2\pi}.
\end{equation}
The quantity $\gamma$ is called the \defn{Zak phase} \cite{Zak1989}. It has geometric origin, but it is not equivalent to the Berry phase, as we will see in the next section.

\subsubsection{Geometric interpretation}
To interpret the Zak phase geometrically, we first need to identify the relevant vector bundles. For each $k$ in the Brillouin zone, the  Bloch Hamiltonian $H_k$ acts on the Hilbert space $\mathcal H_k$. The dimension of $\mathcal H_k$ is determined by the number of degrees of freedom in each unit cell (for our diatomic chain, $\dim \mathcal H_k =2$). Taken together, these spaces form the vector bundle $\{\mathcal H_k\}$ over the Brillouin zone. The lower-energy, or valence band eigenstate $|\psi_k\rangle$ of $H_k$ defines the valence-band eigenspace $V^v_k\subset\mathcal H_k$. We are interested in the geometry of the eigenspace bundle of the valence band $V^v$. The base space is the one-dimensional circle of the Brillouin zone.

Let us find the Berry phase acquired by the eigenstate $|\psi_k\rangle$ after going around the Brillouin zone. Note that this means taking the contour integral of the potential rather than a physical process of the adiabatic evolution. For the Berry potential, one could write
\begin{equation}
A(k) = i\langle \psi_k|\partial_k \psi_k\rangle,
\end{equation}
but this product is not well-defined. As discussed in Sec.~\ref{Sec:BerryPT}, it depends on the choice of the ``constant basis'' in ${\mathcal H_k}$ that one uses to write $|\psi_k\rangle$ in components. Here, we have  two standard bases of Bloch waves over the Brillouin zone, $|\alpha_k\rangle$ and $|\widetilde{\alpha_k}\rangle$. As they are related to each other by the $k$-dependent transformation, the conditions $\partial_k|\alpha_k\rangle=0$ and $\partial_k |\widetilde{\alpha_k}\rangle=0$ clearly cannot be satisfied at the same time. Which one is to be declared constant? As we mentioned earlier, there is no ``conservation law'' in the bundle of Hilbert spaces that would give us a preferred basis.

Formally, this issue can be resolved as follows. Introduce two differential operators in the bundle of Hilbert spaces $\{\mathcal H_k\}$:
\begin{equation}\label{Eq:Dk}
D_k |\alpha_k\rangle =0 \qquad \widetilde{D_k} |\widetilde{\alpha_k}\rangle=0.
\end{equation}
Both operators act on scalar functions as an ordinary derivative $\partial_k$. These equations are to be understood as definitions of the operators rather than statements about the basis vectors. Note that $\widetilde{D_k}$ is well-defined despite the discontinuity of $|\widetilde{\alpha_k}\rangle$, since $|\widetilde{\alpha_k}\rangle$ and $|\widetilde{\alpha_{k+2\pi}}\rangle$ differ by a $k$-independent factor.

Each of the operators, when projected to the eigenspace $|\psi_k\rangle$, gives rise to the corresponding Berry potential:
\begin{equation}
A(k)=i \langle \psi_k| D_k \psi_k\rangle = i \sum_\alpha \overline{\psi_{\alpha k}}\partial_k \psi_{\alpha k},
\end{equation}
and
\begin{equation}
\widetilde A (k)= i \langle \psi_k|\widetilde{ D_k} \psi_k\rangle = i\sum_\alpha \overline{u_{\alpha k}}\,\, \partial_k u_{\alpha k}.
\end{equation}
We conclude that the Zak phase is given by
\begin{equation}
\gamma = \int_{BZ} \widetilde A (k) dk,
\end{equation}
where $\widetilde A (k)$ is the Berry potential obtained by the projection of the derivative $\widetilde{ D_k}$ in $\{\mathcal H_k\}$ defined by the condition $\widetilde{D_k} |\widetilde{\alpha_k}\rangle=0$.

Note that these operators can be thought of as covariant derivatives in the bundle of Hilbert spaces, defined in the spirit of Eq.~\eqref{Eq:CDfromPT}.  Each basis defines its own parallel transport in $\{\mathcal H_k\}$, according to Eq.~\eqref{Eq:Dk}. The fact that the basis $|\widetilde{\alpha_k}\rangle$ is not periodic in the Brillouin zone does not harm this geometric picture. Indeed, the parallel transported vector need not return to its initial state after going around a closed loop.

\subsubsection{Two Berry potentials, classical currents, and electric polarization}
From the relation between the components of $|\psi_k\rangle$ in two bases, $\psi_{\alpha k}= e^{ik\tau_\alpha}u_{\alpha k}$, one finds that the Berry potentials $\widetilde A(k)$ and $A(k)$ are related as
\begin{equation}
i\sum_\alpha \overline{u_{\alpha k}} \partial_k u_{\alpha k} =\sum_{\alpha} \overline{\psi_{\alpha k}} \partial_k \psi_{\alpha k} + \sum_\alpha \tau_{\alpha} |\psi_{\alpha k}|^2.
\end{equation}
In this section, we will develop the physical interpretation of this equation.

First, we multiply both sides by $\frac{e}{2\pi}$ and integrate each term over the Brillouin zone. Then, using Eq.~\eqref{Eq:wpsi} to transform the first term on the right, we obtain
\begin{equation}\label{Eq:ZakDecomp}
\frac{e\gamma}{2\pi} = e \sum_{\alpha m}m |w^0_{\alpha m}|^2 + \sum_\alpha \tau_\alpha\rho_{\alpha m}.
\end{equation}
The term on the left is proportional to the Zak phase. The last term on the right has the meaning of the average dipole moment of a unit cell $P_{dip}$ (with lattice constant set to unity).

We need to find an interpretation for the first term on the right. Let us denote it $\Delta Q$. Recall that the charge density can be expressed as a sum of contributions from individual Wannier functions:
\begin{equation}
\rho_{\alpha m} = \sum_n e |w^n_{\alpha m}|^2 \equiv
\sum_n \rho^n_{\alpha m}.
\end{equation}
In each contribution, we further sum over orbitals and obtain the coarse-grained version:
\begin{equation}
\rho^n_m = \sum_{\alpha} \rho^n_{\alpha m},
\end{equation}
which measures the amount of charge carried by the $n$-th Wannier function in the $m$-th unit cell. The quantity $\Delta Q$ can be expressed as
\begin{equation}
\Delta Q = \sum_m m\rho^0_m.
\end{equation}

\begin{figure}[h]
\begin{center}
\psfrag{Qp}{$Q_+$}
\psfrag{Qm}{$Q_-$}
\psfrag{gr}{$\rho^n_m$}
\psfrag{n}{$n$}
\psfrag{m}{$m$}
\includegraphics[scale=1.3]{img/28Wannier}
\caption{Contributions $\rho_m^n$ to the charge density. Circles represent electric charge carried by the $n$-th Wannier function in the $m$-th unit cell.\label{Fig:Wannier} }
\end{center}
\end{figure}
To understand the meaning of this sum, we plot the part of the graph of $\rho^n_m$  near the origin of the $(m, n)$ torus, as shown in Fig.\ref{Fig:Wannier}. Each column of circles represents the charge distribution for some Wannier function.
A little reflection shows that the quantity $\Delta Q$ measures the difference between total charges contained in the selected triangles,
\begin{equation}
\Delta Q = Q_+- Q_-.
\end{equation}
To see this, note that Wannier functions are related by lattice translations, and thus the values of $\rho^n_m$ are constant along diagonals.

Now suppose that we perform the experiment described above, and change the state of the crystal in a way that results in redistribution of charge. In this case Wannier functions will change their shape and shift along vertical lines. Observe that the increase rate $\dot Q_+$ and decrease rate $-\dot Q_-$, taken together, measure exactly the current that flows through the boundary of the zeroth unit cell. Amusingly, in spite of equality \eqref{Eq:ChargeDens}, the decomposition of charge density $\rho_{\alpha m}$ in terms of Wannier functions contains more information than that based on Bloch waves.

Now note  that the terms on the right in
\begin{equation}
\partial_t\biggl(\frac{e\gamma}{2\pi}\biggr) =
\Delta\dot Q + \partial_t \bigl(\sum_\alpha \tau_\alpha\rho_{\alpha m} \bigr)
\end{equation}
correspond to the current through the unit cell boundary $j(a)$ and time derivative of the dipole moment of a unit cell, $\partial_t P_{dip}$. Thus we have just obtained the quantum versions of the terms on the right in the defining property of the electric polarization, Eq.~\eqref{Eq:ClassPol}. We conclude that the changes in the electronic contribution to the  polarization are determined by the Zak phase:
\begin{equation}\label{Eq:DPZak}
\Delta P = \frac{e}{2\pi}\int_{t_i}^{t_f} \dot \gamma dt.
\end{equation}
With some caution in mind, this can be rewritten as
\begin{equation}\label{Eq:P2p}
\Delta P = \frac{e}{2\pi} \biggl(\gamma(t_f) - \gamma(t_i) \biggr).
\end{equation}
From the last equation, one defines the electronic contribution to the polarization as
\begin{equation}\label{Eq:PolZak}
P = \frac{e \gamma}{2\pi}, \qquad  \gamma = \int_{BZ} \widetilde A(k) dk.
\end{equation}
In other words, polarization is defined as a dipole moment of zeroth Wannier function charge distribution:
\begin{equation}
P = e\langle w^0| \hat x|w^0\rangle =
e\sum_{\alpha m} (m+\tau_\alpha) |w^0_{\alpha m}|^2 =
\sum_{\alpha m} (m+\tau_\alpha) \rho^0_{\alpha m}
\end{equation}
and equals the dipole moment of an elementary charge placed in the zeroth Wannier center.

\begin{figure}[h]
\begin{center}
\psfrag{gam2p}{$\frac{\gamma}{2\pi}$}
\psfrag{eQ}{$\frac{1}{e}\Delta Q$}
\psfrag{ePdip}{$\frac{1}{e} P_{dip}$}
\psfrag{pi}{$\pi$}
\psfrag{2pi}{$2\pi$}
\psfrag{0}{$0$}
\psfrag{1}{$1$}
\psfrag{-1}{\hspace{-0.1cm}$-1$}
\psfrag{=}{\hspace{-0.1cm}$=$}
\psfrag{+}{$+$}
\includegraphics[scale=0.8]{img/ZakTerms}
\caption{Two components of the Zak phase as functions of pumping parameter $p$. The first (second) row corresponds to the extension of the  protocol shown in Fig.~\ref{Fig:Currents} on the left (right). For calculations, see Supplementary material of Ref.~\protect\cite{Sergeev2018}.
\label{Fig:ZakTerms} }
\end{center}
\end{figure}

To illustrate the roles of terms in Eq.~\eqref{Eq:ZakDecomp}, we consider the adiabatic current that flows through the dimerized chain shown in Fig.~\ref{Fig:Currents} on the left. Suppose that charge pumping is then continued in a similar manner between unit cells for $p\in[\pi, 2\pi]$, while internal hopping is set to zero.  Finally, at $p=2\pi$ the system returns to its initial state with each Wannier center shifted by one unit cell to the right. The values of terms of  Eq.~\eqref{Eq:ZakDecomp} as functions of the pumping parameter $p$ are shown in the first row in Fig.~\ref{Fig:ZakTerms}.
The second row shows similar graphs for the analogous extension of the charge-pumping protocol shown in Fig.~\ref{Fig:Currents} on the right. In both cases, the dipole moment $P_{dip}$ is determined solely by the charge density distribution, and returns to its initial value $P_{dip}=0$ after the complete cycle. The current through the boundary, shown in the middle, makes the whole difference. It determines the direction of the shift of Wannier centers, and thus the total current.

Equations \eqref{Eq:ClassPol} and \eqref{Eq:PolZak} form the basis of the \defn{modern theory of electric polarization}. It solves the problem of finding boundary charges form the bulk microscopic theory.
This interplay between bulk and boundary physics of polarization, as well as its geometric interpretation, lies at the heart of the theory of topological insulators.  As we will see, tracking of coordinates of Wannier centers provides an effective way to detect and characterize various topological phases.

\subsection{Magnetic flux and polarization quantum}
Finally, we need to discuss is why it takes caution to switch from Eq.~\eqref{Eq:DPZak} to the two-point formula \eqref{Eq:P2p}. To see this, we compare changes in polarization with insertion of magnetic flux into a circular crystal.

\subsubsection{Magnetic flux and Peierls substitution}\label{Sec:Peierls}
Consider a periodic diatomic chain as a ring, and let it be threaded by the magnetic flux $\Phi$, as shown in Fig.~\ref{Fig:RSFlux}. We choose the vector potential to be constant along the ring,  $\Phi = \int \emf A dl$. We set the coordinate of site $a$ to zero, and that of $b$ site to $\tau$. Since the integral of $\emf A$ over unit cell is $\frac{\Phi}{N}$, we have
\begin{equation}
\int_a^b \emf Adl=\tau \frac{\Phi}{N} \qquad \int_b^{a'} \emf Adl = (1-\tau)\frac{\Phi}{N},
\end{equation}
where the lattice constant is set to unity. The first integration is performed from the $a$ site to the $b$ site inside a unit cell, and the second one from $b$ site to the $a$ site of the next unit cell.

In the tight-binding approximation, one includes magnetic field by altering hopping amplitudes, procedure known as \defn{Peierls substitution}:
\begin{equation}
t_{ab}\quad  \to  \quad t_{ab}\exp\biggl(i\frac{q}{\hbar}\int_a^b\emf Adl\biggr)
\end{equation}
For example, the internal hopping amplitude is changed as
\begin{equation}
t_{in} \quad \to \quad  t_{in} \exp\biggl( -i \Delta k\frac{\Phi}{\Phi_0}\tau\biggr),
\end{equation}
where $\Delta k = \frac{2\pi}{N}$.

\begin{figure}[h]
\begin{center}
\psfrag{P}{$\Phi$}
\psfrag{RS}{$RS$}
\psfrag{Dk}{$\Delta k$}
\psfrag{DkP}{$\Delta k\frac{\Phi}{\Phi_0}$}
\psfrag{kp}{$k'$}
\psfrag{BZ}{$BZ$}
\includegraphics[scale=1]{img/29RSFlux}
\caption{Magnetic flux insertion in the real space ($RS$) leads to the shift of momentum eigenstates along the graph of the spectrum in the Brillouin zone ($BZ$), which is parameterized by $k' = k + \Delta k\frac{\Phi}{\Phi_0}$.
\label{Fig:RSFlux} }
\end{center}
\end{figure}
Since we are concerned with real-space positions of the orbitals, we should use the basis $|\widetilde{\alpha_k}\rangle$. The corresponding Bloch Hamiltonian matrix $\widetilde{H}_k$   transforms as
\begin{equation}
\widetilde H_k = \begin{pmatrix}
U_a & \overline{t_{in}}e^{ik\tau} + t_{ex}e^{-ik(1-\tau)} \\
t_{in}e^{-ik\tau} + \overline{t_{ex}}e^{ik(1-\tau)} & U_b \\
\end{pmatrix}\quad  \to \quad
\widetilde{H}_{k+\Delta k\frac{\Phi}{\Phi_0}}
\end{equation}
The effect of the magnetic field here is similar to that in the case of particle on a ring (see Sec.~\ref{Sec:ParticleSpec}). Points that represent momentum eigenvalues are shifted along the curve of the spectrum by the amount proportional to the flux. The difference with the particle on a ring is that periodicity of the crystal in the reals space makes the momentum space compact, so the spectrum is also periodic. At the half-filling, all the states in the valence band are occupied. Thus, after the insertion of the flux quantum $\Phi_0$, the spectrum coincides with the initial one.

\subsubsection{Berry flux and polarization quantum}
Consider the Brillouin zone circle as a boundary of some fictitious disc (see Fig.~\ref{Fig:BZFlux}). Now the Zak phase can be thought as a Berry flux ``threading'' this disc: $\gamma = \int_{BZ} \widetilde A dk = \Phi_B$. The position of the Wannier charge center inside unit cell is given by  $\frac{\Phi_B}{2\pi}$. Thus we have an interesting symmetry between real and momentum space. Geometric phase acquired in one space leads to the shift of ``particles'' in the reciprocal space.

\begin{figure}[h]
\begin{center}
\psfrag{PB}{$\Phi_B$}
\psfrag{RS}{$RS$}
\psfrag{PB2p}{$\frac{\Phi_B}{2\pi}$}
\psfrag{DkP}{$\Delta k\frac{\Phi}{\Phi_0}$}
\psfrag{kp}{$k'$}
\psfrag{BZ}{$BZ$}
\includegraphics[scale=1]{img/30BZFlux}
\caption{The shift of the Wannier charge center in the real space ($RS$) is described by the Berry flux threading a fictitious disc bounded by the Brillouin zone ($BZ$).
\label{Fig:BZFlux} }
\end{center}
\end{figure}
The insertion of the magnetic flux quantum into the ring in the real space shifts the momenta of eigenstates by $\Delta k$. There is a similar process for polarization: one can change crystal parameters in such a way that the Wannier centers shift by one lattice constant. We will consider such process, known as charge pumping, in the next section.

The difference between real-space magnetic flux and momentum-space Berry flux is that the latter is defined only modulo $2\pi$. Indeed, since the points of the disc bounded by the Brillouin zone do not physically exist, one can always make on the boundary a large gauge transformation  with non-zero winding number (see Sec.~\ref{Sec:FluxAmb} and Eq.~\eqref{Eq:LargeGT}). In terms of polarization, this ambiguity corresponds to the possible relabeling of unit cells that shift all Wannier center coordinates by an integer number of lattice constants\footnote{We stress that the flux insertion (a tedious process of charge pumping) is not the same as the large gauge transformation (an instant mathematical operation).}. Thus, the polarization $P$ given by Eq.~\eqref{Eq:PolZak} is defined modulo \defn{polarization quantum}, which coincides with the elementary charge $e$ for one-dimensional systems.

When the state of the crystal is changing in time, we have a set of branches of the Zak phase $\gamma(t) + 2\pi n$, where $n\in \mathbb{Z}$. It is important to check that two phases in Eq.~\eqref{Eq:P2p} belong to the same branch. Geometrically, this means that the basis section $|\psi_k(t)\rangle$ must be smooth for all $k$ and $t$ during the process.

The polarization quantum is present even in the classical picture described in Sec.~\ref{Sec:ClassPol}, once we assume that the electric charge comes in portions of $e$. Indeed, for a given state of the crystal, imagine taking a single electron from one surface and moving it to the opposite side. As a result, polarization of the crystal is changed by the quantum, while the bulk remained intact.

\section{Charge pumping and topology}\label{Chap:Pump}
Suppose that the Hamiltonian of 1D crystal undergoes a periodic adiabatic evolution, controlled by the parameter $p$, so that $\hat H(p+2\pi) = \hat H(p)$. Consider the trajectory of some Wannier center as a function of $p$. Since the process is periodic in $p$, the final position of the Wannier center must coincide with the initial one (up to a lattice constant). Thus the Wannier center can either return to its initial position, or shift into similar position inside another unit cell. In the latter case, the variation $\hat H(p)$ is said to describe an \defn{adiabatic charge pump}, also known as \defn{Thouless pump} \cite{Thouless1983}.

From the discussion in the previous section, we know that the electric polarization can be equivalently described by putting a point charge in the position of a Wannier center. Thus the shift of a Wannier center between two cells corresponds to the current through the boundary between the cells. Moreover, the point-charge model of polarization suggests that the total charge pumped during the cycle must be quantized in units of $e$.

\subsection{Examples of charge pumps}
We start by describing two charge pumping protocols, which lead to shift of the Wannier center by one lattice constant.
\newpar{Dimerized charge pump.}\label{Sec:DimerPump}
One can easily construct a pump out of charge-pumping dimers. Consider a chain of decoupled dimers. In the first half of the cycle, charge center shifts inside each unit cell, as described in Sec.~\ref{Sec:QuantPol}. In the second half, we perform a similar process between pairs of neighboring atoms of the different unit cells, while the internal hopping amplitude is set to zero. In this protocol, all electrons are always localized, either inside a unit cell or between two unit cells.

\newpar{Charge pump with delocalized electrons.}
A more realistic model would include non-zero values of both hopping amplitudes for all $p$. Let the Hamiltonian parameters vary as follows:
\begin{equation}\label{Eq:DelocPump}
\Delta(p) = -\Delta_0 \cos p, \qquad
t_{in / ex} (p) = t_0  \pm \delta \sin p.
\end{equation}
The electrons are delocalized in this case, but the shift of Wannier centers still gives us information about the current, as shown in Fig.~\ref{Fig:DelocPump}.
\begin{figure}[h]
\begin{center}
\psfrag{p}{$p$}
\psfrag{0}{$0$}
\psfrag{ti}{$t_{in}(p)$}
\psfrag{pi}{$\pi$}
\psfrag{2pi}{$2\pi$}
\psfrag{gte}{$\textcolor{half-gray}{t_{ex}(p)}$}
\psfrag{Dp}{$\Delta(p)$}
\includegraphics{img/31DelocPump}
\caption{Charge pumping protocol with delocalized electrons and the trajectory of the Wannier center (marked by triangle) during the cycle. Double lines indicate stronger bonds. Color of circles shows on-site potential.
\label{Fig:DelocPump} }
\end{center}
\end{figure}
Note that since Wannier centers represent an observable quantity, they must respect symmetries of the crystal. For $p=0$, the chain is inversion-symmetric, with  symmetry center on an atomic site. Thus the Wannier center can also be placed only on the atom, either of $a$ or $b$ type. Since the electron will be present mostly on the atom with lower on-site potential, we find that at $p=0$ the Wannier center has coordinate $0$.  At $p=\frac{\pi}{2}$, the inversion center lies between atoms, and one can expect that the Wannier center lies at the middle of the stronger bond, that is, inside the unit cell. Similar reasoning allows one to find the Wannier center for other symmetric states and then interpolate to all $p$ by continuity.

\subsection{Topological invariant of a charge pump}
Now we wish to calculate the shift of the charge center during one cycle for a general charge pump with a given two-band Hamiltonian $\hat H(k, p)$. Since the Hamiltonian is periodic in both variables, the parameter space has the form of the torus $T^2$. The collection of valence band eigenspaces form a smooth vector bundle $V^v$ over $T^2$. We orient the torus in such a way that coordinates $(k, p)$ are positively-oriented.

\begin{figure}[h]
\begin{center}
\psfrag{k}{$k$}
\psfrag{p}{$p$}
\psfrag{Sig}{$\Sigma$}
\psfrag{C0}{$\mathcal{C}_0$}
\psfrag{C1}{$\mathcal{C}_p$}
\psfrag{Vv}{$V^v$}
\includegraphics[scale=1]{img/32PumpChern}
\caption{Vector bundle of valence band eigenstates over parameter space of a charge-pumping Hamiltonian.
\label{Fig:PumpChern} }
\end{center}
\end{figure}

Let us find the change in the Zak phase as a function of the pumping parameter $p$:
\begin{equation}
\Delta \gamma (p)= \gamma(p)- \gamma(0) =
\int_{\mathcal{C}_p} \widetilde A dk - \int_{\mathcal{C}_0} \widetilde A dk,
\end{equation}
where $\mathcal{C}_p$ is the Brillouin zone circle at $p$. Two circles form the boundary of the cylindrical segment $\Sigma$ of the torus.  Taking orientation into account, we have  $\partial\Sigma = \mathcal{C}_0-\mathcal{C}_p$. Suppose that the states $|\psi_k(p)\rangle$, or the basis section, are smooth on $\Sigma$. Then, according to the Stokes' theorem,
\begin{equation}
\Delta \gamma (p)= -	\int _{\partial \Sigma} \widetilde A dk = -\int_\Sigma \widetilde f_{kp} dkdp.
\end{equation}
Thus the shift of the Wannier center during the full cycle is given by
\begin{equation}
\Delta x^0 =\frac{\gamma(p)}{2\pi}\bigg\rvert_0^{2\pi} = -\frac{1}{2\pi} \int_{T^2} \widetilde f_{kp} dkdp = -c(V^v),
\end{equation}
and is determined by the topology of $V^v$. This gives a ``physical'' way to understand why Chern number takes integer values: after periodic evolution, Wannier center can shift only by integer number of lattice constants.

Recall that for a two-band system, the Bloch Hamiltonian can be expressed as $H_k(p) = \sum_\alpha \sigma_\alpha h_\alpha(k,p)$, where $h_\alpha$ are components of a smooth vector field $\vec h(k, p)$. The vector field $\vec h(k, p)$ defines a map $h: T^2\to S^2$ from the torus to the Bloch sphere. Thus the vector bundle $V^v$ can be thought of as a pullback under this map,  $V^v = h^*D^-$, and
the shift of Wannier center reads
\begin{equation}
\Delta x^0 = - c(h^*D^-) =  -\deg h.
\end{equation}
Note however that the curvature $\widetilde f_{kp}$ is not the pullback of the curvature $f_{\theta\varphi}$ of $D^-$.

\exr{For the dimerized charge pump, plot the vector field $\vec h (k, p)$ and check that corresponding map has non-zero degree.}

\subsection{End states of a finite pump}
In this section, we will consider the spectrum of the system as a function of the pumping parameter $p$. If the chain is finite, the process of charge pumping in the bulk leads to the appearance of unusual states at the ends of the chain. First, we will focus on the two pumping protocols described above, and then discuss the general case.

\subsubsection{Dimers and avoided level crossing}\label{Sec:EndStatesDimer}
We start with a single charge-pumping dimer with Hamiltonian matrix $H=t \sigma_x + \Delta \sigma_z$ with respect to the $\{|a\rangle, |b\rangle\}$ basis. Let parameters vary as shown in Fig.~\ref{Fig:DimerPump}. Then the spectrum is given by
\begin{equation}
\varepsilon(p) = \pm \sqrt{\Delta^2(p) + t^2(p)},
\end{equation}
and has a typical form of avoided level crossing. During the process, the charge shifts from the $a$ site to the $b$ site. Now imagine that one tries to pump the charge back to the $a$ site, but forgets to turn the hopping on. Then electron remains at the $b$ site, the energy levels cross, and the electron ends up in the high-energy state. We will refer to this situation as a ``disconnected dimer''. The corresponding spectra  are shown in Fig.~\ref{Fig:DimerFin} on the left.
\begin{figure}[h]
\begin{center}
\psfrag{0}{$0$}
\psfrag{p}{$p$}
\psfrag{pi}{$\pi$}
\psfrag{2pi}{$2\pi$}
\psfrag{e}{$\varepsilon$}
\includegraphics{img/33DimerFin}
\caption{Dimerized charge pump of finite length. Left: spectra of the charge-pumping and disconnected dimers. Middle: charge pumping in the real space. Right: spectrum of the finite chain as a function of pumping parameter.
\label{Fig:DimerFin} }
\end{center}
\end{figure}

Now consider a dimerized charge pump that consists of $N$ unit cells. Such pump for $N=3$ is shown in the middle panel of Fig.~\ref{Fig:DimerFin}. In the first half of the cycle, we have $N$ charge-pumping dimers.  However, in the second half, there are $N-1$ dimers that pump charge between unit cells, and a single disconnected dimer with zero hopping due to broken periodic boundary conditions. It follows that the spectrum of the whole system is highly degenerate and takes the form shown in Fig.~\ref{Fig:DimerFin} on the right. In the second half of the cycle, there are two branches of the spectrum that connect lower and higher bands. Importantly, the states that correspond to these branches are localized at the opposite ends of the chain.

What will happen to the electron that occupies the rightmost orbital, in the next cycle? The electron starts in the high-energy state. Recall that in discussion of the single charge-pumping dimer, we focused on the occupied state with low energy. Now note that the high-energy state moves in opposite direction during the process. It follows that Wannier centers of the high-energy band of a charge pump move backwards. This can be expressed in terms of Chern numbers for valence and conduction bands as
\begin{equation}
c(V^v) + c(V^c) = 0,
\end{equation}
similar to Eq.~\eqref{Eq:DChernSum}. Thus, after $N$ cycles, all electrons will be in the high-energy state, and will start to return to the low-energy state via the end state on the left.

\subsubsection{Destructive interference of delocalized electrons}
The end states might seem to be an artifact of the dimerized chain, where electrons are locked in the orbitals at the ends. What if they will be able to escape once both hopping amplitudes are non-zero? Let us consider the second charge-pumping protocol for a finite chain. There is a problem: in the absence of periodicity, the Bloch theorem does not apply and one has to diagonalize the $2N\times 2N$ real-space Hamiltonian matrix. First, note that the majority of the states will be extended sates that resemble Bloch waves in the bulk (but behave differently near the ends). For large $N$, breaking a single bond should not significantly alter the energies of these bulk states. In particular, we will assume that the finite chain has the same bulk gap as the periodic crystal.

We are interested in the states that are  localized at the ends of the chain and have energies that lie inside the  bulk gap. Following Ref.~\cite{Bergholtz2017}, we make a slight modification: we remove a single orbital from the right end, so that the chain ends with orbitals of the $a$ type at the both sides. Then we try the following ansatz:
\begin{equation}
|\psi_a\rangle  = \sum_n r^n |a_n\rangle,
\end{equation}
where $r\in\mathbb{C}$ is a fixed complex number. Consider the part of the Hamiltonian that consists only of the on-site terms, $\hat H_{site} = \sum_{\alpha m} U_a |\alpha_m\rangle \langle \alpha_m |$. Then $|\psi_a\rangle$ is its eigenstate:
\begin{equation}
\hat H_{site} |\psi_a\rangle = U_\alpha |\psi_a\rangle.
\end{equation}
To make $|\psi_a\rangle$ an eigenstate of the full Hamiltonian, demand that hopping part acts on it as $\hat H_{hop}|\psi_a\rangle=0$. One finds that
\begin{equation}
 r=- \frac{t_{in}}{\overline{t_{ex}}},
\end{equation}
thus with this value of $r$ the state $|\psi_a\rangle$ is an eigenstate of the full Hamiltonian, $\hat H |\psi_a\rangle = U_a|\psi_a\rangle$. By construction, the state is supported only on the orbitals of type $a$. On can say that the components of $|\psi_a\rangle$  vanish on $b$ sites as a result of the destructive interference of the hopping amplitudes from the adjacent $a$ sites. This is the reason why one needs two similar atoms at the ends of the chain.

For the protocol given by Eq.~\eqref{Eq:DelocPump},  we have $|r|>1$ for $p\in(0, \pi)$ and $|r|<1$ for $p \in (\pi, 2\pi)$. Thus the state $|\psi_a\rangle$ is exponentially localized at the right end for the first half of the cycle and resides at the left end for the second half. Its energy is $U_a(p) = -\Delta_0 \cos p$. The Fig.~\ref{Fig:DelocPumpSpec} shows the spectrum of the chain computed numerically with the \textsc{PythTB} package \cite{PythTB}. One can check that the spectrum of the bulk states is very close to that of the periodic system given by Eq.~\eqref{Eq:BlochBands}. We also have an edge state with the spectrum given exactly by the function $U_a(p)$.
\begin{figure}[h]
\begin{center}
\psfrag{0}{$0$}\psfrag{5}{$5$}\psfrag{10}{$10$}
\psfrag{-5}{$-5$}\psfrag{-10}{$-10$}
\psfrag{p}{$p$}
\psfrag{pi}{$\pi$}
\psfrag{2pi}{$2\pi$}
\psfrag{e}{$\varepsilon$}
\includegraphics{img/DelocPump}
\caption{Spectrum of the finite pump with $a$ type orbitals at both ends. The chain consists of $N=15$ full unit cells, and the other parameters are:  $\Delta_0= 5$, $t_0 = 6$, $\delta=5$.
\label{Fig:DelocPumpSpec} }
\end{center}
\end{figure}

Note that for $p=0, \pi$ the state $|\psi_a\rangle$ is delocalized, since $|r|=1$. At these points, the corresponding branch of the spectrum touches the bulk bands. In a sharp contrast with dimerized charge pump, here the current is generated by a single state that bounces back and forth through the whole crystal.

\subsection{General case}\label{Sec:PumpsTop}
In this section, we argue that the end states must  appear in any finite chain that pumps the charge in the bulk.
\subsubsection{Topological stability}
We start with the Hamiltonian $H_1(k,p)$ that describes the dimerized charge pump. The corresponding vector field $\vec h_1(k,p)$ defines the map $h_1: T^2 \to S^2$ with $\deg h_1 = -1$. Let us consider what will happen to the system under smooth deformations of the Hamiltonian, $\vec h_1 \to \vec h_1'$, that preserve the bulk gap, $|\vec h_1'| \neq 0$.

Suppose first that the crystal is periodic. In the momentum space, smooth deformations preserve the degree of the map, thus $\deg h_1' = \deg h_1$. The Chern number of the valence bundle remains the same, as well as the shift of the Wannier center per cycle in the real space.

For the finite chain, we consider smooth deformations of the corresponding real-space Hamiltonian $H_1(m, p)$. The bulk states shown in Fig.~\ref{Fig:DimerFin} will probably become non-degenerate, but the bulk gap will not close. The end states will remain exponentially localized at the ends, since their energies lie in the bulk gap.

\begin{figure}[h]
\begin{center}
\psfrag{1}{$1$}\psfrag{2}{$2$}\psfrag{3}{$3$}
\psfrag{4}{$4$}
\includegraphics{img/EdgeStab}
\caption{Topological stability of the end state spectrum branches.
\label{Fig:EdgeStab} }
\end{center}
\end{figure}

Consider the number $m_c$ of states transported from the low-energy band to the high-energy band per one pumping cycle. This number is invariant under smooth deformations of $H_1(m, p)$, which can be seen as follows. During the smooth variation of the parameters, eigenvalues must also vary smoothly. We are interested in the branch of the spectrum that corresponds to the end state. Let us try to get rid of this branch by smooth deformations of the spectrum, as shown in Fig.\ref{Fig:EdgeStab}.
\begin{enumerate}
\item The branch can be transformed into several branches, but this does not change the total number of the transported states.
\item The coupling with the bulk states is not an option, since the bulk gap must remain open.
\item We can try to couple it with an ``ordinary'' end state pulled out from the bulk bands, and this does not change the situation.
\item The only way to destroy the branch without closing the bulk gap is to couple it with the branch that goes down from the higher band to the lower one. We do have such a branch, but it resides at the opposite end of the chain.  These two branches cannot be coupled due to exponential localization of the corresponding eigenstates.
\end{enumerate}
We conclude that the number $m_c$ of states transported between bands by end state branches of the spectrum remains invariant under smooth deformation of the Hamiltonian.

\subsubsection{Classification of two-band charge pumps}
Finally, we move to the global perspective. Consider the space of all possible smooth vector fields $\vec h(k, p)$, shown schematically in Fig~\ref{Fig:ClassPump}. Each point in this space encodes a specific field configuration. We assume that nearby points correspond to ``nearby vector fields''. In other words, a curve connecting two points represents a smooth deformation. There is a set of disconnected regions formed by nowhere-vanishing vector fields. Such fields define maps $h: T^2 \to S^2$, which are classified by the degree $\deg h$. Each region consists of maps with a fixed degree. Regions are disconnected since the vector field must pass through zero at some $(k, p)$ in order to change degree of the corresponding map. There is a fact that any two maps of the same degree are connected by a smooth deformation.
\begin{figure}[h]
\begin{center}
\psfrag{h1}{$\vec h_1$}
\psfrag{h-1}{$\vec h_{-1}$}
\psfrag{h0}{$\vec h_0$}
\psfrag{h1p}{$\vec h_1'$}
\psfrag{deg-1}{\small $\deg h =-1 $}
\psfrag{deg0}{\small $\deg h =0 $}
\psfrag{deg1}{\small $\deg h =1 $}
\includegraphics{img/ClassPump}
\caption{Classification of vector fields that represent charge-pumping Hamiltonians. White regions consist of nowhere-vanishing fields.
\label{Fig:ClassPump} }
\end{center}
\end{figure}

Now let us describe this picture in physical terms. Each point in the large space describes a family of Hamiltonians $H_k(p)$. Each of the regions consists of gapped Hamiltonians with a fixed Chern number of the valence eigenspace bundle. To go from one region to another, one needs to close the bulk gap.

\begin{figure}[h]
\begin{center}
\psfrag{h}{$\vec h$}
\psfrag{-cV=-dh}{\hspace{-0.5cm}$-c(V^v)=-\deg h$}
\psfrag{=}{$=$}
\psfrag{Dx0}{$\Delta x^0$}
\psfrag{mc}{$m_c$}
\psfrag{p}{$p$}
\psfrag{x}{$x$}
\psfrag{k}{$k$}
\psfrag{e}{$\varepsilon$}
\psfrag{Lx}{$L_x$}
\psfrag{periodic}{\textcolor{half-gray}{periodic system}}
\psfrag{finite}{\textcolor{half-gray}{finite system}}
\psfrag{real}{\textcolor{half-gray}{real space}}
\psfrag{momentum}{\textcolor{half-gray}{momentum space}}
\includegraphics{img/PumpBBc}
\caption{Equivalence between Chern number, shift of Wannier centers per one pumping cycle, and number of end states in a finite pump.
\label{Fig:PumpBBC} }
\end{center}
\end{figure}

In each region, we pick a representative model. The dimerized pump $\vec h_1$ gives such a model with $\deg h_1 =-1$. One can easily construct other representatives $\vec h_n$ by redefinition of the pumping parameter, $p' = np$, where $n\in \mathbb{Z}$. For each representative, the number of branches of the end states in the finite chain is determined by the Chern number of its periodic version.
Now note that any model with a given Chern number is connected to a representative by a smooth deformation. Since the number of end states in the finite chain is invariant under such deformation, we conclude that any finite charge-pumping chain has end states. The number $m_c$ of their branches in the spectrum is determined by the corresponding Chern number, or by the shift of the Wannier centers in one pumping cycle, as summarized in Fig.~\ref{Fig:PumpBBC}.

\section{Chern insulators}\label{Chap:Chern}
Consider a two-dimensional insulator described by the Bloch Hamiltonian $H(k_x, k_y)$. The eigenspaces that correspond to the valence band form together the \defn{valence bundle} $V^v$ over the Brillouin zone torus. If this vector bundle is non-trivial, such crystal is called a \defn{Chern insulator}.

\subsection{Two-band Chern insulators}
We begin with the case of two bands, one of which is occupied. How can one define a Hamiltonian with a non-trivial eigenspace bundle? Note that we have in our disposal the Hamiltonian of a charge pump $H(k_x, p)$ which is defined on a torus and has  the required property. This suggests the following recipe for constructing a Chern insulator: start with a Hamiltonian of a charge pump and interpret the pumping parameter $p$ as a crystal momentum $k_y$. This substitution might seem artificial, but it is legitimate since both parameters are periodic.

\subsubsection{Consequences of non-trivial topology}
Let $H(k_x, k_y)$ be a Hamiltonian of a Chern insulator obtained from a charge pump under replacement $p\to k_y$. Suppose that the valence bundle has the Chern number $c(V^v)=-1$.  We assume that the coordinates  $(k_x, k_y)$ are positively-oriented. Below we consider effects of the non-trivial topology in different settings.

\newpar{Momentum space $(k_x, k_y)$.} In terms of the Bloch Hamiltonian, two-band Chern insulator is characterized by non-zero degree of the map  $h: T^2 \to S^2$ defined by the vector field $\vec h(k_x, k_y)$. Non-trivial topology also implies that
one cannot choose the phase of the eigenstate $|\psi_{k_x k_y}\rangle$ in a smooth way over the whole Brillouin zone torus. In other words, there must be a singularity in the phase at some point.

\newpar{Mixed space $(m_x, k_y)$.} Let us perform the \defn{partial Fourier transform $\mathcal F^{-1}_x$} of basis vectors
\begin{equation}
|\alpha_{k_x k_y}\rangle \xrightarrow{\mathcal F^{-1}_x} |\alpha_{m_x k_y}\rangle
\end{equation}
and obtain the family of real-space Hamiltonians $H_{m_x}(k_y)$ of the charge-pumping chain. The partial Fourier transform of Bloch eigenstates gives us the \defn{hybrid} Wannier functions
\begin{equation}
|w_{m_x} (k_y)\rangle = \mathcal F^{-1} |\psi_{k_xk_y}\rangle.
\end{equation}
The corresponding Wannier charge centers shift by one unit cell along $x$ direction when $k_y$ increases from $0$ to $2\pi$.

\begin{figure}[h]
\begin{center}
\psfrag{kx}{$k_y$}
\psfrag{e}{$\varepsilon$}
\psfrag{ef}{$\varepsilon_F$}
\includegraphics{img/EdgeMode}
\caption{In a Chern insulator, Fermi level intersects spectral branches of edge states.
\label{Fig:EdgeMode} }
\end{center}
\end{figure}

Now suppose that the crystal has finite size in $x$ direction, but retains periodicity in $y$ direction. From discussion of charge pumps, we know that such crystal must have edge states. However, these states acquire now a new physical meaning. The energy of the state with $x>0$ increases as function of $k_y$, and the state connects the valence and conduction bands. Since the Fermi level lies in the bulk gap, it intersects the dispersion relation $\varepsilon(k_y)$ of the state. This intersection describes an \defn{edge mode} moving along the right edge with Fermi velocity $v_F = \frac{d\varepsilon}{dk_y}$.  We define its \defn{chirality} as
\begin{equation}
m_i = \sign (v_F).
\end{equation}
More generally, there might be several intersections. Then the \defn{total number of edge modes} is defined as a sum of chiralities $m = \sum_i m_i$.

\newpar{Mixed space $(k_x, m_y)$.} Consider now a crystal which is periodic in $x$ direction, but has finite size in $y$ direction. What will happen on the edges? In momentum space, the Hamiltonian is described by a vector field $\vec h(k_x, k_y)$. Let us introduce new coordinates $(k_y, -k_x)$. Then the same vector field can be expressed as a new function: $\vec h'(k_y, -k_x) = \vec h (k_x, k_y)$. Since the pairs of coordinates have the same orientation, degrees of corresponding maps coincide: $\deg h' = \deg h$. It follows that the chain described by $H_{m_y}(-k_x)$ pumps charge in positive direction for decreasing $k_x$. Thus, on the side with $y>0$ there will be an edge mode with increasing energy in $-k_x$ direction.

\newpar{Real space $(m_x, m_y)$.} Consider now a real-space picture of a crystal with periodic boundary conditions. Define two-dimensional Wannier functions $|w_{m_xm_y}\rangle$ as inverse Fourier images of Bloch eigenstates $|\psi_{k_xk_y}\rangle$. In the Chern insulator, these Wannier functions cannot be exponentially localized in both directions. The reason is that the phase singularity in $|\psi_{k_xk_y}\rangle$ leads to appearance of the power-law tails of $|w_{m_xm_y}\rangle$. Sometimes such impossibility of localization of Wannier functions is taken as a basis for  definition of a topological phase.

Finally, if the crystal is finite in both directions, each side supports an edge mode. For the case of $c(V^v)=-1$, they form together a single chiral mode that goes around the sample in positive direction. Since this mode carries electric current,  Chern insulators break time-reversal symmetry.

\exr{Let $H(m_x, k_y)$ be the Hamiltonian of the Chern insulator obtained from the charge pump described by Eq.~\eqref{Eq:DelocPump}. Perform the partial Fourier transform $\mathcal F^{-1}_y$ and find the real-space tight-binding model for this Chern insulator.}

\subsubsection{Topological classification}
Now we translate the discussion of topology of charge pumps in Sec.~\ref{Sec:PumpsTop} into the context of Chen insulators. The total number $m$ of chiral edge modes is invariant under smooth deformations of the Hamiltonian that does not close the bulk gap. Moreover, this number is determined by the Chern number of the valence bundle, the result known as \defn{bulk-boundary correspondence}:
\begin{equation}\label{Eq:BBC}
-c(V^v) = m.
\end{equation}
Note that this equality connects quantities from entirely different settings. Chern number is defined for the valence bundle over the Brillouin zone, and requires periodic boundary conditions. On the other hand, edge modes live in a finite crystal, for which momentum space does not exist.

\begin{figure}[h]
\begin{center}
\psfrag{kx}{$k_x$}
\psfrag{ky}{$k_y$}
\psfrag{x}{$x$}
\psfrag{y}{$y$}
\psfrag{e}{$\varepsilon$}
\psfrag{h}{$\vec h$}
\includegraphics{img/ChernBBC}
\caption{Bulk-boundary correspondence for two-band Chern insulators. Topology of Bloch Hamiltonian determines the number of chiral edge modes.
\label{Fig:ChernBBC} }
\end{center}
\end{figure}

There is another argument for the bulk-boundary correspondence based on charge conservation \cite{Taherinejad2014}.  Suppose that the band structure has a single chiral edge mode, as in Fig.~\ref{Fig:ChernBBC}. Only lower part of the chiral branch of the spectrum is occupied. Now consider the local charge density on the end of the effective 1D crystal along $x$ direction as a function of $k_y$. At some value of $k_y$, the chiral branch crosses the Fermi level, and charge density drops down. But the charge density must be a periodic function of $k_y$, so this change must be compensated by the charge inflow form the bulk. This is described by the shift of the Wannier charge centers, or the Chern number of $V^v$. More generally, $m$ chiral edge modes require shift of Wannier center by $m$ unit cells.

Thus we can characterize any two-band Chern insulator by an integer, which gives \defn{$\mathbb{Z}$-classification} of these phases. Any two phases in the same class are related by a smooth, gap-preserving deformation of the Hamiltonian. Switching between classes changes the number of edge modes and requires closing of the bulk gap.

\subsubsection{Topological transition}
Let us consider a process, which transforms a trivial two-band  insulator into the Chern insulator with $c(V^v)=-1$. Suppose that trivial insulator is described by a constant vector field, say, $\vec h = (0, 0, 1)^T$. The Chern number is determined by the integral of the Berry curvature  $f_v$ of the valence band bundle over the Brillouin zone, or Berry flux:
\begin{equation}
c(V^v) = \frac{1}{2\pi}\int_{BZ} f_v dk_x dk_y = \frac{\Phi_B}{2\pi}.
\end{equation}
Then, in order to change the Chern number we need some sort of flux insertion. Recall from Sec.~\ref{Sec:FluxDef} that increase of curvature inside the solenoid requires rotation of the vector field on the boundary and in the surrounding space, which is the basis of the Aharonov-Bohm effect. But here we have a closed surface as a base space of the bundle, and any smooth, gap-preserving deformation of the field $\vec h$ will leave the Chern number invariant.

\begin{figure}[h]
\begin{center}
\psfrag{e}{$\varepsilon$}
\psfrag{ef}{$\varepsilon_F$}
\psfrag{h}{$\vec h$}
\includegraphics{img/TopTransition}
\caption{Changing of the Chern number by piercing of the Brillouin zone by a flux tube in fictitious parameter space. The transition requires closing the bulk gap and creates a chiral edge mode. The final configuration of $\vec h$ is a skyrmion-like structure similar to the one shown in Fig.~\ref{Fig:ChernBBC}.
\label{Fig:TopTransition} }
\end{center}
\end{figure}

To insert the Berry flux into the Brillouin torus, one needs to close the gap at some point. One can imagine this process as piercing the torus with a flux tube that ends with a monopole. This picture requires fictitious parameter space in additional dimension (which will be useful in discussion of Weyl semimetals). As shown in Fig.~\ref{Fig:TopTransition}, we start from the constant $\vec h$ field and flat bands. When the monopole moves closer to the surface of the torus, the gap narrows. Once the monopole touches the surface, $\vec h$ vanishes at this point, and the bulk gap closes. Note that the Chern number is not defined at this point, since the eigenstate bundles are not well-defined. Indeed, any vector in $\mathcal H$ is an eigenstate of a zero Hamiltonian. Then the monopole moves further, and gap reopens. Now the Brillouin zone carries Berry flux, and two bands have opposite Chern numbers. We will call such a process a \defn{topological transition}. It is the only way to create or destroy a chiral edge mode between two bands.

\subsection{Quantum anomalous  Hall effect}
The state at the edge of the Chern insulator is a unique chiral metal that cannot be realized as a genuine one-dimensional system.  To see this, we consider the typical response of one-dimensional crystals to the applied electric field. Since Bloch theory requires periodic boundary conditions, it is natural to consider such crystal as a ring. Then electric field can be applied by varying the magnetic flux that threads the ring. As discussed in Sec.~\ref{Sec:Peierls}, the insertion of one flux quantum $\Phi_0$  leads to the shift of the momentum values by $\Delta k = \frac{2\pi}{N}$.

\begin{figure}[h]
\begin{center}
\psfrag{e}{$\varepsilon$}
\psfrag{k}{$k_x$}
\psfrag{x}{$x$}
\psfrag{y}{$y$}
\psfrag{P}{$\Phi$}
\psfrag{Pez}{$\Phi=0$}
\psfrag{PeP0}{$\Phi=\Phi_0$}
\psfrag{P}{$\Phi$}
\includegraphics{img/QAH}
\caption{Spectra of a metallic circular crystal and edges of a cylindrical Chern insulator before and after insertion of the magnetic flux quantum. Filled circles indicate occupied states.
\label{Fig:QAH} }
\end{center}
\end{figure}

If the crystal is an insulator, the Fermi level lies inside the gap, so the valence band is fully occupied. After the shift of the states by $\Delta k$, nothing changes. The insulator does not react to the electric field.

The metallic crystal shown in Fig.~\ref{Fig:QAH} has a Fermi surface that consists of two points. They correspond to the modes with the group velocities $v_F$ of opposite signs. The shift of momentum states due to the electric field increases the occupation near one of the points, and reduces the number of electrons at the other. Thus, the electric field induces current in the  metallic crystal.

Consider now the edge state at the top of the cylindrical Chern insulator. It is a metal, but the Fermi surface consists only of a single point. The shift of the momentum values increases the number of occupied states, which seems to violate the charge conservation. This suggests that one should take both edges into account. Since the other edge mode has opposite chirality, its number of occupied states decreases. Thus the flux insertion leads to the transfer of a single electron between the edges. In other words, the electric field in $x$ direction results in the difference in potential along the $y$ direction, similar to the Hall effect.

In fact, Chern insulators are also known as \defn{quantum anomalous Hall} insulators. ``Quantum'' stands for the quantized Hall conductivity, which is determined by the Chern number. The Hall effect is called anomalous if it exists without external magnetic field. The first theoretical example of a Chern insulator, the Haldane model, was designed to produce the quantum Hall effect with zero net magnetic field. Later we will discuss this model in detail.

\subsection{Chern insulators with multiple occupied bands}
The vector field $\vec h$ gives a visual way to understand topology of two-band Chern insulators. However, when the number of bands is larger, no simple picture is available. Here, we discuss the ways to compute the Chern number for the multiband case, and its relation to the number of edge modes.

\subsubsection{Isolated bands}
First, we consider the case of isolated, or non-degenerate bands. In this case, each band has a well-defined complex line bundle $V^i$ of eigenspaces. Number of singularities in any section of this bundle determines its Chern number $c(V^i)$. More generally, one can define the Chern number for any group of isolated bands, according to Eq.~\eqref{Eq:ChernSum}. One such group of particular interest is the set $v$ of valence bands.

Note that one cannot change the Chern number $c(V^i)$ of a single band without affecting the other bands. Indeed, the Hilbert space bundle is trivial, since it has global non-vanishing sections given by the Bloch basis $\{|\alpha_k\rangle\}$. On the other hand, this bundle is the direct sum of eigenstate bundles of individual bands. It follows from Eq.~\eqref{Eq:ChernSum} that the sum of individual Chern numbers over all bands must vanish:
\begin{equation}
\sum_i c(V^i)=0.
\end{equation}
Thus, a process that changes Chern numbers $c(V^i)$ must involve at least two bands. An example of such process is the topological transition considered above, which changes Chern numbers by $\pm 1$.

For instance, consider an insulator with four trivial isolated bands, two of which are occupied. Let us first perform the topological transition between two occupied valence bands. The result will look similar to Fig.~\ref{Fig:TopTransition}, but with both bands lying now below the Fermi level. The branch of the spectrum that corresponds to the edge states will not give rise to the edge \emph{mode}, as it does not intersect Fermi energy. The bundle of valence band eigenstates $V^v$ will remain trivial, since
\begin{equation}
c(V^v) = \sum_{i\in v} c(V^i) = 0.
\end{equation}
Note also that this process does not close the bulk gap.

Now imagine that the topological transition takes place between two bands, one of which is occupied and the other is not. In this case, bulk gap closes, and we end up with a non-trivial Chern insulator, which has chiral edge mode and non-zero Chern number of the valence bundle $c(V^v)$. In a general situation, the number of edge modes is determined by the number of topological transitions required to create them, so the bulk-boundary correspondence \eqref{Eq:BBC} holds in the same form.

\begin{figure}[h]
\begin{center}
\psfrag{1}{$1$}
\psfrag{-1}{$-1$}
\psfrag{2}{$2$}
\psfrag{-2}{$-2$}
\psfrag{0}{$0$}
\includegraphics{img/Chern4}
\caption{Transformation of a four-band Chern insulator with $c(V^v)=-2$. Numbers show individual Chern numbers of bands.
\label{Fig:Chern4} }
\end{center}
\end{figure}
Consider the sequence of transformations shown in Fig.~\ref{Fig:Chern4}. We start with the ``direct sum'' of two decoupled two-band Chern insulators, which we will denote as $C_{-1}\oplus C_{-1}$. Then we use a topological transition among valence bands, and the other one among conduction bands to change the individual Chern  numbers without closing the bulk gap. Finally, the intersecting edge states are gapped out, and as a result we have the direct sum of Chern insulator with a trivial one, $C_{-2}\oplus I$. In this way, one can redistribute Chern numbers inside the valence or conduction band subspaces, while staying inside the same equivalence class, which is  determined  here by $c(V^v) =-2$.

Any classification is based on some equivalence relation that connects objects inside each class. Here we declare two systems to be equivalent if they are connected by a smooth deformation of the Hamiltonian that preserves
the Chern number $c(V^v)$ associated with occupied bands. Note that the systems shown in Fig.~\ref{Fig:Chern4} have the same Chern number $c(V^v)$ as a \emph{two-band} Chern insulator $C_{-2}$. Are they equivalent? Changing the dimension of Hilbert space, or addition of orbitals into the unit cell, does not look like a smooth deformation, so these systems must belong to different classes. For Chern insulators,  such distinction  does not have physical significance, and one can include addition of trivial bands into definition of the equivalence relation. But in some cases the difference becomes crucial: there are \defn{fragile} phases, which loose their topological properties when trivial bands are added to the band structure \cite{Po2018}.

\subsubsection{Degenerate bands and Wilson loop}
In the presence of degeneracy between valence bands, the corresponding eigenstate bundles are not well-defined. Indeed, any linear combination of the eigenstates of two bands at the degeneracy point
is also an eigenstate with the same energy. To characterize the topology of a set of degenerate bands, one considers them as a single ``composite band''. Geometrically, the vector bundle in question is the collection of subspaces of the Hilbert space $\mathcal H$ spanned by the eigenstates of bands from the set. Here we briefly outline how the concepts of covariant derivative, Chern  number, and electric polarization are generalized to the case when two or more bands are degenerate.

Let $\{|\chi_k^n\rangle\}$ be a set of states that span the occupied band subspace, or a set of basis sections for the valence band bundle $V^v$. Note that these basis vectors need not be Hamiltonian eigenstates, so the index $n$ does not in general enumerate bands. We define covariant derivative in $V^v$ as
\begin{equation}
\nabla_k |\varphi_k\rangle  = \Proj^v \widetilde D_k |\varphi_k\rangle,
\end{equation}
where $\Proj^v = \sum_n |\chi_k^n\rangle\langle \chi_k^n|$ is the projection onto the space of valence bands, $\widetilde D_k$ is the derivative in the Hilbert space bundle $\{\mathcal H\}$ that vanishes on $|\widetilde \alpha_k\rangle$, and $|\varphi_k\rangle$ is a section of $V^v$. In components, we have
\begin{equation}
\nabla_k |\varphi_k\rangle = \nabla_k \sum_m \varphi_k^m |\chi_k^m\rangle =
\sum_m \bigl(\partial_k \varphi^m_k - i \sum_n \varphi^n_k \widetilde A^{mn}_k \bigr) |\chi^m_k\rangle,
\end{equation}
where $\widetilde A^{mn}_k = i \langle \chi_k^m|\widetilde D_k|\chi_k^n\rangle$ is an element of the connection \emph{matrix}. The parallel transport around a closed loop results now in a unitary transformation, instead of a scalar phase factor in the single-state case.

The curvature can be defined as in Exercise \ref{Ex:CurvComm}, via the commutator of covariant derivatives in directions $k_x$ and $k_y$:
\begin{equation}
[\nabla_x, \nabla_y]|\varphi_k\rangle =
\frac{1}{i} \biggl(\partial_x \widetilde A_y - \partial_y \widetilde A_x  -i [\widetilde A_x, \widetilde A_y]\biggr)|\varphi_k\rangle \equiv \frac{1}{i} \widetilde F_{xy} |\varphi_k\rangle,
\end{equation}
as the reader may check. Note the appearance of the commutator term due to the non-Abelian nature of connection.

The Chern number of the valence band bundle is defined as the integral of the trace of curvature matrix:
\begin{equation}
c(V^v) = \frac{1}{2\pi} \int_{BZ} \tr (\widetilde F_{xy}) dk_x dk_y,
\end{equation}
which reduces to the sum of individual Chern  numbers for the case of isolated bands.

The electric polarization of a one-dimensional crystal in is given by
\begin{equation}
P_x  = \frac{e}{2\pi} \sum_n \theta_x^n,
\end{equation}
where $\theta_x^n$ are the phases of eigenvalues $e^{i\theta_x^n}$ of the parallel transport matrix $U$.  The matrix $U$ transforms the basis vectors $|\chi_k^n\rangle$ at $k_x=0$ into their parallel transported counterparts, obtained after completing the the loop $k_x = 0\to 2\pi$. In contrast with the single-band case, here we have additional freedom to choose the basis states $|\chi^n_k\rangle$. Different choices are related by multi-band gauge transformations.  Under smooth gauge transformations, the values of individual phases $\theta_x^n$ may vary, while their sum stays invariant. The change of polarization during a charge pumping cycle is given by the sum of ``winding numbers'' of the phases $\theta_x^n$. In the context of Chern insulators, this can be expressed as
\begin{equation}
c(V^v) = \frac{1}{2\pi} \sum_n \int_0^{2\pi}
(\partial_{k_y} \theta_x^n) dk_y,
\end{equation}
where one calculates phases $\theta_x^n$ for an effective one-dimensional crystal in $x$ direction and interprets $k_y$ as a pumping parameter. Thus, the Chern number again determines the number of chiral edge modes, by the charge conservation argument.

One basis choice for $V^v$ is of particular importance. The inverse Fourier transform of such basis vectors defines the \defn{maximally localized Wannier functions} \cite{Marzari2012}, which have minimal spread in the real space. It turns out that this basis also respects parallel transport, in the sense that the vectors do not mix in the process, and each vector undergoes only a phase rotation at a constant rate. The connection matrix with respect to this basis is diagonal and $k$-independent:
\begin{equation}
\widetilde A_k^{mn} = \frac{\lambda^n}{2\pi} \delta_{mn},
\end{equation}
so that $e^{i\lambda^n}$ are the eigenvalues of the parallel transport matrix $U  = \exp(i\int A dk)$. Note that in general, the last formula requires path-ordering.
The set of phases $\{\lambda^n\}$ describes the coordinates of charge centers of maximally-localized Wannier orbitals.  In somewhat different terminology, the set $\{\lambda^n\}$ is known as the \defn{Wilson loop spectrum}. The name comes from high-energy physics, where Wilson loop operator is also related to the non-Abelian parallel transport. Analysis of this spectrum provides a powerful tool for characterization of topological phases \cite{Alexandradinata2014, Benalcazar2017}.

We conclude that Chern insulators can be classified by the number of the chiral edge modes. For two-band insulators, one finds this number as a degree of a map $h: T^2 \to S^2$. In the multi-band case, Chern number can be found as a sum of winding numbers of eigenvalues $\theta^n$, which also gives an integer.
Thus, we have established the $\mathbb{Z}$-classification for all Chern insulators.

\section{Role of symmetry}\label{Chap:Sym}
The existence of Chern insulator does not require any symmetry, except the discrete translational symmetry of the crystal lattice. In fact, such crystal must break time-reversal symmetry, since otherwise there would be no chiral edge modes. Many materials are invariant  under time reversal and can have spatial crystalline symmetries as well. Such materials can also be topological, and in this case symmetry strongly affects their topological classification.

In this section, we will see how a simple spatial symmetry restricts possible values of the Zak phase and  will discuss topological classification of time-reversal invariant insulators.

\subsection{Reflection-symmetric 1D chain}
\subsubsection{Symmetry of Bloch Hamiltonian}
Consider a one-dimensional crystal which is invariant under reflection in a plane perpendicular to the chain. Assume that the unit cell center lies at the symmetry plane. Then the reflection is represented by a linear operator $\hat R$ that maps orbitals in the $m$-th unit cell to those in the cell with index $-m$. The matrix of $\hat R$ is defined as follows:
\begin{equation}
\hat R |\alpha_m\rangle = \sum_\beta |\beta_{-m}\rangle\langle \beta_{-m}| \hat R |\alpha_m\rangle \equiv \sum_\beta |\beta_{-m}\rangle R^{\beta \alpha}.
\end{equation}
Then the action of $\hat R$ on the Bloch basis states reads
\begin{equation}
\hat R |\alpha_k\rangle = \frac{1}{\sqrt{N}} \sum_{m\beta} e^{ikm} |\beta_{-m}\rangle R^{\beta\alpha} =
\sum_\beta |\beta_{-k}\rangle R^{\beta\alpha}.
\end{equation}
The crystal is reflection-symmetric if the symmetry maps hopping amplitudes and on-site potentials to those of equal strength. In such case, the symmetry operator  commutes with the Hamiltonian: $\hat R\hat H = \hat H\hat R$. To translate this into momentum space, we evaluate each product on a Bloch basis state:
\begin{equation}
\hat R\hat H |\alpha_k\rangle =
\hat R\sum_\beta |\beta_k\rangle H^{\beta\alpha}_k =
\sum_\gamma |\gamma_{-k}\rangle R^{\gamma\beta}H^{\beta\alpha}_k,
\end{equation}
\begin{equation}
\hat H\hat R |\alpha_k\rangle =
\hat H\sum_\beta |\beta_k\rangle R^{\beta\alpha} =
\sum_\gamma |\gamma_{-k}\rangle H^{\gamma\beta}_{-k}R^{\beta\alpha}.
\end{equation}
It follows that the symmetry requires the Bloch Hamiltonian to satisfy
\begin{equation}
RH_k = H_{-k}R \quad \Rightarrow \quad RH_kR^{-1}=H_{-k}.
\end{equation}

\subsubsection{Quantization of Zak phase}\label{Sec:QuantZak}

\begin{figure}[h]
\begin{center}
\psfrag{am}{$|a_m\rangle$}
\psfrag{gbmm}{\textcolor{half-gray}{$|b_{-m}\rangle$}}
\includegraphics{img/ReflectionOrb}
\caption{Reflection-symmetric chain with two $p$-orbitals in unit cell. Dashed line indicates reflection plane.
\label{Fig:ReflectionOrb}}
\end{center}
\end{figure}

For example, consider a chain with two atomic $p$-orbitals in the unit cell, shown in Fig.~\ref{Fig:ReflectionOrb}. Assume that both orbitals are placed at the center of the cell, and that the reflection symmetry exchanges the orbitals. The reflection operator acts as $\hat R|a_m\rangle =|b_{-m}\rangle$, and similarly $\hat R|b_m\rangle =|a_{-m}\rangle$. Thus its matrix is $R=\sigma_x$. As usual, we express Bloch Hamiltonian matrix as $H_k = \sum_i h_i(k)\sigma_i$.  From conjugation properties of Pauli matrices\footnote{Recall that the conjugation by $\sigma_x$ describes $\pi$ rotation about $x$ axis in the space of Pauli matrices.} we deduce that
\begin{equation}\label{Eq:Symm_h}
h_x(k)=h_x(-k) \qquad h_y(k)=-h_y(-k) \qquad h_z(k)=-h_z(-k).
\end{equation}
Let us find the points $ {k^\star} $ in momentum space that are fixed under reflection. These are solutions of the equation $k=-k$. Since the Brillouin zone is a circle, we have two such points: $ {k^\star} =0, \pi$. Then $h_y( {k^\star} )$ and  $h_z( {k^\star} )$ are forced to vanish, and $H( {k^\star} )\sim \sigma_x$. There are two essentially different situations: the signs of $h_x$ at $ {k^\star} $ may be the same or different.

First, consider the case when $\sign[h_x(0)] = \sign[h_x(\pi)]$. The simplest Hamiltonian in this class is given by $H^0_k= \sigma_x$. One can deform it away from  $\sigma_x$ at intermediate points $k\neq  {k^\star} $. However, the condition (\ref{Eq:Symm_h}) requires that the images of the two halves of the Brillouin zone on the Bloch sphere are two symmetric loops with opposite orientation, as shown in Fig.~\ref{Fig:ReflectionBloch} on the left.

\begin{figure}[h]
\begin{center}
\psfrag{O}{$\Omega$}
\psfrag{O1}{$\Omega_1$}
\psfrag{sx}{$\sigma_x$}
\psfrag{sy}{$\sigma_y$}
\psfrag{sz}{$\sigma_z$}
\includegraphics{img/ReflectionBloch}
\caption{Two classes of reflection-symmetric two-band Hamiltonians as curves in the space of Pauli matrices.
\label{Fig:ReflectionBloch}}
\end{center}
\end{figure}

For the second case, $\sign[h_x(0)] = - \sign[h_x(\pi)]$, a representative Hamiltonian $H_k^1$ can be described as a uniform $2\pi$ rotation in the $xy$ plane in the space of Pauli matrices. The symmetry constraints require that any possible deformation is symmetric, and the loop always encloses half of the Bloch sphere (Fig.~\ref{Fig:ReflectionBloch}, right).

Next, let us find the Zak phase in both cases. Since the orbitals have the same spatial position, we set $\tau_a= \tau_b =0$, thus $|\alpha_k\rangle = |\widetilde \alpha_k\rangle$, so the subtleties of basis choice can be safely ignored. The Berry connections coincide $A=\widetilde A$, and the Zak phase equals the Berry phase of a two-level system. It can be easily found from the solid angle spanned by the image of the Brillouin zone on the Bloch sphere.

In the first case, if a single loop encloses a solid angle $\Omega_1$, we have for the polarization:
\begin{equation}
P = \frac{e}{2\pi} \biggl(\frac{\Omega_1}{2} - \frac{\Omega_1}{2} \biggr)=0.
\end{equation}
In the second case the solid angle  always equals $2\pi$, so that
\begin{equation}
P =\frac{e}{2\pi} \frac{\Omega}{2} = \frac{e}{2}.
\end{equation}
We conclude that under the symmetry constraints given by Eq.~\eqref{Eq:Symm_h}, the Zak phase can be either $0$ or $\pi$. In other words, the reflection symmetry requires that Wannier charge centers are placed either at a reflection plane or in the middle between two planes. Thus, the quantization of the Zak phase is a consequence of the classical symmetry of charge distribution, and can be expected in a general situation.

Now we use two representative Hamiltonians to illustrate several important points. We will call the state described by $H^0$ \defn{trivial}, and that described by $H^1$ \defn{polarized} state of a reflection-symmetric chain.

\begin{figure}[h]
\begin{center}
\includegraphics{img/ChargeCenters}
\caption{Charge center distributions for trivial and polarized chains of finite length.
\label{Fig:ChargeCenters}}
\end{center}
\end{figure}

\newpar{Polarization of a finite chain.} Fig.~\ref{Fig:ChargeCenters} shows charge center distributions for two finite chains described by the representative Hamiltonians. At half-filling, there will always be extra charge at one of the ends of the second chain. It breaks the symmetry of the finite crystal  and  gives it macroscopic polarization.

This also illustrates the fact that the bulk polarization is naturally defined as a geometric phase rather than a vector. Indeed, any reflection-odd vector must vanish in a symmetric system:
\begin{equation*}
\vec v = - \vec v \quad \Rightarrow \quad \vec v=0.
\end{equation*}
In contrast, the phase is defined modulo $2\pi$, and
\begin{equation*}
\gamma = -\gamma \quad \Rightarrow \quad \gamma = 0, \pi.
\end{equation*}
Note also that the bulk states with $\gamma=\pi$ and $\gamma=-\pi$ are equivalent. In a finite crystal, these sates have extra charge on one of the ends, which results in the macroscopic polarization $P = \pm \frac{e}{2}$ that differs by the polarization quantum.

\newpar{End states of the polarized chain.} Our experience with simple models of charge pumps suggests that any finite symmetric crystal in the polarized state will have the mid-gap edge states in the spectrum. However, these states are not protected, as they can be pushed away by some surface potential. Such a potential need not respect the symmetry, since at the end of the crystal the symmetry is already broken.

\newpar{Symmetry and topological classification.} Now suppose that we want to find an interpolation between the trivial and polarized states. One way is to shift charge centers form $0$ to $\frac{1}{2}$ by charge pumping. This will break the symmetry, since only these positions are invariant. In terms of Hamiltonians on the Bloch sphere, a smooth interpolation between a trivial loop and a diameter requires breaking the symmetry constraints.

The other way is to preserve the symmetry during the deformation. Since $H( {k^\star} )\sim \sigma_x$ and has different signs in two states, this will lead to the gap closing at some point. For example, consider
\begin{equation}\label{Eq:HamTheta}
H^\theta = (1-\theta) H^0 + \theta H^1
\end{equation}
with $\theta \in [0; 1]$. The gap at $k=\pi$ closes at $\theta=\frac{1}{2}$, and the chain becomes metallic. In this state Wannier functions are completely delocalized, and positions of the charge centers are ill-defined.

This suggests the general scheme of classification of topological matter with symmetries: two phases are defined to be inequivalent, if any smooth interpolation between them either closes the bulk gap or breaks the symmetry. In our case, we have two classes of reflection-symmetric models, with $H^0$ and $H^1$ serving as representatives.

\subsubsection{Eigenvalues of symmetry operator}
Now we switch back to the diatomic chain with two atomic $s$-orbitals, shown in Fig.~\ref{Fig:TB}. Suppose that the chain is inversion-symmetric, and inversion center lies in the middle of the unit cell. Inversion $\mathcal{I}$ exchanges the orbitals and is represented by $\sigma_x$ in the standard basis. Mathematically, this system is analogous to the one considered above. The only difference is that Berry potentials do not coincide $A\neq \widetilde A$, but the quantization of the Zak phase still holds for the classical reasons.

At fixed points $ {k^\star} $ in the Brillouin zone, the Hamiltonian commutes with inversion operator:
\begin{equation}
\hat H_{ {k^\star} } \mathcal{I} = \mathcal{I} \hat H_ {k^\star} ,
\end{equation}
and one can choose eigenstates of $\hat H_ {k^\star} $ such that they are also eigenstates of $\mathcal{I}$. In our case, the matrix of inversion happens to coincide, up to a scalar factor,  with the Hamiltonian matrix $H_ {k^\star} $. Thus  the eigenstates $|\psi_ {k^\star} \rangle$ have \defn{inversion eigenvalues} $\lambda( {k^\star} )=\pm 1$, depending on the band and the sign of $h_x( {k^\star} )$.

\begin{figure}[h]
\begin{center}
\psfrag{e}{$\varepsilon$}
\psfrag{pi}{$\pi$}
\psfrag{k}{$k$}
\includegraphics{img/InvEigval}
\caption{Spectra of three Hamiltonians described by Eq.~\eqref{Eq:HamTheta} for $\theta = 0, \frac{1}{2}, 1$. Filled (empty) circles indicate $+1 (-1)$ values of the inversion eigenvalue $\lambda( {k^\star} )$ for each eigenstate.
\label{Fig:InvEigval}}
\end{center}
\end{figure}

One can label the bands at the high-symmetry points by inversion eigenvalues. Band structures of three typical Hamiltonians with such labels are shown in Fig.~\ref{Fig:InvEigval}. This illustrates that the only way to change the symmetry eigenvalues, while preserving the symmetry, is to close the bulk gap. In the figure, gap closes at $ {k^\star} =\pi$, which allows bands to swap inversion eigenvalues. Thus we can use the sign  of the product of symmetry eigenvalues for the occupied band, $\sign[\lambda(0) \lambda(\pi)]$, to distinguish between trivial and polarized states of the crystal. For the diatomic chain, this is just a rephrasing of our previous conclusions. However, this approach can be applied to much more general cases.

\exr{Plot eigenstates $|\psi_ {k^\star} \rangle$ in the real space for both bands of $H^0_k$ and $H^1_k$. Convince yourself that they are inversion-symmetric and have eigenvalues as shown in Fig~\ref{Fig:InvEigval}.}

\subsection{Time reversal}
Our next goal is to construct a time-reversal invariant topological insulator. We start with review of how time-reversal symmetry act on spinless and spinful particles. Then we will give qualitative description of the classification of $T$-invariant topological insulators. Finally, we will derive a topological invariant for systems that have both inversion and time-reversal symmetries.

\subsubsection{Spinless particles}
One can guess the form of time reversal operator $T$ from the action on the plane wave $\psi_p=e^{ipx}$. We demand that the coordinate is left invariant under $T$, while momentum must be reversed. Then
\begin{equation}
 T \psi_p = \psi_{-p} = e^{-ipx} = \overline{\psi_p},
 \end{equation}
which suggests that time-reversal acts as complex conjugation, $T=K$. Indeed, one checks that for any solution $\psi(x,t)$ of the the Schr\"odinger equation, the conjugate wave function $\overline{\psi}(x,t)$ gives the solution for the time-reversed problem.

The time reversal operator is \defn{anti-linear}:
\begin{equation}
K (a|\psi\rangle)= \overline a K|\psi\rangle
\end{equation}
for a complex scalar $a$. Such an operator cannot be represented by a matrix and depends on the basis choice:
\begin{equation}
 K|\alpha\rangle = |\alpha\rangle \quad \Rightarrow \quad
 K|\psi\rangle = \sum_\alpha \overline{\psi_\alpha} |\alpha\rangle.
 \end{equation}
In the plane wave example above it is natural to assume that $T|x\rangle = |x\rangle$. Thus $T=K$ holds in the basis of position operator eigenstates.

Let us see how $T$ interacts with inner product. Since $\langle K\psi| = \sum_\alpha \langle\alpha| \psi_\alpha|$, we have
\begin{equation}
 \langle K\psi|K\varphi\rangle = \overline{\langle \psi|\varphi\rangle} =
 \langle\varphi|\psi\rangle.
 \end{equation}
The operators with this property are called \defn{anti-unitary}.

Consider now Bloch waves in a crystal. From $T|\alpha_m\rangle=|\alpha_m\rangle$ one finds that $T|\alpha_k\rangle = |\alpha_{-k}\rangle$. Thus $T$-invariant Bloch Hamiltonian must satisfy:
\begin{equation}\label{Eq:TRHam}
TH_kT^{-1}=H_{-k} \qquad \Rightarrow \qquad \overline{H_k} = H_{-k}.
\end{equation}

\subsubsection{Spinful particles}
When acting on a spinful particle, time reversal must also  reverse spin. For spin-$\frac{1}{2}$ particles, this is realized by
\begin{equation}
T=\exp\biggl(\frac{\pi}{2}\frac{\sigma_y}{i}\biggr)K = -i\sigma_y K = \begin{pmatrix}
0&-1\\ 1&0\\
\end{pmatrix}K.
\end{equation}
Now $T$ is a product of a unitary and anti-unitary operators, and is again anti-unitary.

Consider the action of $T$ on a state $\psi_{\vec s}$ on the Bloch sphere. For the parametrization given by Eq.~(\ref{Eq:PsiPlus}), the complex conjugation $K:\varphi\mapsto -\varphi$ acts as a reflection in the plane $\varphi=0$, and the exponential acts as a $\pi$ rotation about $y$ axis. The resulting state is proportional to the $\psi_{-\vec s}$ state for eigenstates  both of high and low energy.

\exr{Act with $T$ on spin eigenstates $\psi_\pm(\theta, \varphi)$ and check that time reversal operator flips the spin direction.  }

A crucial property of $T$ for spin-$\frac{1}{2}$ particles is that
\begin{equation}
T^2 = \begin{pmatrix}
0&-1\\ 1&0\\
\end{pmatrix}^2 = -\mathbb{I}.
\end{equation}
This leads to the following \defn{Kramers theorem}.
Suppose that the Hamiltonian $\hat H$ commutes with some operator $T$, which is anti-unitary and squares to $-\mathbb{I}$.  Then each energy level is at least two-fold degenerate. To prove this, first note that
\begin{equation}
\hat H|\psi\rangle = \varepsilon |\psi\rangle \quad\Rightarrow\quad
\hat HT|\psi\rangle = T\hat H|\psi\rangle  = T\varepsilon |\psi\rangle = \varepsilon T|\psi\rangle,
\end{equation}
so $T|\psi\rangle$ is an eigenstate with the same eigenvalue as $|\psi\rangle$. From the properties of $T$ one finds that these states are orthogonal:
\begin{equation}
\langle \psi |T\psi\rangle = \langle T^2 \psi|T\psi\rangle =
-\langle \psi |T\psi\rangle\quad  \Rightarrow \quad |\psi\rangle\perp |T\psi\rangle,
\end{equation}
and degeneracy follows from their orthogonality.

In order to describe spinful electrons in a crystal, one replaces each real-space orbital $|\alpha_m\rangle$ with a basis of a two-level spin system $|\alpha_{m\sigma}\rangle$, where $\sigma = \spinup, \spindown$. In this way, Bloch eigenstates become $|\psi_k\rangle = \sum_{\alpha \sigma} \psi_{\alpha k \sigma} |\alpha_{k \sigma}\rangle$. Time reversal sends spin eigenstate $\psi_{\vec s}$ at crystal momentum $k$ to the state proportional to $\psi_{-\vec s}$ at $-k$. As in the case of spatial symmetries, fixed points $ {k^\star}  = 0,\pi$ play an important role. In this context they are known as \defn{time-reversal invariant momenta}, or TRIM.

\subsection{$T$-invariant topological insulators and $\mathbb{Z}_2$ classification}
\subsubsection{Two dimensions: classification by edge sates}
Now we can construct a time-reversal invariant topological insulator. To this end, we put two Chern insulators with opposite Chern numbers into spin-up and spin-down subspaces of a generic spinful four-band model. On one hand, time reversal changes the sign of the Chern number and direction of the chiral edge mode. On the other hand, it flips spin. Hence these Chern insulators are mapped to each other by $T$, and the whole system remains time-reversal invariant.  This construction is somewhat artificial, since any realistic Hamiltonian would contain terms that do not conserve $s_z$, so in general, one cannot define spin-up and spin-down subspaces. Nevertheless, we can use such a model as a starting point and look what will happen under general $T$-invariant deformations of the Hamiltonian.

\begin{figure}[h]
\begin{center}
\psfrag{e}{$\varepsilon$}
\psfrag{pi}{$\pi$}
\psfrag{k}{$k$}
\psfrag{0}{$0$}
\psfrag{cpm1}{$c_{\spinup\spindown}=\pm 1$}
\psfrag{cpm2}{$c_{\spinup\spindown}=\pm 2$}
\psfrag{cpm0}{$c_{\spinup\spindown}= 0$}
\includegraphics[scale  = 1.2]{img/QSH}
\caption{Band crossings at TRIM (filled circles) are protected by Kramers degeneracy against any $T$-invariant perturbations. The spectrum corresponding to $c_{\spinup\spindown}=\pm 2$ can be smoothly deformed into the spectrum of a trivial insulator.
\label{Fig:QSH}}
\end{center}
\end{figure}

Fig.~\ref{Fig:QSH} shows typical spectra of edge states of spin-polarized pairs of Chern insulators.  Suppose first that Chern numbers for the spin subsystems are $c_{\spinup \spindown} =\pm 1$. Such an insulator has two spin-polarized counter-propagating edge modes related by time reversal. The corresponding branches intersect at one of the TRIM, say, $ {k^\star} =\pi$. Note that at TRIM $H_{ {k^\star} }T = TH_{ {k^\star} }$, and level crossings at these points are protected by Kramers degeneracy. It follows that the pair of the edge states cannot be removed by a $T$-invariant perturbation. This model gives an example of the \defn{quantum spin Hall insulator} \cite{Kane2005}.

However, if we start from the Chern insulators with $c_{\spinup \spindown} =\pm 2$, we can easily gap out the crossings at $k\neq {k^\star} $, and push what remains of the edge modes into bulk bands. In general, we will end up with a single protected crossing when initial  $c_{\spinup \spindown}$ are an odd,  and with a trivial phase otherwise.

\begin{figure}[h]
\begin{center}
\psfrag{e}{$\varepsilon$}
\psfrag{pi}{$\pi$}
\psfrag{k}{$k$}
\psfrag{0}{$0$}
\includegraphics[scale  = 1.2]{img/ZigZag}
\caption{Two ways of connecting Kramers pairs in edge states spectrum: topological (left) and trivial (right). Arrows indicate possible deformations of the spectrum. Only half of the Brillouin zone is shown.
\label{Fig:ZigZag}}
\end{center}
\end{figure}

One can also obtain this classification starting from the spectrum of the edge states. Suppose we first select positions of mid-gap Kramers pairs at TRIM, as shown in Fig.~\ref{Fig:ZigZag}. What are compatible configurations of the edge modes? There are two inequivalent shapes the edge modes can take. First,
Kramers pairs can be connected by a single mode that goes between them in a zigzag fashion. Such a mode cannot be destroyed by any smooth symmetry-preserving transformation, and indicates a non-trivial phase. In the other scenario,  Kramers pairs are connected by corresponding pairs of edge modes. Such modes can be pushed into the bulk bands, so the insulator is trivial.

In contrast to Chern insulators, which are classified by integers, here we have only two classes: either an insulator is topological, or it is trivial. One says that $T$-invariant insulators have \defn{$\mathbb{Z}_2$-classification}. In group theory, $\mathbb{Z}_2$ denotes the additive group of integers modulo $2$. It has two elements, which correspond to even and odd integers, and addition is defined as $[e]+[o] = [o]$, etc. As we will see, in some situations $\mathbb{Z}_2$ class of an insulator can be related to the parity of the Chern number.

There is a number of ways to determine whether a given Hamiltonian $\hat H(k_x, k_y)$ describes a non-trivial $T$-invariant insulator. A minimal model of a gapped spinful system contains four bands. Thus, geometrically, we have a vector bundle of valence bands $V^v$ with a fiber isomorphic to $\mathbb{C}^2$. The bundle $V^v$ is necessarily trivial, since otherwise we would have a Chern insulator, which breaks time-reversal symmetry. However, there is a symmetry constraint on its possible sections, or  eigenstates of $\hat H(k_x, k_y)$ over the Brillouin zone. If there is no non-vanishing section that satisfies symmetry constraints, the insulator is non-trivial. This can be detected by several equivalent topological invariants. For detailed description of these invariants and proof of their equivalence, see Ref.~\cite{Fruchart2013}. Another way is to track coordinates of Wannier centers in one direction as functions of crystal momentum along the other. If there are counter-propagating flows of Wannier centers related by time reversal, we have a non-trivial phase \cite{Gresch2016}.

\subsubsection{Two dimensions: classification by inversion eigenvalues}
The classification of $T$-invariant insulators takes the simplest form for crystals that also have inversion symmetry $\mathcal{I}$. Recall how the  symmetries act on real space coordinate, crystal momentum, and spin state:
\begin{center}
$\begin{tabu}{l  c c   r}
   & x & k & s \\
T & + & - & - \\
\mathcal{I} & - & - & +
\end{tabu}$
\end{center}

First, let us consider effect of the inversion symmetry on a spinless Chern insulator. Because of the symmetry, Wannier center trajectory in the mixed position-momentum  space must be invariant under inversion $(x, k_y) \to (-x, -k_y)$. In particular, Wannier centers at $ k^\star _y=0, \pi$ are pinned to inversion centers of the crystal in the real space, $x = 0, \frac{1}{2}$, as shown in Fig.~\ref{Fig:2DInv} on the left.

\begin{figure}[h]
\begin{center}
\psfrag{pi}{$\pi$}
\psfrag{kx}{$k_x$}
\psfrag{ky}{$k_y$}
\psfrag{x}{$x$}
\psfrag{12}{$0.5$}
\psfrag{0}{$0$}
\includegraphics[scale=1.2]{img/2DInv}
\caption{Wanier center trajectory (left) and inversion eigenvalues (right) for inversion-symmetric Chern insulator. Filled (empty) circles denote eigenvalues $\lambda = +1 (-1)$.
\label{Fig:2DInv}}
\end{center}
\end{figure}

Note that in the momentum space $(k_x, k_y)$, the lines $ k^\star _y = 0, \pi$ are invariant under inversion (although are not pointwise fixed). The restrictions of the Bloch Hamiltonian to these lines
\begin{equation}
H(k_x, 0), \qquad H(k_x, \pi)
\end{equation}
describe inversion-symmetric chains. Depending on the inversion eigenvalues, each of these chains can be either in trivial or polarized states. If the states of two chains are different, symmetry forces Wannier center to shift by odd number of unit cells per ``pumping cycle'' in $k_y$. Otherwise, when two chains are in the same state, charge centers shift by even number of cells. Two situations are distinguished by the value of product $\nu$ of inversion eigenvalues at four inversion-symmetric points. We conclude that the number $\nu=\pm 1$ detects the \emph{parity} of the Chern number $c(V^v)$:
\begin{equation}
\nu = \prod_{k\in\{( k^\star _x,  k^\star _y)\}} \lambda(k) = (-1)^{c(V^v)}.
\end{equation}

With this in mind, we can use two inversion-symmetric Chern insulators with $c_{\spinup \spindown}=\pm n$ to construct a $T$-invariant insulator. As above, odd value of $n$ corresponds to the topological phase, while even $n$ gives a trivial insulator. But now we can determine the parity of $n$ by calculating the bulk invariant, $\nu$. In doing so, it is enough to use inversion eigenvalues for any of two occupied bands. Indeed, since time reversal commutes with inversion, $[T, \mathcal{I}]=0$, both states that form a Kramers pair at $ {k^\star} $ must have the same inversion eigenvalue.

\subsubsection{Three dimensions}
Finally, we apply this approach to a 3D $T$-invariant crystal with inversion symmetry. There are eight TRIM in the Brillouin zone, and three pairs of parallel planes that are invariant under the  symmetries. Restriction of the Hamiltonian $ H(k_x, k_y, k_z)$ to any of these planes gives a Hamiltonian of 2D $T$-invariant insulator. We define $\nu$ as a product of eight inversion eigenvalues of one of the occupied bands:
\begin{equation}
\nu = \prod_{k \in \{( k^\star _x,  k^\star _y,  k^\star _z)\}} \lambda(k).
\end{equation}
Depending on the $\lambda( {k^\star} )$, some of the six effective 2D insulators can be trivial or topological. If all 2D insulators are trivial, we have a trivial 3D insulator. Otherwise, there are two cases:
\begin{itemize}
\item $\nu=1$: in each pair of parallel $T$-invariant planes,  either both 2D insulators are trivial, or both are topological.
\item $\nu=-1$: in each pair of parallel $T$-invariant planes, one insulator is trivial and the other is topological.
\end{itemize}
In both cases, the resulting 3D insulators are non-trivial, but they have distinct character. An insulator with $\nu=1$ can be obtained by stacking layers of non-trivial 2D insulators along some crystallographic direction, and is known as a \defn{weak} topological insulator. On the contrary, the \defn{strong} topological insulator with $\nu=-1$ is a truly three-dimensional phase with no 2D counterparts.

\begin{figure}[h]
\begin{center}
\psfrag{e}{$\varepsilon$}
\psfrag{pi}{$\pi$}
\psfrag{kx}{$k_x$}
\psfrag{ky}{\hspace{-0.1cm}$k_y$}
\psfrag{kz}{$k_z$}
\includegraphics[scale=1.2]{img/3DCone}
\caption{Inversion eigenvalues at TRIM of the 3D Brillouin zone and spectrum of surface states at the surface orthogonal to the $x$ direction. Filled (empty) circles denote eigenvalues $\lambda = +1 (-1)$.
\label{Fig:3DCone}}
\end{center}
\end{figure}

A possible combination of inversion eigenvalues of a strong topological insulator is shown in  Fig.~\ref{Fig:3DCone}. From inversion eigenvalues, one can build a compatible spectrum of surface states. For example, consider the surface orthogonal to $x$ direction. There are four $T$-invariant planes in the Brillouin zone that are parallel to $x$ direction, given by
\begin{equation}
k_y=0 \qquad k_z =0 \qquad k_y=\pi \qquad k_z=\pi.
\end{equation}
According to the eigenvalues, the first pair of planes describe trivial 2D insulators, while the second pair correspond to topological ones. We put typical spectra of edge states at the corresponding lines on the surface Brillouin zone, connect them inside the quarter and finally extend to the full range of the crystal momenta. In this way, we obtain the \defn{Dirac cone}, the hallmark of 3D $T$-invariant strong topological insulators. Each surface of the crystal will have a similar cone in the spectrum. Note that the states at the opposite points on the circumference of the cone must have opposite spin directions, as they are related by time-reversal symmetry. Such spin polarization of the Dirac cone has indeed been observed in experiments.

\subsection{Symmetries in general}
Classification of topological matter with symmetries is a vast topic, and we have barely scratched its surface. We conclude this section by a brief discussion of the general situation.

Time-reversal symmetry (TR) is one of fundamental quantum mechanical symmetries, along with particle-hole symmetry (PH) and chiral symmetry (C), which is the combination of TR and PH. Each of the symmetries may be  present or absent, and the first two, when present, can square to plus or minus identity operator. There are various real-space dimensions of the systems and ten possible combinations of symmetries. The information of classification scheme for each dimension and  symmetry class fits into \defn{periodic table} of topological insulators and superconductors \cite{Kitaev2009}. The part of the table that describes phases discussed above reads
\begin{center}
\begin{tabular}{c c c | c c c}
TR & PH & C & $1$ & $2$ &$3$ \\ \hline
$0$ & $0$& $0$&$0$&$\mathbb{Z}$&$0$\\
$-1$& $0$& $0$& $0$  &$\mathbb{Z}_2$&$\mathbb{Z}_2$\\
\end{tabular}
\end{center}
The first row corresponds to Chern insulators, which break all symmetries, exist only in two-dimensional systems and have $\mathbb{Z}$-classification. The time-reversal invariant topological insulators with spinful electrons have $\mathbb{Z}_2$-classification in 2D and 3D. The classification given by the table appears in several physical and mathematical contexts, as discussed in Refs.~\cite{Stone2010} and \cite{Zirnbauer2015}.

The situation becomes much more complicated when one takes into account spatial crystalline symmetries. In mathematics, the classification of vector bundles is the subject of K-theory. The symmetries can be thought of as additional constraints, which relate fibers of the vector bundle at different points of the base space. Classification of vector bundles with symmetries by K-theoretic methods is a profoundly difficult problem that is to be solved on case-by-case basis.

Fortunately, there is a ``physical way'' to obtain such classification in physically relevant cases, known as \defn{topological quantum chemistry} (see Ref.~\cite{Cano2020} for a recent overview). In a nutshell, it has two key ingredients: band structure combinatorics and induction of bands from the real-space atomic orbitals. The idea is to enumerate first all possible band structures compatible with a given symmetry, by means of the representation theory. There are special subspaces of the Brillouin zone (points, lines, planes), which are fixed under the action of the symmetry group or one of its subgroups. Then the eigenstates at the special points must transform according to the irreducible representations of the corresponding (sub)group. Furthermore, representation theory provides the rules of how such representations can split and combine along the various paths in the Brillouin zone. As a result, the list of all possible band structures comes as a solution to the combinatorial problem of assigning irreducible representations to the special points in a way compatible with the rules.

But which of the band structures are topological? Since we have already listed all possibilities, this is equivalent to asking which of them are trivial. The answer comes form the localized Wannier orbitals. It turns out that one can take a real-space distribution of atomic orbitals compatible with the symmetry and produce the corresponding band structure in momentum space, via the procedure known as induced representation. This gives the list of all band structures that correspond to the localized Wannier functions, and thus are trivial by construction. Then all band structures that are in the first list and do not appear in the second one must represent topologically non-trivial phases compatible with the symmetry. Remarkably, the classification based on band structure combinatorics agrees with that given by K-theoretic methods whenever the latter is available (see, for example, Ref.~\cite{Kruthoff2017}).

\section{Topological semimetals}\label{Chap:Semi}
\defn{Semimetals} are gapless materials with vanishing density of states at the Fermi level. For Weyl  semimetals, this happens because valence and conduction bands have conical intersections at some points. If all such crossings occur at the Fermi energy, the Fermi surface will consist of isolated points. Weyl semimetals are closely related to Chern insulators: the band crossings act as sources and sinks of Berry curvature and lead to the appearance of surface states.

In this section, we discuss topological properties of Weyl semimetals and crystals with symmetry-enforced band crossings. But first, we consider another famous two-dimensional semimetal.

\subsection{Graphene}
\subsubsection{Brillouin zone and band structure}
The simplest tight-binding model of graphene includes two orbitals per unit cell, as shown in Fig.~\ref{Fig:Graphene}. The Bravais lattice is hexagonal, while the atomic sites form the honeycomb lattice. We label unit cells by vectors $\vec m = m_1 \vec a_1 + m_2 \vec a_2$, where $\{\vec a_i\}$ is a  real space basis. In a similar way, any vector in reciprocal space is decomposed in terms of the basis $\{\vec b_i\}$ as $\vec k = k_1\vec b_1+k_2\vec b_2$. The Bloch states  are defined as
\begin{equation}
|\alpha_{\vec k}\rangle = \frac{1}{\sqrt{N}} \sum_{\vec m} e^{i\vec k\cdot \vec m} |\alpha_{\vec m}\rangle,
\end{equation}
where the dot product means
\begin{equation}\label{Eq:DotRR}
\vec k \cdot \vec m = \sum_i k_i m_i.
\end{equation}

\begin{figure}[h]
\begin{center}
\psfrag{ex}{$\vec e_x$}
\psfrag{ey}{$\vec e_y$}
\psfrag{a1}{$\vec a_1$}
\psfrag{a2}{$\vec a_2$}
\psfrag{b1}{$\vec b_1$}
\psfrag{b2}{$\vec b_2$}
\psfrag{k1}{$k_1$}
\psfrag{k2}{$k_2$}
\psfrag{a}{$a$}
\psfrag{b}{$b$}
\psfrag{sx}{$\sigma_x$}
\psfrag{sy}{$\sigma_y$}
\psfrag{Dp}{$\vec D_+$}
\psfrag{Dm}{$\vec D_-$}
\psfrag{2p3}{$\frac{2\pi}{3}$}
\psfrag{4p3}{$\frac{4\pi}{3}$}
\includegraphics{img/Graphene}
\caption{Left: crystal lattice of graphene. Black lines show hexagonal Bravais lattice. Heavy gray lines form the honeycomb lattice and indicate hopping amplitudes $t$. Right: Brillouin zone of graphene with Dirac points $\vec D_\pm$. Gray arrows show linearized Bloch Hamiltonian near Dirac points.
\label{Fig:Graphene}}
\end{center}
\end{figure}

Let us express $\vec a_i$ in terms of orthonormal Cartesian basis $\{\vec e_x, \vec e_y\}$:
\begin{equation}
 \vec a_{1,2} = \begin{pmatrix}
 \pm \tfrac{1}{2}\\  \tfrac{\sqrt 3}{2}     \\
 \end{pmatrix},
\end{equation}
where lattice constant is set to unity. One can also associate $\vec b_i$ with real-space vectors, and interpret the dot product above as the standard inner product on the plane. For example, let us find the direction of $\vec b_1$. The corresponding Bloch wave $e^{i\vec b_1 \cdot \vec m} = e^{ik_1m_1}$ has the wave front along $\vec a_2$. Since $\vec b_1$ plays the role of wave vector, we find that $\vec b_1 \perp \vec a_2$, and similarly $\vec b_2 \perp \vec a_1$. This can also  be understood algebraically: Eq.~\eqref{Eq:DotRR} holds only if $\vec b_i\cdot \vec a_j = \delta_{ij}$. We have
\begin{equation}
\vec b_{1,2} = \begin{pmatrix}
 \pm 1\\
 \tfrac{1}{\sqrt 3}\\
 \end{pmatrix},
\end{equation}
and the Brillouin zone assumes the shape shown in Fig.~\ref{Fig:Graphene}.

Now we construct the tight-binding model. Suppose that only nearest-neighbor hopping is present, with hopping amplitude $t\in \mathbb{R}$. For example, hopping from $b$ site of a cell with coordinate $\vec m$ to the neighboring $a$ sites is described by three terms:
\begin{equation}
t\biggl(|a_{\vec m}\rangle\langle b_{\vec m}| +  |a_{\vec m+\vec a_1}\rangle\langle b_{\vec m}|+  |a_{\vec m+\vec a_2}\rangle\langle b_{\vec m}| \biggr).
\end{equation}
In momentum space this becomes
\begin{equation}
H^{ba}_{\vec k} = t\bigl(1+e^{-i\vec k\cdot \vec a_1} +e^{-i\vec k\cdot \vec a_2} \bigr) =
t\bigl(1+e^{-i k_1} +e^{-i k_2} \bigr) \equiv f_{\vec k}.
\end{equation}
Thus the Bloch Hamiltonian in $\{|a_{\vec k}\rangle, |b_{\vec k}\rangle\}$ basis reads:
\begin{equation}\label{Eq:GrapheneH}
H_{\vec k} =\begin{pmatrix}
0& f_{\vec k}\\
\overline{f_{\vec k}}& 0\\
\end{pmatrix}.
\end{equation}

There are two special points $\vec D_\pm$ in the Brillouin zone, in which the Hamiltonian is gapless,  $H_{\vec D_\pm}=0$. Their coordinates are easily found to be
\begin{equation}
\vec D_\pm =\pm \frac{2\pi}{3}\begin{pmatrix}
1\\-1\\
\end{pmatrix}.
\end{equation}
These points are known as \defn{Dirac points} because of the conical form of the band touchings, reminiscent of the linear dispersion relation of massless relativistic particles.

Let us find the expansion of the Hamiltonian around Dirac points $\vec D_\pm$ in terms of a long-wavelength parameter $|\vec q|\ll 2\pi$. We start with the expression
\begin{equation}
f_{\vec D_\pm +\vec q} = t\biggl( 1 + \sum_j e^{-i(\vec D_\pm +\vec q)\cdot\vec a_j}\biggr)
\end{equation}

\exr{Show that expansion of $f_{\vec D_\pm +\vec q}$ to linear order in $\vec q$ reads \begin{equation}
f_{\vec D_\pm +\vec q}\approx -\frac{t\sqrt 3}{2}( \pm \vec q\cdot \vec e_x - i \vec q \cdot \vec e_y).
\end{equation}
}

Denoting dot products in the last expression  by $q_x$ and $q_y$, we obtain the Hamiltonian
\begin{equation}
H_{\vec D_\pm +\vec q} \approx  -\frac{t\sqrt 3}{2}(\pm q_x \sigma_x +q_y \sigma_y),
\end{equation}
which is linear in terms of Pauli matrices.

\subsubsection{Tuning graphene into Chern insulator}
While graphene itself is not a topological phase, there is a property of its Hamiltonian, which is topological in nature. Since this model of graphene includes only hopping between different sublattices, the Hamiltonian does not have terms proportional to $\sigma_z$. The Hamiltonian is thus described by a planar vector field with periodic boundary conditions. We can think about this vector field as a section of a trivial real vector bundle with two-dimensional fiber. Then topology requires that the sum of indices of singularities vanishes, which is indeed the case. In this context, one calls the index of singularity at the Dirac point a \defn{winding number} and its sign is known as \defn{chirality}.

Now let us see what happens when we add $\sigma_z$ terms, which open the gap at the Dirac points and turn graphene into an insulator. The easiest way to do so is to add staggered on-site potentials $\pm \Delta$ with opposite signs at different sublattices. Physically this is realized in boron nitride, which has the same honeycomb lattice as graphene. The Hamiltonian becomes
\begin{equation}
H_{\vec k}^{BN} = \begin{pmatrix}
\Delta & f_{\vec k}\\
\overline{f_{\vec k}} & -\Delta\\
\end{pmatrix}.
\end{equation}
To determine the Chern number, consider the Hamiltonian as a map $h: T^2 \to S^2$ from the Brillouin zone torus to the Bloch sphere. For graphene, the map is defined everywhere except Dirac points, and the image coincides with the equator. Once the constant term $\Delta \sigma_z$ is added, the image becomes a cap above the equator. Clearly, it does not cover the full sphere, so the Chern number must be zero. One can also deduce this from the Eq.~(\ref{Eq:ChernPBSum}): pre-images of the north pole are two Dirac points, and Jacobian determinant has opposite signs at $\vec D_+$ and $\vec D_-$. To understand the triviality in physical terms, note that in the limit $\Delta \gg t$ the system becomes an atomic insulator that consists of decoupled orbitals.

\begin{figure}[h]
\begin{center}
\includegraphics{img/HaldaneSym}
\caption{Haldane model. Black arrows indicate hopping amplitudes $it_2$. Inversion center is marked by cross. Dashed line connects two orbitals related by inversion.
\label{Fig:HaldaneSym}}
\end{center}
\end{figure}

To open the gap in topologically non-trivial way, we follow Haldane \cite{Haldane1988} and add imaginary hopping $it_2$ between next-nearest neighbors of the same type, as shown in Fig.~\ref{Fig:HaldaneSym}. Phase accumulation after hopping around each loop indicates the presence of magnetic field. However, there is no phase associated with going around the hexagon of the honeycomb lattice. Thus the unit-cell average of magnetic field must be zero.

The hopping elements between $a$ sites are:
\begin{equation}
it_2\biggl(|a_{\vec m+\vec a_1}\rangle\langle a_{\vec m}| +
|a_{\vec m-\vec a_2}\rangle\langle a_{\vec m}|+
|a_{\vec m+\vec a_2 -\vec a_1}\rangle\langle a_{\vec m}|\biggr) + h.c.
\end{equation}
In momentum space one has
\begin{equation}
H^{aa}_{\vec k}  =  2 t_2 (\sin k_1 -\sin k_2 + \sin (k_1-k_2) ) \equiv \Delta_{\vec k},
\end{equation}
so the Hamiltonian reads
\begin{equation}
H_{\vec k}^{HM} = \begin{pmatrix}
\Delta_{\vec k}&  f_{\vec k}\\
\overline{f_{\vec k}} & -\Delta_{\vec k}\\
\end{pmatrix}.
\end{equation}
Crucially, the $\sigma_z$ component is now $k$-dependent. It is anti-symmetric function of momentum, and in particular  $\Delta_{\vec D_+} = -\Delta_{\vec D_-}$. Thus the north pole of the Bloch sphere has now a single pre-image, and $H_{\vec k}^{HM}$ describes a Chern insulator. This is the celebrated \defn{Haldane model} of the quantum Hall effect with zero net magnetic field.

\subsubsection{Symmetry considerations}
Among the three models considered above, graphene Hamiltonian \eqref{Eq:GrapheneH} is the most symmetric one. Since $f_{\vec k} = \overline{f_{-\vec k}}$, the Hamiltonian satisfies Eq.~\eqref{Eq:TRHam} and has time-reversal symmetry $T$. It is also invariant under inversion with the center in the middle of a unit cell, as shown in  Fig.~\ref{Fig:HaldaneSym}. We put the inversion center into the unit cell with $\vec m=0$, so the symmetry is represented by $\mathcal{I} |a_{\vec m}\rangle = |b_{-\vec m}\rangle$.

In momentum space, the Hamiltonian satisfies
\begin{equation}
\sigma_x H_{\vec k} \sigma_x = H_{-\vec k}.
\end{equation}
The combined $\mathcal{I} T$ symmetry acts in momentum space locally:
\begin{equation}
\sigma_x \overline{H_{\vec k}} \sigma_x  = H_{\vec k}.
\end{equation}
In the space of Pauli matrices, the transformation on the left is the combination of reflection in the $xz$ plane with $\pi$ rotation around $x$ axis. The result is reflection in the $xy$ plane. Thus the $\mathcal{I} T$ symmetry requires $h_z$ component to vanish everywhere.

The Hamiltonian of boron nitride retains the time-reversal symmetry, but inversion symmetry is broken by the staggered on-site potential. Hamiltonians at $\vec k$ and $-\vec k$ are related by reflection in the $xz$ plane, so the $h_z$ components at the Dirac points must be the same. In contrast, the Haldane model has the inversion symmetry, which can be seen from Fig.~\ref{Fig:HaldaneSym}:  inversion preserves direction of arrows, or the sign of complex hopping amplitudes. However, it breaks time reversal since $\Delta_{\vec k} \neq \overline{\Delta_{-\vec k}}$. The Hamiltonians at opposite momenta are related by $\pi$ rotation about $x$ axis, and $h_z$ components at Dirac points have opposite signs.

\subsection{Weyl semimetals}
Weyl semimetal can be thought of as a three-dimensional generalization of graphene. Here we start with the dimension-counting argument underlying this analogy, and then discuss the surface states of Weyl semimetals.

\subsubsection{Dimension of zero locus of a section}
\begin{figure}[h]
\begin{center}
\psfrag{g1}{\textcolor{half-gray}{1}}
\psfrag{g2}{\textcolor{half-gray}{2}}
\psfrag{g0}{\textcolor{half-gray}{0}}
\psfrag{ge}{\textcolor{half-gray}{$\varnothing$}}

\psfrag{1}{1}
\psfrag{2}{2}
\psfrag{3}{3}
\psfrag{n}{$n$}
\psfrag{d}{$d$}

\includegraphics[scale=1.1]{img/ZeroLoci}
\caption{Zero loci of sections of real vector bundles. Numbers indicate dimension of base space $d$, dimension of fiber $n$ and dimension of zero locus (gray numbers).
\label{Fig:ZeroLoci}}
\end{center}
\end{figure}

Let $\vec v$ be a section of a real vector bundle over a base $\mathcal{B}$. We are interested in a generic shape of zero locus of $\vec v$, that is, the set of points $p\in \mathcal{B}$ such that $\vec v(p)=0$. By ``generic'' one means that such zeroes cannot be eliminated by a smooth local deformation of $\vec v$. The dimension of this space depends on the combination of dimension of the base $d = \dim \mathcal{B}$ and dimension of the fiber $n$. Fig.~\ref{Fig:ZeroLoci} shows typical shapes of zero loci for various pairs $(n,d)$. One can see that dimension of zero locus is given by $d-n$ if $d \geqslant n$ and there are no stable zeroes for $d < n$. We have encountered some of this situations before. In examples below we treat complex line bundles as real plane bundles.
\begin{center}
\begin{tabular}{cl}
$(n, d)$ & Examples: \\
$(1, 1)$ & singularity on a M\"obius band in Fig.~\ref{Fig:RTopoBun}\\
$(2, 1)$ & complex function of one variable does not have stable zeroes \\
$(1, 2)$ & real function of two variables vanishes along one-dimensional curves \\
$(2, 2)$ & Dirac points in graphene, singularities on $TS^2$\\
$(3, 2)$ & opening the gap in graphene by breaking symmetries \\
$(2, 3)$ & nodal line of zeroes of wave function near magnetic monopole \\
$(3, 3)$ & Berry monopoles in Fig.~\ref{Fig:TopTransition}.
\end{tabular}
\end{center}

Suppose now that we have a crystal with a pair of isolated bands separated by Fermi energy. If no symmetry constraints are present, these bands are described by a general two-level Bloch Hamiltonian, or a three-component vector field $\vec h_{\vec k}$. The gap between bands closes whenever $|\vec h_{\vec k}|=0$. It follows from the discussion above that one needs three dimensions of parameter space to have stable point-like zeroes of $|\vec h_{\vec k}|$. This is the idea behind \defn{Weyl semimetal}: in a three-dimensional crystal, crossing points of a an isolated pair of bands cannot be removed by a local perturbation. In this way, Weyl semimetal can be thought of as a three-dimensional generalization of graphene, where symmetry requires that $\vec h_{\vec k}$ lies in the plane and thus stabilizes Dirac points.

\subsubsection{Weyl points and Fermi arcs}
 Let us construct a hypothetical two-band Hamiltonian of a Weyl semimetal. We start form the planar vector field $\vec h_{\vec  k}$ with two singularities, shown in Fig.~\ref{Fig:Weyl} on the left. Imagine that this is a three-dimensional vector field with cylindrical symmetry about the horizontal axis that goes through the singular points. Now the singularities become \defn{Weyl points} of opposite chirality.

Note that the Hamiltonian component proportional to the unit matrix $\mathbb I$, which controls overall shift of energy levels is also important for semimetallic behavior. For example, high anisotropy of such term near a Weyl point can tilt the conical band intersection so strongly that the cone generator will cross Fermi level, resulting in \defn{type-II} Weyl point \cite{Soluyanov2015}.

Each Weyl point is characterized by the Chern number of the eigenstate bundle over the sphere that contains the point. In Fig.~\ref{Fig:Weyl}, we have $c=1$ for the left point, while near the right one the map $h$ reverses orientation, and the Chern number is $c=-1$ (compare with the pullback of $TS^2$ under reflection in $yz$ plane). Thus the Weyl points act as sources and sinks of the Berry flux. In our picture, two points are connected by the flux tube with skyrmion-like field structure in the cross section.

\begin{figure}[h]
\begin{center}
\psfrag{Wp}{$\vec W_+$}
\psfrag{Wm}{$\vec W_-$}
\psfrag{kx}{$k_x$}
\psfrag{ky}{\hspace{-0.1cm}$k_y$}
\psfrag{kz}{$k_z$}
\psfrag{ep}{$\varepsilon$}
\includegraphics[scale=1.2]{img/Weyl}
\caption{Left: vector field $\vec h_{\vec k}$ that contains two Weyl points. Right: formation of the Fermi arcs in a Weyl semimetal.
\label{Fig:Weyl}}
\end{center}
\end{figure}

Now we place this $\vec h_{\vec k}$ configuration inside the torus $T^3$ of the Brillouin zone of a three-dimensional crystal, as shown in Fig.~\ref{Fig:Weyl} on the right. We further select two cylinders that contain Weyl points, and a plane between them. Because of periodic boundary conditions on the faces of Brillouin zone, each of these surfaces is a two-dimensional torus $T^2$. Note that they are pierced by the line of the Berry fulx. It follows that restriction $\vec h_{\vec k}\rvert_{T^2}$ to any of these tori describes a Hamiltonian of a Chern insulator. Thus, the crystal must support edge modes on the top and bottom surfaces. Above the top surface, we plot the possible spectra $\varepsilon (k_x, k_y)$ of the edge modes of these effective Chern insulators. Similar modes will be present at the boundary circles of any cylinder that surrounds a Weyl point. Each mode intersects the Fermi level in one point. By continuity, these points form together a \defn{Fermi arc}, or a Fermi surface that consists of an open curve that connects projections of Weyl points. There will be another Fermi arc on the bottom surface, and also a pair of arcs on the surfaces orthogonal to $y$ direction. However, there need not  be any arcs on the left and right surfaces, since projections of Weyl points to the $yz$ plane coincide.

\subsection{Topology of local symmetry constraints}\label{Sec:TopoSymCons}
For a two-band model, Bloch Hamiltonian is represented by a vector field $\vec h_{\vec k}$ with values in the space of Pauli matrices. More generally, $\vec h_{\vec k}$ can be thought of as a section of a real vector bundle: plane bundle for graphene, or bundle of three-dimensional spaces for a Weyl semimetal. In both cases, vector bundles are trivial, since there are global basis sections, $\{\sigma_x, \sigma_y\}$ and $\{\sigma_x, \sigma_y, \sigma_z\}$, respectively. However, for certain symmetries, the space of all compatible Hamiltonians assumes the form of a non-trivial vector bundle.

\subsubsection{Non-symmorphic symmetry in 1D}

\begin{figure}[h]
\begin{center}
\psfrag{a}{$a$}
\psfrag{b}{$b$}
\psfrag{t}{$t$}
\psfrag{k}{$k$}
\psfrag{0}{$0$}
\psfrag{pi}{$\pi$}
\psfrag{sx}{$\sigma_x$}
\psfrag{sy}{$\sigma_y$}
\includegraphics{img/NS}
\caption{Non-symmorphic diatomic chain and the space of compatible Hamiltonians.
\label{Fig:NS}}
\end{center}
\end{figure}

A symmetry operation is called \defn{non-symmorphic} if it includes a translation by a fraction of the lattice constant. For example, the chain shown in Fig.~\ref{Fig:NS} on the left is invariant under reflection with respect to the dashed line, combined with translation by a half of the unit cell. In real space, the symmetry acts as
\begin{equation}
\hat G |a_m\rangle = |b_m\rangle \qquad \qquad
\hat G |b_m\rangle = |a_{m+1}\rangle
\end{equation}
In momentum space, the matrix of symmetry operator becomes
\begin{equation}
G_k = \begin{pmatrix}
0 & e^{-ik} \\ 1 & 0
\end{pmatrix}.
\end{equation}
We rewrite it as
\begin{equation}
G_k = e^{-ik/2} \begin{pmatrix}
0& e^{-ik/2} \\ e^{ik/2} & 0 \\
\end{pmatrix} = i e^{-ik/2} \biggl(
\cos \frac{\pi}{2} \mathbb{I} +
\sin \frac{\pi}{2} \, \vec w_k\cdot  \frac{\vec \sigma}{i}\biggr),
\end{equation}
where $\vec w_{k} = (\cos \frac{k}{2} , \sin \frac{k}{2}, 0)^T$.

Unlike inversion symmetry, $\hat G$ acts in momentum space locally:
\begin{equation}
G_k H_k G_k^{-1}=H_k.
\end{equation}
The conjugation of Hamiltonian by $G_k$ describes $\pi$ rotation about the axis defined by $\vec w_k$. Thus the vector $\vec h_k$ must belong to this line.
Now observe that when $k$ varies from $0$ to $2\pi$, the line defined by $\vec w_k$ rotates by $\pi$. Thus, the space of compatible Hamiltonians has the shape of the M\"obius band, as shown in Fig.~\ref{Fig:NS} on the right. One can consider $\vec w_k$ as a (discontinuous) basis section, and express the Hamiltonian as
\begin{equation}
\vec h_k = f_k \vec w_k,
\end{equation}
where $f_k$ satisfies $f_{k+2\pi} = -f_k$.

\exr{Find the function $f_k$ for the chain shown in Fig.~\ref{Fig:NS} with complex nearest-neighbor hopping amplitude $t = t_1 + i t_2$.}

Since the vector field that corresponds to the Hamiltonian is a section of a M\"obius band, the gap must close at least in one point. In this way, non-symmorphic symmetry enforces semimetallic behavior. Alternatively, one can see this by observing the eigenvalues of $G_k$ \cite{Zhao2016}. From $G_k^2 = \mathbb I e^{-ik}$ and $\tr G_k=0$, one finds that eigenvalues are $\lambda_{1,2} = \pm e^{-i\frac{k}{2}}$. As $k$ varies through the Brillouin zone, the eigenvalues switch places. Since one can label Hamiltonian eigenstates with $\lambda_i$, the bands  must also change their order in energy, and gap must close somewhere.

\subsubsection{Space-time inversion symmetry in 2D}
Consider now a two-band model with \defn{space-time inversion} symmetry: the system is invariant under combined action of inversion $\mathcal{I}$ and time reversal $T$, while both symmetries can be broken individually.
Note that both symmetries send $k\to -k$, so their combination acts in momentum space locally.

The lattice is shown in Fig.~\ref{Fig:IT} on the left.
We place inversion center on the $a$ site of the  unit cell with coordinates $(m_x, m_y) = (0, 0)$. Inversion acts on the orbitals differently:
\begin{equation}
\mathcal{I} |a_{m_x m_y}\rangle = |a_{-m_x \,\,-m_y}\rangle \qquad\qquad
\mathcal{I} |b_{m_x m_y}\rangle = |b_{-m_x-1\,\, -m_y-1}\rangle.
\end{equation}
For site of type $a$ we have, as usual, $\mathcal{I}: \vec m \to -\vec m$, while for $b$ sites there is an additional shift.

\begin{figure}[h]
\begin{center}
\psfrag{a}{$a$}
\psfrag{b}{$b$}
\psfrag{kx}{$k_x$}
\psfrag{ky}{$k_y$}
\psfrag{mx}{$m_x$}
\psfrag{my}{$m_y$}
\psfrag{0}{$0$}
\psfrag{pi}{$\pi$}
\psfrag{sx}{$\sigma_x$}
\psfrag{sy}{$\sigma_y$}
\includegraphics{img/IT}
\caption{Two-dimensional crystal with space-time inversion symmetry $\mathcal{I} T$. Left: checkerboard lattice. Dotted lines connect sites related by inversion symmetry. Arrows show pure imaginary hopping amplitudes, which break $\mathcal{I}$ and $T$ individually, but are invariant under $\mathcal{I} T$. Right: space of symmetry-compatible terms of Hamiltonian proportional to $\sigma_x$ and $\sigma_y$.
\label{Fig:IT}}
\end{center}
\end{figure}

The matrix of inversion in Bloch basis reads
\begin{equation}
I_k = \begin{pmatrix}
1 & 0 \\
0& e^{i(k_x+k_y)}
\end{pmatrix} =
e^{i\frac{k_x+k_y}{2}} \begin{pmatrix}
e^{-i\frac{k_x+k_y}{2}} & 0 \\
0& e^{i\frac{k_x+k_y}{2}}
\end{pmatrix}.
\end{equation}
Bloch Hamiltonian compatible with space-time inversion must satisfy
\begin{equation}
I_k \overline{H_k} I_k^{-1}  = H_k.
\end{equation}
In the space of Pauli matrices, transformation on the left is reflection in the $xy$ plane (complex conjugation) followed by $(k_x+k_y)$ rotation about $z$ axis (conjugation by $I_k$). It follows that at each $\vec k$, the vector $\vec h_{\vec k}$ must lie in the plane that contains $z$ axis and a line in the $xy$ plane from the family shown in Fig.~\ref{Fig:IT} on the right.

The Hamiltonian vector field $\vec h_{\vec k}$ is thus a section of a real plane bundle over a two-dimensional torus, so one can expect appearance of point-like zeros, or Dirac points. To trace their origin, we consider separately $(h_x, h_y)$ and $h_z$ components.

In-plane components $(h_x, h_y)$ correspond to the off-diagonal terms, which describe hopping between orbitals that belong to different sublattices. Real line bundle shown in Fig.~\ref{Fig:IT} has the shape of the M\"obius band in both directions. Hence, there must be at least one point of band crossing on each line given by $k_x = k_x^0$ or $k_y = k_y^0$. By continuity, these points form together nodal lines that cannot be removed by any symmetry-preserving perturbation of the $(h_x, h_y)$ components.

On-site potentials and hopping amplitudes inside sublattice give rise to diagonal terms, proportional to the identity matrix and $\sigma_z$. For instance, nearest-neighbor real horizontal hoppings $a\leftrightarrow a$ and $b\leftrightarrow b$ result in
\begin{equation}
h_z(k) = (t_a - t_b) \cos k_x.
\end{equation}
Unlike in-plane component, $h_z(k)$ need not vanish. But if it does, this happens again along some lines. We conclude that Dirac points $\vec h_{\vec k}=0$ will appear at the intersections of nodal lines of $(h_x, h_y)$ and those of $h_z$.

Recall that a vector bundle is  non-trivial, if one cannot define $n$ nowhere-vanishing linearly independent sections, where $n$ is the dimension of the fiber. Here, M\"obius-like twists of the in-plane component prevent one from defining the global basis, so the bundle is non-trivial. But in contrast with line bundles, in this case non-triviality does not require existence of singularities in any section. Indeed, $\vec h_{\vec k} = (0, 0, 1)^T$ is an example of gapped Hamiltonian compatible with symmetry. Smooth deformations of the Hamiltonian can create and annihilate pairs of Dirac points, similar to the situation in graphene and Weyl semimetals. The difference is that here it is impossible to define the chirality of Dirac points, or signs of winding numbers, in a global self-consistent way. The reason is that the bundle is \defn{non-orientable}: one cannot coherently choose orientation of planes.

The absence of the well-defined chirality can lead to apparent paradoxes of non-conservation of winding numbers in creation and annihilation processes for Dirac points. The authors of Ref.~\cite{Montambaux2018} consider such situation and resolve the problem by introducing ``winding vector'' as additional parameter of a Dirac point. In our terms, this vector is the plane normal, which allows one to define orientation of fibers locally.