# Driven quantum many-body systems and out-of-equilibrium topology

Souvik Bandyopadhyay(1), Sourav Bhattacharjee(1) and Diptiman Sen(2)
(1)Department of Physics, Indian Institute of Technology Kanpur, Kanpur
208016, India
(2)Center for High Energy Physics and Department of Physics, Indian Institute
of Science, Bengaluru 560012, India

###### Abstract

In this review we present some of the work done in India in the area of driven
and out-of-equilibrium systems with topological phases. After presenting some
well-known examples of topological systems in one and two dimensions, we
discuss the effects of periodic driving in some of them. We discuss the
unitary as well as the non-unitary dynamical preparation of topologically non-
trivial states in one and two dimensional systems. We then discuss the effects
of Majorana end modes on transport through a Kitaev chain and a junction of
three Kitaev chains. Following this, transport through the surface states of a
three-dimensional topological insulator has also been reviewed. The effects of
hybridization between the top and bottom surfaces in such systems and the
application of electromagnetic radiation on a strip-like region on the top
surface are described. Two unusual topological systems are mentioned briefly,
namely, a spin system on a kagome lattice and a Josephson junction of three
superconducting wires. We have also included a pedagogical discussion on
topology and topological invariants in the appendices, where the connection
between topological properties and the intrinsic geometry of many-body quantum
states is also elucidated.

###### Contents

  1. I Introduction
  2. II Topology of fermionic systems in equilibrium
    1. II.1 Su-Schrieffer-Heeger Model
    2. II.2 Kitaev chain
    3. II.3 Kitaev honeycomb model
    4. II.4 Haldane model of graphene
  3. III Periodically driven topological systems
    1. III.1 Periodically driven Kitaev chain
    2. III.2 Topological invariants for periodic driving
    3. III.3 Periodically driven Kitaev honeycomb model
  4. IV Generically driven one-dimensional topological systems
    1. IV.1 Explicit symmetry breaking: unitary preparation of topological states
      1. IV.1.1 Preparation of Floquet topological state
      2. IV.1.2 Stroboscopic emergence of edge states
    2. IV.2 Dynamical symmetry breaking: A symmetry classification of non-equilibrium systems
    3. IV.3 Non-unitary dynamics of topological states
      1. IV.3.1 Lindblad master equation and steady state topology
      2. IV.3.2 Slow quenches in a dissipative environment
      3. IV.3.3 Majorana correlations in localized phases
      4. IV.3.4 Mixed state topological phase transitions
  5. V Generically driven two-dimensional Chern topological systems
    1. V.1 The unitary no-go theorem
    2. V.2 The non-unitary no-go theorem
    3. V.3 Dissipative preparation of mixed states
      1. V.3.1 Dissipators having a finite range
      2. V.3.2 Macroscopic electric polarization and mixed state Chern number
      3. V.3.3 Long-range dissipators and Floquet Chern insulators
    4. V.4 Going around the unitary no-go theorem: bilayer Chern systems
  6. VI Majorana modes in one-dimensional systems and transport
  7. VII Transport on surfaces of topological insulators
  8. VIII Other topological systems
  9. IX Summary and future directions
  10. A A pedagogical introduction to the geometry and topology of quantum states
    1. A.1 Adiabatic transformations
      1. A.1.1 Topology of adiabatically connected states
    2. A.2 The underlying geometry of quantum states
    3. A.3 Topological invariants
      1. A.3.1 Chern invariants
      2. A.3.2 Bott invariant
  11. B A perturbative expansion of the steady state Majorana correlations

##  I Introduction

Topological phases of quantum matter have been investigated extensively for
many years hasan10 ; qi11 . Typically, systems in such phases have bulk states
which are gapped and states at the boundaries which are gapless; the boundary
states lie within the bulk gap and they contribute to electronic transport at
low temperatures. In addition, the number of species of boundary states (for,
say, a given value of the boundary momentum) is given by a topological
invariant which can be calculated from the bulk properties of the system. The
properties of the boundary states (such as their energies) are often robust
against small amounts of disorder. The form of the topological invariant
depends on the dimensionality of the system and the symmetries of its
Hamiltonian; for example, it may be a winding number, a Chern number, a Bott
invariant, and so on. The transport properties of the boundary states are of
great interest as they often have universal features, such as quantized values
of the differential conductance or peaks in the differential conductance at
special values of the voltage bias.

The dynamics of quantum systems is another subject of considerable interest
from the point of view of topology. For instance, periodic driving or the
quenching of some parameters in the Hamiltonian can sometimes turn a non-
topological system which has no boundary states into a topological system
which has boundary states, and vice versa. Further, the value of the
topological invariant for such systems can be controlled by parameters such as
the amplitude and frequency of the periodic driving or the choice of initial
and final Hamiltonians of the quenching protocol. The effects of dissipation
on the dynamics is also of great interest.

In this review, we plan to provide a glimpse of various systems in which
topology and dynamics combine to play an important role. The plan of this
review is as follows. In Sec. II, we introduce some well-studied topological
systems such as the Su-Schrieffer-Heeger (SSH) model and the Kitaev chain in
one dimension, and the Kitaev honeycomb model and the Haldane model of
graphene in two dimensions, and the different phases that these models have
depending on the system parameters. In Sec. III, we discuss the rich structure
of phases which can arise when the Kitaev chain and the Kitaev honeycomb model
are driven periodically in time. In Sec. IV, we discuss some aspects of one-
dimensional topological systems which are out-of-equilibrium, in particular,
their topological classification and their non-unitary dynamics in the
presence of dissipation. We have primarily focused on the dynamical
preparation of topological states which also exhibit the dynamical emergence
of a topological bulk-boundary correspondence. In Sec. V, we study out-of-
equilibrium Chern topological systems in two dimensions, and the possibility
of using dynamics and dissipation to experimentally produce a desired
topological phase. In Sec. VI, we study the Kitaev chain and the effect that
the boundary states (which are zero energy Majorana modes) of such a system
have on the transport through such a system. Transport through a junction of
three Kitaev chains is also discussed.

Apart from one-dimensional symmetry protected topological insulators and Chern
insulators in two dimensions, there also exist various three-dimensional
systems which exhibit an intricate connection between topology and the
discrete symmetries, such as time-reversal symmetric three-dimensional
topological insulators hasan10 ; qi11 . Although the equilibrium topological
properties of three-dimensional systems are beyond the scope of this review,
we will nevertheless briefly discuss a few recent studies into the non-
equilibrium transport through topological surface states in such systems in
Sec. VII. For example, we discuss the effects of hybridization between the
opposite surfaces of a thin topological insulator and the effects of periodic
driving (which can be produced by applying electromagnetic radiation) on a
strip-like region of the top surface. In Sec. VIII, two unusual examples of
topological systems are presented. These include a spin system on a kagome
lattice in two dimensions where the magnon bands can have Berry curvature, and
a Josephson junction of three superconducting wires with different pairing
phases in which the bands of Andreev bound state at the junction have a Berry
curvature as a function of the pairing phases. We summarize our review in Sec.
IX. We end with some Appendices. Appendix A presents a pedagogical review of
the geometry of topological systems and the various topological invariants.
Appendix B discusses a perturbative expansion for the edge correlations which
arise in the dissipative dynamics of a Kitaev chain.

##  II Topology of fermionic systems in equilibrium

In this section we introduce a few minimalistic systems which manifest
topologically non-trivial phases of matter, particularly topological
insulators (we refer to Appendix A.1 for a mathematical interlude). In this
section, we introduce various quantum many-body topological systems and also
discuss the physical manifestations of topology in their properties. We will
see that non-trivial topological states of thermodynamically large systems
hold robust non-local signatures which are protected against local
perturbations and can therefore find extensive use in processing of quantum
information. Depending on the dimension and symmetries respected by the
respective Hamiltonians, every topological system can be classified into a
generic periodic table of topological insulators. To begin with, we consider
the simplest topological systems in one dimension and then move on to discuss
more complex models in two dimensions.

###  II.1 Su-Schrieffer-Heeger Model

Figure 1: The SSH chain with two sublattices (shown in red and blue and
labeled as sublattices 1 and 2 respectively). A typical unit cell containing a
site from each sublattice has been circled with dotted lines.

The SSH modelshen12 ; bernevig13  belongs to the BDI symmetry class of
topological insulators and is the simplest one-dimensional model exhibiting an
underlying topological structure and end statesasboth16 . Physically, it
describes a lattice with a two atom sublattice structure in which the intra-
cell hopping amplitude is in general different from the inter-cell hopping
amplitude. The Hamiltonian for the SSH model can be written in terms of the
(spin polarized) fermion creation and annihilation operators as,

$H=\sum_{n=1}^{\mathcal{N}}(vc^{\dagger}_{n,1}c_{n,2}+wc^{\dagger}_{n,2}c_{n+1,1}+{\rm
H.c.}),$ (1)

where H.c. denotes the Hermitian conjugate, $v$ and $w$ are the intra-cell
(see Fig. 1) and inter-cell hopping amplitudes respectively and $\mathcal{N}$
is the total number of unit cells in the chain. The complex fermionic operator
$c^{\dagger}_{n,i}$ ($c_{n,i}$) creates (destroys) a fermion in the sublattice
position $i\leavevmode\nobreak\ (i=1,2)$ of the $n^{th}$ unit cell and
satisfies the fermionic anti-commutation rules,

$\\{c^{\dagger}_{p},c_{q}\\}=\delta_{pq}\text{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
}\\{c_{p},c_{q}\\}=\\{c^{\dagger}_{p},c^{\dagger}_{q}\\}=0.$ (2)

Note that with periodic boundary conditions, the lattice has a discrete
translational symmetry for a translation by a complete unit cell, i.e., the
system remains invariant under a translation by two consecutive lattice
points. Utilizing this translational invariance, one can rewrite the
Hamiltonian in Fourier space in terms of the spinors $\Psi_{k}$, defined as

$\Psi_{k}=\left(\begin{array}[]{c}c_{k,1}\\\ c_{k,2}\\\
\end{array}\right)=\frac{1}{\mathcal{\sqrt{N}}}\sum\limits_{l=1}^{\mathcal{N}}e^{-ikl}\left(\begin{array}[]{c}c_{l,1}\\\
c_{l,2}\\\ \end{array}\right),$ (3)

such that the Hamiltonian in Eq. (1) is decoupled for each conserved
quasimomentum mode $k$ as,

$\begin{split}H=\bigoplus_{k}\mathcal{H}_{k}=\bigoplus_{k}\Psi^{\dagger}_{k}H(k)\Psi_{k},\\\
H(k)=\vec{h}(k)\cdot\vec{\sigma},\end{split}$ (4)

where $\vec{\sigma}$ is the vector of Pauli matrices denoting pseudo-spin
operators and $k=2\pi m/\mathcal{N}$, where $m=0,1,\dots,\mathcal{N}$. For the
SSH model, the functions $\vec{h}(k)$ and hence, the single-particle
Hamiltonian $H(k)$ are found to be

$\displaystyle h_{x}(k)$ $\displaystyle=$ $\displaystyle
Re(v)+|w|\cos(k+arg(w)),$ $\displaystyle h_{y}(k)$ $\displaystyle=$
$\displaystyle-Im(v)+|w|\sin(k+arg(w)),$ $\displaystyle\noindent h_{z}(k)$
$\displaystyle=$ $\displaystyle 0,$ (5)

with the lattice parameter set equal to unity.

The single-particle Hamiltonian $H(k)$ has the following eigenvalue spectrum,

$E(k)=\pm|\vec{h(k)}|,$ (6)

and the respective eigenvectors are

$|\pm\rangle=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\pm e^{-i\phi(k)}\\\
1\\\ \end{array}\right),$ (7)

where $\phi=\tan^{-1}\left(\frac{h_{y}}{h_{x}}\right)$. It is clear from the
Eq. (II.1), that $\vec{h}(k)$ is periodic in $k$, $0\leq k\leq 2\pi$ (the
Brillouin zone (BZ) has the geometry of $\mathcal{S}^{1}$). Sweeping over the
BZ is thus equivalent to tracing out a closed loop in the parameter space,
with $k$ parametrizing the curve. Further, it follows from Eq. (6) that the
system is in the ground state at half-filling, i.e., when only the lowest band
is completely occupied. The different phases of the SSH model are therefore
classified through the bulk topological winding number $\nu$, which is given
as (see Eq. (244) of Appendix A),

$\nu=\frac{i}{2\pi}\oint_{BZ}\langle\psi_{0}^{k}|\partial_{k}|\psi_{0}^{k}\rangle
dk=\frac{i}{2\pi}\int^{\pi}_{-\pi}dk\frac{d}{dk}\ln(h_{x}+ih_{y}).$ (8)

The second equality in the above equation can be understood as follows. We
note that $\vec{h}(k)$ traces out a closed curve over a plane as $k$ varies
over the first BZ. The winding number is proportional to the change in the
argument of $\vec{h}(k)$ as $k$ varies over the first BZ and is thus integer
quantized. If the curve does not enclose the origin, $\nu$ is zero (i.e., if
$|v|>|w|$). On the other hand, $\nu$ is one if the curve encloses the origin
(i.e., if $|v|<|w|$). It is also evident from Eq. (6) that the energy gap
between the two bands vanishes at $|v|=|w|$ and the winding number becomes
undefined. Thus if the gap is not closed, $\nu$ is well-defined and is robust
to external changes in the Hamiltonian and hence is a topological invariant
clearly demarcating the topologically trivial and non-trivial phases.

It is also known that a topological non-triviality of the bulk system induces
conducting boundary states under open boundary conditions. This bulk-boundary
correspondenceshen12  is simple to understand in the particular situation of
the SSH Hamiltonian when $v,w\in\mathbb{R}$. Consider the two different
topological phases characterized by the sign of the quantity $m=v-w$ separated
by the QCP at $v=w$ ($m=0$). At the QCP, the energy gap vanishes at
$k^{*}=\pi$. We therefore focus on the low energy spectrum of the SSH model
near the QCP ($m\sim 0$) such that $k=k^{*}+\delta k$, where $\delta k<<1$.
Plugging in these approximations in the Hamiltonian in Eq. (II.1) one sees
that

$H_{c}(k)\simeq-\delta k\sigma^{y}+m\sigma^{x}.$ (9)

The linear dispersion at the QCP at $m=0$ is why it is called the Dirac point.
Also, the low-energy Hamiltonian $H_{c}$ in real space translates to the one-
dimensional Dirac Hamiltonian with $m$ being the mass,

$H_{c}=i\partial_{x}\sigma^{y}+m\sigma^{x}.$ (10)

This equation can then be utilized to model a junction of two systems along
the x-axis such that there is a mass discontinuity. Specifically, we take the
mass to be of different signs on the positive and negative x-axis respectively
signifying that the two systems across the junction are topologically
inequivalent. It is straightforward to check that under such a mass
discontinuity, the Dirac Hamiltonian hosts zero energy eigenstates which are
exponentially localized at the junction. This exemplifies the existence of
edge states localized at the boundary of topologically non-trivial systems.
These edge states are protected by the bulk topology of the thermodynamically
large system. In other words, they are robust against all such local
perturbations which are incapable of changing the bulk topology. The
thermodynamic limit is necessary as otherwise the exponentially localized
states at the opposite edges might hybridize to produce finite tunneling into
the bulk.

The topological classification of non-interacting many-body quantum systems is
performed by considering three different discrete symmetries viz., the time-
reversal symmetry ($\mathcal{T}$), the particle-hole symmetry ($\mathcal{P}$)
and the sublattice (chiral) symmetry ($\mathcal{S}$). The constraints imposed
upon the Hamiltonian of a system possessing the above symmetries in the quasi
momentum basis are expressed as,

$\begin{split}\mathcal{T}^{-1}H(k)\mathcal{T}=H(-k),\\\
\mathcal{P}^{-1}H(k)\mathcal{P}=-H(-k),\\\
\mathcal{S}^{-1}H(k)\mathcal{S}=-H(k),\end{split}$ (11)

where $\mathcal{T}$ and $\mathcal{P}$ are anti-unitary operators such that
$\mathcal{T}^{2}=\pm\mathbb{I}$ and $\mathcal{P}^{2}=\pm\mathbb{I}$, whereas
$\mathcal{S}$ is an unitary operator satisfying $\mathcal{S}^{2}=\mathbb{I}$
and $\mathbb{I}$ is the $2\times 2$ identity operator. We also note that the
sublattice symmetry is a combined effect of the time-reversal symmetry and the
particle hole symmetry as,

$\mathcal{S}=\mathcal{T}\mathcal{P}.$ (12)

It is now evident from the Hamiltonian of the SSH model in Eqs. (4) and (II.1)
and the symmetry transformations in Eq. (11) that the SSH model is symmetric
under the sublattice transformation $\mathcal{S}=\bigotimes_{k}\sigma^{z}$
which results in the vanishing of $h_{z}(k)$. Also, if the hopping
coefficients $v$ and $w$ are real, the Hamiltonian possesses time-reversal
symmetry $\mathcal{T}=\bigotimes_{k}\mathcal{K}$, $\mathcal{K}$ being the
complex conjugation operator. Hence, it is clear from Eq. (12) that the system
is also symmetric under the particle-hole/charge conjugation operation with
$\mathcal{P}=\bigotimes_{k}\mathcal{K}\sigma^{z}$. As such the SSH model is
classified in the BDI class of Hamiltonians within the topological
classification scheme. It must also be noted that the Hamiltonian $H(k)$ spans
a $\mathcal{S}^{1}$ on a plane because $h_{z}(k)=0$ is satisfied for all $k$.
This in turn is guaranteed when the system lies in the BDI class or, in
particular, respects the chiral/sublattice symmetry $\mathcal{S}$. In other
words, the chiral symmetry protects the topological phases of the system and
hence, these phases are also known as symmetry protected topological (SPT)
phases.

