A journey through the SSH model

# Exploring ideas in topological quantum phenomena

Anantha Hegde  Adarsh Kumar  Adhip Agarwala and Bhaskaran Muralidharan
111hegdeanantha@gmail.com,kumar.adarsh@iitb.ac.in,adhip<EMAIL_ADDRESS>
###### Abstract

Geared as an invitation for undergraduates, beginning graduate students, we
present a pedagogical introduction to one-dimensional topological phases – in
particular the Su-Schrieffer-Heeger model. In the process, we delve upon ideas
of entanglement using the correlator method and the von-Neumann density-matrix
method, geometric phase, polarization, transport signatures and the role of
electron-electron interactions. Through hands-on numerical experiments, whose
codes are shared, we try to drive home the message why a program of simulating
quantum electronics with topological toy models is the store house for
discovering fantastic physics ideas.

Anantha and Adarsh just graduated from IIT Bombay (India) with a B.Tech and an
M.Tech in Electrical Engineering. Adhip is a postdoctoral fellow at Max Planck
Institute for the Physics of Complex Systems, Dresden (Germany). Bhaskaran is
a Professor in the Department of Electrical Engineering at IIT Bombay (India).

July 2021 GENERAL ARTICLE

##  1 Phases of matter that surround us

The understanding of the various phases of matter which surround us is often
characterised by the way we react to them. For instance, it would hardly be
the best idea to touch a live wire which is attached to the electric socket.
That a material is a metal is often defined by its ability to conduct electric
current. Therefore while we move around, we intuitively classify phases around
us, things which are metallic (copper, the steel dining plates, the spoon) and
things which are insulating (that abhor current) such as a piece of wood,
glass etc. Things which stick to each other (magnets) and things that would
rather burn than to conduct heat (plastics or polymers, say over a metal).
However intuitive and obvious while these may look, the microscopic mechanisms
behind their behavior is rather intricate and has a rich display of quantum
mechanics which we often fail to appreciate. The following pages will
introduce you, the reader, to some of these fascinating stories which are
hidden behind seemingly plain phenomena. The write-up will assume you are an
undergraduate student who has had a first course in quantum mechanics at the
level of the introductory book by D.J. Griffiths and some exposure to solid-
state physics. The latter is not an essential pre-requisite though.

##  2 Taking the first shot

While every material has a unique set of atoms, bond-strengths, ionization
energies, stiffness etc. and therefore solving any material using the rule of
quantum mechanics is essentially impossible. However using the most simplest
of assumptions and keeping the minimal but essential ingredients one can still
understand a variety of phenomena. For instance a class of materials which we
understand in a relatively straightforward manner are metals and insulators.

First three momentum states for a particle in a one dimensional box

In quantum mechanics 101, we have solved for an electron in a potential box.
There, the laws of quantum mechanics dictate that each electron, instead of
being thought of as a particle should rather be thought of as a wave
characterized by a wavelength $\lambda$ (and the momentum vector
$\vec{k}=\frac{2\pi}{\lambda}\hat{\bm{r}}$). When we wish to fill in many
electrons we invoke Pauli exclusion principle – and demand that each of these
states (labelled by a distinct $\bm{k}=\frac{\pi n}{L}$ ($n=1,2,\ldots$))
should be occupied by just one electron. If one wishes to add more electrons
one needs to add them in a different momentum state. If our universe was just
a box, and all we had were non-interacting electrons – our universe could be
exactly solvable but also boring. Thankfully, even in a usual metal – that is
not the case. Electrons as they move in a material are constantly pulled and
pushed by the ions which surround them. In reality electrons also push and
pull each other – because they have the same electronic charge and repel.
Often their interactions cannot be ignored, one of the reasons many of us are
still doing research.

But, as is true for many of us when we appear for exams, we answer what we
know best – and see how good or bad it really is.

So let us assume electrons don’t interact – what does that mean? We can solve
for just one electron in a potential landscape which gives us all the one-
particle states it can occupy – and then when we have more electrons to worry
about, we just fill them one by one keeping in mind the Pauli exclusion
principle. If electrons did interact – addition of another electron could have
disturbed the existing ones and therefore we would have had to solve for them
again.

###  2.1 Modelling a solid

To model a solid, what potential landscape should one then use? Schrodinger’s
equation is a differential equation in space and time. And again solving this
is not very easy – specially when potentials are complicated. Now given a
material has a host of atoms each with different number of protons, one could
imagine that the electrostatic potential they generate for an electron can be
really weird. One could of course try to solve it numerically and obtain
stationary solutions, by discretizing the Schrodinger equation as follows:

$\displaystyle E\psi(\bm{r})$ $\displaystyle=$
$\displaystyle\Big{(}-\frac{\partial^{2}}{\partial^{2}\bm{r}}+V(\bm{r})\Big{)}\psi(\bm{r})$
(1) $\displaystyle E\psi(\bm{r}_{i})$ $\displaystyle=$
$\displaystyle-\frac{\psi(\bm{r}_{i+1})+\psi(\bm{r}_{i-1})-2\psi(\bm{r}_{i})}{\bm{r}_{i+1}-\bm{r}_{i}}+V(\bm{r}_{i})\psi(\bm{r}_{i})$
(2)

This leads to a matrix form of this equation which can be numerically
diagonalized to find the eigenvalues.

Even while the above method may just seem like a numerical convenience –
diagonalizing a matrix form of a Hamiltonian is a very standard practice in
quantum condensed matter and has very concrete physical motivations to do so.
This framework of solving Hamiltonians is called the tight-binding framework.
It introduces some assumptions in order to make our job of analysis simpler
for e.g. that the atoms which are placed in real space are sufficiently far
such that their atomic orbital wavefunctions don’t overlap.

More concretely, consider two atoms $A,B$ governed with a Hamiltonian $H$ with
their orbital wavefunctions $|\phi_{A}\rangle$ and $|\phi_{B}\rangle$. Even
while we assume $\langle\phi_{A}|\phi_{B}\rangle=0$, this does not imply that
$\langle\phi_{A}|H|\phi_{B}\rangle$ is zero. In fact, the term
$\langle\phi_{A}|H|\phi_{B}\rangle$ is crucially non-zero and is called the
hopping amplitude for an electron to move between site $A$ and site $B$. This
reflects the idea that electrons can freely move in a lattice potential
allowing for such hopping processes. The above way of defining our system
already assumes that a general solution for the wavefunction for this system
is a linear combination of the orbital wavefunctions on the two atoms
separately. Are such assumptions physically meaningful in a a realistic
setting? In order to make sense of this assumptions we consider two finite
potential wells separated by a barrier and analyze what happens to the ground
state and the first excited state as the width of the barrier is tuned (see
Example 1).

Example 1  Example: For an electron being shared between two atoms, is it fair
to assume (i) that the atomic orbitals of the atoms don’t overlap? and (ii)
that the complete solution is just a linear combination of individual
wavefunctions on the two atoms? Consider a one dimensional (D) system of two
finite potential wells (of depth $V_{0}$) each of width $a$ and separated by a
distance $b$ (see Fig. 1). We will call these atoms $A$ and $B$. A single
potential well of width $a$ centered at $x=x_{o}$ with a depth of $V_{0}$ has
a lowest energy state $\phi(x)$ which can be numerically obtained. Now our
orbital wavefunctions on atoms $A$ and $B$ would be shifted versions of
$\phi(x)$. We’d have $\phi_{A}(x)=\phi(x+(a+b)/2)$ and
$\phi_{B}(x)=\phi(x-(a+b)/2)$. Since $\phi(x)$ would have an exponentially
decreasing tail, the overlap $\langle\phi_{A}|\phi_{B}\rangle$ would decrease
as $b$ increases. What about the ground state for the double well system?
According to our assumption, a trial solution is
$\psi_{trial}=\frac{1}{\sqrt{2}}\Big{(}\phi_{A}(x)+\phi_{B}(x)\Big{)}$ (3)
Let’s try and calculate the actual solution. What will the ground state
$\psi_{0}$ (with energy $E$) of such a system be?  Figure 1: A double
potential well of depth $V_{0}$ in one dimension. Width of each well is $a$
and they are separated by a distance $b$. The system is centered at $x=0$.
The potential landscape divides the 1D space into five distinct regions.
Looking for solutions with $E<0$ such that in regions of potential $V=-V_{o}$
we have a standing wave, while for regions with $V=0$ we have a decaying
solution we assume the following form of the wavefunction $\displaystyle
Ae^{-kx}$ $\displaystyle\ \ \ in\ region\ I$ $\displaystyle
B\cos(px)-C\sin(px)$ $\displaystyle\ \ \ in\ region\ II$ $\displaystyle
D(e^{kx}+e^{-kx})$ $\displaystyle\ \ \ in\ region\ III$ $\displaystyle
B\cos(px)+C\sin(px)$ $\displaystyle\ \ \ in\ region\ IV$ $\displaystyle
Ae^{-kx}$ $\displaystyle\ \ \ in\ region\ V$ where
$k=\sqrt{\frac{-2mE}{\hbar^{2}}}$ and
$p=\sqrt{\frac{2m(V_{0}+E)}{\hbar^{2}}}$. Matching the wavefunction and its
derivative at the boundaries of the various regions, we arrive at the
following quantization condition
$k^{2}(1-e^{-kb})\sin(pa)-p^{2}(1+e^{-kb})\sin(pa)+2kp\cos(pa)=0$ (4) Looks
daunting, but let’s see what happens when $b\gg a$. In that case we can drop
the $e^{-kb}$ term as it would be much smaller compared to 1. Thus, we get
$k^{2}\sin(pa)-p^{2}\sin(pa)+2kp\cos(pa)=0$ Trying to solve this as a
quadratic equation in $k$, we get two possible conditions as $k=p\tan(pa/2)\ \
or\ \ k=-p\cot(pa/2)$ which are the quantization conditions for the even and
odd solutions of a single potential well! One can numerically solve eqn. (4)
for the ground state $\Psi_{0}$ and plot the wavefunctions for various values
of $b$ to see what is happening. Figure 2: The ground state $\Psi_{0}$ for
$b=0.1$nm and $b=4$nm. Width of each well is $a=1$nm and $V_{0}=1.14$eV.  The
code for fig. 2 shows the numerical calculations. An interesting exercise is
to check what happens to the overlap $\langle\psi_{trial}|\psi_{0}\rangle$ as
a function of $b/a$. We find that for $b\gg a$ the wavefunction for a double
potential well problem can be very approximated by a linear combination of two
separate potential wells when the distance between these ’atoms’ is larger
than the size of the atoms.

The above example shows that instead of solving for a differential equation in
all of space we can just worry about two coefficients in the tight-binding
limit and the Hamiltonian just becomes a $2\times 2$ matrix. In general, for a
system with $N$ sites, each having $m$ atomic orbitals, the Hamiltonian
becomes a $Nm\times Nm$ matrix. This is now a matrix problem which can be
diagonalized to get the energy eigenvalues for the electrons and the
corresponding wave functions. The hopping parameters which enter this
Hamiltonian matrix are often motivated from the chemistry of the material or
theoreticians often choose them to model a particular phenomena. In the
following examples (Example 2 and 3) we will visit this framework for some toy
problems.

Example 2  Example: Sweet spot for the electron
Let’s model an electron on a molecule of four sites forming a letter $T$.
Consider fermions hopping on four sites labelled by $1,2,\ldots,4$ and let the
hopping amplitude be $t$. The Hamiltonian then is determined by the hopping
process on the three bonds $H=t\Big{(}|1\rangle\langle 2|+|2\rangle\langle
1|\Big{)}+t\Big{(}|2\rangle\langle 3|+|3\rangle\langle
2|\Big{)}+t\Big{(}|2\rangle\langle 4|+|4\rangle\langle 2|\Big{)}$ (5) The
Hamiltonian matrix is then be given by:
$|1\rangle$ $|2\rangle$ $|3\rangle$ $|4\rangle$ $|1\rangle$ 0 t 0 0
$|2\rangle$ t 0 t t $|3\rangle$ 0 t 0 0 $|4\rangle$ 0 t 0 0 which can be
straightforwardly diagonalized. Before you diagonalize the matrix, where do
you think the electron will like to reside the most?

