2.1.Discretizing-Hamiltonians-computing-spectra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discretizing Hamiltonians and computing a spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tight-binding models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kwant is designd to work with what commonly is called a \"tight-binding\" model. In practice, what we mean by that is that the Hamiltonian has a set of *discrete* degrees of freedom, i.e. that it can be written as\n",
    "\n",
    "$$\n",
    "H = \\sum_{ij} \\, t_{ij}\\, \\left|i\\rangle\\langle j\\right|\n",
    "$$\n",
    "\n",
    "The states $\\left|i \\right>$ will typically correspond to sites on a lattice and additional\n",
    "degrees of freedom, such as orbitals, spin, ...\n",
    "\n",
    "Let's look at a few examples!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graphene"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A tight-binding model originates naturally when describing materials in terms of atomic orbitals on a lattice. Electrons can then only hop to neighboring orbitals, forexample in graphene between $p_z$-orbitals:\n",
    "\n",
    "<img src=\"images/graphene.svg\" style=\"width: 50%\">\n",
    "\n",
    "The value of the hopping parameters must then be taken from experiment/DFT calculations ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finite differences for continuum Hamiltonians"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Discretizing a 1D Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quite often, the starting point of a study is not a tight-binding Hamiltonian, but a continuum Hamiltonian such as\n",
    "\n",
    "$$H = -\\frac{\\hbar^2}{2 m^*} \\frac{d^2}{dx^2}$$\n",
    "\n",
    "(for simplicity, we start with a 1D example, and generalize to higher dimensions later)\n",
    "\n",
    "So how can we solve this in kwant?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The usual trick is to use the method of *finite differences*. We approximate:\n",
    "\n",
    "$$ \\frac{d \\psi(x)}{d x} \\approx \\frac{\\psi(x+a) - \\psi(x-a)}{2a}$$\n",
    "\n",
    "$$ \\frac{d \\psi(x)}{d x^2} \\approx \\frac{2 \\psi(x) - \\psi(x-a) - \\psi(x+a)}{a^2}$$\n",
    "\n",
    "We see that now only the wave function at points that are spaced by $a$ are entering here. Hence we can replace the continuous coordinate space by a discrete lattice with lattice spacing $a$:\n",
    "\n",
    "<img src=\"images/discretization.svg\" style=\"width: 70%\">\n",
    "\n",
    "where $x_i = i a$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "In the finite difference approximation, the Schrödinger equation then reads\n",
    "\n",
    "$$ H \\psi(x_i) = \\frac{\\hbar^2}{2 m^* a^2} \\left[2 \\psi(x_i) - \\psi(x_i-a) - \\psi(x_i+a) \\right] = E \\psi(x_i) $$\n",
    "\n",
    "It is useful to define the quantity $t = \\frac{\\hbar^2}{2 m^* a^2}$ which has units of energy (and depends on the lattice constant $a$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now everything to turn our continuum Hamiltonian into a tight-binding version. We use the lattice points $x_i$ as a basis: $\\left|x_i \\right>$. How can we now read off the values of the onsite and hopping enerrgies from the finite difference equation?\n",
    "\n",
    "When the Hamiltonian acts on a wave function, it \"brings\" the value of the wave function from some neighboring lattice point to the lattice point $x_i$ (see the right-hand side $E \\psi(x_i)$). Hence, if we find a term proportional to $t \\psi(x_j)$ in the finite difference\n",
    "equation, this corresponds to a hopping energy *from* $x_j$ *to* $x_i$, i.e. a tight-binding term $t \\left|x_i \\rangle \\langle x_j\\right|$.\n",
    "\n",
    "In our simple example we can thus read off:\n",
    "- onsite energy: $2 t$\n",
    "- hopping between nearest neighbors $-t$\n",
    "\n",
    "We thus finally end up with a tight-binding system that looks like this:\n",
    "\n",
    "<img src=\"images/tb_1d.svg\" style=\"width: 40%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Discretizing in higher dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The generalization to higher dimensions is straight-forward: For\n",
    "\n",
    "$$H = -\\frac{\\hbar^2}{2 m^*} \\nabla^2 $$\n",
    "\n",
    "we find\n",
    "- onsite energy $4 t$ (in 2D) or $6 t$ (in 3D)\n",
    "- hopping energy $-t$ for hopping between neighbors in $x$, $y$ and $z$-direction\n",
    "\n",
    "(*Control question: Why does the onsite energy depend on dimensionality?*)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Validity of the finite difference approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously, the discretization only works if it is fine enough. In particular, it must be fine enough to capture the changes in the wave function, i.e. $a \\ll \\lambda_F$. In practice, a good rool of thumb is that the approximation is good for energies $E \\lesssim t$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The first kwant system!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now put our knowledge into practice! Let's first activate kwant:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import kwant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes it's as easy as that!\n",
    "\n",
    "We'll do some more import of things we need later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%run matplotlib_setup.ipy\n",
    "\n",
    "import scipy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will describe a particle in 2D with the Hamiltonian $H= -\\frac{\\hbar^2}{2m^*} \\nabla^2$. The finite difference approach thus yields a *square lattice* with lattice constant $a$ (which we will set to 1). So in kwant we do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lat = kwant.lattice.square(a=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to confine the electrons in a finite region. Here, we choose a circle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t = 1\n",
    "r = 15\n",
    "\n",
    "def circle(pos):\n",
    "    x, y = pos\n",
    "    return x**2 + y**2 < r**2\n",
    "\n",
    "sys = kwant.Builder()\n",
    "sys[lat.shape(circle, (0,0))] = 4 * t\n",
    "sys[lat.neighbors()] = -t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kwant.plot(sys);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we are finished building the kwant system, we \"finalize\" it: this means the system will be brought into a format better suited for doing calculations, and it's shape cannot\n",
    "be changed any more (its values still can - see later tutorial)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sys = sys.finalized()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Building a kwant system means defining the Hamiltonian! We can get it from the finalized kwant system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ham = sys.hamiltonian_submatrix()\n",
    "plt.matshow(ham==0, cmap=\"gray\", interpolation=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But now we can start doing nice things! For example, we can now compute the eigenvalues and eigefunctions of this system very easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "eval, evec = la.eigh(ham)\n",
    "\n",
    "kwant.plotter.map(sys, abs(evec[:,0])**2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use some ipython magic and make an interactive plot showing the first 30 wave functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "def plot_wf(i=0):\n",
    "    print(\"Plotting wave function with index\", i)\n",
    "    print(\"energy:\", eval[i],\"x t\") \n",
    "    kwant.plotter.map(sys, abs(evec[:, i])**2)\n",
    "    \n",
    "interact(plot_wf, i=(0, 30))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your turn!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modify the code shown above and play with some possibilities. For example, you could:\n",
    "\n",
    "- Change the shape of the confined region. The circle has a rotational symmetry, which leads\n",
    "to degenerate eigenvalues (well, due to the discretization rather nearly degenerate). Take a shape that breaks this symmetry! For example, an ellipse can be defined through the equation $(x/r_1)^2+(y/r_2)^2 < 1$\n",
    "- Above, we made a homogeneous circle. With ``kwant.Builder()``, you can also make things position-dependent. For example, try to add an additional potential at site `(i,j)` using, say,\n",
    "``sys[lat(i, j)] = 8 * t``\n",
    "- Try to cut a hole in the above example. You could do this either by modifying the circle\n",
    "function, or by deleting points using ``del``.\n",
    "- Plot also the higher wave functions! At which point do they stop looking like you would expect it? Is this in agreement with the criterion for validity of the finite difference approximation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "IPython 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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}