###  II.2 Kitaev chain

Similar to the SSH chain, the Kitaev chain is a one dimensional lattice system
of spin-polarized fermions that represents a prototype model of $p$-wave
superconductors (Since the system is spin-polarized we will not write the spin
label below). As we will see below, it has topological phases in which a long
system has a pair of zero energy Majorana modes, with one localized at each
end kitaev01 ; lutchyn10 ; oreg10 ; gott1 ; gott2 . For a chain of linear
dimension $N$, with the lattice parameter once again set to unity, the many-
body Hamiltonian is given by

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N-1}[\gamma(c_{n}^{\dagger}c_{n+1}+c_{n+1}^{\dagger}c_{n})+\Delta(c_{n}c_{n+1}+c_{n+1}^{\dagger}c_{n}^{\dagger})]$
(13) $\displaystyle-\sum_{n=1}^{N}\mu(2c_{n}^{\dagger}c_{n}-1),$

where $\gamma$ and $\Delta$ are the hopping and superconducting pairing
amplitudes between neighboring sites, respectively, while $\mu$ is the on-site
chemical potential. All of these parameters are considered to be real and one
may assume that $\gamma>0$ without loss of generality.

The topological features of the Kitaev chain, particularly the existence of
localized end states, are more explicitly manifested in the Majorana
representation. Defining the Majorana operators,

$a_{2n-1}=c_{n}+c_{n}^{\dagger}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
a_{2n}=i(c_{n}-c_{n}^{\dagger}),$ (14)

where $n=1,2,\cdots,N$, Eq. (13) assumes the form

$\displaystyle H$ $\displaystyle=$ $\displaystyle i\leavevmode\nobreak\
\sum_{n=1}^{N-1}\leavevmode\nobreak\ [\leavevmode\nobreak\
J_{x}a_{2n}a_{2n+1}\leavevmode\nobreak\ -\leavevmode\nobreak\
J_{y}a_{2n-1}a_{2n+2}]$ $\displaystyle+\leavevmode\nobreak\
i\leavevmode\nobreak\ \sum_{n=1}^{N}\leavevmode\nobreak\ \mu a_{2n-1}a_{2n},$
$\displaystyle J_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\gamma-\Delta)\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ J_{y}=\frac{1}{2}(\gamma+\Delta).$ (15)

The Majorana operators are Hermitian operators satisfying
$\\{a_{m},a_{n}\\}=2\delta_{mn}$. A careful look into the above Hamiltonian
uncovers an interesting connection with the SSH Hamiltonian. To be precise,
the Kitaev chain in the Majorana representation is exactly the SSH chain (see
Eq. (9)) with Majorana fermions replacing Dirac fermionic operators.

Furthermore, a Jordan-Wigner transformation of the Majorana operators lieb ,
given as,

$\displaystyle a_{2n-1}$ $\displaystyle=$
$\displaystyle\left(\prod_{j=1}^{n-1}\sigma_{j}^{z}\right)\leavevmode\nobreak\
\sigma_{n}^{x},$ $\displaystyle a_{2n}$ $\displaystyle=$
$\displaystyle\left(\prod_{j=1}^{n-1}\sigma_{j}^{z}\right)\leavevmode\nobreak\
\sigma_{n}^{y},$ (16)

maps the above system to a spin-1/2 $XY$ chain placed in a magnetic field
pointing in the $\hat{z}$ direction, with the Hamiltonian modified as,

$H\leavevmode\nobreak\ =\leavevmode\nobreak\
-\sum_{n=1}^{N-1}\leavevmode\nobreak\ [\leavevmode\nobreak\
J_{x}\sigma_{n}^{x}\sigma_{n+1}^{x}\leavevmode\nobreak\ +\leavevmode\nobreak\
J_{y}\sigma_{n}^{y}\sigma_{n+1}^{y}]\leavevmode\nobreak\
-\sum_{n=1}^{N}\leavevmode\nobreak\ \mu\sigma_{n}^{z},$ (17)

where the $\sigma_{n}^{a}$ denote the Pauli matrices at site $n$. In all the
numerical calculations that follow, we set $\gamma=-\Delta$ which implies
$J_{y}=0$ and $J_{x}=\gamma$; the system is then equivalent to an Ising model
(with interaction $J_{x}$) in a transverse magnetic field $\mu$.

Figure 2: Phase diagram of the model in Eq. (13) as a function of $\mu/\gamma$
and $\Delta/\gamma$. Phases I and II are topological while III and IV are non-
topological.

To elucidate the bulk properties, let us first consider a chain with periodic
boundary conditions. Substituting the Fourier transform
$c_{k}=\frac{1}{\sqrt{N}}\leavevmode\nobreak\ \sum_{n=1}^{N}c_{n}e^{ikn}$ in
Eq. (13), we obtain,

$\displaystyle H$ $\displaystyle=$ $\displaystyle
2(\gamma-\mu)c_{0}^{\dagger}c_{0}\leavevmode\nobreak\ +\leavevmode\nobreak\
2(-\gamma-\mu)c_{\pi}^{\dagger}c_{\pi}$ (21)
$\displaystyle+\leavevmode\nobreak\ \sum_{0<k<\pi}\leavevmode\nobreak\
\left(\begin{array}[]{cc}c_{k}^{\dagger}&c_{-k}\end{array}\right)\leavevmode\nobreak\
h_{k}\leavevmode\nobreak\ \left(\begin{array}[]{c}c_{k}\\\
c_{-k}^{\dagger}\end{array}\right),$ $\displaystyle h_{k}$ $\displaystyle=$
$\displaystyle{\vec{h}}(k)\cdot{\vec{\tau}}=2(\gamma\cos
k-\mu)\leavevmode\nobreak\ \tau^{z}\leavevmode\nobreak\ +\leavevmode\nobreak\
2\Delta\sin k\leavevmode\nobreak\ \tau^{y},$ (22)

where ${\vec{h}}(k)=(h^{x}(k),h^{y}(k),h^{z}(k))=(0,2\Delta\sin k,2(\gamma\cos
k-\mu))$ and ${\vec{\tau}}=(\tau^{x},\tau^{y},\tau^{z})$ is the vector of
Pauli matrices denoting pseudo-spin. Note that the effective BZ is halved
($0<k<\pi$) as the pairing term in Eq. (13) couples the momenta $k$ and $-k$.
The dispersion relation follows from Eq. (22) and is given by gott1 ; gott2

$E_{k}\leavevmode\nobreak\ =\leavevmode\nobreak\ \pm\sqrt{4(\gamma\cos
k-\mu)^{2}\leavevmode\nobreak\ +\leavevmode\nobreak\ 4\Delta^{2}\sin^{2}k}.$
(23)

As shown in Fig. 2, the Kitaev chain, or equivalently the $XY$ spin chain, can
exist in four different phases gott1 ; gott2 . In the spin chain
representation, phases I and II correspond to long-range ferromagnetic orders
of $\sigma^{x}$ and $\sigma^{y}$, respectively. On the other hand, phases III
and IV correspond to paramagnetic phases with no long-range order. The four
phases are separated from each other by quantum critical lines where the
energy $E_{k}$ vanishes for some values of $k$. The critical lines are given
by $\mu/\gamma=\pm 1$ for all values of $\Delta$, and $-1\leq\mu/\gamma\leq 1$
for $\Delta=0$. Interestingly, these phases are also topologically distinct
and can be characterized by a winding number.

To define the winding number, we proceed in a similar way as we did in the
case of the SSH chain. First, note that the vector ${\vec{h}}(k)$ is confined
to the $y-z$ plane. Let $\phi_{k}=\tan^{-1}(h^{z}(k)/h^{y}(k))$ denote the
angle sustained by the vector $\vec{h}(k)$ with respect to the $\hat{z}$ axis.
The winding number can be defined as in Eq. (8) by following the change in
$\phi_{k}$ as we go around the full BZ ($-\pi\leq k\leq\pi$) niu ; tong ,
i.e.,

$W\leavevmode\nobreak\ =\leavevmode\nobreak\
\int_{-\pi}^{\pi}\leavevmode\nobreak\ \frac{dk}{2\pi}\leavevmode\nobreak\
\frac{d\phi_{k}}{dk}.$ (24)

This can take any integer value and is a topological invariant, namely, it
does not change under small changes in $h_{k}$ unless $h_{k}$ happens to pass
through zero for some value of $k$ in which case the winding number becomes
ill-defined; this can only happen if the energy $E_{k}=0$ at some value of $k$
which means that the bulk gap is zero. In a gapped phase, therefore, Eq. (24)
defines a $Z$-valued topological invariant.

In the ferromagnetic phases, it is easy to check that $W=\pm 1$; the system is
said to exist in a topologically non-trivial phase. As in the case of the SSH
chain, in the topologically non-trivial phase, a long and open Kitaev chain
hosts a pair of zero energy Majorana modes at each end gott2  and
correspondingly, the two-point correlation function of the Majorana end modes
in the ground state of the Hamiltonian, defined as $\bra{GS}\theta\ket{GS}$,
where $\theta=ia_{1}a_{2N}$, also remains finite and approaches unity as
$\mu\to 0$. This Majorana end correlation in a thermodynamically large system
is a remarkable manifestation of topologically protected long range
correlations. This can be verified by considering the extreme case $J_{x}>0$
and $J_{y}=\mu=0$ in Eq. (15), in which case the Hamiltonian is independent of
the terms $a_{1}$ and $a_{2N}$; hence the corresponding modes have zero
energies. In the paramagnetic phase, one finds $W=0$ and hence the system is
topologically trivial. Also, no zero energy Majorana modes are present at
either end of an open chain. The winding number thus exhibits a bulk-boundary
correspondence.

It is important to analyze the role played by symmetry in the existence of the
topological structure discussed above. The crucial ingredient is once again
the confinement of the vector $\vec{h}(k)$ to a 2d parameter space, which
results from the underlying chiral symmetry of the Kitaev chain. This can be
verified from the symmetry transformation
$\tau_{x}\vec{h}(k)\tau_{x}=-\vec{h}(k)$; which is satisfied in the absence of
any $x-$ component of the vector $\vec{h}(k)$. The Kitaev chain also possesses
the particle-hole, $\mathcal{P}\vec{h}(k)\mathcal{P}^{\dagger}=-\vec{h}(-k)$,
and the time-reversal symmetries,
$\mathcal{T}\vec{h}(k)\mathcal{T}^{\dagger}=\vec{h}(-k)$. The symmetry
operators can be identified as $\mathcal{P}=\mathcal{K}\tau_{x}$ and
$\mathcal{T}=\mathcal{K}$, where $\mathcal{K}$ is the complex conjugation
operator. In general, the presence of any two of these symmetries guarantees
that the third must also be present. This establishes that the completely
symmetric Kitaev chain also lies in the BDI symmetry class.

###  II.3 Kitaev honeycomb model

The Kitaev model is a two-dimensional model of spin-$1/2$’s placed on the
sites of a honeycomb lattice, with a Hamiltonian of the form

$H\leavevmode\nobreak\ =\leavevmode\nobreak\
\sum_{j+l=even}(J_{1}\sigma^{x}_{j,l}\sigma^{x}_{j+1,l}+J_{2}\sigma^{y}_{j-1,l}\sigma^{y}_{j,l}+J_{3}\sigma^{z}_{j,l}\sigma^{z}_{j,l+1}),$
(25)

where $j,l$ are the column and row indices, respectively. For simplicity, we
assume that $J_{1},J_{2},J_{3}\geq 0$.

Figure 3: Kitaev model on the honeycomb lattice with $xx$ coupling $J_{1}$,
$yy$ coupling $J_{2}$ and $zz$ coupling $J_{3}$. $\vec{M_{1}}$ and
$\vec{M_{2}}$ are the spanning vectors of the lattice, and $A$ and $B$ denote
the two sites of a unit cell. (See Ref. thakurathi14, ).

Fig. 3 shows a schematic representation of the honeycomb lattice. The unit
cells are chosen as the vertical bonds, with the sublattices labeled as $A$
and $B$; these have $j+l$ equal to odd and even integers, respectively. It is
convenient to set the nearest-neighbor distance to be $1/\sqrt{3}$, which
allows us to label each unit cell by the vector,
$\vec{n}=\hat{i}n_{1}+(\frac{1}{2}\hat{i}+\frac{\sqrt{3}}{2}\hat{j})n_{2}$,
where $n_{1},\leavevmode\nobreak\ n_{2}$ are integers that can be calculated
from the coordinates of the $B$ site in the corresponding unit cell as
$n_{1}=(j-l)/2$ and $n_{2}=l$. Fig. 3 shows the spanning vectors
$\vec{M_{1}}=\frac{1}{2}\hat{i}+\frac{\sqrt{3}}{2}\hat{j}$ and
$\vec{M_{2}}=\frac{1}{2}\hat{i}-\frac{\sqrt{3}}{2}\hat{j}$ which join some
neighboring unit cells.

As with the Kitaev chain, the topological aspects of the Kitaev model are also
best understood in the Majorana representation. We define the Majorana
operators kitaev06 ; feng07 ; chen07 ; nussinov08 ; baskaran07 ; lee07

$\displaystyle\hat{a}_{j,l}$ $\displaystyle=$
$\displaystyle\left(\prod_{i,k<l}\sigma^{z}_{i,k}\right)\left(\prod_{i<j}\sigma^{z}_{i,l}\right)\sigma^{y}_{j,l},$
$\displaystyle\hat{a}^{\prime}_{j,l}$ $\displaystyle=$
$\displaystyle\left(\prod_{i,k<l}\sigma^{z}_{i,k}\right)\left(\prod_{i<j}\sigma^{z}_{i,l}\right)\sigma^{x}_{j,l}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ j+l=\text{even},$
$\displaystyle\hat{b}_{j,l}$ $\displaystyle=$
$\displaystyle\left(\prod_{i,k<l}\sigma^{z}_{i,k}\right)\left(\prod_{i<j}\sigma^{z}_{i,l}\right)\sigma^{x}_{j,l},$
$\displaystyle\hat{b}^{\prime}_{j,l}$ $\displaystyle=$
$\displaystyle\left(\prod_{i,k<l}\sigma^{z}_{i,k}\right)\left(\prod_{i<j}\sigma^{z}_{i,l}\right)\sigma^{y}_{j,l}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ j+l=\text{odd},$

These are Hermitian operators satisfying the anticommutation relations
$\\{\hat{a}_{m,n},\hat{a}_{m^{\prime},n^{\prime}}\\}=2\delta_{mm^{\prime}}\delta_{nn^{\prime}}$,
$\\{\hat{a}^{\prime}_{m,n},\hat{a}^{\prime}_{m^{\prime},n^{\prime}}\\}=2\delta_{mm^{\prime}}\delta_{nn^{\prime}}$,
$\\{\hat{b}_{m,n},\hat{b}_{m^{\prime},n^{\prime}}\\}=2\delta_{mm^{\prime}}\delta_{nn^{\prime}}$,
$\\{\hat{b}^{\prime}_{m,n},\hat{b}^{\prime}_{m^{\prime},n^{\prime}}\\}=2\delta_{mm^{\prime}}\delta_{nn^{\prime}}$,
and
$\\{\hat{a}_{m,n},\hat{a}^{\prime}_{m^{\prime},n^{\prime}}\\}=\\{\hat{a}_{m,n},\hat{b}_{m^{\prime},n^{\prime}}\\}=\\{\hat{a}_{m,n},\hat{b}^{\prime}_{m^{\prime},n^{\prime}}\\}=\\{\hat{a}^{\prime}_{m,n},\hat{b}_{m^{\prime},n^{\prime}}\\}=\\{\hat{a}^{\prime}_{m,n},\hat{b}^{\prime}_{m^{\prime},n^{\prime}}\\}=\\{\hat{b}_{m,n},\hat{b}^{\prime}_{m^{\prime},n^{\prime}}\\}=0$.
The Hamiltonian then takes the form

$H\leavevmode\nobreak\ =\leavevmode\nobreak\ i\leavevmode\nobreak\
\sum_{\vec{n}}(J_{1}\hat{b}_{\vec{n}}\hat{a}_{\vec{n}-\vec{M_{1}}}+J_{2}\hat{b}_{\vec{n}}\hat{a}_{\vec{n}+\vec{M_{2}}}+J_{3}\hat{D}_{\vec{n}}\hat{b}_{\vec{n}}\hat{a}_{\vec{n}}),$
(28)

where
$\hat{D}_{\vec{n}}=i\hat{b}^{\prime}_{\vec{n}}\hat{a}^{\prime}_{\vec{n}}$. It
is straightforward to see that the operators $\hat{D}_{\vec{n}}$ commute with
each other and also with the Hamiltonian; hence they represent conserved
quantities. In particular, the eigenvalues of $\hat{D}_{\vec{n}}$ can assume
the values $\pm 1$, independently, for each $\vec{n}$, thereby decomposing the
$2^{N}$-dimensional Hilbert space ($N$ being the total number of sites) into
$2^{N/2}$ sectors. Further, it is known that the ground state of the model
lies in the sector in which the eigenvalue of $\hat{D}_{\vec{n}}$ equals $1$
for all $\vec{n}$ and as such, we will restrict ourselves to this sector. We
note that the operators $\hat{D}_{\vec{n}}$, when expressed in terms of spin
operators using Eqs. (II.3) and (LABEL:eq_2dmajob), are highly non-local
string operators.