Example 3  Example: Electron in a ’one-dimensional’ solid. Consider a lattice
of 10 sites labelled by $1,2,\ldots,10$ with a distance of 1 unit between
them. We shall consider a hopping term $t$ only between neighbouring sites. In
order to make our life simpler, we will assume a periodic boundary condition
(PBC) and also consider a hopping term between site 1 and site 10. The
reasoning behind this assumption is model an infinite chain where any electron
will never encounter a boundary. The Hamiltonian then becomes
$H=t\Big{(}\sum_{j=1}^{9}|j\rangle\langle j+1|+|j+1\rangle\langle
j|\Big{)}+t\Big{(}|1\rangle\langle 10|+|10\rangle\langle 1|\Big{)}$ (6) The
Hamiltonian matrix is then given by: $|1\rangle$ $|2\rangle$ $|3\rangle$
$|4\rangle$ $|5\rangle$ $|6\rangle$ $|7\rangle$ $|8\rangle$ $|9\rangle$
$|10\rangle$ $|1\rangle$ 0 t 0 0 0 0 0 0 0 t $|2\rangle$ t 0 t 0 0 0 0 0 0 0
$|3\rangle$ 0 t 0 t 0 0 0 0 0 0 $|4\rangle$ 0 0 t 0 t 0 0 0 0 0 $|5\rangle$ 0
0 0 t 0 t 0 0 0 0 $|6\rangle$ 0 0 0 0 t 0 t 0 0 0 $|7\rangle$ 0 0 0 0 0 t 0 t
0 0 $|8\rangle$ 0 0 0 0 0 0 t 0 t 0 $|9\rangle$ 0 0 0 0 0 0 0 t 0 t
$|10\rangle$ t 0 0 0 0 0 0 0 t 0 How would we diagonalize this? Borrowing from
math (the $n_{th}$ roots of unity), you can check that one can write the
solution as
$|\psi\rangle=\sum_{j=1}^{10}\psi_{j}|j\rangle=\frac{1}{\sqrt{10}}\sum_{j=1}^{10}e^{ikj}|j\rangle$
(7) where $i=\sqrt{-1}$, and $k$ is a parameter as in the particle in a box
case. Now we need to find a condition on $k$ that would lead to a consistent
solution. In the eigenvalue equation, observe that for rows 2 to 9 of the
matrix, we get the following equation $te^{ik(j-1)}+te^{ik(j+1)}=Ee^{ikj}$
which suggests $E=2t\cos(k)$. However for the first row,
$te^{ik(10)}+te^{ik(2)}=Ee^{ik(1)}$ If we just had $e^{i10k}=1$ it would solve
our problem. Check that this takes care of the last row as well. This leads to
quantization of $k$ where $k$ can be set to $\frac{2\pi m}{10}$ where $m$
ranges from 1 to 10, and each $k$ gives a state with energy $2t\cos(k)$ . An
interesting observation to make is that, all the coefficients have the same
magnitude and hence the electron is equally likely to reside at any of these
10 sites. Since, we were initially looking for 10 eigenvalues and we have
found 10 values of $k$, we are done! Usually we often define $k$ within the
range of $-\pi$ to $\pi$ also referred to as the momentum space or the
Brillouin zone.

If you have worked out the examples, they may seem too simple to be of any
practical relevance – but there are good reasons to discuss them.

Let’s move to a bit more complicated compound - copper or silicon. But before,
we need a small help with notation. We are going to introduce a second
quantization notation [2ndq, altland_simons_2006]. where each site ($i$) of
the system can contribute two states – (i) vacuum $|0\rangle$ and (ii) one
with an electron $|1\rangle$. Operators called creation operators
$c^{\dagger}_{i}$ creates an electron on a vacuum state and annihilation
operators ($c_{i}$) kill the same.

$\displaystyle c^{\dagger}_{i}|0\rangle\rightarrow|1\rangle$ (8)
$\displaystyle c_{i}|1\rangle\rightarrow|0\rangle$ (9)

Thus another way of writing the Hamiltonian is where every row and column of
the matrix can be associated with creation and annihilation operators on a
particular position $i$. When an electron moves from position $i$ to $j$ one
could write this as a term.

$\text{electron moves from site $i$ to $j$}\rightarrow c^{\dagger}_{j}c_{i}$
(10)

So the Hamiltonians we discussed just now in Examples 2 and 3 can also be
written in second quantized form as

Second quantized form  for Example 2:
$H=t(c^{\dagger}_{1}c_{2}+c^{\dagger}_{2}c_{1})+t(c^{\dagger}_{2}c_{3}+c^{\dagger}_{3}c_{2})+t(c^{\dagger}_{2}c_{4}+c^{\dagger}_{4}c_{2})$
(11) for Example 3:
$H=t\Big{(}\sum_{j=1}^{9}c^{\dagger}_{j}c_{j+1}+c^{\dagger}_{j+1}c_{j}\Big{)}+t(c^{\dagger}_{1}c_{10}+c^{\dagger}_{10}c_{1})$
(12)

Let us move on to some real life materials. Consider copper, usually the
material that forms the wires of our home, or salt which we eat. They are
crystals and electrons move around in them.

Question: Is brass, or steel a crystal? If not, what are they? How does one
draw band diagrams for them?

Then, should just the lattice and saying electrons move from one site to
another – tell us everything about them? The answer is clearly no. How many
atoms we have, how different are they, what are the different hopping
parameters – all forms the first basic skeletal structure to help us start
making the first guesses. Real life systems have extremely complicated band
structures, which can be found to a great extent of accuracy using modern
computational tools where the chemistry of the atoms are taken into account to
a great deal. Since the momentum space is in three dimensions – now we have
three momentum points $(k_{x},k_{y},k_{z})$ which are associated with various
labels such as $\Gamma,X,W$ etc. A typical band diagram is shown in Fig. 3.
Notice the ways various momentum eigenstates move around in energies as we
move around in momentum space. Interestingly at zero energy – notice in copper
we have a line with cuts through, saying that energy required to shift an
electron to a new state is $zero$ while for silicon there is a ‘gap’! Which
means one needs a finite amount of energy to add a particle to the next band
(see Fig. 3). This is the crucial feature which distinguishes a metal from an
insulator.

Let’s therefore ask which is the simplest way of modeling a metal or an
insulator in one dimension?

(a)

(b)

Figure 3: (a) Band structure of copper [Burdick_PR_1963] which is a metal
commonly used in electrical wires. Energies are measured in Rydbergs
($=13.6eV$). Fermi energy (see $E_{f}$) crosses a band signifying a metal. (b)
Band structure of silicon [silicon_bs] which is a semiconductor used
extensively to form the substrate of electronic chips. Fermi energy sits at
$E=zero$ which lies in the band gap.

Example 4 Example: On-site energies makes things interesting!
So far, we have considered the onsite energy to be $zero$. That is if you
check Example 2, $\langle j|H|j\rangle$ is taken to be $zero$ for every site
$j$. It could have been some constant for every site, but that would just be
equivalent to adding a scaled identity matrix to the Hamiltonian. This would
not change the eigenvectors and the eigenvalues would just have a constant
shift. But now, let’s consider the case where we have a onsite energy
$\epsilon_{0}$ for the odd numbered sites and different onsite energy
$\epsilon_{1}$ for the even numbered sites. The Hamiltonian then becomes
$\displaystyle H=$ $\displaystyle t\Big{(}\sum_{j=1}^{9}|j\rangle\langle
j+1|+|j+1\rangle\langle j|\Big{)}+t\Big{(}|1\rangle\langle
10|+|10\rangle\langle 1|\Big{)}$ $\displaystyle+$
$\displaystyle\Big{(}\sum_{j=1}^{5}\epsilon_{0}|2j-1\rangle\langle
2j-1|+\epsilon_{1}|2j\rangle\langle 2j|\Big{)}$ (13) The Hamiltonian matrix is
then given by: $|1\rangle$ $|2\rangle$ $|3\rangle$ $|4\rangle$ $|5\rangle$
$|6\rangle$ $|7\rangle$ $|8\rangle$ $|9\rangle$ $|10\rangle$ $|1\rangle$
$\epsilon_{0}$ t 0 0 0 0 0 0 0 t $|2\rangle$ t $\epsilon_{1}$ t 0 0 0 0 0 0 0
$|3\rangle$ 0 t $\epsilon_{0}$ t 0 0 0 0 0 0 $|4\rangle$ 0 0 t $\epsilon_{1}$
t 0 0 0 0 0 $|5\rangle$ 0 0 0 t $\epsilon_{0}$ t 0 0 0 0 $|6\rangle$ 0 0 0 0 t
$\epsilon_{1}$ t 0 0 0 $|7\rangle$ 0 0 0 0 0 t $\epsilon_{0}$ t 0 0
$|8\rangle$ 0 0 0 0 0 0 t $\epsilon_{1}$ t 0 $|9\rangle$ 0 0 0 0 0 0 0 t
$\epsilon_{0}$ t $|10\rangle$ t 0 0 0 0 0 0 0 t $\epsilon_{1}$ How do we
diagonalize this? We need to be able to exploit some periodic pattern.
Consider the following solution
$|\psi\rangle=\frac{1}{\sqrt{10}}\sum_{j=1}^{5}e^{ikj}\Big{(}a_{k}|2j-1\rangle+b_{k}|2j\rangle\Big{)}$
(14) From our previous experience let us also put a constraint on $e^{i5k}$ to
be 1. Now, if we observe the eigenvalue equation in the $3^{rd}$ and $4^{th}$
row, we get $\begin{bmatrix}\epsilon_{0}&t(1+e^{-ik})\\\
t(1+e^{ik})&\epsilon_{1}\end{bmatrix}\begin{bmatrix}a_{k}\\\
b_{k}\end{bmatrix}=E_{k}\begin{bmatrix}a_{k}\\\ b_{k}\end{bmatrix}$ (15) You
can check that the same equation is obtained for other pairs of rows as well.
For a given $k$ the energy eigenvalues are obtained by solving eqn. (15) . It
is a quadratic equation and yields two values of energy. From the condition of
$e^{i5k}=1$, $k$ is $\frac{2\pi m}{5}$ where $m$ ranges from -2 to 2 (keeping
$k$ between $-\pi$ and $\pi$). For each value of $k$ we have two values of
energy giving 10 eigenvalues as is expected. What’s interesting to note now is
that $a_{k}$ and $b_{k}$ need not have the same magnitude and the electron can
now have a preference for one of the sites based on whether it is odd or even
numbered! Another takeaway is that if we plot energy as a function of $k$, we
get 2 bands as opposed to one band in the previous example. If $\epsilon_{0}$
and $\epsilon_{1}$ are unequal, then there exists a gap between these bands.
Now if the system is half-filled i.e., the Fermi energy lies in the middle of
the two bands (Fermi energy is the energy below which all states are filled),
then the gapped system is an insulator while the gapless one is a metal! For
$\epsilon_{0}=\epsilon_{1}$ there is no gap and for
$\epsilon_{0}=\epsilon_{1}=0$ this problem is identical to the one in Example
3. But we have two bands here and we only had one band previously. Work out
the wavefunctions for both cases and convince yourself that the solutions
actually turn out to be the same. Labelling the parameter space in different
ways is just an effective tool, to interpret the eigenvalues and eigenvectors
we obtain.

While the simplistic examples may miss everything about what copper or silicon
does, we find that just these two parameters allow us to model one essential
phenomena which we saw in the real material – of having a gapless and gapped
energy spectra.

##  3 The curious case of a polymer

Polymers are compounds where a ‘large molecule’ repeats itself periodically.
For instance in the usual polythene which we often use at home to carry fruits
and vegetables is a chain made of a molecular unit which looks like this

-C(-[2]H)(-[6]H)-C(-[2]H)(-[6]H)-

Here each carbon with a valency of four is strongly covalently bonded with two
other carbons and two hydrogen atoms. It is an insulating system with quite a
large band gap and low melting point. Let’s focus on another polymer now which
looks like this:

-C(-[6]H)=C(-[2]H)-C(-[6]H)=C(-[2]H)-C(-[6]H)=C(-[2]H)-C(-[6]H)=C(-[2]H)-

Unlike the previous example, here every carbon just has one associated
hydrogen atom – and carbons now share three electrons. Notice that even bonds
share two electrons each! But the fact we drew the double bonds on even
numbered bonds is just a choice. We could have instead put them on the odd
ones. For instance this,

=C(-[6]H)-C(-[2]H)=C(-[6]H)-C(-[2]H)=C(-[6]H)-C(-[2]H)=C(-[6]H)-C(-[2]H)=

In the actual molecules the bond angles are not $90^{\circ}$ but we are going
to ignore this specificity in our discussion here. Both these configurations
have the same energies, and therefore in a system they would be equally
likely. Let’s call them version PA-$\alpha$, and PA-$\beta$. Now, polymers are
long chains where such units just attach to each other. So it may seem that we
can have two possible ways of attaching them:

$\text{Junction 1:}\ \ \boxed{PA-\alpha}--\boxed{PA-\alpha}\ or\
\boxed{PA-\beta}--\boxed{PA-\beta}$ (16) $\text{Junction 2:}\ \
\boxed{PA-\alpha}--\boxed{PA-\beta}\ or\ \boxed{PA-\beta}--\boxed{PA-\alpha}$
(17)

Are these two ways really the same?

History of this polymer is quite amazing. The physics of attaching these
various chains and that it was far from trivial is tied to some remarkable
work which showed an anomalous signal in a magnetic susceptibility experiment.

Benzene too has a similar bond structure with even and odd bonds, just that it
has 6 atoms of carbon. Do they have two different versions of the compound –
like polyacetelyne? If yes, why? If not, why not? Answer: It doesn’t. Now
think about why not. W.P. Su who was then a graduate student in UPenn, John
Schrieffer, who had received the Nobel prize a few years earlier for his work
on superconductivity and was a theorist in the same department, and A.J.
Heeger – who was an expert on experiments on polymers, again working there –
came together to understand this and wrote a series of papers on just trying
to understand this simple looking polymer that had nothing but carbon and
hydrogens. While these set of papers founded a new sub-field and ideas that
stemmed topological physics – Heeger continued his investigations on polymers.
He eventually won the Nobel prize in Chemistry in 2000 for his work on
polymers, in particular polyacetelyne – the polymer you see above. In case,
you haven’t appreciated – just sit back and think for a minute the impact
these polymers have had on physics, chemistry, and our whole planet! And all
of it had people scribbling funny bond diagrams.

##  4 Welcome to the SSH chain

In the last few sections we looked at a few tight binding Hamiltonians. We
also noticed on how to think about metals and insulators.

