import sys
sys.path.append('../code')
from init_mooc_nb import *
init_notebook()

from scipy import linalg as la
from functools import reduce

pi_ticks = [(-np.pi, r"$-\pi$"),
            (-np.pi / 2, r"$-\pi/2$"),
            (0, r"$0$"),
            (np.pi / 2, r"$\pi/2$"),
            (np.pi, r"$\pi$")]


def checkerboard(W=None):
    lat = kwant.lattice.general([[2, 0], [1, 1]], [(0, 0), (1, 0)])
    a, b = lat.sublattices
    if W:
        def lead_shape(pos):
            (x, y) = pos
            return (0 <= y < W and 0 <= x < W)
        syst = kwant.Builder(kwant.TranslationalSymmetry((1, 1)))
        syst[a.shape(lead_shape, (0, 0))] = 0
        syst[b.shape(lead_shape, (1, 0))] = 0
    else:
        syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))    
        syst[lat.shape(lambda pos: True, (0, 0))] = 0
    syst[kwant.HoppingKind((0, 0), b, a)] = lambda s1, s2, p: -p.t1
    syst[kwant.HoppingKind((-1, 1), b, a)] = lambda s1, s2, p: -p.t2
    syst[kwant.HoppingKind((1, 0), a, b)] = lambda s1, s2, p: -p.t3
    syst[kwant.HoppingKind((0, 1), a, b)] = lambda s1, s2, p: -p.t4
    return syst


def evolution_operator(hamiltonians, T):
    n = len(hamiltonians)
    exps = [la.expm(-1j * h * T / n) for h in hamiltonians]
    return reduce(np.dot, exps)


def get_h_k(lead, p):
    bands = kwant.physics.Bands(lead, args=[p])
    h, t = bands.ham, bands.hop
    return lambda k: h + t * np.exp(-1j * k) + t.T.conj() * np.exp(1j * k)

Populated the namespace with:
np, matplotlib, kwant, holoviews, init_notebook, SimpleNamespace, pprint_matrix, scientific_number, pretty_fmt_complex, plt, pf, display_html
from code/edx_components:
MoocVideo, MoocDiscussion, MoocCheckboxesAssessment, MoocMultipleChoiceAssessment, MoocSelfAssessment
from code/functions:
spectrum, hamiltonian_array, h_k, pauli
Using kwant 1.3.2 and holoviews 1.7.0-x-gf3a9f4f
Executed on 2019-01-31 at 20:43:59.970421.

Introduction

Today's topic, Floquet topological insulators, is introduced by Mark Rudner from the Niels Bohr Institute at Copenhagen.

MoocVideo("1peVp_IZ7Ts", src_location="11.1-intro")

Periodically driven systems

We will now learn about a new generalization of topology, namely how it applies to the quantum evolution of systems with a time-dependent Hamiltonian. As you may recall, we've already encountered time dependence, back when we considered quantum pumps. However, back then we assumed that the time evolution was very slow, such that the system stayed in the ground state at all times, i.e. it was adiabatic. Can we relax the adiabaticity constraint? Can we find an analog of topology in systems that are driven so fast that energy isn't conserved?

For the same reasons as before, we'll consider periodic driving

$$ H(t + T) = H(t). $$

Once again, this is necessary because otherwise, any system can be continuously deformed into any other, and there is no way to define a gap.

Before we get to topology, let's refresh our knowledge of time-dependent systems.

The Schrödinger equation is:

$$ i\frac{d \psi}{dt} = H(t) \psi. $$

It is a linear equation, so we can write its solution as

$$ \psi(t_2) = U(t_2, t_1) \psi(t_1), $$

where $U$ is a unitary time evolution operator. It solves the same Schrödinger equation as the wave function, and is equal to the identity matrix at the initial time. It is commonly written as

$$ U(t_2, t_1) = \mathcal{T} \exp\,\left[-i\int_{t_1}^{t_2} H(t) dt\right], $$