Assuming periodic boundary conditions, the Majorana operators can be Fourier
transformed as

$\displaystyle\hat{a}_{\vec{n}}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{4}{N}}\sum_{\vec{k}\in\frac{1}{2}\text{BZ}}{}(\hat{a}_{\vec{k}}e^{i\vec{k}\cdot\vec{n}}+\hat{a}_{\vec{k}}^{\dagger}e^{-i\vec{k}\cdot\vec{n}}),$
$\displaystyle\hat{b}_{\vec{n}}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{4}{N}}\sum_{\vec{k}\in\frac{1}{2}\text{BZ}}{}(\hat{b}_{\vec{k}}e^{i\vec{k}\cdot\vec{n}}+\hat{b}_{\vec{k}}^{\dagger}e^{-i\vec{k}\cdot\vec{n}}),$
(29)

which satisfy the anticommutation relations
$\\{\hat{a}_{\vec{k}},\hat{a}_{\vec{k}^{\prime}}^{\dagger}\\}=\\{\hat{b}_{\vec{k}},\hat{b}_{\vec{k}^{\prime}}^{\dagger}\\}=\delta_{\vec{k},\vec{k}^{\prime}}$.
Note that the sums over $\vec{k}$ in Eq. (29) only go over half the BZ; a
convenient choice of the BZ is given by a rhombus whose vertices lie at
$(k_{x},k_{y})=(\pm 2\pi,0)$ and $(0,\pm 2\pi/\sqrt{3})$. The Hamiltonian in
Eq. (28) can then be written in a form similar to the Kitaev chain Hamiltonian
in Eq. (22)

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{\vec{k}\in\frac{1}{2}\text{BZ}}\leavevmode\nobreak\
\left(\begin{array}[]{cc}\hat{a}_{\vec{k}}^{\dagger}&\hat{b}_{\vec{k}}^{\dagger}\end{array}\right)\leavevmode\nobreak\
H_{k}\leavevmode\nobreak\ \left(\begin{array}[]{c}\hat{a}_{\vec{k}}\\\
\hat{b}_{\vec{k}}\end{array}\right),$ (33) $\displaystyle H_{\vec{k}}$
$\displaystyle=$ $\displaystyle
2[J_{1}\sin(\vec{k}.\vec{M}_{1})-J_{2}\sin(\vec{k}.\vec{M}_{2})]\tau^{x}$ (34)
$\displaystyle+2[J_{3}+J_{1}\cos(\vec{k}.\vec{M}_{1})+J_{2}\cos(\vec{k}.\vec{M}_{2})]\tau^{y}.$

The dispersion relation can be derived from Eq. (34); it consists of two bands
with energies

$\displaystyle E_{\vec{k}}^{\pm}$ $\displaystyle=$ $\displaystyle\pm
2[\\{J_{1}\sin(\vec{k}.\vec{M}_{1})-J_{2}\sin(\vec{k}.\vec{M}_{2})\\}^{2}$
$\displaystyle+\\{J_{3}+J_{1}\cos(\vec{k}.\vec{M}_{1})+J_{2}\cos(\vec{k}.\vec{M}_{2})\\}^{2}]^{1/2}.$

To construct a phase diagram of the model, it is convenient to impose the
condition $J_{1}+J_{2}+J_{3}=1$, which describes points lying within (or on)
an equilateral triangle. This triangle can be divided into four smaller
equilateral triangles as shown in Fig. 4, namely, $A_{x}$, $A_{y}$, $A_{z}$
and $B$. It turns out kitaev06  that the system is gapped in the three $A$
phases, with $E_{\vec{k}}$ being non-zero for all $\vec{k}$, and is gapless in
the $B$ phase, with $E_{\vec{k}}=0$ for some value of $\vec{k}$ whose value
depends on the location of the point in that phase. The four phases are
separated from each other by quantum critical lines where one of the couplings
is equal to the sum of the other two.

Figure 4: (a) Phase diagram for zigzag edge states in the triangle with
$J_{1}+J_{2}+J_{3}=1$. Majorana modes exist in the regions $A_{z}$ and $B$.
(b) Zigzag and (c) armchair edge configuration of the Kitaev model. Majorana
fermions with a momentum $k$ along the edge are indicated. (See Ref.
thakurathi14, ).

Figure 5: (a) Phase diagram for armchair edge states in the triangle with
$J_{1}+J_{2}+J_{3}=1$. Majorana modes of type $\hat{b}$ ($\hat{a}$) exist in
the regions $A_{y}$ ($A_{x}$) and the left (right) half of $B$. These regions
are indicated by dark (light) shades. (b) Armchair edge states and (c) corner
states for a system with $N_{x}\times N_{y}=27\times 14$ sites, with
$J_{1}=0.7,\leavevmode\nobreak\ J_{2}=0.15$ and $J_{3}=0.15$. (See Ref.
thakurathi14, ).

The nature of the edge modes which can appear in the Kitaev model is dependent
on the geometry of the edge themselves. Depending on the orientation, the
edges in a semi-infinite lattice can have different structures, two of such
structures studied commonly are the zigzag and armchair nakada ; kohmoto
edges. Assuming that the edges are infinitely long, translational invariance
permits us to label the edge states by their momentum $k$. To find the edge
states, we first write the Hamiltonian in Eq. (28) in the form

$H\leavevmode\nobreak\ =\leavevmode\nobreak\ 2i\leavevmode\nobreak\
\sum_{\alpha\beta}\leavevmode\nobreak\ \hat{b}_{\beta}\leavevmode\nobreak\
L_{\beta\alpha}\leavevmode\nobreak\ \hat{a}_{\alpha},$ (36)

where $\alpha,\leavevmode\nobreak\ \beta$ label the sites, and
$L_{\alpha\beta}$ is a real matrix. We now use the Heisenberg equations of
motion $d\hat{a}_{\alpha}/dt\leavevmode\nobreak\ =\leavevmode\nobreak\
i\leavevmode\nobreak\ [H,\hat{a}_{\alpha}]$ and similarly for
$\hat{b}_{\beta}$. We then obtain

$\displaystyle\frac{d\hat{a}_{\alpha}}{dt}$ $\displaystyle=$
$\displaystyle-4\leavevmode\nobreak\ \sum_{\beta}\leavevmode\nobreak\
\hat{b}_{\beta}\leavevmode\nobreak\ L_{\beta\alpha},$
$\displaystyle\frac{d\hat{b}_{\beta}}{dt}$ $\displaystyle=$ $\displaystyle
4\leavevmode\nobreak\ \sum_{\alpha}\leavevmode\nobreak\
L_{\beta\alpha}\leavevmode\nobreak\ \hat{a}_{\alpha}.$ (37)

Zero energy modes, if they exist, can only appear in the form of localized
edge states and therefore can be found by equating the above Heisenberg
equations to zero. We will henceforth denote wave functions by alphabet
letters without hats (such as $a$ and $b$) to distinguish them from operators
which are denoted by $\hat{a}$ and $\hat{b}$.

In Fig. 4, the wave functions for the Majorana operators of type $\hat{b}$ are
denoted by $b_{m,n}$, where $n$ goes from $-\infty$ to $\infty$ and increases
towards the right along the edge, and $m=1,2,3,...$ increases as we go down
away from the top edge and into the bulk of the system. Further, they can be
decomposed into Fourier components as $b_{m,n}=b_{m}e^{ikn}$ or
$b_{m}e^{ik(n+1/2)}$ depending on whether $m$ is odd or even. Similarly, the
wave functions for Majorana operators of type $\hat{a}$ (not shown in Fig. 4)
are denoted by $a_{m,n}=a_{m}e^{ikn}$ or $a_{m}e^{ik(n+1/2)}$. Substituting in
Eq. (37) equating with zero reveals that there exists zero energy solutions
with $a_{m}=0$ for all $m$ and

$J_{1}b_{m}e^{ikn}+J_{2}b_{m}e^{ik(n+1)}+J_{3}b_{m+1}e^{ik(n+\frac{1}{2})}=0$
(38)

for all $m\geq 1$. From Fig. 4, it is clear that the zero energy solution
corresponds to localized edge states at the top where the Majorana operators
are of type $\hat{b}$. Further, assuming $b_{m}=(\lambda_{k})^{m}$ leads to
$\lambda_{k}=-(J_{1}e^{-ik/2}+J_{2}e^{ik/2})/J_{3}$. For a normalizable edge
state, we require $|\lambda_{k}|<1$; this occurs if

$\cos k\leavevmode\nobreak\ <\leavevmode\nobreak\
\frac{J_{3}^{2}-J_{1}^{2}-J_{2}^{2}}{2J_{1}J_{2}}.$ (39)

The above condition holds true for all values of $k$ in region $A_{z}$ and for
a finite range of values of $k$ in region $B$. On the contrary, no solution
exist for Eq. (39) in the regions $A_{x}$ and $A_{y}$, which implies the
absence of edge states in these regions. A similar analysis for zero energy
solutions with $b_{m}=0$ for all $m$ shows the presence of localized edge
states at the bottom of the semi-infinite lattice (see Fig. 4) where the
Majorana operators are of type $\hat{a}$.

For armchair edges (see Fig. 4), the $\hat{a}$ and $\hat{b}$ Majorana
operators have wave functions given by $a_{m}e^{ikn}$ and $b_{m}e^{ikn}$,
respectively. In this case, the Heisenberg equations of motion have zero
energy solutions with $a_{m}=0$ for all $m$, provided that

$J_{1}b_{m}+J_{2}b_{m+2}+J_{3}b_{m+1}e^{-ik}\leavevmode\nobreak\
=\leavevmode\nobreak\ 0$ (40)

for all $m\geq 1$, and

$J_{2}b_{2}+J_{3}b_{1}e^{-ik}\leavevmode\nobreak\ =\leavevmode\nobreak\ 0.$
(41)

Assuming $b_{m}=(\lambda_{k})^{m}$, Eq. (40) has two roots given by
$\lambda_{k\pm}=(-J_{3}e^{-ik}\pm\sqrt{J_{3}^{2}e^{-2ik}-4J_{1}J_{2}})/(2J_{2}).$
Eqs. (40-41) imply that a normalizable edge state will exist if both
$|\lambda_{k\pm}|<1$. We find that this occurs if $J_{1}<J_{2}$ and $J_{3}\leq
J_{1}+J_{2}$. The corresponding region in the phase diagram is shown as a dark
shaded region on the left side of Fig. 5. A similar analysis shows that zero
energy edge modes of type $\hat{a}$ (i.e., with $b_{m}=0$ for all $m$) also
exist, the corresponding region is shown as light shaded in Fig. 5.

Finally, it is important to consider a finite system with armchair edges along
one direction and zigzag edges along the other. Interestingly, as shown in
Figs. 5 and 5, while some edge modes are found to be localized along the
armchair edges, the rest appear only at the corners.

We remark here that although the Kitaev honeycomb model is capable of hosting
Majorana edge modes, the latter do not arise as a consequence of an underlying
bulk topological structure of the complete two-dimensional Hamiltonian. This
can be argued from a symmetry analysis of the system. From Eq. (34), it is
straightforward to see that the Hamiltonian possesses chiral, time-reversal as
well as particle-hole symmetries, with the symmetry operators identified as
$\mathcal{S}=\tau^{z}$, $\mathcal{T}=i\tau^{y}$ and $\mathcal{P}=\tau^{x}$,
respectively. The model thus belongs to the BDI symmetry class which lacks any
topological structure for two-dimensional systems. However, as shown in Ref.
thakurathi14, , the two-dimensional semi-infinite Kitaev ribbon can be mapped
to an aggregate of one-dimensional Kitaev chains for each quasimomentum along
the infinitely long direction (see also Ref. chitra21, ) which might host
topological end Majoranas for some particular quasimomenta. It is these
Majorana modes which appear at the edges of the Kitaev Honeycomb lattice.

###  II.4 Haldane model of graphene

The Haldane model of graphene is a paradigmatic system which encompasses the
basic physics of all Chern insulators and have tremendous application in the
design of topologically protected quantum states. The bare Hamiltonian for the
Haldane model haldane83  is obtained by breaking the time-reversal and
sublattice symmetries of graphene,

Figure 6: (a) Hexagonal graphene lattice showing the nearest-neighbor (N1) and
next-nearest-neighbor (N2) hopping vectors $\vec{\Delta}_{1i}$ and
$\vec{\Delta}_{2i}$, respectively, where the lattice constant is set to be
$a=1$. The hollow and the filled atoms represent the B and A sublattices
respectively. (b) The Brillouin zone of graphene spanned by the reciprocal
lattice vectors $\vec{b}_{1}$ and $\vec{b}_{2}$ containing two inequivalent
Dirac points $K$ and $K^{\prime}$ (the cartesian directions has been labeled
by $k_{x}$ and $k_{y}$ respectively). The color density shows the absolute
value of the band gap $\Delta_{k}$ of the reciprocal space graphene
Hamiltonian showing vanishing gaps at the Dirac points for a $600\times 600$
lattice size having the N1 hopping strength $t_{1}=1.0$ and the N2 hopping
$t_{2}=0$.

$H_{\alpha,\beta,n,m}=-t_{1}\sum\limits_{\left<m\alpha,n\beta\right>}a_{m,\alpha}^{\dagger}a_{n,\beta}+M\sum\limits_{n}a_{n,A}^{\dagger}a_{n,A}-M\sum\limits_{n}a_{n,B}^{\dagger}a_{n,B}-t_{2}{\sum\limits_{\left<\left<m\alpha,n\alpha\right>\right>}[e^{i\phi_{m,n}}a_{m,\alpha}^{\dagger}a_{n,\alpha}}+{\rm
H.c.}],$ (42)

where the real nearest-neighbor (N1) hopping $t_{1}$ (with $t_{2}=0,M=0$)
comprises the bare graphene Hamiltonian; the indices $n$ and $\alpha$
represent site and sublattice respectively. The diagonal staggered mass
(Semenoff mass) $M$ explicitly breaks the sublattice symmetry of the model.
Further, the complex phase $\phi_{m,n}$ is chosen such that within each
hexagonal plaquette, $\phi_{m,n}=\phi$ if the hopping is in the clockwise
sense while $\phi_{m,n}=-\phi$ if the corresponding hopping is in the anti-
clockwise sense. The complex phase thus breaks the time-reversal symmetry in
the next nearest-neighbor hopping. However, the net flux through each
hexagonal plaquette remains zero. The Haldane model does not respect any of
the $\mathcal{S}$, $\mathcal{P}$ or $\mathcal{T}$ symmetries. Depending on the
parameters $M$, $t_{1}$, $t_{2}$ and $\phi$, it is known to exhibit non-
trivial Chern topology in its ground state when the lower energy band is
completely filled. This is somewhat different from the topology manifested in
the SSH and Kitaev models considered earlier, where the topology of those
models were protected by discrete symmetries. However, the Haldane Chern
insulator explicitly breaks all symmetries yet being a two-dimensional system,
can still host a non-trivial topology.

For a half-filled Haldane system, the topology of the Hamiltonian with
periodic boundary conditions is precisely characterized by a gauge invariant
Chern topological invariant (see Eq. (289) of Sec. A.3),

$C=\frac{1}{2\pi}\oint_{BZ}dk_{1}dk_{2}\mathcal{F}_{12}(\ket{\psi_{k}}),$ (43)

where $\mathcal{F}_{12}(\ket{\psi_{k}})$ is the $U(1)$ curvature defined over
the ground state $\ket{\psi_{k}}$ of the Hamiltonian $H^{k}$, i.e.,

$\begin{split}\mathcal{F}_{12}(\ket{\psi_{k}})=i\partial_{k_{2}}\langle{\psi_{k}|\partial_{k_{1}}|\psi_{k}}\rangle-i\partial_{k_{1}}\langle{\psi_{k}|\partial_{k_{2}}|\psi_{k}\rangle}.\end{split}$
(44)

The Chern invariant is integer quantized as long as the system does not
approach a QCP where the Chern number becomes ill-defined. Different integer
values of the Chern number characterize distinct topological phases separated
by QCPs (see Fig. 7).

Each point on the underlying Bravais lattice of the graphene honeycomb lattice
can be referenced in terms of the vector,

$\vec{r}=n_{1}\vec{a}_{1}+n_{2}\vec{a}_{2},$ (45)

where $\vec{a}_{1}$ and $\vec{a}_{2}$ are the primitive vectors and
$n_{1},n_{2}$ are integers. We choose the vectors $\vec{a}_{1}$ and
$\vec{a}_{2}$ to be the next-nearest-neighbor (N2) hopping vectors such that

$\displaystyle\vec{a}_{1}$ $\displaystyle=$ $\displaystyle\vec{\Delta}_{22},$
$\displaystyle\vec{a}_{2}$ $\displaystyle=$ $\displaystyle-\vec{\Delta}_{21},$
(46)

where $\vec{\Delta}_{2i}$ are the $N2$ vectors as shown in Fig. 6. Their
values are found to be,