Let’s recap. Basically given a Hamiltonian we look at its energy spectrum.
This is often done assuming periodic boundary conditions and invoking the
Bloch wave functions. The energy levels which represent allowed states, when
closely spaced, can form bands.

Consider polythene – whose band gap was $\sim 8eV$ – how much voltage would
one need to apply for the electron to conduct?

Often when the number of electrons in a lattice is half of the number of atoms
– we have a system which we call half filling. The energy upto which the
electrons gets filled is called Fermi energy or the chemical potential. The
lower (filled) bands are called the valence bands, and the upper bands (empty)
are termed the conduction bands. For insulators the chemical potential or the
Fermi energy lies in the band gap and hence need sufficient (often humongous
amounts of) energy to excite electrons and get the material to conduct. For
metals such as copper, there is no gap, these bands criss-cross each other and
therefore any little energy allows for electrons to move.

So what about polyacetelyne?

Su-Schrieffer-Heeger (SSH) modelled the chain saying we have two kinds of
atoms – A and B (see Fig. 4). Now even though both atoms are really carbon,
their bond-strengths aren’t the same. We know the system likes to either have
stronger even or odd bonds and can have two configurations. So consider the
Hamiltonian has parameters $w$ and $v$ which models the bond strengths.

Figure 4: Schematic for the SSH Model. A and B within the same unit cell are
coupled with the intracellular hopping parameter $v$, and adjacent unit cells
are coupled through the intercellular hopping parameter $w$.

The Hamiltonian is then given by

$H=v\sum_{j=1}^{N/2}|2j-1\rangle\langle
2j|+w\sum_{j=1}^{N/2}|2j+1\rangle\langle 2j|+h.c.$ (18)

In second quantized form this becomes

$H=v\sum_{j=1}^{N/2}c_{2j-1}^{\dagger}c_{2j}+w\sum_{j=1}^{N/2}c_{2j+1}^{\dagger}c_{2j}+h.c.$
(19)

Given two distinct atoms it is helpful to consider a unit cell comprising of
two atoms. The first atom is then referred to as the $A$ site of the unit
cell, and the second atom is referred to as the $B$ site of the unit cell. The
Hamiltonian in eqn. (18), is written in such a way that the odd numbered sites
correspond to the $A$ sites and the even numbered sites correspond to the $B$
sites in Fig. 4. The Hamiltonian after relabelling this numbering to consider
unit cells with their $A$ and $B$ sites is

$H=v\sum_{j=1}^{N/2}|j,A\rangle\langle
j,B|+w\sum_{j=1}^{N/2}|j+1,A\rangle\langle j,B|+h.c.$ (20)

$H=v\sum_{j=1}^{L}c_{j,A}^{\dagger}c_{j,B}+w\sum_{j=1}^{L}c_{j+1,A}^{\dagger}c_{j,B}+h.c.$
(21)

where $L$ is the number of unit cells. The number of atoms are $N$ and this
gives us the relation $N=2L$ and $j$ labels the unit cells.

###  4.1 Band Diagram

Given the Hamiltonian for the SSH chain, we can Fourier transform and
diagonalize the Hamiltonian. Solving in a similar fashion as in Example 4 we
get the following $2\times 2$ matrix eigenvalue problem

$\begin{bmatrix}0&(v+we^{-ik})\\\
(v+we^{ik})&0\end{bmatrix}\begin{bmatrix}a_{k}\\\
b_{k}\end{bmatrix}=E_{k}\begin{bmatrix}a_{k}\\\ b_{k}\end{bmatrix}$ (22)

The energy eigenvalues turn out to be

$E_{k}=\pm\sqrt{v^{2}+w^{2}+2vw\cos(k)\ }$ (23)

One thing that might strike as odd at this place is that the energy expression
seems to be symmetric with respect to $v$ and $w$ implying that PA-$\alpha$
and PA-$\beta$ are the same when viewed from this lens. The band diagram for
some values of $v$ and $w$ is illustrated in Fig. 5

Figure 5: Band diagram of the SSH chain for various parameter values

We now track the band diagram as $w$ and $v$ is altered. We keep in mind that
we have a system which is half-filled which means the Fermi energy is pinned
to $zero$. We find the system is always an insulator when $v\neq w$ but
becomes a metal when $v$ equals $w$!

###  4.2 Modelling the junction of polyacetelyne

Now for polyacetelyne – we had two types of phases – the one with even and odd
bonds stronger and then two ways of making junctions out of them. We are now
at the stage of asking whether all of these junctions behave similarly.

(a)

(b)

Figure 6: Modeling the polyactelyne chain junctions (a) Connecting the last
atom of a $PA\mathchar 45\relax\alpha$ chain to the first atom of another
$PA\mathchar 45\relax\alpha$ chain. (b) Connecting the last atom of a
$PA\mathchar 45\relax\alpha$ chain to the first atom of a $PA\mathchar
45\relax\beta$ chain.

Let’s model an odd and even numbered bonds between carbon atoms via parameters
$v$ and $w$ respectively. Let’s have $PA\mathchar 45\relax\alpha$ with
$v_{\alpha}=1,w_{\alpha}=0.5$ and $PA\mathchar 45\relax\beta$ with
$v_{\beta}=0.5,w_{\beta}=1.0$ and each having 10 carbon atoms. As seen in fig.
6 trying to connect them with periodic boundary conditions, we have two
possible cases. Connecting two $PA\mathchar 45\relax\alpha$s or two
$PA\mathchar 45\relax\beta$s would be equivalent just giving us a regular
chain with the first junction type from eqn. (16). Connecting two different
versions of PA as seen in Fig. 6(b) would give us the second junction type
from eqn. (17). We now have two consecutive double bonds at one end and two
consecutive single bonds at the other end. The energy eigenvalues for the two
cases have been plotted in Fig. 7.

Figure 7: Energy eigenvalue spectrum for the two setups formed in Fig. 6 has
been shown here. Edge states occur when $PA-\alpha$ and $PA-\beta$ are joined
end-to-end.

Notice in fig. 7 for the second type of chain the eigenvalues show some energy
states close to zero energy. In fact if we look at the probability density of
the electron it resides on one of these junctions. So these electrons just
like to reside on the junction and not worry about other things!

This is a fantastic result – if for a large polymer which is just made of one
version of the PA - there are no such localized electrons, but say when such
junctions are present at a finite density – that means a bunch of electrons
just sit on the junctions. It is these electrons which gives an anomalous
signal in the experiments and was the discovery SSH gang made.

###  4.3 A single chain with open boundaries

While the polyacetelyne gave us the motivation to model the junctions – to
think why bond strengths can be different. It really gave us a toy model to
play with. We can now forget about the compound, and start finding out what
does a single chain tells us. What happens if we now think about $w$ and $v$
as some parameters which can be tuned in a lab.

Let us look at the energy spectrum, for $w$ = 1 and $v$ varying from 1 to 4 in
Fig. 8(a). If we focus on small values of $v$, we see that the spectrum
resembles an insulator except for two states with energy very close to zero.
If we keep the Fermi energy close to zero, we can excite these two edge
states. But what do these states correspond to? Looking at the probability
density of these states, we notice that it is highly localized on the edge
sites as shown in Fig. 8(b)! So although the material is insulating in the
bulk, it has some anomalous states at the edges.

(a)

(b)

Figure 8: (a) Energy spectrum of a SSH chain with 12 atoms and open boundary
conditions. We note that the zero energy states exist in the topological
regime and (b) $|\Psi|^{2}$ for two states, in the topological regime
($v=0.1$, $w=1$). The $E=0$ state is highly localized on the edge sites.

To obtain analytic expressions for the wave functions in a arbitrary open
system is in general hard. Let’s see if a zero energy wavefunction exists in
this system [short_topo] –

Consider the following form for the wavefunction

$|\psi\rangle=\sum_{j=1}^{L}\Big{(}a_{j}|j,A\rangle+b_{j}|j,B\rangle\Big{)}$

For a zero energy solution, we need to find coefficients $a_{j}$ and $b_{j}$
such that $H|\psi\rangle=0$, where $H$ is the Hamiltonian for the open SSH
chain (see eqn. (19)). Looking at each row, we get the following set of
equations

$va_{j}+wa_{j+1}=0\ \ \ \ \ for\ j=1\ldots L-1$ $wb_{j}+vb_{j+1}=0\ \ \ \ \
for\ \ j=1\ldots L-1$ $va_{L}=0\ \ \ \ \ vb_{1}=0$

This would give us the following recurrence relations and a condition

$\displaystyle a_{j}=\frac{-w}{v}a_{j+1}\ \ \ \ \ for\ j=1\ldots L-1$ (24)
$\displaystyle b_{j+1}=\frac{-w}{v}b_{j}\ \ \ \ \ for\ j=1\ldots L-1$ (25)
$\displaystyle a_{L}=b_{1}=0$ (26)

Looking at these equations one might think that having $b_{1}=0$ would force
all other $b_{j}s$ to 0 (similarly for $a_{j}s$) and that no such solution
having energy to be exactly zero could exist!

But for large $L$ we could try and cook up exponentially decaying coefficients
to get an approximate zero energy solution. Let us define a parameter $\xi$ as
follows

$\xi=\frac{1}{\log{|w|}-\log{|v|}}$

From the recurrence relations we now get

$\displaystyle|a_{N/2}|=e^{-\frac{L-1}{\xi}}|a_{1}|$ (27)
$\displaystyle|b_{1}|=e^{-\frac{L-1}{\xi}}|b_{L}|$ (28)

If $\xi$ were negative, as $L\to\infty$ we could satisfy eqn. (26)! And $\xi$
is indeed negative if $v<w$!

Great, so now we can have two possible solutions, one by setting $a_{1}$ to 0
and $b_{L}$ as non-zero and vice-versa. What do these solutions correspond to?
Having $a_{1}$ as non-zero corresponds to the electron having the highest
probability of residing in site $|j=1,A\rangle$ and then decaying
exponentially as we move to the right. While having $b_{L}$ as non-zero
corresponds to the electron having the highest probability of residing in site
$|j=L,B\rangle$ and then decaying exponentially as we move the left. For this
reason $\xi$ is also termed as the localization length. The fun part is that
both these solutions reside at the edges!

##  5 It’s not just a phase

In the above analysis we noticed something quite strange. The two Hamiltonians
for the polyacetalyne versions seemed quite the same. The eigenvalues in the
periodic system were in fact exactly the same – both describing insulators at
half-filling. Yet, in an open system one of them hosts boundary modes and the
other doesn’t. Are these phases of electrons really the same? Or are these
insulators different?

Turns out they are very subtly different. And this is what we discuss next.

###  5.1 Polarization

Insulators, while electrically may not conduct – they may still be different.
This is because, insulators polarize in presence of electric fields – and the
way they polarize can be different. This is the electrostatic response to an
electric field, and to understand its effect we first need to understand how
to think about polarization within a tight-binding framework.

Notice we have a charge neutral system so the number of electrons and number
of positive charges should be the same. This means for the SSH system, we
should effectively consider that every two atoms has one positive charge to
balance one shared electron. This of course assumes we have ignored all other
(protons’) charges and deeper electron orbitals.

Consider two systems (see Fig. 9) where electrons (represented by a shaded
orange ellipse) and positive charges shown in squared rectangles have the same
centre of mass. Notice this has no dipole moment, however in the other case
where the electron cloud is slight shifted – here the system has a dipole
moment.

Figure 9: Two insulators, but with different dipole moments.

Both these systems are insulators, given electrons are localized. And their
wavefunctions must have an observable that distinguishes them. Clearly,
energies to excitations etc. cannot. One thing which potentially can is the
centre of mass of the electrons. Consider operator $n_{jA}$ which measures the
electron density on the $A^{th}$ site of the $j^{th}$ unit cell (see Fig. 4) -
the centre of mass is measured by

$\hat{M}=\sum_{j=1,\ldots,L}j(n_{jA}+n_{jB})$ (29)

However we run into trouble, how can one think about this in a periodic
system? After all, the site we label $1$ or $L$ is just a choice we made. And
physically meaningful quantities, such as observables shouldn’t depend on such
choices. We therefore place the chain on a ring – and calculate the center of
mass of the electrons on a chain using this operator instead –

$\hat{O}=\exp(i\frac{2\pi\hat{M}}{L})=\exp\Big{(}i\frac{2\pi}{L}\sum_{j=1,\ldots,L}j(n_{jA}+n_{jB})\Big{)}$
(30)

Figure 10: Placing the chain on a ring

On a circle, this operator is well defined and doesn’t change its eigenvalues
for a half-filled system. This operator, as we will see below will play a
crucial role in discovering the physics of the SSH chain. However before we
calculate these explicitly let us develop the understanding of this operator a
bit more.

###  5.2 Applying a small electric field

In the last section we looked at $\hat{O}$ as a centre of mass for electrons,
but there is another way of thinking about this. Imagine we wish to apply a
small electric field to an electronic system, in general we could do this by
adding a potential of the form to a discrete Hamiltonian

$V=v\sum_{j}jn_{j}$ (31)

where the electric field in dimensionless units is just $v$. But then we run
into the same trouble – where we do not know how to keep periodicity intact.
So we remind ourselves that electric fields can be applied in another way – by
threading a magnetic flux since

$E=-\frac{\partial A}{\partial t}$ (32)