where $\mathcal{T}$ represents time-ordering (not time-reversal symmetry). The time-ordering is just a short-hand notation for the need to solve the full differential equation, and it is necessary if the Hamiltonians evaluated at different times in the integral do not commute.

The time evolution operator satisfies a very simple multiplication rule:

$$ U(t_3, t_1) = U(t_3, t_2) U(t_2, t_1), $$

which just says that time evolution from $t_1$ to $t_3$ is a product of time evolutions from $t_1$ to $t_2$ and then from $t_2$ to $t_3$. Of course an immediate consequence of this is the equality $U(t_2, t_1)^\dagger = U(t_2, t_1)^{-1} = U(t_1, t_2)$.

Floquet theory

The central object for the study of driven systems is the evolution operator over one period of the driving,

$$ U(t + T, t) \equiv U, $$

which is called the Floquet time evolution operator. It is important because it allows us to identify the wave functions that are the same if an integer number of drive periods passes. These are the stationary states of a driven system, and they are given by the eigenvalues of the Floquet operator:

$$ U \psi = e^{i \alpha} \psi. $$

The stationary states are very similar to the eigenstates of a stationary Hamiltonian, except that they are only stationary if we look at fixed times $t + nT$. That's why the Floquet time evolution operator is also called a stroboscopic time evolution operator.

We can very easily construct a Hermitian matrix from $U$, the Floquet Hamiltonian:

$$ H_\textrm{eff} = i T^{-1} \,\ln U. $$

Its eigenvalues $\varepsilon = \alpha / T$ are called quasi-energies, and they always belong to the interval $-\pi < \alpha \leq \pi$.

If the system is translationally invariant, we can study the effective band structure of $H_\textrm{eff}(\mathbf{k})$, find an energy in which the bulk Hamiltonian has no states, and study the topological properties of such a Hamiltonian: most of the things we already know still apply.

Of course, selecting a single quasi-energy as the Fermi level is arbitrary, since the equilibrium state of driven systems doesn't correspond to a Fermi distribution of filling factors, but at least it seems close enough for us to try to apply topological ideas.

question = ("But wait, we arbitrarily chose the starting point $t$ in time for calculating the "
            "Floquet operator. What if we chose a different one?")

answers = ["The starting time is just an extra parameter of our system, and topology depends on it.",
           "It doesn't matter, the wave function evolution within one period "
           "can be neglected, since we are interested in many periods.",
           "There's only one correct starting point in time.",
           "It doesn't matter since the quasienergies are independent of the starting point."]

explanation = ("Choosing a different starting point applies a unitary transformation "
               "to the Floquet evolution operator, and so it keeps the quasienergies the same.")

MoocMultipleChoiceAssessment(question=question, answers=answers, correct_answer=3, explanation=explanation)

But wait, we arbitrarily chose the starting point $t$ in time for calculating the Floquet operator. What if we chose a different one?

The starting time is just an extra parameter of our system, and topology depends on it.
It doesn't matter, the wave function evolution within one period can be neglected, since we are interested in many periods.
There's only one correct starting point in time.
It doesn't matter since the quasienergies are independent of the starting point.

Show answer

The correct answer is:
It doesn't matter since the quasienergies are independent of the starting point.
Choosing a different starting point applies a unitary transformation to the Floquet evolution operator, and so it keeps the quasienergies the same.

Driven Majorana wire

Let us start by considering something we know very well, namely the superconducting Majorana nanowire model from week 2. This model has three important parameters which determine whether the wire is in the topological Majorana phase or not: the chemical potential $\mu$, the superconducting gap $\Delta$, and the magnetic field $B$. The topological phase with unpaired Majorana modes at zero energy is realized for $B > \sqrt{\mu^2 + \Delta^2}$.

Now, imagine that we can periodically drive some of these parameters. For instance, consider the simple example when