$\begin{split}\vec{\Delta}_{11}=\frac{a}{2}\\{1,\sqrt{3}\\},\leavevmode\nobreak\
\leavevmode\nobreak\ \vec{\Delta}_{12}=\\{-a,0\\},\leavevmode\nobreak\
\leavevmode\nobreak\ \vec{\Delta}_{13}=\frac{a}{2}\\{1,-\sqrt{3}\\},\\\
\vec{\Delta}_{21}=\frac{a}{2}\\{-3,\sqrt{3}\\},\leavevmode\nobreak\
\leavevmode\nobreak\
\vec{\Delta}_{22}=\frac{a}{2}\\{3,\sqrt{3}\\},\leavevmode\nobreak\
\leavevmode\nobreak\ \vec{\Delta}_{23}=\\{0,-a\sqrt{3}\\},\end{split}$ (47)

where $a$ is the distance between neighboring sites.

Invoking the discrete translational invariance of the Hamiltonian, one can
employ a discrete Fourier transform to decouple the Hamiltonian $H$ in
momentum space. The reciprocal space is spanned by the vector
$\vec{b}=k_{1}\vec{b}_{1}+k_{2}\vec{b}_{2}$, where the primitive vectors are
given by

$\vec{b}_{1}=\frac{2\pi}{3a}\\{1,\sqrt{3}\\}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\vec{b}_{2}=\frac{2\pi}{3a}\\{1,-\sqrt{3}\\},$ (48)

where $k_{1}$ and $k_{2}$ assume integer values. We choose a rhomboidal BZ
spanned by reciprocal lattice vectors $\vec{b}_{1}$ and $\vec{b}_{2}$ (see
Fig. 6). It contains two independent Dirac points $K$ and $K^{\prime}$ see
(Fig. 6),

$K=\frac{2\pi}{3}\left(1,\frac{1}{\sqrt{3}}\right)\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
K^{\prime}=\frac{2\pi}{3}\left(1,-\frac{1}{\sqrt{3}}\right),$ (49)

where we have set $a=1$.

The bare Haldane Hamiltonian gets decoupled in the momentum space where $H(k)$
can be written in the basis $\ket{k,A}$ and $\ket{k,B}$ as,

$H(k)=\vec{h}(k)\cdot\vec{\sigma}=h_{x}(k)\sigma^{x}+h_{y}(k)\sigma^{y}+h_{z}(k)\sigma^{z},$
(50)

such that

$\begin{split}h_{x}(k)=-t_{1}\sum\limits_{i=1}^{3}\cos{\left(\vec{k}\cdot\vec{\Delta}_{1i}\right)},\\\
h_{y}(k)=-t_{1}\sum\limits_{i=1}^{3}\sin{\left(\vec{k}\cdot\vec{\Delta}_{1i}\right)},\\\
{h_{z}(k)=M-2t_{2}\sin{\phi}\sum\limits_{i=1}^{3}\sin{\left(\vec{k}\cdot\vec{\Delta}_{2i}\right)},}\end{split}$
(51)

Similar to the one-dimensional models, the Chern insulator also exhibits a
topological bulk-boundary correspondence – the topological non-triviality of
the bulk system with periodic boundary conditions is reflected as dissipation-
less conducting states on the boundary surface under open boundary conditions.
This is seen from the Kubo formula which explicitly connects the bulk Chern
number with the transverse Hall conductivity in clean systems. Using linear
response theory, it can be rigorously established that at sufficiently low
temperatures, the transverse Hall conductivity is given byhajdu87 ;
vasilopoulos11 ,

$\sigma_{xy}=-\frac{e^{2}}{2\pi\hbar}C,$ (52)

where $e$ is the electronic charge, $\hbar$ is the Planck’s constant and $C$
is the integer-quantized topological Chern number defined in Eq. (43). This
connection lies at the root of the highly stable integer quantized anomalous
quantum Hall effect manifested in such systems without the presence of any
external magnetic fields. We have thus found a robust physical manifestation
of the complete geometric and topological structure of quantum states which
can be experimentally probed. Indeed, the anomalous quantum Hall effect and
the quantum Hall effect are some of the benchmark experiments in physics till
date which show near perfect quantization of the transverse Hall conductivity
thus firmly establishing the role of topology in quantum many-body systems.
Interestingly, the mere semblance between the electromagnetic theory and the
geometry of quantum states that we discussed extensively earlier has given
rise to a completely new field called Berry electrodynamics where the action
of the gauge field $F_{\mu\nu}$ on moving charges is studied. In simple terms,
it concerns the motion of charges on the curved parameter-manifold and hence,
one must account for its curvature which itself appears as the gauge field,
similar to the "magnetic field" in electrodynamics.

Figure 7: Topological phase diagram of the Haldane model of graphene where the
inequivalent phases (with different Chern numbers) are separated by gapless
critical lines (solid lines in the figure).

We have reached the end of the first part of our review whose primary goal was
to introduce the rich concept of topology in quantum many-body systems. We
next move on to the second part in which we will be focusing on the
topological physics of dynamical quantum many-body systems and their
experimental realization.

##  III Periodically driven topological systems

It has been seen that irradiated quantum many-body systems develop novel out-
of equilibrium properties which are absent in the equilibrium system. This has
initiated a plethora of studies in periodically driven quantum many-body
systems to explore exciting new phases of matter oka08 ; kitagawa11 ;
lindner11 ; moessner13 . For, example, it has been seen that bare graphene
which does not have a well-defined topological classification in equilibrium,
may exhibit non-trivial topological properties when irradiated with
lightkitagawa11 . In this section, we will discuss two such interesting
systems, i.e., periodically driven Kitaev models both as a Kitaev
chainthakurathi13 ; smitha13  and the Kitaev honeycomb latticethakurathi14 .

###  III.1 Periodically driven Kitaev chain

Let us first consider the case of a periodically driven Kitaev chain
thakurathi13 . In particular, we explore the simple situation where the
chemical potential $\mu$ in the Kitaev chain Hamiltonian, given in Eq. (13),
is periodically kicked in time,

$\mu(t)\leavevmode\nobreak\ =\leavevmode\nobreak\ c_{0}\leavevmode\nobreak\
+\leavevmode\nobreak\ c_{1}\sum_{n=-\infty}^{\infty}\delta(t-nT).$ (53)

Here, $T=2\pi/\omega$ is the time interval between subsequent kicks. In the
Majorana representation, the Hamiltonian can be represented in the general
form

$H(t)\leavevmode\nobreak\ =\leavevmode\nobreak\
i\sum_{m,n=1}^{2N}\leavevmode\nobreak\ a_{m}M_{mn}(t)a_{n},$ (54)

where $M(t)$ is a real antisymmetric matrix satisfying $M(t+T)=M(t)$. Noting
that the Heisenberg operators $a_{n}(t)$ satisfy,

$\frac{da_{m}(t)}{dt}\leavevmode\nobreak\ =\leavevmode\nobreak\
4\sum_{n=1}^{2N}\leavevmode\nobreak\ M_{mn}(t)\leavevmode\nobreak\ a_{n}(t),$
(55)

the time-evolution operator is found to be,

$\displaystyle U(t_{2},t_{1})$ $\displaystyle=$ $\displaystyle{\cal
T}e^{4\int_{t_{1}}^{t_{2}}dtM(t)},$ (56)

where $\cal T$ denotes the time-ordering symbol.

We will restrict ourselves to observing the system only at stroboscopic
instants, separated by time interval $T$. Floquet theory dictates that the
evolution of a quantum system at such instants is generated by the Floquet
operator $U(T,0)$ floquet_sens_47 ; shirley65 ; bukov_ap_139 ;
eckardt_rmp_011004 . It is important to note that the eigenvalues of $U(T,0)$
are given by phases, $e^{i\theta_{j}}$, and they come in complex conjugate
pairs if $e^{i\theta_{j}}\neq 1$. This is because
$U(T,0)\psi_{j}=e^{i\theta_{j}}\psi_{j}$ implies that
$U(T,0)\psi_{j}^{*}=e^{-i\theta_{j}}\psi_{j}^{*}$. For eigenvalues
$e^{i\theta_{j}}=\pm 1$ (these eigenvalues may, in principle, appear with no
degeneracy), the eigenvectors can be chosen to be real; this can be seen by
noting that $U(T,0)\psi_{j}=\pm\psi_{j}$ also implies
$U(T,0)\psi_{j}^{*}=\pm\psi_{j}^{*}$ (as $U(T,0)$ is a real orthogonal matrix)
and one can choose the eigenvectors to be the real combinations
$\psi_{j}+\psi^{*}_{j}$ and $i(\psi_{j}-\psi^{*}_{j})$.

For the periodic $\delta$-function kick in Eq. (53), the Floquet operator
assumes the form

$U(T,0)\leavevmode\nobreak\ =\leavevmode\nobreak\
e^{4M_{1}}\leavevmode\nobreak\ e^{4M_{0}T},$ (57)

where $M_{0/1}$ are $(2N)$-dimensional antisymmetric matrices whose non-zero
matrix elements can be found using Eqs. (15), (53) and (54):

$\displaystyle(M_{0})_{2n+1,2n}$ $\displaystyle=$
$\displaystyle-\leavevmode\nobreak\ (M_{0})_{2n,2n+1}\leavevmode\nobreak\
=\leavevmode\nobreak\ -\frac{1}{4}(\gamma-\Delta),$
$\displaystyle(M_{0})_{2n-1,2n+2}$ $\displaystyle=$
$\displaystyle-\leavevmode\nobreak\ (M_{0})_{2n+2,2n-1}\leavevmode\nobreak\
=\leavevmode\nobreak\ -\frac{1}{4}(\gamma+\Delta),$
$\displaystyle(M_{0})_{2n-1,2n}$ $\displaystyle=$
$\displaystyle-\leavevmode\nobreak\ (M_{0})_{2n,2n-1}\leavevmode\nobreak\
=\leavevmode\nobreak\ \frac{c_{0}}{2},$ $\displaystyle(M_{1})_{2n-1,2n}$
$\displaystyle=$ $\displaystyle-\leavevmode\nobreak\
(M_{1})_{2n,2n-1}\leavevmode\nobreak\ =\leavevmode\nobreak\ \frac{c_{1}}{2},$
(58)

for an appropriate range of values of $n$. A more convenient form is to use
the symmetrized expression

$U(T,0)\leavevmode\nobreak\ =\leavevmode\nobreak\
e^{2M_{1}}\leavevmode\nobreak\ e^{4M_{0}T}\leavevmode\nobreak\ e^{2M_{1}}.$
(59)

It is easy to show that the Floquet operators in Eqs. (57) and (59) have the
same eigenvalues, while their eigenvectors are related by a unitary
transformation. We will see below that the symmetrized form in Eq. (59) leads
to some simplifications when we derive an effective Hamiltonian and a
topological invariant.

The inverse participation ratio (IPR) serves as a useful probe to find the
eigenvectors of $U(T,0)$ which are localized at the ends. Assuming that the
eigenvectors are normalized so that $\sum_{m=1}^{2N}|\psi_{j}(m)|^{2}=1$ for
each value of $j$, where $m=1,2,\cdots,2N$ labels the components of the
eigenvector, the IPR of an eigenvector is defined as
$I_{j}=\sum_{m=1}^{2N}|\psi_{j}(m)|^{4}$. If $\psi_{j}$ is extended in real
space, so that $|\psi_{j}(m)|^{2}\approx 1/(2N)$ for each $m$, then
$I_{j}=1/(2N)$; which vanishes as $N\to\infty$. On the contrary, if $\psi_{j}$
is localized over a distance $\xi$ (which is of the order of the decay length
of the eigenvector and remains constant as $N\to\infty$), then
$|\psi_{j}(m)|^{2}\sim 1/\xi$ in a region of length $\xi$ and $\sim 0$
elsewhere; the IPR thus remains finite as $N\to\infty$. For sufficiently large
$N$, one can hence identify localized states, if any. A subsequent analysis of
the plot of the probabilities $|\psi_{j}(m)|^{2}$ versus $m$ reveals whether
the localized state is indeed an end state.

Figure 8: (a) IPRs of different eigenvectors of the Floquet operator for a
200-site system with a periodic $\delta$-function kick with $\gamma=1$,
$\Delta=-1$, $c_{0}=2.5$, $c_{1}=0.2$ and $\omega=12$. The two eigenvectors
with the largest IPRs both have an IPR equal to $0.142$ and Floquet eigenvalue
equal to $-1$. (b) The corresponding Majorana end states of the two
eigenvectors with the largest IPRs. (c) Plot of the number of end states
versus $\omega$ for a 200-site system, with Floquet eigenvalues $+1$ ($N_{+}$,
blue squares) and $-1$ ($N_{-}$, red stars). Values of other relevant
parameters are same as that in (a) (See Ref. thakurathi13, ).

Fig. 8 shows the IPRs of the different eigenvectors for a finite system with
$N=200$ with parameters chosen as $\gamma=1$, $\Delta=-1$, $c_{0}=2.5$,
$c_{1}=0.2$ and $\omega=12$. Two of the IPRs clearly stand out with a value of
$0.142$ each. The corresponding eigenvectors are found to have Floquet
eigenvalues $e^{i\theta}=-1$; importantly, the value of $\theta=\pi$ is
separated by a gap of $0.148$ from the values of $\theta$ for all the other
eigenvalues. These eigenvectors are localized at the two ends of the system
and are real; the corresponding probabilities are shown in Fig. 8. The state
at the left end has non-zero $a_{m}$ only if $m$ is even, while the state at
the right end has non-zero $a_{m}$ only for $m$ odd. Thus, the periodic
driving produces Majorana end modes even though for the parameter values given
above, the ground state of the instantaneous Hamiltonian is at all times in
topologically trivial phase. Remarkably, the number of Majorana end modes is
crucially dependent on the driving frequency $\omega$. This is illustrated in
Fig. 8, where the number of eigenvalues lying near $+1$ and $-1$, denoted by
$N_{+}$ and $N_{-}$, respectively, are plotted as a function of $\omega$ for a
finite system of size $N=200$ with $\gamma=1$, $\Delta=-1$, $c_{0}=2.5$ and
$c_{1}=0.2$. Although the number of end modes is not a monotonic function of
$\omega$, the number generally increases as $\omega$ decreases. Recently, it
has also been shown that the number and locality of these Majorana modes can
be controlled if the system is driven by multiple periodic drives having
commensurate frequenciespablo20 .

###  III.2 Topological invariants for periodic driving

As in the equilibrium situation, to define a topological invariant, let us
consider a system with periodic boundary conditions. Once again, Fourier
transforming to the quasimomentum basis decouples each momentum mode; for each
value of $k$, one can define a Floquet operator $U_{k}(T,0)$ which is a
$2\times 2$ unitary matrix. Using Eqs. (22), (53) and (59), it follows that

$U_{k}(T,0)=e^{ic_{1}\tau^{z}}e^{-i2T[(\gamma\cos k-c_{0})\tau^{z}+\Delta\sin
k\tau^{y}]}e^{ic_{1}\tau^{z}},$ (60)

where $k$ is taken to lie in the full range $-\pi\leq k\leq\pi$. A couple of
remarks about the Floquet operator $U_{k}(T,0)$ are in order. First, assuming
that $\Delta\neq 0$, it is easy to show that $U_{k}(T,0)\neq\pm\mathbb{I}$ for
any value of $k\neq 0$ or $\pi$, where $\mathbb{I}$ is the identity operator.
Second, for $k=0$ or $\pi$, $U_{k}(T,0)\neq\pm\mathbb{I}$ unless
$2T(c_{0}\pm\gamma)+2c_{1}=n\pi$, or equivalently unless

$\omega\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{4\pi(c_{0}\pm\gamma)}{n\pi-2c_{1}}$ (61)

for some integer value of $n$. The $\pm$ sign in Eq. (61) corresponds to
$k=\pi$ and $0$ respectively. In other words, if the condition set by Eq. (61)
is not satisfied, the eigenvalues of the Floquet operator are non-degenerate
for all values of $k$ and thus the ‘Floquet eigenspectrum’ is gapped. When Eq.
(61) is satisfied, the spectrum is gapless at modes $k=0$ or $\pi$, which can
be associated with a non-equilibrium phase transition.

Utilizing the properties of the Floquet operator discussed above, one can now
construct a topological invariant as follows. First, let us define an
effective Hamiltonian $h_{\mathrm{eff},k}$ as

$U_{k}(T,0)\leavevmode\nobreak\ =\leavevmode\nobreak\
e^{-ih_{\mathrm{eff},k}}.$ (62)

The symmetric structure of Eq. (60) (as chosen in Eq. (59)) ensures that
$h_{\mathrm{eff},k}$ assumes the form

$h_{\mathrm{eff},k}\leavevmode\nobreak\ =\leavevmode\nobreak\
h^{y}(k)\leavevmode\nobreak\ \tau^{y}\leavevmode\nobreak\
+\leavevmode\nobreak\ h^{z}(k)\leavevmode\nobreak\ \tau^{z}$ (63)

as in Eq. (22). Further, to define $h_{\mathrm{eff},k}$ uniquely, we impose
the condition that the coefficients in Eq. (63) satisfy
$0<\sqrt{h^{y}(k)^{2}+h^{z}(k)^{2}}<\pi$, for $U_{k}(T,0)\neq\pm\mathbb{I}$. A
winding number $W$, analogously to the equilibrium situation, can thus be
computed as described in Eq. (24).