We need to therefore understand in a tight-binding framework how should one
electron sense a magnetic field. As electron moves from a site $i$ to $j$ it
should sense a phase
$\int_{\bm{r}_{i}}^{\bm{r}_{j}}\mbox{\boldmath$A$}.d\bm{r}$. This often goes
by the name of Peierl’s substitution [Sakurai] [peierl]

$t_{ij}:=t_{ij}e^{-i\frac{e}{\hbar}\int_{i}^{j}\mathbf{A}\cdot d\mathbf{l}}$
(33)

For instance consider a tight-binding chain with a one of kind of atom as we
saw in Example 4 whose Hamiltonian is given by

$\hat{H}=\sum_{i=1}^{L}-\left[tc_{i+1}^{\dagger}c_{i}+tc_{i}^{\dagger}c_{i+1}\right]$
(34)

An operator given by

$\hat{O}=\exp\Big{(}i\frac{\phi}{L}\sum_{j=1,\ldots,L}jn_{j}\Big{)}$ (35)

changes

$c_{l}^{\dagger}\rightarrow
e^{i\frac{\phi}{\Phi_{0}}l}c_{l}^{\dagger}\quad\quad
c_{l+1}^{\dagger}c_{l}\rightarrow e^{i(\phi/L)}$ (36)

implementing a Ahranov-Bohm phase [Griffiths2004Introduction] arising from a
flux $\phi$ in a ring. Therefore, for $\phi=2\pi$ this really becomes the
operator for dipole moment that we were looking for. This implies that as an
electron moves through the ring it gains a phase of $2\pi$.

Example 5 Example: Let us consider a 4 atom square as shown in Fig. 11.
Figure 11: 4 atom periodic chain in a magnetic flux The Hamiltonian of this
setup would be: $H=\begin{bmatrix}0&t&0&t^{\ast}\\\ t^{\ast}&0&t&0\\\
0&t^{\ast}&0&t\\\ t&0&t^{\ast}&0\\\ \end{bmatrix}$ (37) Here, $t$ is the
hopping parameter between nearest neighbours of the chain which can be set to
1 without loss of generality. $t^{\ast}$ denotes the complex conjugate of $t$.
Now if a magnetic flux $\phi$ is passed through this setup then the hopping
parameters would change as given by Peierl’s substitution. Consider $\phi$ =
$\phi_{0}$ where $\phi_{0}$ is a constant. If the length of the entire chain
is $L$ (i.e., distance between adjacent sites is $L/4$), then the magnetic
vector potential is $A$ = $\phi_{0}/L$. Now, using eqn. (33), we get that the
hopping parameter modifies to: $t:=te^{i\frac{\phi_{0}}{4}}$ (38) Let us now
plot the energy eigenvalues of the Hamiltonian under consideration as
$\phi_{0}$ changes from 0 to $2\pi$. Figure 12: Energy eigenvalues of the
4-atom Hamiltonian given in eqn. (37) as a function of magnetic field
$\phi_{o}$.

Therefore for a state $|\Psi\rangle$ which comprises of single particle states
given by say the momenta labels $k_{1},k_{2},k_{3}..k_{N_{p}}$ where $N_{p}$
is the number of particles

$|\Psi\rangle=c^{\dagger}_{k_{1}}c^{\dagger}_{k_{2}}\ldots
c^{\dagger}_{k_{N_{p}}}|\Omega\rangle$ (39)

and $|\Omega\rangle$ is the vacuum state. We could interpret
$\langle\Psi|\hat{O}|\Psi\rangle$ as an expectation of the dipole moment (when
placed on a ring) or, how the state has changed after an action of a small
electric field (or a flux of $2\pi$). This has a particularly neat
interpretation for what we know of metals and insulators. For instance, let us
ask how a single $k$ state behaves under $\hat{O}$

$Oc_{k}^{\dagger}O^{\dagger}=c_{k+2\pi/L}^{\dagger}$ (40)

So the action of operator $\hat{O}$ on a half-filled band, for instance where
all the $k$ states aren’t filled – shifts each of the $k$ points by $2\pi/L$
thereby changing the many body state leading to $\langle\hat{O}\rangle=0$,
while for a filled band where all the $k$ states are already filled – it can’t
change it since $k=\pi+\frac{2\pi}{L}\rightarrow k=-\pi+\frac{2\pi}{L}$
leading to the fact that $\langle\hat{O}\rangle\neq 0$ (see Fig. 13).

Figure 13: Action of operator $\hat{O}$ on a metal and an insulator. While on
the metal it takes it to a new state, on an insulator it brings it to itself
albeit with a phase.

In fact for a multi-band system (such as SSH, where we have two bands for
every $k$ point) one can fill one of the bands – and as the state comes back
to itself it can carry a non-trivial phase. The general value of the
expectation on the operator is given by

$\langle\Psi|\hat{O}|\Psi\rangle=re^{i\gamma}$ (41)

where $\gamma$, the phase turns out to play a crucial role, something we
discuss in the next section.

###  5.3 Geometric phase

In the above discussion we found that for an insulating state, action of
operator $\hat{O}$ on an insulator just shifts each $k$ point to the next $k$
point. Therefore the value of $\gamma$ could in principle be calculated from
just the single particle wavefunctions when they are calculated at each value
of $k$

$\gamma=-\arg\left[\prod_{k=-\pi}^{\pi}\left\langle\psi_{k}\mid\psi_{k+2\pi/L}\right\rangle\right]$
(42)

For a multiband system where we have a band index one can calculate this for
each band

$\gamma_{n}=-\arg\left[\prod_{k=-\pi}^{\pi}\left\langle\psi^{n}_{k}\mid\psi^{n}_{k+2\pi/L}\right\rangle\right]$
(43)

Therefore if $|\psi^{n}_{k}\rangle$ is thought of as a vector, and given $k$
goes from $-\pi$ to $\pi$ (i.e a ring) the phase we are worrying about is in
fact a geometric phase which is captured when a vector twists over a ring.

###  5.4 Geometric phase in discrete SSH chain

We now go back to the SSH model we were discussing before, and see how
$\hat{O}$ acts. Diagonalizing the Hamiltonian (see eqn. (22)) gives us the
following eigenstates

$|\pm k\rangle=\frac{1}{\sqrt{2}}\begin{bmatrix}\pm e^{-i\phi(k)}\\\ 1\\\
\end{bmatrix}\ \ ,\ \
\phi(k)=\tan^{-1}\left(\frac{w\sin(k)}{v+w\cos(k)}\right)$ (44)

The creation operator for the two bands can be written as:

$\alpha_{k\pm}^{\dagger}=\frac{1}{\sqrt{2}}\left[e^{-i\phi(k)}c_{k+}^{\dagger}\pm
c_{k-}^{\dagger}\right]$ (45)

where for the lower band

$\alpha_{k-}^{\dagger}\equiv\frac{1}{\sqrt{2}}\frac{1}{r_{(k)}}\left[\left(v+we^{-ik}\right)c_{k+}^{\dagger}-r_{(k)}c_{k-}^{\dagger}\right]$
(46)

Here, $r_{(k)}$ is the absolute value of the eigenvalue. The ground state
wavefunction is given by $|GS\rangle=\prod_{k\in
BZ}\alpha_{k-}^{\dagger}|\Omega\rangle$. Therefore,

$\hat{O}|GS\rangle=\prod_{k\in
BZ}\hat{O}\alpha_{k-}^{\dagger}\hat{O}^{-1}\hat{O}|\Omega\rangle$ (47)

This simplifies to

$\hat{O}|GS\rangle=\prod_{k\in
BZ}\frac{1}{\sqrt{2}}\frac{1}{r_{(k+2\pi/L)}}\left[\left(v+we^{-i(k+2\pi/L)a}\right)c_{k+}^{\dagger}-r_{(k+2\pi/L)}c_{k-}^{\dagger}\right]$
(48)

This leads to

$\gamma=-\sum_{k\in
BZ}\arg\left[\frac{\left(v+we^{ik}\right)\left(v+we^{-i(k+2\pi/L)}\right)}{r_{(k)}r_{(k+2\pi/L)}}+1\right]$
(49)

To do a sanity check on the above equation, we can put $v$ = 1 and $w$ = 0 to
get $\gamma$ = 0 and putting $v$ = 0 and $w$ = 1 gives us $\gamma$ = $\pi$.
Interestingly, remember that at $v=0$ or $w=0$ we have a set of decoupled
dimers. So these two insulators are different – they polarize differently
under a electric field – something we had anticipated before. Even more
surprisingly, this feature just doesn’t happen at extreme limits. For the
entire regime when $v<w$, we have $\gamma=0$ and then it suddenly changes to
$\pi$ at $v=w$ (see Fig. 15). These transitions are sensitive to finite size –
as the system size increases the transition happens exactly at $v=w$.

(a)

(b)

(c)

(d)

Figure 15: (a) Geometric phase ($\gamma$) switches from $0$ to $\pi$ at
$\delta=0.5$ albeit with finite size effects. (b) Behavior of the overlap
angle (eqn. (49)) as a function of $k$ for $N=16$ for $v/w<1$, i.e. the
topological regime and for $v/w>1$ in (c). The behavior of overlap angle when
very close to gap closing $(v\sim w)$ is shown in (d) for $N=200$ system.

The geometric phase which gets captured between one point of $k$ and the next
$k$ point can tell us where in the band is the non-trivial geometric phase
located. We calculate this for each individual $k$ and show as a function of
$k$ in Fig. 15. An interesting thing to note is while the qualitative shape of
the curve does not change for $v>w$ and $v<w$ (see (b) and (c) of Fig. 15)
they integrate to different quantities: zero or $\pi$ depending on whether we
are in a trivial or a topological phase. Another observation is that when we
are close to the phase transition point ($v=w$) the entire contribution to the
geometric phase comes from the points where the band gap closes in the
Brillouin zone! (see Fig. 15 (d)).

For a fun introduction to topology, do watch YouTube lectures by Tadashi
Tokieda.

The reason we call the $v<w$ regime as a topological phase, is because it is
characterized by a nontrivial integer $\gamma/\pi$ which by construction
cannot change smoothly but only suddenly by an integer. Moreover, this was
found by integrating a geometric phase over the complete Brillouin zone. This
is remarkably different the way we think about conventional condensed matter
phases where we can identify the nature of the phase by quantities which
change smoothly and can be observed or measured locally. For example, a
conventional magnet has magnetization everywhere in the sample which can
smoothly go to zero as temperature is increased.

Another feature of a topological phase such as $v<w$ regime, is that the
boundary physics is stable to disorder and perturbations – something we won’t
delve into here. This property that there are physical systems which are
characterized by quantities such as integers which remain stable to smooth
changes shares a deep similarity with the subject of topology in mathematics,
for instance such as the number of punctures in any material can only be an
integer and cannot change smoothly. These connections can be made very
concrete and goes under the name of topological quantum field theory –
something we won’t discuss at all.

Example 6 Example: Dipoles everywhere Does the energy spectrum tell us all we
need to understand about any system? Consider the two chains given below. The
first one is just our usual SSH chain. The second one seems a little odd, the
$w$ hopping is now connected to the next nearest unit cell instead of the
nearest unit cell. Let’s take a case where $v=0$ for both the chains. Can you
write the real space Hamiltonian for both the chains?  Figure 16: Chain 1
Figure 17: Chain 2 It may seem that the second system is just a rearrangement
of the first one, given both of them are a bunch of disconnected dimers with
the same strength of hopping. In fact elongating the dimers in the first chain
can smoothly take us to the second chain. The surprise is that these two
systems are in fact not smoothly deformable, and one needs to necessarily go
via a metallic phase. To see this one has to invoke the polarization framework
we had discussed before. The two Hamiltonians in the $k$ space are
$H_{1}=\begin{bmatrix}0&e^{-ika}\\\ e^{ika}&0\\\ \end{bmatrix}$ (50)
$H_{2}=\begin{bmatrix}0&e^{-i2ka}\\\ e^{i2ka}&0\\\ \end{bmatrix}$ (51) where
$a$ is the inter-unit cell distance. As is expected for disconnected dimers,
the spectrum for both $H_{1}$ and $H_{2}$ are identical and have flat bands.
To interpolate between them let us construct a new Hamiltonian,
$H_{\lambda}=\lambda H_{1}\ +\ (1-\lambda)H_{2}$ (52) Does our system remain
an insulator while we vary $\lambda$ from 0 to 1?

(a) $\lambda$ = 0, $\lambda$ = 1

(b) $\lambda$ = 0.3

(c) $\lambda$ = 0.5

(d) $\lambda$ = 0.8

Figure 18: Energy spectrum for the Hamiltonian in eqn. (52), for different
values of $\lambda$. The gap closes for $\lambda$ = 0.5 Turns out, one can’t
do that because we encounter a metallic phase at $\lambda=0.5$ (see Fig. 18).
This suggests that $\lambda<0.5$ and $\lambda>0.5$ are two different phases.
Are they topologically different? Can you calculate the geometric phase for
$H_{1}$ and $H_{2}$?

##  6 Entanglement entropy

Up until last section we found polarization is one way that can distinguish
the two phases of the SSH system. We now discuss, quantum entanglement, an
idea which has no classical analogue - and how it also plays a crucial role in
these topological phases.

###  6.1 What is entanglement entropy?

Given any system we can partition the complete Hilbert space as the tensor
product of two Hilbert sub-spaces.

$H_{AB}=H_{A}\otimes H_{B}$

where $A$ and $B$ represent the sub-spaces. In this context one can pose: for
a wavefunction defined on the complete Hilbert space, can it also be written
in a separable form between these two sub-spaces? One way to answer this is to
Schmidt decompose a wave function