$$ \mu = \left\{ \begin{matrix} \mu_1 \quad \text{for } 0 < t < T/2 \\ \mu_2 \quad \text{for } T/2 < t < T \end{matrix}\right. $$

Then, the integral to find the time evolution operator is easy to evaluate, and we simply have

$$ U = \exp(i T H_2 / 2) \exp(i T H_1 / 2) $$

with $H_1$ and $H_2$ the nanowire Hamiltonians with chemical potential $\mu_1$ and $\mu_2$. A peculiar property of driven systems is that as the period becomes large, the band structure 'folds': if the driving is very weak, and the original Hamiltonian has energy $E$, the Floquet Hamiltonian has a much smaller quasienergy $(E\bmod 2\pi /T)$. This means that even when $H_1$ and $H_2$ correspond to trivial systems, we can still obtain nontrivial topology if we make the period large enough, as you can see for yourself:

%%opts Path {+axiswise}

def nanowire_chain(L=None):
    lat = kwant.lattice.chain()

    def onsite(site, p):
        return (2*p.t - p.mu) * pauli.szs0 + p.B * pauli.s0sz + p.delta * pauli.sxs0

    def hopping(site1, site2, p):
        return -p.t * pauli.szs0 + 0.5 * 1j * p.alpha * pauli.szsx

    if L:
        syst = kwant.Builder()
    else:
        syst = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
        L = 1

    syst[(lat(x) for x in range(L))] = onsite
    syst[kwant.HoppingKind((1,), lat)] = hopping

    return syst


def calculate_finite_spectrum(periods, hamiltonians):
    energies = []
    for T in periods:
        U = evolution_operator(hamiltonians, T)
        phases = np.angle(la.eigvals(U))
        phases = np.sort(np.abs(phases))
        ev = np.sort([(-1)**n * val for n, val in enumerate(phases)])
        energies.append(ev)
    return np.array(energies).real


def calculate_bands(momenta, hamiltonians_k, T):
    energies = []
    for k in momenta:
        hamiltonians = [h_k(k) for h_k in hamiltonians_k]
        U = evolution_operator(hamiltonians, T)
        phases = np.angle(la.eigvals(U))
        phases = np.sort(np.abs(phases))
        ev = np.sort([(-1)**n * val for n, val in enumerate(phases)])
        energies.append(ev)
    return np.array(energies).real


J = 2.0
p1 = SimpleNamespace(t=J/2, mu=-1*J, B=J, delta=2*J, alpha=J)
p2 = SimpleNamespace(t=J/2, mu=-3*J, B=J, delta=2*J, alpha=J)

syst = nanowire_chain(L=20).finalized()
H1 = syst.hamiltonian_submatrix(args=[p1])
H2 = syst.hamiltonian_submatrix(args=[p2])

lead = kwant.wraparound.wraparound(nanowire_chain(L=None)).finalized()
h1_k = lambda kx: lead.hamiltonian_submatrix(args=[p1, kx])
h2_k = lambda kx: lead.hamiltonian_submatrix(args=[p2, kx])

periods = np.linspace(0.2 / J, 1.6 / J, 100)
momenta = np.linspace(-np.pi, np.pi)

energies = calculate_finite_spectrum(periods, [H1, H2])
spectrum = np.array([calculate_bands(momenta, [h1_k, h2_k], T) for T in periods])


def plot(n):
    T = J * periods[n]

    plot_1 = holoviews.Path((J * periods, energies),
                         kdims=[r'Driving period $(JT)$', r'Quasi-energy $(ET)$'],
                         label='Finite system')(plot={'xticks': 5, 'yticks': pi_ticks})

    VLine = holoviews.VLine(T)(style={'color': 'b', 'linestyle': '--'})

    plot_2 = holoviews.Path((momenta, spectrum[n]),
                         kdims=['$k$', '$E_kT$'],
                         label='Floquet bands')(plot={'xticks': pi_ticks,
                                                      'yticks': pi_ticks,
                                                      'aspect': 'equal'})
    return plot_1 * VLine + plot_2

holoviews.HoloMap({n: plot(n) for n in np.arange(0, 100, 10)}, kdims=['n']).collate()

n: