{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graphene and Kane-Mele model\n",
    "\n",
    "We are going to:\n",
    "* Deal with 2D band structures\n",
    "* Use a more general lattice (honeycomb lattice of graphene)\n",
    "* Construct the very first topological insulator\n",
    "* Learn about topological protection in presence of time-reversal symmetry\n",
    "\n",
    "*Parts of this tutorial are based on the [online course on topology in condensed matter](http://topocondmat.org)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We'll have 3D plotting and 2D band structure, so we need a handful of helper functions.\n",
    "\n",
    "%run matplotlib_setup.ipy\n",
    "\n",
    "from types import SimpleNamespace\n",
    "\n",
    "from ipywidgets import interact\n",
    "import matplotlib\n",
    "from matplotlib import pyplot\n",
    "from mpl_toolkits import mplot3d\n",
    "import numpy as np\n",
    "\n",
    "import kwant\n",
    "from wraparound import wraparound\n",
    "\n",
    "\n",
    "def momentum_to_lattice(k):\n",
    "    \"\"\"Transform momentum to the basis of reciprocal lattice vectors.\n",
    "    \n",
    "    See https://en.wikipedia.org/wiki/Reciprocal_lattice#Generalization_of_a_dual_lattice\n",
    "    \"\"\"\n",
    "    B = np.array(graphene.prim_vecs).T\n",
    "    A = B.dot(np.linalg.inv(B.T.dot(B)))\n",
    "    return np.linalg.solve(A, k)\n",
    "\n",
    "\n",
    "def dispersion_2D(syst, args=None, lim=1.5*np.pi, num_points=200):\n",
    "    \"\"\"A simple plot of 2D band structure.\"\"\"\n",
    "    if args is None:\n",
    "        args = []\n",
    "    momenta = np.linspace(-lim, lim, num_points)\n",
    "    energies = []\n",
    "    for kx in momenta:\n",
    "        for ky in momenta:\n",
    "            lattice_k = momentum_to_lattice([kx, ky])\n",
    "            h = syst.hamiltonian_submatrix(args=(list(args) + list(lattice_k)))\n",
    "            energies.append(np.linalg.eigvalsh(h))\n",
    "    \n",
    "    energies = np.array(energies).reshape(num_points, num_points, -1)\n",
    "    emin, emax = np.min(energies), np.max(energies)\n",
    "    kx, ky = np.meshgrid(momenta, momenta)\n",
    "    fig = pyplot.figure()\n",
    "    axes = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "    for band in range(energies.shape[-1]):\n",
    "        axes.plot_surface(kx, ky, energies[:, :, band], cstride=2, rstride=2,\n",
    "                          cmap=matplotlib.cm.RdBu_r, vmin=emin, vmax=emax,\n",
    "                          linewidth=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Graphene\n",
    "\n",
    "Quantum Hall effect creates protected edge states using a strong magnetic field. Another way to create those is to start from a system with Dirac cones, and open gaps in those.\n",
    "\n",
    "There is a real (and a very important) two-dimensional system which has Dirac cones: [graphene](http://en.wikipedia.org/wiki/Graphene). So in this chapter we will take graphene and make it into a topological system with chiral edge states.\n",
    "\n",
    "Graphene is a single layer of carbon atoms arranged in a honeycomb lattice. It is a triangular lattice with two atoms per unit cell, type $A$ and type $B$, represented by red and blue sites in the figure:\n",
    "\n",
    "![](images/graphene_lattice.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "graphene = kwant.lattice.general([[1, 0], [1/2, np.sqrt(3)/2]],  # lattice vectors\n",
    "                                 [[0, 0], [0, 1/np.sqrt(3)]])  # Coordinates of the sites\n",
    "a, b = graphene.sublattices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now create a `Builder` with the translational symmetries of graphene, and calculate the bulk dispersion of graphene."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, the wave function in a unit cell can be written as a vector $(\\Psi_A, \\Psi_B)^T$ of amplitudes on the two sites $A$ and $B$. Taking a simple tight-binding model where electrons can hop between neighboring sites with hopping strength $t$, one obtains the Bloch Hamiltonian:\n",
    "\n",
    "$$\n",
    "H_0(\\mathbf{k})= \\begin{pmatrix} 0 & h(\\mathbf{k}) \\\\ h^\\dagger(\\mathbf{k}) & 0 \\end{pmatrix}\\,,\n",
    "$$\n",
    "\n",
    "with $\\mathbf{k}=(k_x, k_y)$ and\n",
    "\n",
    "$$h(\\mathbf{k}) = t_1\\,\\sum_i\\,\\exp\\,\\left(i\\,\\mathbf{k}\\cdot\\mathbf{a}_i\\right)\\,.$$\n",
    "\n",
    "Here $\\mathbf{a}_i$ are the three vectors in the figure, connecting nearest neighbors of the lattice [we set the lattice spacing to one, so that for instance $\\mathbf{a}_1=(1,0)$].  Introducing a set of Pauli matrices $\\sigma$ which act on the sublattice degree of freedom, we can write the Hamiltonian in a compact form as\n",
    "\n",
    "$$H_0(\\mathbf{k}) = t_1\\,\\sum_i\\,\\left[\\sigma_x\\,\\cos(\\mathbf{k}\\cdot\\mathbf{a}_i)-\\sigma_y \\,\\sin(\\mathbf{k}\\cdot\\mathbf{a}_i)\\right]\\,.$$\n",
    "\n",
    "The energy spectrum $E(\\mathbf{k}) = \\pm \\,\\left|h(\\mathbf{k})\\right|$ gives rise to the famous band structure of graphene, with the two bands touching at the six corners of the Brillouin zone:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bulk_graphene = kwant.Builder(kwant.TranslationalSymmetry(*graphene.prim_vecs))\n",
    "bulk_graphene[graphene.shape((lambda pos: True), (0, 0))] = 0\n",
    "bulk_graphene[graphene.neighbors(1)] = 1\n",
    "\n",
    "dispersion_2D(wraparound(bulk_graphene).finalized())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also create 1D ribbons of graphene.\n",
    "\n",
    "There are two nontrivial directions: armchair and zigzag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zigzag_ribbon = kwant.Builder(kwant.TranslationalSymmetry([1, 0]))\n",
    "zigzag_ribbon[graphene.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = 0\n",
    "zigzag_ribbon[graphene.neighbors(1)] = 1\n",
    "\n",
    "kwant.plotter.bands(zigzag_ribbon.finalized());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Your turn!\n",
    "\n",
    "*Calculate a dispersion of an armchair nanoribbon. You'll need to figure out what is its period.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "armchair_ribbon = kwant.Builder(kwant.TranslationalSymmetry([0, np.sqrt(3)]))\n",
    "armchair_ribbon[graphene.shape((lambda pos: abs(pos[0]) < 9), (0, 0))] = 0\n",
    "armchair_ribbon[graphene.neighbors(1)] = 1\n",
    "\n",
    "kwant.plotter.bands(armchair_ribbon.finalized(), fig_size=(12, 8));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Your turn!\n",
    "\n",
    "Add potentials of opposite sign to the zigzag nanoribbon, and see what happens to the dispersion relation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "zigzag_ribbon = kwant.Builder(kwant.TranslationalSymmetry([1, 0]))\n",
    "zigzag_ribbon[a.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = 0.2\n",
    "zigzag_ribbon[b.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = -0.2\n",
    "zigzag_ribbon[graphene.neighbors(1)] = 1\n",
    "\n",
    "kwant.plotter.bands(zigzag_ribbon.finalized(), fig_size=(12, 8));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now opened a gap, but there are no protected states inside it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Haldane model of anomalous quantum Hall effect\n",
    "\n",
    "The more interesting way to open the gap in graphene dispersion is introduced by Duncan Haldane, [Phys. Rev. Lett. **61**, 2015 (1988)](http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.61.2015)\n",
    "\n",
    "The idea of this model is to break inversion symmetry that protects the Dirac points by adding next-nearest neighbor hoppings\n",
    "\n",
    "![](images/haldane_hoppings.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a))\n",
    "nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b))\n",
    "nnn_hoppings = nnn_hoppings_a + nnn_hoppings_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def nnn_hopping(site1, site2, params):\n",
    "    return 1j * params.t_2\n",
    "\n",
    "def onsite(site, params):\n",
    "    return params.m * (1 if site.family == a else -1)\n",
    "\n",
    "def add_hoppings(syst):\n",
    "    syst[graphene.neighbors(1)] = 1\n",
    "    syst[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings]] = nnn_hopping\n",
    "\n",
    "haldane = kwant.Builder(kwant.TranslationalSymmetry(*graphene.prim_vecs))\n",
    "haldane[graphene.shape((lambda pos: True), (0, 0))] = onsite\n",
    "haldane[graphene.neighbors(1)] = 1\n",
    "haldane[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings]] = nnn_hopping\n",
    "\n",
    "@interact(t_2=(0, .08, .01))\n",
    "def qshe_dispersion(t_2=0, m=.2):\n",
    "    dispersion_2D(wraparound(haldane).finalized(), [SimpleNamespace(t_2=t_2, m=m)], num_points=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we see that the gap closes in one of the Dirac cones, and does not close in the other half. Let's see what this means for the dispersion relation in a ribbon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Your turn!\n",
    "\n",
    "*Plot a dispersion of either nanoribbon, and see what happens to the edge states*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def nnn_hopping(site1, site2, params):\n",
    "    return 1j * params.t_2\n",
    "\n",
    "def onsite(site, params):\n",
    "    return params.m * (1 if site.family == a else -1)\n",
    "\n",
    "def add_hoppings(syst):\n",
    "    syst[graphene.neighbors(1)] = 1\n",
    "    syst[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings]] = nnn_hopping\n",
    "\n",
    "zigzag_haldane = kwant.Builder(kwant.TranslationalSymmetry([1, 0]))\n",
    "zigzag_haldane[graphene.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = onsite\n",
    "zigzag_haldane[graphene.neighbors(1)] = 1\n",
    "zigzag_haldane[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings]] = nnn_hopping\n",
    "\n",
    "@interact(t_2=(0, .08, .01))\n",
    "def qshe_dispersion(t_2=0, m=.2):\n",
    "    kwant.plotter.bands(zigzag_haldane.finalized(), [SimpleNamespace(t_2=t_2, m=m)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum spin Hall effect in Kane-Mele model\n",
    "\n",
    "(Following: C.L. Kane and E.J. Mele, Phys. Rev. Lett. **95**, 226801 (2005))\n",
    "\n",
    "Haldane model breaks time-reversal symmetry and inversion symmetry. Lattice-scale hoppings that break time-reversal symmetry do not appear in non-magnetic materials. We can make the Hamiltonian invariant under inversion and time-reversal by making the next-nearest neighbor hoppings spin-dependent.\n",
    "\n",
    "So if we take those hoppings equal to $i t_2 \\sigma_z$, we get teh "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Pauli matrices \n",
    "s0 = np.identity(2)\n",
    "sx = np.array([[0, 1], [1, 0]])\n",
    "sy = np.array([[0, -1j], [1j, 0]])\n",
    "sz = np.diag([1, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def spin_orbit(site1, site2, params):\n",
    "    return 1j * params.t_2 * sz\n",
    "\n",
    "def onsite(site, params):\n",
    "    return s0 * params.m * (1 if site.family == a else -1)\n",
    "\n",
    "def add_hoppings(syst):\n",
    "    syst[graphene.neighbors(1)] = s0\n",
    "    syst[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings]] = spin_orbit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bulk_kane_mele = kwant.Builder(kwant.TranslationalSymmetry(*graphene.prim_vecs))\n",
    "bulk_kane_mele[graphene.shape((lambda pos: True), (0, 0))] = onsite\n",
    "add_hoppings(bulk_kane_mele)\n",
    "\n",
    "@interact(t_2=(0, .3, .01))\n",
    "def qshe_dispersion(t_2=0, m=.2):\n",
    "    dispersion_2D(wraparound(bulk_kane_mele).finalized(), [SimpleNamespace(t_2=t_2, m=m)], num_points=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zigzag_kane_mele = kwant.Builder(kwant.TranslationalSymmetry([1, 0]))\n",
    "zigzag_kane_mele[graphene.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = onsite\n",
    "add_hoppings(zigzag_kane_mele)\n",
    "\n",
    "@interact(t_2=(0, .12, .01))\n",
    "def qshe_zigzag_dispersion(t_2=0, m=.2):\n",
    "    kwant.plotter.bands(zigzag_kane_mele.finalized(), [SimpleNamespace(t_2=t_2, m=m)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Robustness of quantum spin Hall effect\n",
    "\n",
    "We have an open important question: what protects these new edge states? Is it spin conservation? If we don't break the conservation of $\\sigma_z$, then it is obvious that the edge states don't disappear (we have two copies of Haldane model after all).\n",
    "\n",
    "The most interesting property of the quantum spin Hall effect is that it does not rely on any conservation law. The reason for this is [Kramers degeneracy](https://en.wikipedia.org/wiki/Kramers_theorem), that prevents two states that are time-reversal partners of each other from coupling to each other.\n",
    "\n",
    "As a final test of topological protection, let's add an extra parameter that breaks spin conservation and check what happens to the edge states."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Your turn!\n",
    "\n",
    "Add a small perturbation proportional to $i \\sigma_x$ or $i \\sigma_y$ to some hopping. Check how the dispersion changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "zigzag_kane_mele = kwant.Builder(kwant.TranslationalSymmetry([1, 0]))\n",
    "zigzag_kane_mele[graphene.shape((lambda pos: abs(pos[1]) < 9), (0, 0))] = onsite\n",
    "zigzag_kane_mele[graphene.neighbors(1)] = s0\n",
    "zigzag_kane_mele[kwant.builder.HoppingKind((0, 0), b, a)] = s0 + 0.3j * sx\n",
    "zigzag_kane_mele[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings]] = spin_orbit\n",
    "\n",
    "@interact(t_2=(0, .12, .01))\n",
    "def qshe_zigzag_dispersion(t_2=0, m=.2):\n",
    "    kwant.plotter.bands(zigzag_kane_mele.finalized(), [SimpleNamespace(t_2=t_2, m=m)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Self-study (takes more time)\n",
    "\n",
    "Using the multiterminal conductance calculation, calculate what happens to $\\sigma_{xx}$ and $\\sigma_{xy}$ if you turn magnetic field on. What role does spin conservation play?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}