$|\Psi\rangle\
=\sum_{i=1}^{n}\alpha_{i}|u_{i}\rangle_{A}\otimes|v_{i}\rangle_{B}$ (53)

where $|u_{i}\rangle_{A}$ and $|v_{i}\rangle_{B}$ are orthonormal states in
the subsystem A and B respectively. The state is called a product (un-
entangled) state if only one of the coefficients $\alpha_{i}$ is non-zero
where such a separation is indeed possible; more generally such a separation
is not possible and the state can carry a finite entanglement.

Consider a system with two spinless electrons on two sites. Are these two
electrons entangled? Given the density matrix of the system

$\rho^{AB}=|\Psi\rangle\langle\Psi|$ (54)

the partial trace of this density matrix over the subsystem A is called
$\rho^{A}$ while over the subsystem B is called $\rho^{B}$.

The entanglement entropy $S$ is defined as

$S=-\text{Tr}[\rho^{A}\log(\rho^{A})]=-\text{Tr}[\rho^{B}\log(\rho^{B})]$ (55)

This equals

$S=-\sum_{i}e_{i}\log(e_{i})$ (56)

where the $e_{i}$ are the eigenvalues of the partial traced density matrix
$\rho^{A}$ or $\rho^{B}$.

It is fairly straightforward to see that in terms of the Schmidt decomposition
coefficients (see eqn. (53)) and definitions of the partially traced density
matrix, $S$ can be represented in terms of the Schmidt coefficients
$\alpha_{i}$ as

$S=-\sum_{i=1}^{n}|\alpha_{i}|^{2}\log(|\alpha_{i}|^{2})$

The reader may heard that two spin 1/2 particles, when in a singlet state, has
an entanglement entropy of $\log(2)\sim 0.693$. At this point it may be
tempting to think that in order to discuss entanglement entropy, we need at
least two or more particles. But rather than the number of particles – it is
the wavefunction and its separability in real space that determines the value
of $S$ (see Example 7).

Example 7 Example: Let’s consider two sites $A$ and $B$ with a hopping term
between them and with one electron. The ground state wavefunction is given by
$|\Psi\rangle=\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\ -1\\\ \end{bmatrix}$ Or in
second quantization notation
$|\Psi\rangle=\frac{1}{\sqrt{2}}(c_{A}^{\dagger}-c_{B}^{\dagger})|\Omega\rangle$
where $|\Omega\rangle$ is the vacuum state. With two sites the total Hilbert
space dimension is $2^{2}=4$. The basis states being, the $zero$ particle
state $|\Omega\rangle$, the 1 particle states $c_{A}^{\dagger}|\Omega\rangle$
and $c_{B}^{\dagger}|\Omega\rangle$ and the two particle state
$c_{A}^{\dagger}c_{B}^{\dagger}|\Omega\rangle$. The state $|\Psi\rangle$ in
the complete Hilbert space is then given by
$|\Psi\rangle=\frac{1}{\sqrt{2}}\begin{bmatrix}0\\\ 1\\\ -1\\\ 0\\\
\end{bmatrix}$ The subsystems we are considering are the site $A$ with states
$|\Omega\rangle$ and $|c_{A}^{\dagger}|\Omega\rangle$ and similarly the single
site $B$ with the states $|\Omega\rangle$ and $c_{B}^{\dagger}|\Omega\rangle$.
The density matrix is given by $\rho^{AB}=\frac{1}{2}\begin{bmatrix}0&0&0&0\\\
0&1&-1&0\\\ 0&-1&1&0\\\ 0&0&0&0\\\ \end{bmatrix}$ The partial trace over
subsystem A is $\rho^{A}=\frac{1}{2}\begin{bmatrix}1&0\\\ 0&1\\\
\end{bmatrix}$ with the entanglement entropy given by
$S=-\Big{(}\frac{1}{2}\log(\frac{1}{2})+\frac{1}{2}\log(\frac{1}{2})\Big{)}$
$S=\log(2)=0.693$ Therefore, even for a simple system with just one electron
hopping between two sites, there is a non-zero entanglement entropy associated
with the partitioning the system into two halves. This toy problem also gives
an intuition that there is an entropy of $\log(2)$ attached to breaking of a
bond, something we will later find useful.

Figure 19: Entanglement entropy as a function of $\delta$ for halfway
partition of a SSH chain with periodic boundary conditions for different
lengths $N$.

Therefore given a wavefunction, one can evaluate a density matrix and trace it
over partially to obtain a reduced density matrix. This when done for the
ground state wavefunction in the SSH model leads us to Fig. 19. While for
$v<w$ the entanglement entropy is small and close to zero, it saturates to
$\log(2)$ when $v>w$! Therefore even entanglement entropy can provide us
indications of having a topological phase!

While a density matrix calculation is easily done for small system sizes, as
the system size increases, the dimensions required for our state and the
density matrices would blow up exponentially (for a system of size $N$, the
density matrix size is $2^{N}\times 2^{N}$). Is there a way then to extract
the same information in another, but in a simpler way? Turns out that for a
non-interacting fermion model, as is our case, a two-point correlator matrix
could tell us all that we need to know [Calabrese].

The number of bonds being cut is representative of the entanglement entropy
between the two halves We had a look at the following dipoles (see below)
earlier in Example 6. Would the entanglement entropy of the left half with the
right half be different for the cases considered here?

###  6.2 Correlators

In condensed matter physics, the “correlator” plays a very important role in
the understanding of ground state properties as well as the dynamics of the
system. The single particle correlator, which will be invoked a lot in this
work is defined as $C_{\alpha\beta}=\langle
c_{\alpha}^{\dagger}c_{\beta}\rangle$ describes the correlations between
indices $\alpha$ and $\beta$ in any given state, where $\alpha$ and $\beta$
denote any arbitrary basis labels. For instance, consider an insulator where
the ground state wave function is given by
$|\Phi_{G}\rangle=\prod_{k}c_{k}^{\dagger}|\Omega\rangle$ (i.e., all the
momentum states are filled). Here, we simply get
$C_{k,k^{\prime}}=\langle\Phi_{G}|c_{k}^{\dagger}c_{k^{\prime}}|\Phi_{G}\rangle=\delta_{kk^{\prime}}$,
which indicates that we do not have long “distance” correlation. Of course,
here distance simply points to the difference in the label number in an
arbitrary basis.

To put the concept of distance to a real context, we can now calculate the
averaged correlation amplitude of processes that connect two sites $i$ and
$j$\- namely the simultaneous annihilation of an electron at site $j$ and the
creation of an electron at site $i$. One may immediately recognize that if
$i=j$, we do have that the correlator represents the electron density at the
site $i=j$. Let us investigate this a little further. Here we suggest the use
of the word “simultaneous”, literally meaning equal time, a concept that is
meaningful in the unperturbed system. With the introduction of a time
dependent perturbation, as we will shortly see, the correlators have to be
extended to two different times. For now, let us focus on the equal time
static correlators, in connection with unperturbed ground states.

Let us see what happens to the long distance correlation between two sites in
an eigenstate, i.e.,
$C_{ij}=\langle\Phi_{G}|c_{i}^{\dagger}c_{j}|\Phi_{G}\rangle$. In order to
understand this, we ask the question: what happens when we remove an electron
from one site and add to another simultaneously in (a) metallic state and b)
insulating state? Before we get to the math, let us summarise the essentials.
In the insulating state to be completely filled such that
$|\Phi_{G}\rangle=\prod_{k\in BZ}c_{k}^{\dagger}|\Omega\rangle$, whereas in a
metallic state we have
$|\Phi_{G}\rangle=\prod_{|k|<k_{F}}c_{k}^{\dagger}|\Omega\rangle$. In the
former case, $BZ$ represents the Brillouin zone and hence a fully filled
system. The latter case represents a partially filled system, such that at
zero temperature, the states $k<k_{F}$, the Fermi wave vector, are filled. A
simple pictorial representation of the two cases can be gathered from Fig. 13.

We can also quite well surmise what will happen in the former case. In a
filled band, with the states being filled in the eigenspace, quite obviously
the sites comprising the real space are also occupied. Keeping the Pauli
exclusion principle in mind, it will not be possible to “simultaneously”
remove an electron in site $i$ and place in site $j$, since site $j$ will also
be occupied, of course unless $i=j$. Thus it is easy to make the right guess
that $C_{ij}=\delta_{ij}$. In the latter case i.e. for the metal, this might
not be the case since it is a partially filled system in the eigenbasis and
hence in real space also! let us now formalize this calculation.

Using

$c_{i}^{\dagger}=\frac{1}{\sqrt{N}}\sum_{k\in BZ}e^{ikx_{i}}c_{k}^{\dagger}$
(57)

for a fully filled state one finds,

$C_{ij}=\frac{1}{N}\sum_{k,k^{\prime}\in
BZ}\delta_{kk^{\prime}}e^{i(kx_{i}-k^{\prime}x_{j})}=\frac{1}{N}\sum_{k\in
BZ}e^{i\left(k\left(x_{i}-x_{j}\right)\right)}=\delta_{ij}.$ (58)

However, for the partially filled band where

$|\Phi_{G}\rangle=\prod_{|k|<k_{F}}c_{k}^{\dagger}|\Omega\rangle$ (59)

one obtains

$C_{ij}=\sum_{\mid k\mid\mid
k^{\prime}\mid<k_{F}}\delta_{k,k^{\prime}}e^{i(kx_{i}-k^{\prime}x_{j})}=\frac{1}{N}\sum_{\mid
k\mid<k_{F}}e^{i\left(k\left(x_{i}-x_{j}\right)\right)}$ (60)

where one uses the fact that
$\langle\Phi_{G}|c_{k}^{\dagger}c_{k^{\prime}}|\Phi_{G}\rangle=\delta_{k,k^{\prime}}$
when $|k|<k_{F}$, denoting filled states and
$\langle\Phi_{G}|c_{k}^{\dagger}c_{k^{\prime}}|\Phi_{G}\rangle=0$ when
$|k|>k_{F}$, denoting empty states. When $i=j$, the correlator
$C_{ii}=\frac{N_{f}}{N}$, denoting that the probability of occupying a given
site $i$ is smaller than unity and equal to the density!

Interestingly for $i\neq j$, assuming $x_{i}=0$ and a thermodynamically large
system

$C_{j}\sim\frac{2}{2\pi}\int_{0}^{k_{F}}\cos(kx_{j})dk=\frac{\sin(k_{F}x_{j})}{\pi
x_{j}}$ (61)

This shows that for a metal not only the correlator has characteristic
oscillations that depend on the Fermi vector but also a power law decay as a
function of distance $x_{j}$ from $x_{i}=0$. This is remarkably different from
an insulating state, and is an alternate way of thinking about metals and
insulators!

Now that we have seen one interesting example of the use of the correlator, we
will now see some more applications of the correlator, specifically when it
comes to determining the topological aspect of a given state. Here, the
correlator- most importantly the one-particle two point correlator will be
used in lieu of the full density matrix that we encountered earlier.

###  6.3 Entanglement Entropy using the correlator matrix

Here we point out that the two-point correlators and entanglement entropy are
in fact intricately related and one can calculate the latter using the former,
in the non-interacting case. We now define a two point correlator matrix.
Given a $N$ site chain, it is a $N\times N$ matrix CM, such that

$\text{CM}_{ij}=\langle\Psi|c_{i}^{\dagger}c_{j}|\Psi\rangle$

where $c_{i}^{\dagger}$ is the creation operator for site $i$ and $c_{j}$ is
the annihilation operator for site $j$ and $\text{CM}_{ij}$ is the expectation
of this operator taken over the many body ground state we had encountered
before: $|\Psi\rangle$. But wait, does this really help? Even though CM is
$N\times N$ , to calculate it we still need to calculate the $2^{N}\times$1
$|\Psi\rangle$. Fortunately one can also write CM, in terms of the one
particle first quantized $N\times$1 wavefunctions.

$\text{CM}=\sum_{i}|\Phi_{i}\rangle\langle\Phi_{i}|$ (62)

where the summation is over all the occupied states having energy less than
the Fermi energy. The way to see this is given a second quantized Hamiltonian
$\hat{H}$, we can write this in terms of a matrix $\mathcal{H}$ as

$\hat{H}=\mathbf{C}^{\dagger}\mathcal{H}\mathbf{C}$ (63)

where C is a column vector comprising of the annihilation operators $c_{i}$ at
site $i$. We can diagonalize $\mathcal{H}$ to rewrite it as

$\hat{H}=\mathbf{C}^{\dagger}R^{T}DR\mathbf{C}$ (64)

where $D$ is a diagonal matrix of the energy eigenvalues $\epsilon_{n}$ and
C=$R\mathbf{C_{f}}$, where $C_{f}$ is column vector comprising of annihilition
operators $f_{n}$ where $f^{\dagger}_{n}$ creates an eigenstate at the
$n^{th}$ eigenstate. For a ground state

$|\Psi\rangle=\prod_{\epsilon_{n}<E_{f}}f^{\dagger}_{n}|\Omega\rangle$ (65)

such that number of filled states are $N/2$ we get

$\langle\Psi|c_{i}^{\dagger}c_{j}|\Psi\rangle=\sum_{qp}^{N}\langle
0|f_{q}^{\dagger}R^{T}_{qn}R_{mp}f_{p}|0\rangle=\sum_{p}^{N/2}R^{T}_{pn}R_{mp}$
(66)