An useful interpretation of the winding number defined above can be made by
noting that, the condition $0<h^{y}(k)^{2}+h^{z}(k)^{2}<\pi$ implies that
$h_{\mathrm{eff},k}$ can be mapped to a point on the surface of a sphere whose
polar angles $(\alpha,\beta)$ are given by
$\alpha=\sqrt{h^{y}(k)^{2}+h^{z}(k)^{2}}$ and
$\beta=\tan^{-1}(h^{z}(k)/h^{y}(k))$. As $k$ goes from 0 to $2\pi$, we obtain
a closed curve which does not pass through the north and south poles. The
integer $W$ can then be related to the winding number of this curve around
either the north pole or the south pole. Note that the winding numbers around
the north and south pole are given by the same integer.

Remarkably, the winding number exhibits a bulk-boundary correspondence, akin
to its equilibrium counterpart. This correspondence can be seen in Fig. 9
which compares the number of Majorana modes at each end of a chain and the
winding number as a function of $\omega$. Note that in the limit
$\omega\to\infty$, i.e., $T\to 0$, Eq. (60) becomes independent of $k$, and we
therefore find a single point in the $(h^{y}(k),h^{z}(k))$ plane. This
corresponds to a curve with zero winding number; thus there is a maximum value
of $\omega$ beyond which no Majorana end modes exist. On the other hand, for
very small $\omega<0.2$, the number of Majorana end modes appear to increase
drastically and diverge in the limit $\omega\to 0$.

It turns out that it is possible to define an alternate topological invariant
that not only exhibits a direct correspondence with the Majorana end modes but
can also distinguish between end modes with eigenvalues $+1$ and $-1$. To
elucidate, let us first define a set of parameters,

$\displaystyle b_{0}$ $\displaystyle=$
$\displaystyle\frac{4(c_{0}-\gamma)}{\omega}\leavevmode\nobreak\
+\leavevmode\nobreak\ \frac{2c_{1}}{\pi},$ $\displaystyle b_{\pi}$
$\displaystyle=$
$\displaystyle\frac{4(c_{0}+\gamma)}{\omega}\leavevmode\nobreak\
+\leavevmode\nobreak\ \frac{2c_{1}}{\pi},$ (64)

so that the Floquet operators at $k=0,\pi$ can be recast as
$U_{0}(T,0)=e^{i\pi b_{0}\tau^{z}}$ and $U_{\pi}(T,0)=e^{i\pi
b_{\pi}\tau^{z}}$. Next, we define a finite line segment, $L_{\omega}$, which
goes from $b_{0}$ to $b_{\pi}$ in one dimension; this line segment is further
assumed to coincide with the $z$-axis.

For $\omega\to\infty$, the line $L_{\omega}$ collapses to a single point given
by $z=2c_{1}/\pi$. As $\omega$ is decreased, $L_{\omega}$ moves and is also
increased in length. With the system parameters chosen as $\gamma=1$,
$c_{0}=2.5$ and $c_{1}=0.2$, the right end of $L_{\omega}$, given by $b_{\pi}$
in Eq. (64), crosses the point $z=n$ with $n=1$ when $\omega\simeq 16.04$. at
this point, the Floquet eigenvalue at $k=\pi$ thus assumes the value
$e^{in\pi}=-1$ and is also degenerate. On further decreasing $\omega$, the
point $z=1$ enters the segment $L_{\omega}$. Importantly, as can be seen from
Fig. 9, a Majorana end mode first appears within the range $16\leq\omega\leq
17$ and also has a Floquet eigenvalue equal to $-1$. As $\omega$ is decreased
further, the right end of $L_{\omega}$ given by $b_{\pi}$ crosses the point
$z=n$ with $n=2$ when $\omega\simeq 7.48$; Eq. (60) then shows that the
Floquet eigenvalue at $k=\pi$ is equal to $e^{in\pi}=1$. Once again, it can be
seen from Fig. 9 that a Majorana end mode appears in the range
$7\leq\omega\leq 8$ with a Floquet eigenvalue equal to $1$. Similarly, as
$\omega$ is decreased even further, the left end of $L_{\omega}$, given by
$b_{0}$ in Eq. (64), crosses the point $z=n$ with $n=1$ at when $\omega=6.88$,
where the Floquet eigenvalue at $k=0$ equals $e^{in\pi}=-1$. Note that as
$\omega$ is decreased a little more, $L_{\omega}$ no longer includes the point
$z=1$. This is reflected in Fig. 9 in the form of a disappearing Majorana end
mode with Floquet eigenvalue equal to $-1$, within the range $6\leq\omega\leq
7$.

In simpler terms, for any value of $\omega$, the number of points $z=n$ (where
$n$ is an integer) which lie inside the line segment $L_{\omega}$ is equal to
the number of Majorana modes at each end of a chain. Further, the numbers of
points with $n$ odd and even will give the numbers of end modes with Floquet
eigenvalue equal to $-1$ and $1$ respectively. It is therefore clear that the
numbers of odd and even integers lying inside $L_{\omega}$ form a set of
topological invariants which can change only at values of $\omega$ where
either $b_{0}$ or $b_{\pi}$ acquires an integer value. When that happens, the
Floquet spectrum at either $k=0$ or $\pi$ becomes gapless, thus bearing the
hallmark of a topological phase transition.

Figure 9: Comparison of the number of Majorana modes at each end of a 200-site
system (black solid, $y$-axis on left) and the winding number (magenta dashed,
$y$-axis on right) as a function of $\omega$ from 1 to 18, for $\gamma=1$,
$\Delta=-1$, $c_{0}=2.5$ and a periodic $\delta$-function kick with
$c_{1}=0.2$. The inset shows a range of $\omega$ from $0.2$ to 1 where there
is a large number of Majorana modes (See Ref. thakurathi13, ).

The above results can be summarized as follows. Assuming that $b_{0/\pi}$ are
not integers, we consider all the integers lying between $b_{0}$ and
$b_{\pi}$. Of these, let $n_{e}^{>}$ ($n_{o}^{>}$) and $n_{e}^{<}$
($n_{o}^{<}$) respectively denote the numbers of even (odd) integers which are
greater than and less than $2c_{1}/\pi$. Then the numbers $N_{+}$ and $N_{-}$
of modes at each end of a chain with Floquet eigenvalues $+1$ and $-1$ are
given by $N_{+}\leavevmode\nobreak\ =\leavevmode\nobreak\
|n_{e}^{>}-n_{e}^{<}|$ and $N_{-}\leavevmode\nobreak\ =\leavevmode\nobreak\
|n_{o}^{>}-n_{o}^{<}|$. Also, note that the winding number $W$, given by
$|W|=|n_{e}^{>}-n_{e}^{<}+n_{o}^{>}-n_{o}^{<}|$, is generally not equal to the
total number of modes, $N_{+}+N_{-}$, at each end of a chain (although
$|W|-(N_{+}+N_{-})$ is always an even integer). In Table I, we list the values
of $N_{+}$, $N_{-}$ and $|W|$ versus $\omega$ for a 200-site system with
$\gamma=1$, $\Delta=-1$, $c_{0}=0.5$ and $c_{1}=0.2$.

$\omega$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$N_{+}$ | 2 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0
$N_{-}$ | 2 | 2 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
$|W|$ | 4 | 2 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0
Table 1: Values of $N_{+}$, $N_{-}$ and $|W|$ versus $\omega$ for $\gamma=1$,
$\Delta=-1$, $c_{0}=0.5$ and $c_{1}=0.2$. $|W|\neq N_{+}+N_{-}$ for
$4\leq\omega\leq 6$.

Finally, in the limit $\omega\to 0$, it is clear from Eq. (64) that the number
of Majorana end modes diverges asymptotically as $8|\gamma|/\omega$ if
$c_{0}\pm\gamma$ have the same sign and as $8|c_{0}|/\omega$ if
$c_{0}\pm\gamma$ have opposite signs.

Before concluding this subsection, we would also like to note that several
other variants of topological invariants have been proposed over the years
which capture the topological features of periodically driven systems kundu13
; asboth14 ; yao17 ; vega18 . Finally, we would like to mention that the
topological properties of a (integrable) Kitaev chain with long range
interactions has been extensively studied (see Ref. maity_review,  and
references therein). Similarly, the emergent topological features of the long
range chain with periodic drives has also been investigated maity_periodic .

###  III.3 Periodically driven Kitaev honeycomb model

Finally, we consider a periodically driven Kitaev honeycomb model thakurathi14
, the equilibrium situation of which was discussed in Sec. II.3. As with the
periodically driven Kitaev chain, we analyze the situation in which one of the
coupling parameters is repeatedly kicked after equal time intervals. For
simplicity, consider the case in which the driven parameter is $J_{3}$, so
that,

$J_{3}(t)\leavevmode\nobreak\ =\leavevmode\nobreak\ J_{0}\leavevmode\nobreak\
+\leavevmode\nobreak\ J_{p}\leavevmode\nobreak\
\sum_{n=-\infty}^{\infty}\leavevmode\nobreak\ \delta(t-nT),$ (65)

where the $T=2\pi/\omega$ is the time-period.

As in the case of the Kitaev chain, the Floquet operator $U(T,0)$ governs the
evolution of the system at stroboscopic times, and this can be evaluated by
recasting Eq. (37) as follows. Given a system with $N=N_{x}N_{y}$ sites, we
introduce a column vector
$\hat{c}=(\hat{a}_{1},\hat{a}_{2},\cdots,\hat{a}_{N},\hat{b}_{1},\hat{b}_{2},\cdots,\hat{b}_{N})^{T}$.
Eqs. (37) can then be written as $d\hat{c}(t)/dt=4M(t)\hat{c}(t)$, where

$M\leavevmode\nobreak\ =\leavevmode\nobreak\
\left(\begin{array}[]{cc}0&-L^{T}\\\ L&0\end{array}\right).$ (66)

The periodicity of $M(t)$ in time implies that the evolution of $\hat{c}$ over
an interval between subsequent kicks is given by

$\displaystyle\hat{c}(T)$ $\displaystyle=$ $\displaystyle
U(T,0)\leavevmode\nobreak\ \hat{c}(0),$ $\displaystyle{\rm
where}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ U(T,0)$
$\displaystyle=$ $\displaystyle{\cal T}e^{4\int_{0}^{T}dtM(t)},$ (67)

and $\cal T$ denotes the time-ordering symbol. For the case of a periodically
kicked system as in Eq. (65), we obtain, $U(T,0)=e^{4M_{1}}e^{4M_{0}T}$, where
$e^{4M_{0}T}$ is the operator which time evolves from $t=0$ to $t=T$, and
$e^{4M_{1}}$ then evolves across the $\delta$-function at $t=T$.

Figure 10: Zigzag edge states for the Kitaev honeycomb model with $N_{x}\times
N_{y}=27\times 14$, $J_{1}=0.7,\leavevmode\nobreak\
J_{2}=0.15,\leavevmode\nobreak\ J_{0}=0.15,\leavevmode\nobreak\ J_{p}=0.2$ and
$\omega=3$ (See Ref. thakurathi14, ).

Proceeding as in the case of the Kitaev chain, one can now look for the
presence of zero energy Majorana edge modes in the Floquet operator.
Remarkably, edge modes are found to appear on both zigzag and armchair edges.
(An example of a Floquet zigzag edge mode is shown in Fig. 10). This is in
contrast to the time-independent version of the model in which all such modes
appear for the armchair edges.

To explore the emergence of the Floquet Majorana modes in the context of a
bulk-boundary correspondence, let is consider an infinite system, so that the
momentum modes are decoupled. The Floquet operator for each mode then assumes
the form