Given that the columns of $R$ are the first quantized wave functions, we
arrive at eqn. (62)!

So now, how do we get the entanglement entropy from the correlation matrix
[hatsugai, Peschel] and how do we partition it?
For most of the examples we’d consider partitioning the chain into halves, but
in general to partition over sites 1..$j$ and $j$+1…$N$, we can simply
consider the truncated matrix CMt which is a $j\times j$ sub matrix of CM.
Alternatively, we could have also considered the other sub-system which
corresponds to a matrix size $(N-j)\times(N-j)$, but the entanglement entropy
would be the same in both cases. Given a truncated CMt matrix, the
entanglement entropy $S$ is given by

$S=-\sum_{i=1}^{j}\Big{(}e_{i}\log(e_{i})+(1-e_{i})\log(1-e_{i})\Big{)}$ (67)

where $e_{i}$ is an eigenvalue of CMt. For a not so quick derivation on how to
arrive at this see Appendix.A.

Example 8 Example: Let us try to obtain the entanglement entropy for a 4 atom
SSH chain with periodic boundary conditions using the correlator matrix method
we discussed. Using the density matrix method we would require a $16\times 16$
size matrices. Let’s see if we can do better this time. Let’s take the
topological case with $v=0$ and $w=1$. The four single-particle eigenvalues of
this setup are -1,-1,1,1. For the half filled system we look at the two lowest
energy levels and their wavefunctions.
$|\Phi_{1}\rangle=\frac{1}{\sqrt{2}}\begin{bmatrix}0\\\ -1\\\ 1\\\ 0\\\
\end{bmatrix}\ \ \ \ \ \
|\Phi_{2}\rangle=\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\ 0\\\ 0\\\ -1\\\
\end{bmatrix}$ Using eqn. (62), we have the correlator matrix $CM$ to be
$CM=\frac{1}{2}\begin{bmatrix}0&0&0&0\\\ 0&1&-1&0\\\ 0&-1&1&0\\\ 0&0&0&0\\\
\end{bmatrix}+\frac{1}{2}\begin{bmatrix}1&0&0&-1\\\ 0&0&0&0\\\ 0&0&0&0\\\
-1&0&0&1\\\ \end{bmatrix}$ $CM=\frac{1}{2}\begin{bmatrix}1&0&0&-1\\\
0&1&-1&0\\\ 0&-1&1&0\\\ -1&0&0&1\\\ \end{bmatrix}$ Now, since we want the
entropy of entanglement of the first half with the second half, we truncate
the top left $2\times 2$ matrix of CM. Referring to this as CMt, we have
$CMt=\frac{1}{2}\begin{bmatrix}1&0\\\ 0&1\\\ \end{bmatrix}$ Both the
eigenvalues of $CMt$ are 0.5. Using eqn. (67), we have the entanglement
entropy S as,
$S=-2\Big{(}\frac{1}{2}\log(\frac{1}{2})+(1-\frac{1}{2})\log(1-\frac{1}{2})\Big{)}$
$S=2\log(2)=1.3863$ We can see this value of entropy in Fig. 19 and we only
had to work with matrices of size $4\times 4$!

###  6.4 A relook at the SSH Chain

Great, equipped with correlator matrices, we can look at entanglement entropy
in SSH chains.

First we look at a SSH chain with periodic boundary conditions. In the plot
Fig. 19 we show the behavior of entanglement entropy as a function of $\delta$
where $v=1-\delta$ and $w=\delta$, for chains with different lengths $N$. Let
us just look at $N$=20 for the rest of this discussion. For $\delta=0$ we have
$v=1$ and $w=0$ – the trivial region. We choose a partition at $j=N/2$ to ask
how entangled the first half is with the second half.

What would one expect the entanglement entropy to be? Site 10 and 11 are
connected by hopping parameter $w$ which is 0 in this case. Intuitively one
would expect the entropy to be 0 and from the plot, it indeed is! Now let’s
look at the other end of the plot where $\delta=1$ and we have $v=0$ and $w=1$
i.e. the extreme topological region. Now, site 10 and 11 are connected and
from our earlier toy example, we expect the entropy for one broken bond to be
$\log(2)$. But the actual value from the plot is twice of it – $2\log(2)$!
Does this mean that two bonds are being broken? The answer is in affirmative,
since what we have is a periodic chain or a ring, and hence by partitioning at
$j=N/2=10$, we cut the bond between site 10 and 11 as well as the bond between
site 20 and 1!

At $\delta=0.5$ there appears to be a peak whose value increases with $N$. At
this point we have $v=w=0.5$ and therefore the chain is a metal. Metals have
higher sub-system entanglement and it grows with system size albeit slowly.
Next, let us look at the SSH chain with open boundary conditions.

(a)

(b)

Figure 20: (a) Entanglement entropy for halfway partition in a SSH chain with
open boundary conditions for different lengths $N$. (b) Entanglement entropy
for halfway partition in a SSH chain with open boundary conditions, after
removing the highest energy electron (so as to remove the edge state electron
in the topological regime) for different lengths $N$.

In Fig. 20(a) for $\delta=0$ we have the entropy to be zero as earlier.
However for $\delta=0.99$ we have the entropy to be $2\log(2)$. This is
counter intuitive, since we no longer have the second bond due to periodicity
between the first and the $N^{th}$ site. So where does the other $\log(2)$
come from? At $\delta\sim 0$ why does the system with $N/2-1$ electrons has
$\log(2)$ entropy? In the topological regime, we also have two edge states in
an open chain. For $v=0.01$ and $w=0.99$ these states hybridize leading to a
bonding orbital between site 1 and site $N$. It is this bonding orbital which
leads to an additional $\log(2)$ contribution.

To check this, let us keep the number of electrons in the system to be one
less than $N$/2. This would lead to the edge state remaining unoccupied. The
result is shown in Fig. 20(b) and as expected, the entanglement entropy now
remains $\log(2)$ when $\delta>0.5$.

###  6.5 U can make all the difference

So far, we have been looking at models where there is no electron electron
interaction. Let’s look at the Hubbard chain where such an interaction is
considered and electrons are now assumed to be spinful. Therefore, each site
has two states in the singly occupied sector[Hubbard_BM], a spin up state and
a spin down state.

Figure 21: A diagram giving a pictorial depiction of the hopping terms given
in eqn. (68).

The Hamiltonian is given by,

$H=-t\sum_{i,\sigma}\
c_{i,\sigma}^{\dagger}c_{i+1,\sigma}+c_{i+1,\sigma}^{\dagger}+U\sum_{i}\
n_{i\uparrow}n_{i\downarrow}$ (68)

where $n_{i,\sigma}$ is the number operator on that site for a particular spin
$\sigma$. Just like the SSH chain we now have $v$ and $w$ hoppings for every
spin sector separately (see Fig. 21).

(a)

(b)

Figure 22: Entanglement entropy for halfway partition in a Hubbard chain(open
boundary conditions) with 4 sites for two values of $U$, (a) $U$ = 0.5 and (b)
$U$ = 10 CM denotes that the entropy has been calculated using the correlator
matrix and DM denotes that the entropy has been calculated using the density
matrix.

The parameter $U$ models the interaction where two electrons of opposite
spins, at the same site, have an energy strength of $U$. If $U$ is negative,
then there is an added incentive for electrons to occupy the same site and
vice-versa. If $U=0$, then the chain in Fig. 21 is simply 2 parallel SSH
chains, a spin up chain and a spin down chain that are not interacting with
each other. Now, since we’ve introduced the model, we can try to calculate the
entanglement entropy of the left half with the right half in presence of
interactions. Here’s the tricky part though, since the product of two number
operators is no longer quadratic, the correlation matrix method to calculate
the entropy would miss something compared to the density matrix method. But
how bad could it be?

Fig. 22 shows this calculation using both the density matrix method and the
correlator matrix method for a 4 site Hubbard chain with open boundary
conditions. For the $U=0$ (blue curve), the entanglement entropy is just twice
that of a single SSH chain and both methods lead to the same result. For
$U=0.5$, the two methods are still quite similar only deviating from each
other in the topological regime of $\delta>0.5$. But for $U=10$, the two
methods do not agree with each other at all! Do we understand these numbers?
Just like the number of bond cuts or other ideas we could use to interpret our
earlier results?

Well, the upshot of this section is to make you appreciate that there is a lot
we do not understand about topological phases of matter, in particular what
role do the interactions play.

##  7 Taking lead

###  7.1 Looking at electronic transport

Up until now, we have looked at periodic systems and open chains. But real
systems are also affected by externalities of the environment to which it is
connected. Let us now try to model such systems. We will consider that the 1D
open chains that we have been analyzing up until now are connected to the
environment using leads. Such systems are called open-quantum systems and one
way to analyze them is through the Keldysh non-equilibrium Green’s function
(NEGF) technique [LNE_datta] [datta_2005_chapter_8][datta_2005_chapter_9]. The
Green’s functions are themselves multi-point correlators taken at different
times, and for the purposes of steady state calculations, we Fourier transform
the time differences into the energy domain. We will leave the reader to these
advanced considerations from other excellent pedagogical references [Jauho].
Our pedagogy here begins with the consideration of steady state situations via
the energy domain Green’s functions.

The leads are also referred to as contacts. We will assume that these contacts
are large sources or sinks that are maintained at a constant chemical
potential. These open systems are important to understand as they can help us
simulate actual experiments and thus give us an idea as to what to expect.
Employing NEGF, would help us simulate their local density of states,
conductance and many other properties as we will see later. We will try to
explain the basics of NEGF through a simple example.

Example 9 Example: Consider a 1D chain of 8 atoms as shown in Fig. 23.
Consider a potential $V(x)$ over this 1D lattice chain and let us assume that
the chain is attached to contacts on both the ends. These contacts are
themselves semi-infinite 1D tight-binding chains. Take, $\psi_{i}$ to be the
discretized value of the wavefunction $\psi(x)$ at the $i^{th}$ lattice
position. Figure 23: An eight atom 1D chain is considered. The setup is an
open quantum system, i.e. boundaries of the lattice are connected to contacts
(semi-infinite 1D chains) The lattice is connected to two contacts – left
contact and right contact. The left and the right contacts have a constant
potential $V_{LC}$ and $V_{RC}$ respectively. We commonly refer the left
contact as the source and the right contact as the drain. Why do we take
contacts to be semi-infinite 1D chains? The semi-infinite chains help them act
as large sinks which maintain a constant chemical potential. Therefore, we can
write wavefunctions in that contact regions as plane waves. The wave function
in the left contact is $e^{ik_{1}x}$ \+ $re^{-ik_{1}x}$. Similarly, wave
function in the right contact is $te^{ik_{2}x}$. Here, $k_{1}$ is
$\sqrt{\frac{2mV_{LC}}{\hbar^{2}}}$ and $k_{2}$ is
$\sqrt{\frac{2mV_{RC}}{\hbar^{2}}}$. The lattice sites 1 and 8 are connected
to the contacts. Note that, $V_{1}$ is equal to $V_{LC}$ and $V_{8}$ is equal
to $V_{RC}$. Why? This is because the first and the eighth lattice sites are a
part of the contact and hence the potential there should be equal to the
potential of the connected contact. Using eqn. (2), we can write the
discretized Schrodinger’s equation for the lattice region as:
$E\psi_{i}-\left[-t_{0}\psi_{i-1}+\left(2t_{0}+V_{i}\right)\psi_{i}-t_{0}\psi_{i+1}\right]=0$
(69) where, $t_{0}$ = $\hbar^{2}/2ma^{2}$. eqn. (69) holds for all $i$ from 1
to 8. More specifically for $i=1$ we have
$E\psi_{1}-\left[-t_{0}\psi_{0}+\left(2t_{0}+V_{1}\right)\psi_{1}-t_{0}\psi_{2}\right]=0$
(70) Now, we employ the boundary conditions to evaluate $\psi_{0}$ which is
the wavefunction amplitude on the first site of the left contact. Given the
first atom of the device (site 1) is located at $x=0$, $\psi_{1}=1+r$, and
$\psi_{0}=\psi_{x=-a}=e^{-ik_{1}a}+re^{ik_{1}a}$. Placing these values, in
eqn. (70) gives us: $E\psi_{1}-\left[\hbox{\pagecolor{red!50}$\displaystyle-
t_{0}e^{ik_{1}a}\psi_{1}$}+\left(2t_{0}+V_{1}\right)\psi_{1}-t_{0}\psi_{2}\right]=\hbox{\pagecolor{red!50}$\displaystyle
t_{0}(e^{ik_{1}a}-e^{-ik_{1}a})$}$ (71) where we have substituted
$\psi_{0}=\psi_{1}e^{ik_{1}a}-\left(e^{ik_{1}a}-e^{-ik_{1}a}\right)$.
Similarly, if we apply boundary conditions at lattice site 8, we would get
$\psi_{8}$ = $te^{ik_{2}(x=7a)}$ = $te^{i7k_{2}a}$ and $\psi_{9}$ =
$te^{ik_{2}(x=8a)}$ = $te^{i8k_{2}a}$. Using, these substitutions in eqn. (69)
written for lattice site 8 would give us:
$E\psi_{8}-\left[\hbox{\pagecolor{red!50}$\displaystyle-
t_{0}e^{ik_{2}a}\psi_{8}$}+\left(2t_{0}+V_{8}\right)\psi_{8}-t_{0}\psi_{7}\right]=0$
(72) The highlighted terms in both eqn. (71) andeqn. (72) are different from
the bulk sites ($2-7$) shown in eqn. (69). Notice that the complete system
(with the two semi-infinite leads) and the device, is essentially an infinite
system which cannot be solved. However by using the boundary conditions, and
the effect of the leads (see eqn. (71) and eqn. (72)) we are able to capture
the effective physics just within the the sites belonging to the device. If
$H$ is the Hamitlonian for just the device with open boundary conditions, then
the total effective Hamiltoian (including the effect of leads) is given by
$[EI-H-\Sigma_{1}-\Sigma_{2}]\psi=S$ (73) where  Figure 24: Self-energy
matrices for the considered setup are non-Hermitian square matrices of size 8
with one element non-zero. Source matrix is a column matrix of size 8 with the
top element non-zero. $\Sigma_{2}$ as similar to $\Sigma_{1}$ in its form
except that it has the non-zero element at the right-bottom $\Sigma_{1}$ and
$\Sigma_{2}$ are called self-energy matrices and $S$ is called the source
matrix. In the above case, $\Sigma_{1}$ is square matrix with all terms 0
except the term at (1,1) which is $t_{0}e^{ik_{1}a}$, $\Sigma_{2}$ has all
terms 0 except the term at (8,8) which is $t_{0}e^{ik_{2}a}$ and $S$ is a
$8\times 1$ column matrix with all terms 0 except the first term which is
$t_{0}(e^{ik_{1}a}-e^{-ik_{1}a})$. The size of the self energy matrices is the
same as the system size but unlike the Hamiltonian they are not necessarily
Hermitian. We can rearrange eqn. (73) and write $\psi={\cal G}(E)S$ (74) where
$\cal{G}$ is called the Green’s function and is equal to
$[EI-H-\Sigma_{1}-\Sigma_{2}]^{-1}$. Briefly, if you remember, Green’s
functions provide us an alternate of solving differential equations where
instead of finding a solution of a differential equation by either a qualified
guess or explicit expansions (such as Frobenius method), one uses an integral
form where a charge put in the system and one thinks about what response it
creates [ARFKEN2013447]. For example in electromagnetic theory, one can
consider putting a single electric charge in real space and asking what
electric field configuration it creates in its surrounding space. Eqn. (74) is
now of the same form where $S$ can be considered like a charge or a source,
$\Psi$, the wavefunction in the device is the effect it creates. Is easy to
see if $H+\Sigma_{1}+\Sigma_{2}$ has a set of eigenvalues
$\epsilon_{n}=\epsilon^{R}_{n}+i\epsilon^{I}_{n}$ and eigenvectors
$|n\rangle$, then ${\cal G}(E)=\sum_{n}\frac{|n\rangle\langle
n|}{E-\epsilon^{R}_{n}-i\epsilon^{I}_{n}}$ (75)

Note that $H+\Sigma$ has complex eigenvalues in general and thus open systems
gives a finite lifetime to eigenstates. While for a hermitian Hamiltonian,
that has just a real eigenvalue $\epsilon_{n}$ for a stationary eigenstate
$\Psi_{n}$, the time evolution of this state is given by

$\Psi_{n}(t)=\exp(-i\epsilon_{n}t)\Psi_{n}$ (76)

which just adds a phase dependent on time. In the case of a complex eigenvalue
which has a real and an imaginary part
($\epsilon_{n}=\epsilon^{R}_{n}+i\epsilon^{I}_{n}$)

$\Psi_{n}(t)=\exp(-i\epsilon^{R}_{n}t)\exp(\epsilon^{I}_{n}t)\Psi_{n}$ (77)

which leads to a decay of the wavefunction if $\epsilon^{I}<0$ within a time
$\tau\sim\frac{1}{\epsilon^{I}}$. So $\tau\sim\frac{1}{\epsilon^{I}}$ is
called the lifetime of the electron. This essentially means that unlike a
periodic Hamiltonian where an electron cannot vanish, in an open system due to
the presence of leads, an electron can decay into the contacts.

Show that the DOS for a one-dimensional tight-binding model goes as
$D(E)\sim\frac{1}{\sqrt{E}}$ at low energies. How can DOS diverge?

These decay physics also leads to a broadening of density of states (DOS) in
the system. DOS essentially quantifies the number of eigenstates at energy
$\epsilon$. For a set of real eigenvalues at energies $\epsilon_{n}$ which are
fairly separated from each other

$D(E)=\sum_{n}\delta(E-\epsilon_{n})$ (78)

which is a $\delta$ function at every real energy $\epsilon_{n}$. However when
these energies have an imaginary part each of these $\delta$ functions gets
modified to a lorentzian where

$\displaystyle D(E)$ $\displaystyle=$
$\displaystyle\sum_{n}\frac{1}{2\pi}\frac{|\epsilon^{I}_{n}|}{(E-\epsilon^{R}_{n})^{2}+\frac{1}{2}(\epsilon^{I}_{n})^{2}}$
(79) $\displaystyle=$
$\displaystyle\frac{i}{2\pi}\Big{[}\Big{(}\frac{1}{E-\epsilon^{R}_{n}-i\epsilon^{I}_{n}}\Big{)}-\Big{(}\frac{1}{E-\epsilon^{R}_{n}+i\epsilon^{I}_{n}}\Big{)}\Big{]}$
(80)

This is called spectral broadening which is characteristic to systems in which
particles have a finite lifetime.

While such a broadening for a single eigenvalue is easy to capture within
lifetimes of every eigenstate, for a lattice system it is more natural to
define a matrix called $\Gamma$, which is the imaginary part of the self-
energy matrix

$\Gamma=i(\Sigma-\Sigma^{\dagger})$ (81)

that characterises the effect of leads on the device, in terms of creating
finite broadening and corresponding decay of the electrons from within the
device. Therefore corresponding to each contact there is a $\Gamma$ matrix for
e.g., $\Gamma_{1}=i(\Sigma_{1}-\Sigma_{1})$ and
$\Gamma_{2}=i(\Sigma_{2}-\Sigma_{2})$.

The NEGF formalism can be used to calculate important properties of a device
like I-V characteristics, conductance (often called transmission here etc.).
Here we would would not go into the details and refer to an excellent book
[datta_2005_chapter_8] [datta_2005_chapter_9], but present the essential
quantities which one calculates such as

$T(E)=\text{Tr}[\Gamma_{1}{\cal G}(E)\Gamma_{2}{\cal G}(E)^{\dagger}]$ (82)

$T(E)$ is the transmission at an energy $E$. Calculating transmission tells us
‘conductance’ through the device at energy $E$. Naturally it depends on how
the electrons decay from both the contacts ($\Gamma_{1},\Gamma_{2}$) and the
way electrons propagate within the device ${\cal G}$.

Given the transmission it is easy to calculate the current ($I$) in the system
which is

$I=(q/h)\int_{-\infty}^{+\infty}\mathrm{d}ET(E)\left(f_{1}(E-\mu_{LC})-f_{2}(E-\mu_{RC})\right)$
(83)

which depends on the chemical potentials of the right and left contacts, their
corresponding occupancies ($f(\epsilon)$) and their transmission. The DOS
($D(E)$) of the complete system is given by

$D(E)=\text{Tr}[i({\cal G}(E)-{\cal G}(E)^{\dagger})/2\pi]$ (84)

which can be seen by combining eqn. (75) and eqn. (80).

What should be the transmission through a clean one dimensional wire? Can you
evaluate this using NEGF?

Now, for a contact which is semi-infinite and has a clean tight-binding
description we found specific forms for $\Sigma$ and $\Gamma$ matrices as
discussed in the previous example. This also allows us to evaluate the local
DOS (LDOS) which refers to the contribution of the DOS at a real space
position. Given a set of complex eigenvalues, the LDOS at site $j$ should be

$D(E,x=j)=\sum_{n}\frac{1}{2\pi}\frac{|\epsilon^{I}_{n}||\langle
n|x=j\rangle|^{2}}{(E-\epsilon^{R}_{n})^{2}+\frac{1}{2}(\epsilon^{I}_{n})^{2}}$
(85)

Clearly when summed over all sites $j$, the LDOS at energy $E$ equals $D(E)$.
This can be calculated quite simply by evaluating

$D(E,x=j)=\langle j|i({\cal G}(E)-{\cal G}(E)^{\dagger})/2\pi|j\rangle$ (86)

which nothing but the $j^{th}$ diagonal element of $i({\cal G}-{\cal
G}^{\dagger})/2\pi$.

However, more realistically, often we are unaware of the exact nature of the
leads. Even one can use the NEGF equations (see eqn. (73)) to model the
junction by choosing $\Gamma,\Sigma$ as some effective parameters. $\Sigma$
matrices in most simplest assumption can have an imaginary value on the first
and the last site of the device such that both the ($1,1$) element of
$\Sigma_{1}$ and ($N,N$) element of $\Sigma_{2}$ is again
$-\frac{i\gamma}{2}$. This leads to the corresponding $\Gamma$ matrices to
have the corresponding (1,1) and ($N,N$) term to be $\gamma$. This simplifies
the calculations remarkably, where even in absence of any information
regarding the leads, and just the device Hamiltonian ($H$) we can calculate
the current and the conductance using NEGF formalism. We will now apply these
tools to a SSH chain in order to investigate its transport signatures, if any.

(a)

(b)

(c) Figure 25: Local density of states (LDOS) in an SSH chain as a function of
lattice index and energy $E$ for three regimes (a) topological
($w=0.9,v=0.1$), (b) trivial ($w=0.1,v=0.9$) (c) metal ($w=0.5,v=0.5$) for a
32 atom system. A finite coupling to the leads is considered where
$\gamma=0.1w$.

###  7.2 Transport properties of the SSH chain

We have already seen that the SSH chain can exist in two phases – (i) a
topological phase when $v<w$ and (ii) a trivial phase when $v>w$. The
transition point from one phase to another is when $v=w$ and the system is a
metal. Just to remind ourselves (see fig. 8(a)), in the topological regime the
system has $\sim zero$ energy eigenstates at the edges. We will now use the
NEGF formalism to see if the topological regime has any transport signatures.

Intuitively, in a thermodynamically large system, even if the edge states
exist – they cannot lead to a finite conductance because these edge states
decay exponentially (see discussion near eqn. (28)) and therefore with
increasing system size, the overlap between the wavefunctions on the right and
the left edge will continue to fall exponentially, leading to no transmission
at close to $\sim zero$ energies. However in a finite size system, these edge
states can hybridize and lead to some weak signatures in the LDOS.

Fig. 25 plots the local density of states (LDOS) (see eqn. (86)) of a finite
SSH chain in the three different regimes – topological, trivial and a metal as
a function of energy $E$ when their is finite coupling to the leads. The first
figure shows the SSH chain in the topological limit and we can clearly notice
LDOS due to the edge sites hosting $\sim zero$ energy states. The second
figure shows the SSH chain in the trivial limit where there are no edge
states. The third figure shows the SSH chain in the conductance limit. As
expected, LDOS are spread uniformly across the entire spectrum of energy.

What do you think would be the transmission signature (transmission as a
function of energy) for these three phases of the SSH chain? (for a hint, see
the code provided.)

(a)

(b)

Figure 26:  The behavior of transmission ($T(E)$) as a function of $\delta$ at
$E=0$ for (a) different values $\gamma$ given a system size $N=32$ and (b)
different values of $N$ for $\gamma=1$.

We will now observe what happens if we try and apply voltages on the two ends
of the chain and pass current through it. We would try and observe the
difference, if any, in the two regimes. The behavior of transmission
(conductance) as a function of $\delta$ where $v=\delta$ and $w=1-\delta$ is
shown in fig. 26 for $E=0$. The conductance, as is expected shows a peak when
$\delta=0.5$, i.e., we are in a metallic system. In both $v<w$ or in $v>w$ the
transmission is close to $zero$ as is expected for insulating systems.
Interestingly if the coupling to the leads (approximated by a parameter
$\gamma$) is varied one can tune the transmission. For instance at $\gamma=0$
where the system is not coupled then the transmission is zero while with
increasing $\gamma$ the transmission increases when $\delta=0.5$. Some
representative plots are shown in fig. 26 for $\gamma=0.1w,1.0$. In fig. 26
(b) we look at the conductance as a function of increasing system size and
find that except at $\delta=0.5$ the conductance falls in the $\delta<0.5$ and
$\delta>0.5$ regime, consistent with the expectations that the insulators even
in the topological limit show no conductance signature.

We must now remark that the transport discussion presented here is for the
$U=0$ limit (absence of interactions). The NEGF technique in the strict sense
works well when the contact coupling $\gamma\gg U$. In the opposite regime, we
are in a situation of weakly coupled quantum dots and the transport is
typically treated within more complicated methods such as the density matrix
master equation approach in the Fock space. A detailed description of that
technique is carried out in various works including [Muralidharan_2008,
LNE_datta] which we will just refer to here.

You may wonder what surprise did we find for the transport signature of the
SSH chain – and the answer none! In fact the topological regime cannot be
detected with a conventional transport experiment for the SSH model. However a
LDOS measurement (via using a scanning tunneling microcope (STM)) can still
lead to some signatures (see Fig. 25) but we have not been able to
successfully do this in experiments yet. This reflects that even while
topological phases such as SSH model can very well be understood in periodic
systems etc., it is hard to detect them experimentally.