$\displaystyle U_{\vec{k}}(T,0)$ $\displaystyle=$ $\displaystyle
e^{-i2J_{p}\tau^{y}}\leavevmode\nobreak\
e^{-iT(X_{\vec{k}}\tau^{z}\leavevmode\nobreak\ +\leavevmode\nobreak\
Y_{\vec{k}}\tau^{y})},$ $\displaystyle X_{\vec{k}}$ $\displaystyle=$
$\displaystyle 2\leavevmode\nobreak\
[J_{1}\sin({\vec{k}}\cdot{\vec{M}}_{1})\leavevmode\nobreak\
-\leavevmode\nobreak\ J_{2}\sin({\vec{k}}\cdot{\vec{M}}_{2})],$ $\displaystyle
Y_{\vec{k}}$ $\displaystyle=$ $\displaystyle 2\leavevmode\nobreak\
[J_{0}\leavevmode\nobreak\ +\leavevmode\nobreak\
J_{1}\cos({\vec{k}}\cdot{\vec{M}}_{1}\leavevmode\nobreak\
+\leavevmode\nobreak\ J_{2}\cos({\vec{k}}\cdot{\vec{M}}_{2})].$

As Majorana modes are associated with the existence of Floquet eigenvalues
equal to $\pm 1$, one can analytically derive the critical drive frequencies
at which Majorana edge modes appear or disappear,

$\omega_{\vec{k}}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{4\pi\leavevmode\nobreak\ [\leavevmode\nobreak\ J_{0}\leavevmode\nobreak\
\pm\leavevmode\nobreak\
\sqrt{J_{1}^{2}+J_{2}^{2}+2J_{1}J_{2}\cos(k_{x})}\leavevmode\nobreak\
]}{n\pi\leavevmode\nobreak\ -\leavevmode\nobreak\ 2J_{p}}.$ (69)

A deeper insight into the appearance of Majorana modes can be obtained by
mapping the two-dimensional Kitaev model to a chain as follows. Consider a
system which has a finite width in the $y$-direction (with zigzag edges along
the top and bottom as indicated in Fig. 4) and is infinitely long in the
$x$-direction. The momentum $k$ along the $x$-axis is a good quantum number.
We now use the Heisenberg equations of motion

$\displaystyle\frac{d\hat{a}_{m}}{dt}$ $\displaystyle=$
$\displaystyle(J_{1}e^{-ik/2}\leavevmode\nobreak\ +\leavevmode\nobreak\
J_{2}e^{ik/2})\hat{b}_{m}\leavevmode\nobreak\ +\leavevmode\nobreak\
J_{3}\hat{b}_{m+1},$ $\displaystyle\frac{d\hat{b}_{m}}{dt}$ $\displaystyle=$
$\displaystyle-(J_{1}e^{ik/2}+J_{2}e^{-ik/2})\hat{a}_{m}-J_{3}\hat{a}_{m-1},$
(70)

for all $m\geq 1$, with the understanding that $\hat{a}_{0}=0$. Next we
rewrite Eq. (15) as,

$\displaystyle H$ $\displaystyle=$ $\displaystyle
i\sum_{n=1}^{\infty}\leavevmode\nobreak\
[-J_{x}\hat{b}_{n+1}\hat{a}_{n}-J_{y}\hat{b}_{n}\hat{a}_{n+1}+\mu\hat{b}_{n}\hat{a}_{n}],$
(71)

which becomes identical to Eq. (15) on substituting $a_{n}\to a_{2n},b_{n}\to
a_{2n-1}$. For $k=0$, it is straightforward to see that the Heisenberg
equations of motion of the operators $\hat{a}_{n}$ and $\hat{b}_{n}$ in Eq.
(71) agree with Eqs. (70) when,

$J_{x}\leavevmode\nobreak\ =\leavevmode\nobreak\ J_{3}/2,\leavevmode\nobreak\
\leavevmode\nobreak\ J_{y}\leavevmode\nobreak\ =\leavevmode\nobreak\
0,\leavevmode\nobreak\ \leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\
\leavevmode\nobreak\ \mu\leavevmode\nobreak\ =\leavevmode\nobreak\
-(J_{1}+J_{2})/2.$ (72)

We now see the situation of $J_{3}$ being periodically kicked corresponds, for
$k=0$, to the situation in which the parameter $J_{x}$ of the chain is kicked
periodically. As discussed in Sec. III.1, the $\delta$-function kicks can
produce Floquet Majorana modes at the ends of the system, which correspond to
the zigzag edges of the two-dimensional model. This idea can also be extended
for $k\neq 0$, where two of the parameters appearing in Eqs. (70) can be
written as,

$\displaystyle J_{1}e^{\pm ik/2}\leavevmode\nobreak\ +\leavevmode\nobreak\
J_{2}e^{\mp ik/2}\leavevmode\nobreak\ =\leavevmode\nobreak\ J_{k}e^{\pm
i\phi_{k}},$ $\displaystyle{\rm where\leavevmode\nobreak\ }$ $\displaystyle
J_{k}\leavevmode\nobreak\ =\leavevmode\nobreak\
\sqrt{J_{1}^{2}\leavevmode\nobreak\ +\leavevmode\nobreak\
J_{2}^{2}\leavevmode\nobreak\ +\leavevmode\nobreak\ 2J_{1}J_{2}\cos k}.$ (73)

The phase $\phi_{k}$ can be removed from Eqs. (70) by a unitary
transformation; this unitary transformation is independent of $J_{3}$ and is
therefore not affected by the periodic kicks in $J_{3}$. The situation thus
becomes similar to that for $k=0$, except that the parameter $\mu$ in Eq. (72)
is now given by $\mu_{k}=-J_{k}/2$. We thus have a family of one-dimensional
problems which are labeled by the parameter $k$ and the appearance of Majorana
modes can be similarly explained for each of them separately.

Figure 11: Blue regions in the $(\omega,k)$ space where Majorana states appear
on the zigzag edges when the parameter $J_{3}$ is given periodic
$\delta$-function kicks. The system being considered has a width of $100$
sites, and $J_{1}=0.7,\leavevmode\nobreak\ J_{2}=0.15,\leavevmode\nobreak\
J_{0}=0.15$ and $J_{p}=0.3$. The empty region in the middle is bounded on the
right and left by two solid red lines which show the analytical results given
in Eqs. (75) and (76) respectively. Two more red solid lines corresponding to
$n=-1$ and 2 are shown. They almost coincide with each other and appear within
the blue regions on the left; they cross near $\omega=2$ and $k\simeq 1.3$
(See Ref. thakurathi14, ).

The mapping to a chain discussed above allows us to precisely identify the
conditions for the appearance or disappearance of Majorana modes. According to
Eqs. (69) and (73), Floquet edge modes with a given momentum $k$ and FE equal
to $(-1)^{n}$ should appear or disappear when

$\omega_{k}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{4\pi\leavevmode\nobreak\ [J_{0}\leavevmode\nobreak\
\pm\leavevmode\nobreak\ J_{k}]}{n\pi\leavevmode\nobreak\ -\leavevmode\nobreak\
2J_{p}},$ (74)

where $n$ is an integer. As the kicking frequency $\omega$ is decreased, Eq.
(74) gives the red solid line on the right side of the empty region in Fig. 11
where a Floquet edge mode disappears with $n=0$, namely,

$\omega_{k}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{4\pi\leavevmode\nobreak\ [J_{k}\leavevmode\nobreak\
-\leavevmode\nobreak\ J_{0}]}{2J_{p}},$ (75)

and the red solid line on the left side of the empty region in Fig. 11 where a
Floquet edge mode appears with $n=1$, namely,

$\omega_{k}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{4\pi\leavevmode\nobreak\ [J_{0}\leavevmode\nobreak\
+\leavevmode\nobreak\ J_{k}]}{\pi\leavevmode\nobreak\ -\leavevmode\nobreak\
2J_{p}}.$ (76)

In general, for $n\leq 0$, we have a line given by

$\omega_{k}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{4\pi\leavevmode\nobreak\ [J_{k}\leavevmode\nobreak\
-\leavevmode\nobreak\ J_{0}]}{2J_{p}-n\pi},$ (77)

while for $n\geq 1$, we have a line given by

$\omega_{k}\leavevmode\nobreak\ =\leavevmode\nobreak\
\frac{4\pi\leavevmode\nobreak\ [J_{0}\leavevmode\nobreak\
+\leavevmode\nobreak\ J_{k}]}{n\pi\leavevmode\nobreak\ -\leavevmode\nobreak\
2J_{p}},$ (78)

where Majorana modes appear or disappear. Fig. 11 also show the red solid
lines for $n=-1$ and 2. These appear within the blue regions; they cross near
$\omega=2$ and $k\simeq 1.3$ where we see a small gap indicating that there
are no Majorana modes in that region. When $\omega$ is decreased below these
two lines, the Majorana mode with FE equal to $-1$ disappears and a mode with
FE equal to $+1$ appears. When $\omega$ is decreased even further, more modes
start appearing which correspond to $n>2$ and $n<-1$.

Before concluding this section, we remark here that the existence of localized
eigenstates of the Floquet operator does not guarantee that the time-evolved
state of the system would also exhibit localized Majorana edge modes. This is
in contrast with equilibrium systems where the topological features are
manifested in the ground state of the system. In the periodically driven
scenarios discussed above, the time-evolved state, in general, has finite
overlap with a large number of Floquet eigenstates. Hence, the localized
Majorana modes can only manifest if one can ensure that the time-evolved state
maintains a high fidelity with the ground state of the effective Floquet
Hamiltonian throughout the evolution. This difficulty is a part of the broader
challenge of dynamically preparing topological states which we will discuss in
the next section. However, the signature of these non-equilibrium Majorana
modes is expected to manifest in stroboscopic heat transport through the
system and the spin density in the $XY$ spin-$1/2$ chain (see Ref. paolo17, ).

##  IV Generically driven one-dimensional topological systems

In the previous section, we saw that in periodically driven systems, the
Floquet Hamiltonian generating the stroboscopic dynamics may be topologically
non-trivial. In fact, there exist countless scenarios when in a driven system
the effective Hamiltonian or simply a quenched Hamiltonian generating the
dynamics becomes topologically non-trivial despite the triviality of the
equilibrium system. However, to observe such dynamical topological phases, one
must figure out a way to dynamically prepare the non-equilibrium
thermodynamically large system in the ground state of the non-trivial
effective Hamiltonian. As discussed in Sec. II, the thermodynamic limit is
required for the perfect exponential localization of the edge states. Such an
endeavor is not simple, as in a generic situation the out-of-equilibrium
dynamics is bound to generate excitations except when the dynamics can be
considered adiabatic. At the same time, it turns out that such adiabatic
transformations connecting inequivalent topological phases is not possible in
thermodynamically large systems as the adiabaticity necessarily breaks down at
the critical point separating the two phases (see Sec. A.1 for details).
Clearly, the thermodynamic limit and the adiabaticity seems too much to ask
for while attempting to dynamically prepare topologically non-trivial quantum
states from trivial ones. In this section, we will focus on the progress made
in the attempt to address this particular bottleneck. We will also see that it
is possible to topologically classify out-of-equilibrium quantum states
without referring to the Hamiltonian generating the dynamics.

###  IV.1 Explicit symmetry breaking: unitary preparation of topological
states

In the previous sections we have seen that the inequivalent SPT phases cannot
be connected adiabatically in thermodynamically large systems as they are
separated by a gapless quantum critical point. However, it is the protecting
symmetry (such as $S$ in the case of the SSH model) which allows for the
existence of distinct inequivalent topological phases in such models. This
suggests that explicitly breaking the protecting symmetry might let us connect
inequivalent topological phases adiabatically even in thermodynamically large
systems. To visualize this protocol, we first look into the dynamical change
of the winding number when the system is subjected to an external time-
dependent driving.

We consider the SSH model and study the temporal evolutionutso19  of the
equilibrium topological invariant, i.e., the winding number under a generic
unitary drive. We begin with an initial state $|\psi_{k}(0)\rangle$, the
system is allowed to evolve under the driven Hamiltonian $H_{k}(t)$. The state
$|\psi_{k}(0)\rangle$ therefore evolves with time as

$\displaystyle|\psi_{k}(t)\rangle$ $\displaystyle=$
$\displaystyle\mathbb{T}e^{-i\int_{0}^{t}H_{k}(t^{{}^{\prime}})dt^{{}^{\prime}}}|\psi_{k}(0)\rangle$
(79) $\displaystyle\equiv$ $\displaystyle e^{-iH_{k}^{\rm
eff}(t)t}|\psi_{k}(0)\rangle$ $\displaystyle=$ $\displaystyle
U_{k}(t)|\psi_{k}(0)\rangle,$

where $H^{\rm eff}_{k}(t)$ is the time-dependent effective Hamiltonian acting
as a generator of the unitary evolution acting on the driven system and
$\mathbb{T}$ denotes the time ordering operator. We now investigate the fate
of the winding number under such a time-dependent dynamics. To analyze this,
let us recall the time-dependent or dynamical connection as

$A_{k}(t)\equiv\left[\langle\psi_{k}(0)|U_{k}^{\dagger}\right]\partial_{k}\left[U_{k}|\psi_{k}(0)\rangle\right],$
(80)

which evolves in time as,

$\displaystyle A_{k}(t)$ $\displaystyle=$
$\displaystyle\langle\psi_{k}(0)|\partial_{k}|\psi_{k}(0)\rangle+\langle\psi_{k}(0)|U_{k}^{\dagger}(\partial_{k}U_{k})|\psi_{k}(0)\rangle$
(81) $\displaystyle=$ $\displaystyle
A_{k}(0)+\langle\psi_{k}(0)|U_{k}^{\dagger}(\partial_{k}U_{k})|\psi_{k}(0)\rangle.$

Hence, the change in the connection at a later time is given by

$\Delta
A_{k}=A_{k}(t)-A_{k}(0)=\langle\psi_{k}(0)|U_{k}^{\dagger}\partial_{k}U_{k}|\psi_{k}(0)\rangle.$
(82)

Recasting the effective Hamiltonian to the following form, $H_{k}^{\rm
eff}(t)=|m(k,t)|\left(\hat{m}(k,t).\vec{\sigma}\right)$ and also denoting
$|m(k,t)|$ simply as $m$, we obtain

$\displaystyle U_{k}^{\dagger}\partial_{k}U_{k}$
$\displaystyle=\partial_{k}m\leavevmode\nobreak\
\sin^{2}(mt)\leavevmode\nobreak\
\\{-it\hat{m}\cdot{\vec{\sigma}}+i\hat{m}\times\partial_{k}\hat{m}\cdot{\vec{\sigma}}\\}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-i\leavevmode\nobreak\ \sin(mt)\cos(mt)\leavevmode\nobreak\
\partial_{k}\hat{m}\cdot{\vec{\sigma}}.$ (83)

The initial state $|\psi_{k}(0)\rangle$ that we consider happens to be the
ground state of the SSH Hamiltonian (belonging BDI class) which can be chosen
to be of the form of Eq. (7) where $\phi(k)$ is an odd function of $k$.
Interestingly, the terms on the right hand side of the Eq. (IV.1) can be shown
to vanish individually when integrated over the entire BZ, pertaining to
certain conditions imposed upon the effective Hamiltonian $H^{\rm eff}_{k}$ as
discussed below.

Let us now analyze the implications of Eqs. (82) and (IV.1). Taking the
expectation value of the first term of the above equation with respect to the
state $|\psi_{k}(0)\rangle$, one observes that the integral of this quantity
over the full BZ vanishes identically if $m_{x}(k)$ is an even function of $k$
and $m_{y}(k)$ is an odd function of $k$. Similarly, analyzing the integral of
the next two terms over the full BZ, we see that both of them vanish
identically if $m_{z}(k)$ is an odd function of $k$ or zero in addition to the
above constraints imposed on $m_{x}(k)$ and $m_{y}(k)$. If the above
conditions are satisfied by the effective Hamiltonian then the winding number
must remain invariant in time.

It is evident from Eq. (11) that the above constraints on the single particle
Hamiltonian in $k$ space, demand the presence of certain symmetries of the
effective Hamiltonian. Namely, one concludes that the equilibrium winding
number remains invariant under temporal evolution if the effective dynamical
Hamiltonian $(H^{\rm eff}_{k})$ respects either of the symmetry combinations,
$\mathcal{T}$ and $\mathcal{P}$ ($m_{x}(k)\rightarrow$ even,
$m_{y}(k)\rightarrow$ odd, $m_{z}(k)\rightarrow 0$) simultaneously or just
$\mathcal{P}$ ($m_{x}(k)\rightarrow$ even, $m_{y}(k)\rightarrow$ odd,
$m_{z}(k)\rightarrow$ odd). The SSH Hamiltonian described in Eq. (1) with real
hopping amplitudes ($v,w\in\mathbb{R}$) lies in the BDI class of topological
Hamiltonians. This suggests that the topological invariant is temporally
invariant under all unitary dynamics unless the BDI symmetries (here
particularly the particle-hole symmetry) are explicitly broken by the
effective Hamiltonian generating the dynamics.

Nonetheless, there appears to be a more serious issue with the above
procedure. As discussed so far, it is clear that in SPT phases, it is the
protecting symmetry which is behind the quantization of the topological
invariant and thus the quantized numbers of experimentally verifiable edge
states. In fact, if the symmetry is explicitly disrespected in the dynamics,
even if one is able to dynamically tune the value of the topological
invariant, it will no longer remain quantized and the topological
classification is rendered useless. In this regard one needs to connect the
changing winding number with an observable physical phenomena in the non-
equilibrium system. This calls for a physical significance of the topological
winding number in one-dimensional systems which turns out to be intricately
connected with the macroscopic dipole electric polarization of the chain. The
macroscopic electric polarization of the system is defined as,

$P=e\braket{\hat{X}},$ (84)

$e$ being the electronic charge which we set equal to unity and $\hat{X}$
being the position operator summed over all sites of the chain, i.e.,
$\hat{X}=\sum\limits_{i}\hat{x}_{i}$. However, since the momentum space and
real space operators are related by a Fourier transform, we can replace the
position operatorresta93 ; bardyn18  as

$\hat{X}\rightarrow\frac{i}{2\pi}\int dk\partial_{k},$ (85)

which suggests that the winding number is nothing but the macroscopic electric
polarization of the system. The equivalence in Eq. (85) may not be immediately
apparent, as with periodic boundary conditions the operator $\hat{X}$ is not
manifestly periodic and one needs to compactify the defining Eq. (84) before
proceeding. We will discuss the derivation in detail in a later section. Now,
furthering this correspondence, it is straightforward to show that in an
arbitrary time-dependent situation the bulk polarization current density
$j(t)$ of the SSH chain is directly proportional to the rate of change of the
topological winding number $(\nu)$ (see Refs. utso19,  and cooper18, ):

$j(t)=\frac{dP(t)}{dt}=\frac{d\nu}{dt}=\frac{1}{2\pi}\int_{BZ}dk\bra{\psi_{k}(t)}\partial_{k}H_{k}(t)\ket{\psi_{k}(t)},$
(86)

where $\ket{\psi(t)}$ is the time-evolved state for each quasimomenta mode $k$
and $H_{k}(t)$ is the instantaneous time-dependent Hamiltonian. Particularly,
in the case of a time periodic drive with a period $T$, the stroboscopic
(measured after a complete period), the variation of the winding number
denoted as $\Delta\nu_{m}$ for the $m-$th stroboscopic interval, is related to
the average change in the bulk polarization density of the chain within the
$(m-1)-$th and the $m-$th period of evolution i.e.,

$\Delta\nu_{m}=\frac{\nu(mT)-\nu((m-1)T)}{T}=\frac{1}{T}\int_{(m-1)T}^{mT}dtj(t),$
(87)

which is nothing but the net charge transferred through the chain in one cycle
of the drive. Moreover, since it is the Floquet Hamiltonian $H_{F}$ which
generates the stroboscopic time evolution, the above discovered symmetry
constraints must be imposed on $H_{F}$ to keep the topological index invariant
in stroboscopic time. However, it might be noted that the Floquet Hamiltonian
is responsible only for the stroboscopic evolution. The micromotion within
consecutive stroboscopic intervals is governed by the most general time-
dependent Hamiltonian itself. It turns out, that there might arise scenarios
when the instantaneous Hamiltonian explicitly breaks the protecting symmetries
while the Floquet Hamiltonian still respects it. In those situations although
the stroboscopic winding number remains invariant, it might vary in
micromotion.

####  IV.1.1 Preparation of Floquet topological state

To demonstrate the applicability of explicit symmetry breaking in the
dynamics, we will demonstrate a protocol which adiabatically prepares a
Floquet topological state from a trivial one. Starting from the trivial ground
state of a BDI symmetric SSH Hamiltonian, we periodically drive the system.
Choosing a sufficiently high driving frequency $(\omega=2\pi/T)$, the periodic
perturbation is switched on adiabatically, allowing the instantaneous state of
the system to follow the instantaneous ground state of the Floquet
Hamiltonian. We choose the protocoldutta19

$\displaystyle H_{k}(t)$ $\displaystyle=$ $\displaystyle H_{0}+\lambda_{k}(t)$
(88) $\displaystyle=$ $\displaystyle
H_{0}(k)+A_{1}(t)V_{0}\sigma^{x}+A_{2}(t)V_{k}(t),$

with $V_{k}(t+T)=V_{k}(t)$ being periodic in time having a time period of $T$.
The amplitudes $A_{1}(t)$ and $A_{2}(t)$ are adiabatic ramps which slowly
switch on and off the periodic perturbation,

$\displaystyle A_{1}(t)$ $\displaystyle=$ $\displaystyle
A_{2}(t)=-t/\tau\text{\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ for\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ }0\leq t\leq\tau,$ $\displaystyle A_{1}(t)$
$\displaystyle=$ $\displaystyle 1\text{\leavevmode\nobreak\
\leavevmode\nobreak\ and\leavevmode\nobreak\ \leavevmode\nobreak\
}A_{2}(t)=\frac{t}{\tau}-2\text{\leavevmode\nobreak\ \leavevmode\nobreak\
for\leavevmode\nobreak\ \leavevmode\nobreak\ }\tau\leq t\leq 2\tau,$
$\displaystyle A_{1}(t)$ $\displaystyle=$ $\displaystyle
1\text{\leavevmode\nobreak\ \leavevmode\nobreak\ and\leavevmode\nobreak\
\leavevmode\nobreak\ }A_{2}(t)=0\text{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ for\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ }t\geq 2\tau,$ (89)

allowing the Floquet Hamiltonian to get adiabatically modified. The periodic
perturbation $V_{k}(t)$ is chosen to be sinusoidal,

$V_{k}(t)=V_{1}\sin(\omega t)\sigma^{x}+V_{1}\cos(\omega t)\sigma^{y}.$ (90)

We note that in the interval $0\leq t\leq\tau$, the periodic perturbation
$V_{k}(t)$ is switched on adiabatically and in the subsequent interval
$\tau\leq t\leq 2\tau$ a part of the periodic drive is switched off, again
adiabatically, to restore all the symmetries.
Evidently, due to the presence of the linear ramping, the protocol is not
perfectly periodic in the time interval $0\leq t\leq 2\tau$. However, as the
ramping is much slower than the time period of the periodic perturbation, the
complete evolution may be partitioned into the ‘fast’ and the ‘slow’ time
variables kitagawa11 : the former is the time scale of the high frequency
periodic drive and the later represents the time scale of the linearly ramped
amplitude. Therefore, for a fixed value of the amplitude, the drive is
completely periodic and a Floquet picture holds.
If the frequency of the periodic drive is much higher than the characteristic
bandwidth of the system and the amplitude of the drive is small enough
$\left(\mathcal{A}^{2}/\omega\ll 1,\mathcal{A}\leavevmode\nobreak\ \text{being
the {\it effective} amplitude of the drive}\right)$, the drive is off
resonant. The effective Floquet Hamiltonian may be approximated just to
incorporate the single-photon virtual processes or nearest-neighbor hopping in
the Floquet Bloch lattice i.e.,

$H_{F}(k)\simeq\frac{1}{T}\int_{0}^{T}H_{k}(t)dt+\frac{1}{\omega}\left[H_{-1}(k),H_{+1}(k)\right]+\mathcal{O}\left(\frac{\mathcal{A}^{4}}{\omega^{2}}\right),$
(91)

where $H_{\pm 1}(k)=\frac{1}{T}\int_{0}^{T}H(t)e^{\pm i\omega t}dt.$ We now
further assume that $A(t)$ changes very slowly $(\tau\gg T)$ such that its
change over a time period $T$ may be neglected at very high frequencies. Using
Eq. (91), the effective Floquet Hamiltonian takes the following form,

$H_{F}^{t}(k)\simeq
H_{0}(k)-A_{1}(t)V_{0}\sigma^{x}+\frac{\left[A_{2}(t)V_{1}\right]^{2}}{2\omega}\sigma^{z},$
(92)

The slowly varying amplitude can now be interpreted as an adiabatic
deformation of the Floquet Hamiltonian $H_{F}^{t}(k)$. The magnitude of
$V_{0}$ is so chosen that the term proportional to $\sigma^{x}$ in the
effective Floquet Hamiltonian Eq. (92) drives the Floquet Hamiltonian across
the QCP to render the final Floquet Hamiltonian topologically non-trivial at
time $t=2\tau$.

The diagonal term on the other hand serves two purposes. It generates a
staggered mass term which opens up a gap between the Floquet eigenstates and
keeps the floquet spectrum gapped at all times. This gap in the Floquet
spectrum protects the instantaneous ground states against the generation of
excitations particularly near the QCP where the gap otherwise vanishes in the
absence of the staggered mass. The mass term also explicitly breaks the
$\mathcal{P}$ and $\mathcal{S}$ symmetry at the same time in the Floquet
Hamiltonian $H_{F}(k)$ causing the bulk topological invariant to change
stroboscopically. At the end of the switching off (i.e., at $t=2\tau$), the
particle hole symmetry is completely restored in the final Floquet Hamiltonian
$H_{F}^{t}(k)$ which now belongs to a completely symmetric BDI SSH model in
the topologically non-trivial sector.

Figure 12: Temporal variation of the winding number as a function of
stroboscopic time period for a periodic boundary condition. The winding number
assumes a constant value (unity) after the $\mathcal{P}$ symmetry have been
completely restored in the Floquet Hamiltonian at $t=2\tau$, with
$\tau=10^{4}T$. (See Ref. dutta19, ).  Figure 13: Overlap of the evolved state
$\ket{\psi_{k}(2\tau)}$ at $t=2\tau$ with the Floquet eigenstate at the same
time $\ket{\psi_{0}^{F}(k)}$ for all $k\in\left[-\pi,\pi\right]$. For a
sufficiently large $\tau$ ($\tau=10^{4}T$), the instantaneous state of the
system has a significant overlap with the Floquet topological ground state in
the high frequency limit. The initial parameters are $v=1.55$, $w=1.50$ and
$L=1000$; $\omega=100.0$ and amplitudes $V_{0}=0.1$, $V_{1}=7.0$. (See Ref.
dutta19, ).

The winding number calculated from the non-equilibrium state of the driven
chain is also observed to become topologically non-trivial (see Fig. 12). In
response to the explicit breaking of the $\mathcal{P}$ symmetry and the
consequent variation of the winding number, we observe a stroboscopic
generation of polarization current in the bulk which vanishes just after the
BDI symmetries are restored in the Floquet Hamiltonian at $t=2\tau$. Since the
$\mathcal{S}$ symmetry has been restored finally, the winding number also
remains invariant further and the final Floquet Hamiltonian indeed becomes
topologically inequivalent to the trivial Hamiltonian. The mass term
necessarily avoids the crossing of a gapless QCP and therefore under an
adiabatic protocol, the instantaneous state of the system
$\ket{\psi_{k}(2\tau)}$ is seen to dynamically follow the ground state
$\ket{\psi_{0}^{F}(k)}$ of the high frequency Floquet Hamiltonian.

The adiabatic ramping functions $A(t)$ in Eq. (IV.1.1) have been chosen so
that its rate of change is much slower than the time period of the periodic
perturbation (i.e., $\tau\gg T$). This ensures that the dynamical state of the
system at stroboscopic intervals of time have a significant fidelity to the
topological eigenstate of the Floquet Hamiltonian. In Fig. 13, we see that at
$t=2\tau$, the stroboscopic state $\ket{\psi_{k}(2\tau)}$ has a significant
overlap with the eigenstate $\ket{\psi_{0}^{F}(k)}$ of the Floquet Hamiltonian
for all $k\in\left[-\pi,\pi\right]$.

####  IV.1.2 Stroboscopic emergence of edge states

In this section we study the time periodic driving protocol given in Eq.
(IV.1.1) in a SSH chain with open boundary conditions. The static SSH
Hamiltonian $H(v_{0},w)$ as in Eq. (1), with $\mathcal{N}$ unit cells under
open boundary conditions (i.e., $c^{\dagger}_{2\mathcal{N}+1}=0$), is a
$2\mathcal{N}\times 2\mathcal{N}$ Hermitian matrix in the single-particle
basis $\mathcal{B}\equiv\\{\ket{i}\\}_{i=1}^{2\mathcal{N}}$. Here $\ket{i}$
denotes the $i^{th}$ site in the chain. We start with an eigenstate
$\ket{\psi(0)}$ of the initial topologically trivial Hamiltonian $H(v_{0},w)$
having $\mathcal{N}=50$ unit cells and adiabatically switch on an unitary time
dependent periodic perturbation,

$H(t)=H(v_{0},w)+v(t)\left(\sum_{i=1}^{2\mathcal{N}}c^{\dagger}_{n,1}c_{n,2}+{\rm
H.c.}\right),$ (93)

where

$\displaystyle v(t)$ $\displaystyle=$
$\displaystyle-\frac{t}{\tau}\left[V_{0}+V_{1}e^{-i\omega
t}\right]\text{\leavevmode\nobreak\ \leavevmode\nobreak\
for\leavevmode\nobreak\ \leavevmode\nobreak\ }0\leq t\leq\tau,$ $\displaystyle
v(t)$ $\displaystyle=$ $\displaystyle-
V_{0}+\left(\frac{t}{\tau}-2\right)V_{1}e^{-i\omega
t}\text{\leavevmode\nobreak\ \leavevmode\nobreak\ for\leavevmode\nobreak\
\leavevmode\nobreak\ }\tau\leq t\leq 2\tau,$ $\displaystyle v(t)$
$\displaystyle=$ $\displaystyle-V_{0}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{\leavevmode\nobreak\
\leavevmode\nobreak\ for\leavevmode\nobreak\ }t\geq 2\tau,$ (94)

and $\omega=2\pi/T$ is the frequency of the drive while $T$ being it’s time
period. The matrix elements of the dynamical unitary propagator $U(t)$ in the
time-independent basis $\mathcal{B}$ under the protocol Eq. (IV.1.2),
satisfies a linear system of $2\mathcal{N}\times 2\mathcal{N}$ differential
equations,

$i\frac{dU(t)}{dt}=H(t)U(t).$ (96)

To find the stroboscopic propagator $U(nT,(n-1)T)$, which dynamically evolves
the state of the system from the time $t=(n-1)T$ to $t=nT$, we iteratively
solve Eq. (96) with appropriate boundary conditions in each interval
$\left[(n-1)T,nT\right]$. By a repetitive application of the unitary
propagator $U(nT,(n-1)T)$ for each stroboscopic interval, we obtain the time-
evolved state of the system at $t=NT=2\tau$,

$\ket{\psi(2\tau)}=\prod_{n=1}^{N=\frac{2\tau}{T}}U(nT,(n-1)T)\ket{\psi(0)}.$
(97)

Figure 14: Under open boundary conditions, at $t=2\tau$, the stroboscopic
time-evolved state shows localization at the ends of the chain after the BDI
symmetries have been restored in the topological Floquet Hamiltonian. The
energy of the localized state $E=0.003$ for a finite chain with 100 lattice
sites. The initial state (same as in (a)) is subjected to a unitary drive as
in Eq. IV.1.1 with $V_{0}=0.1$, $V_{1}=2.0$, $\tau=500T$ and $\omega=100.0$,
where $i$ denotes the lattice site. (See Ref. dutta19, ).

To maintain adiabaticity of the deformation of the Floquet Hamiltonian which
is essential for the adiabatic preparation, we choose $\tau\gg T$
($\tau=500T$) and a sufficiently high frequency $\omega=100$ for the periodic
perturbation. The initial parameters of the topologically trivial SSH chain
$H(v_{0},w)$ were chosen to be $v_{0}=1.55$ and $w=1.50$. With the choice of
the amplitudes $V_{0}=0.1$ and $V_{1}=2.0$ as specified in the (see Fig. 14),
we study the support of the stroboscopic state $\ket{\psi(2\tau)}$ of the
system at the $i^{th}$ site of the chain through the quantity $|\langle
i\ket{\psi(2\tau)}|^{2}$. We observe localization of the state
$\ket{\psi(2\tau)}$ to the two ends of the finite chain corresponding to the
topologically non-trivial winding number of the bulk Floquet Hamiltonian
depicted in Fig. 12. We note that although the generated stroboscopic end
modes de-localize considerably into the bulk of the chain, in the
thermodynamic limit the hybridization between the end modes will eventually
vanish and they will get strictly localized at the edge.

The stroboscopic state thus prepared at $t\geq 2\tau$, is observed to host
topologically protected zero energy modes localized at the ends of a finite
chain with open boundary conditions (see Fig. 14). As the parameters of the
final Floquet Hamiltonian are close to the critical values, the localization
length of the zero energy edge states are large. The end states therefore
hybridize near the midpoint of the chain which consequently lifts the energy
of the localized states from exactly zero; this is also true for an undriven
SSH chain for the same set of parameters. However, for a sufficiently long
chain, the edge states formed would be highly localized and will be perfectly
localized in the thermodynamic limit. It can also be shown that higher winding
number Floquet states can be prepared adiabatically using similar dynamical
protocols starting from trivial states.

###  IV.2 Dynamical symmetry breaking: A symmetry classification of non-
equilibrium systems

In the previous section we discussed the applications of explicitly breaking
the protective symmetries of a SPT insulator in dynamically preparing non-
trivial states. There is however a different kind of symmetry breaking in such
systems which arise purely out of the non-equilibrium dynamics while the
Hamiltonian generating the dynamics still respects all the symmetries of the
initial systemcooper18 . To be clear, this will be addressed as dynamical
symmetry breaking in contrast to the case when the symmetries are explicitly
broken in the Hamiltonian itself. To precisely discuss these scenarios one
must therefore find a way to develop a symmetry classification of the non-
equilibrium states themselves irrespective of the symmetries of the
Hamiltonian that generates the dynamics. To achieve this, we must first look
into the symmetry classification of equilibrium states and how they are
inherited from the symmetries of the Hamiltonian itself. As demonstrated in
Sec. II, all non-interacting quantum systems can be classified with respect to
two anti-unitary symmetries, i.e., charge conjugation/particle-hole
$\mathcal{P}$, time-reversal $\mathcal{T}$ and an unitary symmetry, i.e.,
chiral/sublattice $\mathcal{S}$. In an equilibrium system, if the ground state
is non-degenerate, it is naturally expected that the ground state will
directly inherit all the symmetries respected by the Hamiltonian,

$\mathcal{H}=\sum\limits_{ij}\psi^{\dagger}_{i}H_{ij}\psi_{j},$ (98)

where $\psi_{j}$ and $\psi_{j}^{\dagger}$ are fermionic annihilation and
creation operators respectively and $H$ is the single-particle Hamiltonian. As
discussed previously, the symmetry operators can be re-written as
$\mathcal{P}=P\mathcal{K}$, $\mathcal{T}=T\mathcal{K}$ and $\mathcal{S}$ where
$P$, $T$ and $S$ are unitary. The symmetry transformations corresponding to is
said to respect the time-reversal, charge conjugation and chiral symmetries
can then be represented to act on the creation and annihilation operators such
that the fermionic anti-commutators are preserved,

$\displaystyle\mathcal{T}\psi_{i}\mathcal{T}^{-1}$ $\displaystyle=$
$\displaystyle T_{ij}\psi_{j}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \mathcal{T}i\mathcal{T}^{-1}=-i,$
$\displaystyle\mathcal{P}\psi_{i}\mathcal{P}^{-1}$ $\displaystyle=$
$\displaystyle P_{ij}\psi_{j}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \mathcal{P}i\mathcal{P}^{-1}=-i,$
$\displaystyle\mathcal{S}$ $\displaystyle=$
$\displaystyle\mathcal{T}\mathcal{P}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \mathcal{S}i\mathcal{S}^{-1}=i.$ (99)

Since the charge conjugation operator transforms a particle into its anti-
particle, it acts on the annihilation operators to produce linear combinations
of creation operators. If the system respect theses discrete symmetries then
one obtains,

$\begin{split}\mathcal{P}\mathcal{H}\mathcal{P}^{-1}=H,\end{split}$ (100)

and equivalently for the other two. Similarly, the single-particle density
matrix $\rho$ can be constructed out of the single-particle fermionic
correlation functions in a many-body state $\ket{\Phi}$ as,

$\rho_{ij}=\braket{\Phi}{\psi_{i}^{\dagger}\psi_{j}}{\Phi}.$ (101)

It is then straight forward to show that if the equilibrium state $\rho$
respects the three discrete symmetries, one must obtain,

$\begin{split}T\rho^{*}T^{\dagger}=\rho,\\\
P\rho^{*}P^{\dagger}=\mathbb{I}-\rho,\\\ S\rho
S^{\dagger}=\mathbb{I}-\rho.\end{split}$ (102)

The density matrix thus evolves in time as,

$\rho(t)=U(t)\rho(0)U^{\dagger}(t),$ (103)

where $U(t)$ is the temporal propagator and is generated by a time-dependent
effective Hamiltonian $H^{eff}(t)$ (see Eq. (79)). Unlike our discussion in
the previous section, we further assume that the symmetries of the initial
Hamiltonian are not explicitly broken in the effective Hamiltonian generating
the dynamics and therefore satisfies the constraints in Eq. (11). It can then
be directly inferred from Eq. (102)-(103) that the dynamical state $\rho(t)$
must satisfy the following constraints under the same symmetry operations,

$\begin{split}T\rho^{*}(t)T^{\dagger}=\rho(-t),\\\
P\rho^{*}(t)P^{\dagger}=\mathbb{I}-\rho(t),\\\
S\rho^{*}(t)S^{\dagger}=\mathbb{I}-\rho(-t).\end{split}$ (104)

However, in generic out-of-equilibrium systems the equivalence of $\rho(t)$
and $\rho(-t)$ does not hold, i.e., in most situations $\rho(t)\neq\rho(-t)$.
A comparison with Eq. (102) reveals that the dynamical state no longer
respects the time-reversal and the chiral symmetry even while the effective
time-dependent Hamiltonian stays symmetric. This shows that although the
charge conjugation symmetry is preserved, the chiral and time-reversal and
symmetries are generically broken in the dynamical state without explicitly
breaking the symmetries of the effective Hamiltonian. This phenomena will be
termed as a dynamical symmetry breaking in all the subsequent discussions.
Consequently, if the topological phases are protected by the particle-hole
symmetry, the time-dependent topological invariant should stay invariant as
the dynamics preserves $\mathcal{P}$ even in the out-of equilibrium state.
However, SPT models lying in other symmetry classes not having particle-hole
symmetry will generically experience a dynamical symmetry breaking while out
of equilibrium and hence, the topological invariant need no longer be
quantized. Remarkably, as we have seen previously, a dynamically varying
winding number is also experimentally relevant and manifests as a flow of
polarization current (see Eq. (86)),

$j(t)=\frac{d\nu(t)}{dt}$ (105)

through the bulk system. Thus, even though the effective Hamiltonian respects
the equilibrium symmetries, a dynamical symmetry breaking can indeed be probed
by a bulk polarization current in such systems.

###  IV.3 Non-unitary dynamics of topological states

####  IV.3.1 Lindblad master equation and steady state topology

In the previous sections we have primarily focused on the fate of symmetry
protected topological phases of matter under unitary dynamical and a possible
topological classification of out-of equilibrium quantum states. We now
proceed to discuss the behavior of topologically non-trivial states under non-
unitary dynamics which naturally arise when the system is in interaction with
the environment which is difficult to avoid in thermodynamically large
systems. Dissipative quantum systems, particularly those undergoing a
Markovian evolution always tends to a steady state at asymptotic long times.
In such open quantum systems, it instructive to observe the reduced system of
interest by tracing out the degrees of freedom of the bath. The Markovian
evolution of quantum systems can described by a Master equation generating the
time evolution of the reduced density matrix describing only the system
degrees of freedombreuer ; sourav ,

$\partial_{t}\rho=\mathcal{L}(\rho)=-i[H,\rho]+\mathcal{D}(\rho),$ (106)

where as usual, the Hamiltonian generates the unitary processes within the
system itself and the dissipator functional $\mathcal{D}(.)$ takes care of the
bath induced processes due to the system-bath coupling. The dissipator can be
shown to take the form

$\mathcal{D}(\rho)=\sum\limits_{i}\kappa_{i}\left(2L_{i}\rho
L_{i}^{\dagger}-\\{L_{i}^{\dagger}L_{i},\rho\\}\right),$ (107)

where the operators $L_{i}$ are known as jump operators containing the action
of the bath on the system through the coupling strengths $\kappa_{i}$. For
example, if the bath acts as a particle sink at a site $i$ of a fermionic
system, one might set $L_{i}=c_{i}$ where $c_{i}$ is the annihilation operator
of the system at that site. We will see the effect of such a bath when coupled
to a topological system. Also, a significant property of a Markovian evolution
is that the steady state is completely independent of the initial conditions
when the steady state is non-degenerate. First, we proceed to discuss the
possibility of the dynamical generation of topological steady states through
careful enough designing of the bath. For any choice of te Lindblad jump
operators, the system reaches an asymptotic steady state given by

$\partial_{t}\rho=\mathcal{L}(\rho)=0.$ (108)

It is then possible to engineer the steady state of the system to host non-
trivial properties by constructing the bath accordingly. The exact route and
time taken by such a Markovian system to reach equilibrium has been
extensively studied and is also an interest of current research (see Ref.
vernier20, ).

In fact, it has been shownzoller11  that choosing the Lindblad operators such
that it destroys the Bogoliubov quasiparticles in a $p$-wave Kitaev chain (see
Sec. II.2), one can reach a steady state subspace in which the topological end
Majorana correlations are preserved against dissipation.

####  IV.3.2 Slow quenches in a dissipative environment

Artificially engineered Lindblad operators has been shownzoller11  to cool a
system to specific targeted pure steady states. However, in such cases the
Lindblad operators are highly non-local and do not generically describe real
dissipative scenarios. This is because mostly environmental interactions has
an opposite effect on quantum states, i.e., the destruction of quantum
coherence. For example, the bath $L_{i}=c_{i}$, where $c_{i}$ are the
annihilation operator on site $i$, accounts for the loss of Dirac fermions
from every site of the system. It can be seen that in a purely dissipative
evolution, such a system will asymptotically reach the featureless and boring
fermionic vacuum. In this section we attempt to understand the time and energy
scales involved when such a Kitaev chain coupled to such a system is driven
externally in time across the quantum critical point. Although we saw that a
topological state with a perfect bulk-boundary correspondence does not exist
in finite-size systems, here we will specifically study finite-size systems,
as we know that a finite system is not gapless even at the critical point.
However, since we are solely interested in uncovering the competition between
unitary and dissipative dynamics, we can simply deal with a finite-size system
which can always be driven adiabatically.

Figure 15: (a) Defects generated at the end of the ramp $\chi(\tau)$ as a
function of $\tau$ for a local bath with Lindblad operators of form Eq. (111).
For $\tau<L\lll\tau_{B}$, the dynamics is dominated primarily by the unitary
evolution thereby resulting in lowering of defects with increasing $\tau$. For
$L\ll\tau<\tau_{B}$, the ramp is adiabatic and induces negligible defects on
its own. The defects generated is then primarily through dissipative effects
resulting in more defect generation with increasing $\tau$. The interplay of
the unitary and dissipative dynamics of the bath results in an optimal
$\tau_{o}$ for which the defect generated is minimum. (b) For
$L\ll\tau\ll\tau_{B}$, the defect scales as $\chi(\tau)\sim\kappa\tau$. The
dashed lines correspond to linear fits of the data points. In both (a) and
(b), the relevant parameters chosen are $L=20$, $J=1$, $\mu_{i}=2$ and
$\mu_{f}=0$. (See Ref. souvik20, ).

The presence of a QCP or vanishing gap in the bulk spectrum of a
thermodynamically large chain, presents a conundrum in the context of
preparing a topological state — starting from a trivial phase of the system,
it is impossible to drive the system into a non-trivial phase through a
unitary dynamics. Specifically, in the problem that we consider heresouvik20 ,
the system is initially in the ground state $\ket{\psi(0)}=\ket{\psi_{i}^{0}}$
of an initial Hamiltonian Kitaev chain $H_{i}:|\mu_{i}|>|J|$, following which
the Hamiltonian is ramped across the QCP at $|\mu|=|J|$ to a final
$H_{f}:|\mu_{f}|<|J|$ using the protocol,

$\displaystyle\mu(t)=\left(\mu_{i}+(\mu_{f}-\mu_{i})\frac{t}{\tau}\right)\Theta(\tau-t)+\mu_{f}\Theta(t-\tau),$
(109)

where $\Theta(x)$ is the Heaviside step function. The protocol therefore
linearly ramps the initial $\mu_{i}$ to a final $\mu_{f}$ during time $\tau$
after which $\mu$ remains frozen at the targeted $\mu_{f}$. In the
thermodynamic limit, the quantum adiabatic theorem breaks down; the system
$\ket{\psi(t)}$ therefore cannot be exclusively prepared in the ground state
$\ket{\psi_{f}^{0}}$ of $H_{f}$. This results in generation of defects in the
correlation of the Majorana end modes which we quantify as

$\displaystyle\chi(t)=\bra{\psi_{f}^{0}}\theta\ket{\psi_{f}^{0}}-\theta(t),$
(110)

where $\theta(t)=\bra{\psi(t)}\theta\ket{\psi(t)}$. We remark that although
the ground state $\ket{\psi_{f}^{0}}$ is doubly degenerate in the
thermodynamic limit, we choose the state $\ket{\psi_{f}^{0}}$ as the
Bogoliubov vacuum corresponding to the final Hamiltonian. The choice of the
ground state however, does not affect the results qualitatively. Note that for
the ideal situation of a perfect unitary adiabatic preparation, the quantity
$\chi(t)$ should vanish.

However, for a finite system of size $L$, the bulk spectrum is not truly
gapless; the gap $\delta$ at the QCP scales as $\delta\sim 1/L$. Consequently,
if the time-scale $\tau$ of the ramp is large enough such that $\tau\gg L$,
the bulk of the chain (closed chain with PBC) can be prepared in the ground
state of the final Hamiltonian. However, unlike the bulk, the splitting
between the edge modes vanishes exponentially with $L$. Therefore, even for a
finite system, the defect $\chi(t)$ truly vanishes only in the limit
$\tau\to\infty$. Nevertheless, it is possible for $\chi(t)$ to come
arbitrarily close to zero for a finite-size system even for finite $\tau\gg
L$. The key point of this section is that perfect Majorana correlations cannot
be generated if the ramp duration is finite with the defects scaling as
$\chi(\tau)=e^{-f(L,\mu_{i})\tau}$, where the coefficient $f(L,\mu_{i})\geq 0$
is non-universal.

Now, consider the same protocol in the presence of dissipation due to a
particle-loss type of bath, which effectively acts as a sink of Dirac fermions
at each site of the chain. We model the dissipative effects through local
Lindblad operators that act locally at each site on the Kitaev chain, barring
the edge sites. Specifically, we choose

$L_{j}=c_{j},\quad j\in\\{2,3...,N-1\\},$ (111)

where $c_{j}$ are the fermionic annihilation operators acting on site $j$ of
the chain. In the presence of such local dissipative channels, the Majorana
edge modes which were initially present in the system, are known to decay
exponentially in time due to the finite overlap of the edge modes with the
bulk for $\mu\neq 0$. However, for $\mu=0$, the edge modes are essentially
disconnected from the bulk and are therefore robust against any dissipation
induced in the bulk. In our protocol, we therefore set $\mu_{f}=0$, so that
the edge modes once created (with some defects) remain localized after the
ramp, i.e., $\chi(t>\tau)=\chi(\tau)$. This allows us to focus solely on the
effects of dissipation on the preparation stage of localized edge modes,
thereby excluding any defects generated after the ramp is complete.

The key observations from the numerical results are two-fold. First, as can be
seen from Fig. 15. The second and more significant result is that there exists
an optimal ramp duration for which the defects generated are minimum, which,
as we will elaborate below, arises due to the competition between the unitary
and the dissipative dynamics.

To comprehend these results, it is instructive to compare the three relevant
length/time scales in the dynamics – $L$, $\tau$ and $\tau_{B}$. Assuming a
weak coupling strength $\kappa$ ($1/\tau_{B}$) amounts to setting $L$, $\tau$
$\ll\tau_{B}$. In the adiabatic limit of the ramping protocol, i.e.,
$L\ll\tau\ll\tau_{B}$, the defects generated scale as
$\chi(\tau)\sim\kappa\tau$ (see Fig. 15). Intuitively, this monotonic rise in
defect can be explained as follows — as the edges interact with the bath
indirectly through the bulk, an increase in ramp duration $\tau$ implies that
the edge modes have proportionally increasing time to decay before the
chemical potential is eventually ramped to $\mu_{f}=0$ at $t=\tau$, following
which the edge modes can no longer decay.

Although being an expectation over localized states, this resonates perfectly
with the scaling of defect density in bulk residual energy under similar
dissipation and driving protocols keck17 . See Appendix B for a perturbative
expansion (with $\kappa\tau\ll 1$ as small parameter) of the solution of the
dynamical equations of motion for the two-point Majorana correlations to show
that any observable, which is a linear function of two-point Majorana
correlations, is indeed expected to follow the same linear scaling
$\sim\kappa\tau$ under the action of linear Lindblad operators.

On the other hand, for very fast quenches $\tau<L\lll\tau_{B}$, the chemical
potential $\mu$ is quickly ramped to zero within a short duration $\tau$;
within this duration the environment fails to induce any substantive decay in
the Majorana edge correlations. The short duration of the quench however,
itself results in the generation of defects. The defects scale with the ramp
duration as $\chi(\tau)\sim e^{-f(L,\mu_{i})}$, where the coefficient
$f(L,\mu_{i})\geq 0$ is a model dependent non-universal function which can not
be derived within any analytical framework. The defect generation in the fast
quench limit is therefore dominated by the unitary dynamics and arises due to
the fast (non-adiabatic) ramping.

It follows that, there exists an optimal ramp duration $\tau_{o}$ at which the
defect generation at the end of the ramp is minimized (see Fig. 15). From the
preceding discussions, one can assume that the defect generated for
$\tau\sim\tau_{o}$, has the form as
$\chi(\tau)=\exp\left(-f(L,\mu_{i})\tau\right)+\kappa\tau$. A generic
expression for $\tau_{o}$ can be derived by minimizing the defect with respect
to $\tau$ as

$\frac{d}{d\tau}\left(\exp\left(-f(L,\mu_{i})\tau\right)+\kappa\tau\right)|_{\tau=\tau_{0}}=0,$
(112)

which gives

$\tau_{o}=-\frac{1}{f(L,\mu_{i})}\log\left(\frac{\kappa}{f(L,\mu_{i})}\right).$
(113)

The positivity of $f(L,\mu_{i})$ along with the condition of weak-coupling
strength $\kappa$ implies that there exist a positive definite $\tau_{o}$ at
which the defects are minimized. The existence of the optimal $\tau_{o}$ is a
consequence of the competition between the unitary dynamics which demands a
large $\tau$ for defect minimization and the dissipative effects (also See
Ref. vicari20, ) which requires short $\tau$ for the same.

We emphasize here that the above expression for the optimal ramp duration is
not universal because the scaling of the defects arising from non-adiabatic
effects is model dependent as reflected in the coefficient $f(L,\mu_{i})$.
This is unlike the universal scaling of the optimal ramp time obtained in the
case of defect generation in residual energy keck17 , where the residual
energy is defined as the excess energy of the time-dependent state over the
instantaneous ground state and is quantitatively obtained by replacing
$\theta$ in Eq. (110) with $H(t)$ so that $H(\tau)=H_{f}$. Consequently, the
residual energy is an extensive (bulk) property of the system. The
universality in the scaling of residual energy stems from the fact that the
contribution to the defects from the non-adiabatic excitations follow a
universal Kibble-Zurek scaling.

####  IV.3.3 Majorana correlations in localized phases

In a recent study, it was shown that the decay of the topological end-Majorana
correlations $\theta(t)$ in static topological wires coupled to the particle-
loss Markovian bath (modeled by Lindblad operators $L_{i}=c_{i}$ as discussed
in the previous section) follow a generic exponential trendheyl15 ,

$\theta(t)\rightarrow e^{-\gamma t}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{as}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ t\rightarrow\infty$ (114)

with a rate $\gamma$ set by the dissipative coupling $\kappa$. This
exponential fallout can however be slowed down to a stretched exponential
decay,

$\theta(t)\rightarrow e^{-\left({\gamma
t}\right)^{\alpha}}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{as}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ t\rightarrow\infty,$ (115)

with an exponent $\alpha$ in many-body localized phases of the system when
subjected to disorder. This has also been verified numerically for strong
disorder strengths such that mostly all eigenstates are highly localized. It
should now be stated that the topological Majorana correlations cannot be
preserved against the action of the bath when considered to act on every site
including the edge Majorana sites even with disorder based localization.
However, if the edge Majorana sites are not coupled to the bath directly as in
Sec. IV.3.1 disorder based localization can indeed induce a transition into a
stretched exponential decay of the end-Majorana correlations. This can be
understood as follows:

As the end sites are not directly coupled to the bath, it can only be affected
by the bath through the interaction of the edge modes with the bulk system
which in turn is directly coupled to the dissipator. Therefore, in a localized
phase under disorder, since all the eigenstates are now localized, the
interaction of the edge Majorana modes with the bulk system is minimized which
in turn reduces the interaction of the edge-Majoranas with the bulk. This
slows down the decay of the Majorana edge correlations and one obtains a
stretched exponential decay instead of a much faster exponential one.

####  IV.3.4 Mixed state topological phase transitions

Although the unitary preparation of a topological state must surpass the
obstacle of passing through a gapless quantum critical point, this problem can
be somewhat avoided if one deals with the dissipative preparation of
topological states as discussed in the previous section. However, very often
in such an endeavor one ends up with mixed states rather than pure topological
states. In this regard, one must identify topologically non-trivial density
matrices along with the topological classification of pure states constructed
so far. There has been many attempts in constructing a topological
classification of mixed state density matrices. In this section we will review
a recent development in mixed state topology of one-dimensional Gaussian
systems which also hosts a physical manifestation through observables in mixed
states. We start with the relation between the macroscopic electric
polarization and the topological invariant described in Eq. (84)-(85). As
mentioned previously, one encounters a difficulty while extending the
correspondence in systems with periodic boundary conditions as the many-body
position operator $\hat{X}$ is not manifestly periodic. To surmount this
obstacle, Resta proposed an alternative form of the relation by using
translation operators instead of the simple $\hat{X}$ asresta93 ; bardyn18 ,

$\displaystyle P$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi}{\rm
Im}\ln{\left<\hat{T}\right>},$ $\displaystyle\hat{T}$ $\displaystyle=$
$\displaystyle e^{i\delta\hat{X}},$ (116)