However the above programme of simulating such toy systems under a theoretical
framework has allowed us to visit a range of fantastic physics ideas which
span from ideas of polarization, entanglement, correlators, Green’s functions,
and many more. In fact it is studies like this which guides both the
experiments and the theoretical work which can lead to discoveries of new and
exciting phases of quantum condensed matter.

##  8 The way ahead

In the last sections we visited many-a-ideas in quantum condensed matter
physics to investigate a particular model Su-Schriefer-Heeger chain. SSH chain
being a prototypical example has been addressed in beautiful pedagogical
articles, we particularly recommend for the uninitiated [Anderson:1997vm,
short_topo, Batra_2020]. Here we have addressed some of the ideas which are
often not discussed in detail – particularly the idea of polarization, the
entanglement entropy, and electronic transport under the non-equilibrium
Green’s function. We hope this will encourage you, the reader, to take a more
serious study of these ideas and explore more!

If this article made you feel a bit excited about topological quantum
phenomena, then it is the right time to take the dive into the world of
quantum materials and devices. With the concepts of geometric phase and
winding numbers discussed here, one can move on to the exciting world of
topological electronics.

Just as we noticed that the topological phase of the 1D chain gives us an
“end-state”, these ideas can be generalized to higher dimensions (called bulk-
edge correspondence) where in two dimensions, we will have an edge state and
in three dimensions we will have a surface state! This brings in the concept
of topologically protected states along the boundary. One can immediately
envision that such topologically protected states along a boundary, if they
can propagate, can give rise to possibly dissipation less currents
[RevModPhys.82.3045]. This is because such states are “unimpeded” by
impurities and defects to a large extent. In order to start exploring along
these lines, we suggest a short road map.

For two dimensional structures (2D), a few more ideas such as the calculation
of Berry curvature and its relation to the Hall conductivity are needed
[shankar2018topological] to take the plunge in the 2D topological systems. A
detailed study of the quantum Hall effect and the deep connection of the Hall
conductance quantization with Berry curvature needs to be grasped. The
extension of the 2D system to the hexagonal lattice finally opens up a world
of possibilities including the quantum spin Hall effect and various related
topological phases. With that, a whole range of limitless device possibilities
[Gilbert] featuring such topologically protected states open up.

In 1972, Anderson told us “more is different” (read this today [Anderson393],
if you haven’t), it is fair to say that quantum condensed matter today is
finding that more and more things are topologically different. Don’t miss this
exciting physics journey if you haven’t already started!

##  9 Acknowledgements

A.A. would like to acknowledge Abhishodh Prakash for initiating the TQFT
lecture series at ICTS, where part of this material was discussed as in a set
of lectures. A.A. gladly acknowledges his wonderful years of learning and
research at IISc and ICTS, particularly with Vijay B. Shenoy and Subhro
Bhattacharjee, which seeded many of the ideas. A.A. acknowledges partial
funding from Max Plank Partner Grant at ICTS, Bangalore. A part of the
material was also developed by B.M. during the 2020 iteration of the course
titled “Topological Electronics” (EE-787) at IIT Bombay. B. M. would like to
acknowledge useful discussions with Sven Rogge, Gilles Buchs, Alestin Mawrie
and Ashwin Tulapurkar. B.M. also acknowledges funding under the Visvesvaraya
Ph.D Scheme of the Ministry of Electronics and Information Technology (MEITY),
Government of India, implemented by Digital India Corporation (formerly Media
Lab Asia). This work is also supported by the Science and Engineering Research
Board (SERB), Government of India, Grant No. STR/2019/000030, the Ministry of
Human Resource Development (MHRD), Government of India, Grant No.
STARS/APR2019/NS/226/FS under the STARS scheme.

##  Appendix A Entanglement entropy from correlator matrix

The density matrix when a system is in a state $|\Psi\rangle$

$\rho=|\Psi\rangle\langle\Psi|$ (87)

So any operator expectation on this state is

$\langle\hat{O}\rangle=\text{Tr}\Big{(}\hat{O}\rho\Big{)}$ (88)

While for a pure state, density matrix is indeed given by eqn. (87), in
general the density matrices can correspond to a mixed state or can be used to
model a system with an effective temperature or an ensemble of systems. There
while eqn. (87) is not valid, however eqn. (88) continues to hold.

For example, let’s a consider a single spin in a magnetic field $B_{0}$ along
the $z$ direction [2020Density]. The Hamiltonian is given by
$H=-\gamma\mathbf{S}.\mathbf{B}=-\gamma S_{z}B_{0}$. We have two possible
eigenstates, the $|\uparrow\rangle$ state along the z direction with energy
$-\gamma\hbar B_{0}/2$ and the $|\downarrow\rangle$ state with energy
$\gamma\hbar B_{0}/2$. Let $E_{0}$ be $\gamma\hbar B_{0}/2$. The thermal
density matrix is then given by

$\rho=\frac{1}{e^{-\frac{E_{0}}{kT}}+e^{\frac{E_{0}}{kT}}}\Big{(}e^{\frac{E_{0}}{kT}}|\uparrow\rangle\langle\uparrow|+e^{-\frac{E_{0}}{kT}}|\downarrow\rangle\langle\downarrow|\Big{)}$
(89)

In order to calculate the expectation of the $S_{x}$ operator using eqn. (88)
we get,

$\langle
S_{x}\rangle=\text{Tr}\Big{(}\frac{1}{e^{-\frac{E_{0}}{kT}}+e^{\frac{E_{0}}{kT}}}\begin{bmatrix}0&1\\\
1&0\\\ \end{bmatrix}\begin{bmatrix}e^{\frac{E_{0}}{kT}}&0\\\
0&e^{-\frac{E_{0}}{kT}}\\\ \end{bmatrix}\Big{)}$ (90)

which leads to $\langle S_{x}\rangle=0$. The expectation for $S_{y}$ would
also be 0. It would be non zero for $S_{z}$ indicating that the spin would
like to align itself along the z direction. This example also illustrates that
the thermal density matrices can be written as $\rho=e^{-\beta H}$ with
appropriate normalization constants.

Let’s go back to the density matrix of a pure state as shown in eqn. (87).
Subsystem density matrices which are obtained after tracing over part of the
system are mixed states and can often be written as

$\rho_{A}=Ke^{-\beta H}$ (91)

where $H$ is now called the entanglement Hamiltonian. $K$ is a normalization
constant. $\beta$ is not the temperature but just appears as a scaling factor
here.

Can every density matrix be written in the form of eqn. (91)? In principle,
yes. But one intriguing thing happens when traced density matrix belongs to a
free fermionic state – here the entanglement Hamiltonian is also a “free
Hamiltonian” i.e. it just has quadratic operators. Therefore expectation of
operators such as $c^{\dagger}_{i}c_{j}$ which really is the correlator we
discussed before can now be written as

$\langle
c^{\dagger}_{i}c_{j}\rangle=\text{Tr}\Big{(}c^{\dagger}_{i}c_{j}e^{-\beta
H}\Big{)}$ (92)

Note given a ground state wavefunction, the truncated correlator matrix (CMt)
is easily obtainable (see Section 6.3). However for density matrix, we were
interested in its eigenvalues to obtain the entanglement entropy. Now we
interpret the above equation as a thermal expectation over a different
Hamiltonian $H$ at a temperature $1/\beta$.

Therefore a set of wavefunctions which diagonalizes $H$ has a thermal
occupancy of every single particle state as

$n_{f}(\epsilon)=\frac{1}{1+\exp(\beta\epsilon)}$ (93)

So any observable, including a correlator, when measured on this thermal state
is given by

$\langle
c^{\dagger}_{i}c_{j}\rangle=\sum_{n}\phi^{*}_{in}\phi_{jn}\frac{1}{1+\exp(\beta\epsilon_{n})}$
(94)

But the above equation also shows that $C_{ij}$ matrix is itself
diagonalizable with the wavefunctions $\phi$. Hence the eigenvalues of CMt
matrix can tell us the eigenvalues of $H$! If CMt has eigenvalues $e_{i}$,
then

$e_{i}=\frac{1}{1+e^{\beta\epsilon_{i}}}$ (95)

or

$e^{\beta\epsilon_{i}}=\frac{1-e_{i}}{e_{i}}$ (96)

Given $\rho=e^{-\beta H}$, the trace of $\rho$

$\text{Tr}(\rho)=K\prod_{i}(1+e^{-\beta\epsilon_{i}})$ (97)

where $K$ is a normalization constant. Here the states are summed over all
possible occupancies of the free fermionic states. Given properties of density
matrix [Sakurai]

$K=\prod_{i}\frac{1}{(1+e^{-\beta\epsilon_{i}})}=\prod_{i}(1-e_{i})$ (98)

Notice that entanglement entropy is given by

$\displaystyle S$ $\displaystyle=-\sum_{i}\rho_{i}\log\rho_{i}$ (99)
$\displaystyle=-K\log(K)-K\sum_{i}e^{-\beta\epsilon_{i}}\log
Ke^{-\beta\epsilon_{i}}$ $\displaystyle\qquad-K\sum_{i,j,i\neq
j}e^{-\beta(\epsilon_{i}+\epsilon_{j})}\log
Ke^{-\beta(\epsilon_{i}+\epsilon_{j})}+\ldots$ (100)

It is easy to see that in terms of eigenvalues of the Hamiltonian we have

$S=-\log
K+\sum_{i}\frac{e^{-\beta\epsilon_{i}}\beta\epsilon_{i}}{1+\exp(-\beta\epsilon_{i})}$
(101)

$S=\sum_{i}-e_{i}\log e_{i}-(1-e_{i})\log(1-e_{i})$ (102)

It is interesting to see that $e_{i}$ is same as a Fermi function which tells
us the occupancy of any state $i$. In this sense the entropy gets
contributions both of filled and empty phase space!

## References

  * [1] V. S. Shchesnovich.  The second quantization method for indistinguishable particles (lecture notes in physics, ufabc 2010), 2013.
  * [2] Alexander Altland and Ben Simons.  Second quantization, page 39–93.  Cambridge University Press, 2006.
  * [3] Glenn A. Burdick.  Energy band structure of copper.  Phys. Rev., 129:138–150, Jan 1963.
  * [4] J. C. Phillips.  Band structure of silicon, germanium, and related semiconductors.  Phys. Rev., 125:1931–1936, Mar 1962.
  * [5] János K. Asbóth, László Oroszlány, and András Pályi.  A short course on topological insulators.  Lecture Notes in Physics, 2016.
  * [6] Jun John Sakurai.  Modern quantum mechanics; rev. ed. Addison-Wesley, Reading, MA, 1994.
  * [7] James G. Analytis, Stephen J. Blundell, and Arzhang Ardavan.  Landau levels, molecular orbitals, and the hofstadter butterfly in finite systems.  American Journal of Physics, 72(5):613–618, 2004.
  * [8] David J. Griffiths.  The Adiabatic Approximation.  Pearson Prentice Hall, 2nd edition, 2004.
  * [9] Pasquale Calabrese.  Entanglement entropy, the reduced density matrices of free fermions, and beyond.  Journal of Physics A: Mathematical and Theoretical, 49(42):421001, sep 2016.
  * [10] S. Ryu and Y. Hatsugai.  Entanglement entropy and the berry phase in the solid state.  Physical Review B, 73(24), Jun 2006.
  * [11] Ingo Peschel.  Calculation of reduced density matrices from correlation functions.  Journal of Physics A: Mathematical and General, 36(14):L205–L208, Mar 2003.
  * [12] B. Muralidharan, A. W. Ghosh, and S. Datta.  Conductance in coulomb blockaded molecules - fingerprints of wave-particle duality?  Molecular Simulation, 32(9):751–758, 2006.
  * [13] Supriyo Datta.  Lessons from nanoelectronics: A new perspective on transport.  08 2012.
  * [14] Supriyo Datta.  Level broadening, page 183–216.  Cambridge University Press, 2005.
  * [15] Supriyo Datta.  Coherent transport, page 217–251.  Cambridge University Press, 2005.
  * [16] H. Haug and A.-P. Jauho.  Quantum Kinetics in Transport and Optics of Semiconductors.  Springer, 1996.
  * [17] George B. Arfken, Hans J. Weber, and Frank E. Harris.  Chapter 10 - green’s functions.  In George B. Arfken, Hans J. Weber, and Frank E. Harris, editors, Mathematical Methods for Physicists (Seventh Edition), pages 447–467. Academic Press, Boston, seventh edition edition, 2013.
  * [18] B Muralidharan, L Siddiqui, and A W Ghosh.  The role of many-particle excitations in coulomb blockaded transport.  Journal of Physics: Condensed Matter, 20(37):374109, aug 2008.
  * [19] Philip W. Anderson.  When the electron falls apart.  Phys. Today, 50N10:42–47, 1997.
  * [20] Navketan Batra and Goutam Sheet.  Physics with coffee and doughnuts.  Resonance, 25(6):765–786, Jun 2020.
  * [21] M. Z. Hasan and C. L. Kane.  Colloquium: Topological insulators.  Rev. Mod. Phys., 82:3045–3067, Nov 2010.
  * [22] R. Shankar.  Topological insulators – a review, 2018.
  * [23] Matthew J. Gilbert.  Topological electronics.  Communications Physics, 4(1):70, 2021.
  * [24] P. W. Anderson.  More is different.  Science, 177(4047):393–396, 1972.
  * [25] The Density Matrix, 10 2020.  [Online; accessed 2021-07-27].