such that $\delta=2\pi/L$, $L$ being the length of the system and the
expectation is to be taken over a many body state, say $\ket{\Psi}$. Note that
the operator $\hat{T}$ is periodic by construction when one assumes periodic
boundary conditions. It is then straightforward to show using the properties
of determinants that for any slater determinant $\ket{\Psi}$, with periodic
boundary conditions the Eq. (IV.3.4) reduces to,

$P=\frac{1}{2\pi}\sum_{n}\int_{BZ}A^{nn}dk,$ (117)

where the sum is over all the occupied single-particle bands. This is nothing
but the many-body version of the Berry phase (see Eq. (244)) of Appendix A). A
natural generalization to the case of mixed density matrices is now visible
as,

$\nu\rightarrow{\rm Im}\ln{{\rm Tr}\left[\hat{T}\rho\right]},$ (118)

$\rho$ being the density matrix of the system. Given the density matrix $\rho$
is of the Gaussian form

$\rho\propto e^{-\sum\limits_{ij}\hat{\psi}_{i}G_{ij}\hat{\psi}_{j}},$ (119)

where $\hat{\psi}_{i}$ are fermionic operators, it can be shown that with
periodic boundary conditions, the Eq. (118) simply reduces to,

$\begin{split}\nu={\rm Im}\sum\limits_{n}\left(1+\mathcal{O}_{n}\right),\\\
\mathcal{O}_{n}=\prod_{k}e^{-g_{n}(k)}e^{iA_{nn}\delta},\end{split}$ (120)

such that the summation is over every occupied band of the matrix $G$ (also
known as purity bands) having eigenvalues $g_{n}(k)$ (purity eigenvalues) and
the product is taken over all momenta mode $k$. Since the number of momenta
mode are extensive, one can further rewrite the quantities $\mathcal{O}_{n}$
as,

$\mathcal{O}_{n}=e^{-N\sum\limits_{k}{N^{-1}g_{n}(k)}}\prod_{k}e^{iA_{nn}\delta}.$
(121)

It is now clear that in the thermodynamic limit ($N\rightarrow\infty$), for a
non-degenerate purity spectrum only the single band $\ket{g_{0}}$ having the