calculate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bandstructures for 1D, 2D, and 3D nanowires\n",
    "For more information see: B. Nijholt, A. R. Akhmerov, *Orbital effect of magnetic field on the Majorana phase diagram*, [arXiv:1509.02675](https://arxiv.org/abs/1509.02675) [[pdf](https://arxiv.org/pdf/1509.02675.pdf)], [Phys. Rev. B 93, 235434 (2016)](http://journals.aps.org/prb/abstract/10.1103/PhysRevB.93.235434)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import copy\n",
    "from functools import partial\n",
    "from operator import itemgetter\n",
    "from time import time\n",
    "\n",
    "import adaptive\n",
    "import holoviews as hv\n",
    "import kwant\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sympy\n",
    "\n",
    "import wires\n",
    "\n",
    "adaptive.notebook_extension()\n",
    "sympy.init_printing()\n",
    "\n",
    "# hv.notebook_extension('bokeh')\n",
    "hv.notebook_extension(\"matplotlib\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wires.get_sympy_hamiltonian({}, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# System example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot\n",
    "\n",
    "syst_pars = dict(a=10, L=100, dim=1)\n",
    "syst = wires.make_wire(**syst_pars)\n",
    "kwant.plot(syst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "syst_pars = dict(a=10, L=100, dim=2, r=30)\n",
    "syst = wires.make_wire(**syst_pars)\n",
    "kwant.plot(syst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "syst_pars = dict(a=10, L=100, dim=3, r=30)\n",
    "syst = wires.make_wire(**syst_pars)\n",
    "kwant.plot(syst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 30\n",
    "triangle = wires.get_triangle(R)\n",
    "syst_pars = dict(a=10, L=100, dim=3, r=None, shape=triangle)\n",
    "syst = wires.make_wire(**syst_pars)\n",
    "kwant.plot(syst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Band structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bands(ylim, alpha, mu, Delta, B_x, B_y, B_z, params, lead_SC, lead_N):\n",
    "    params = dict(params, alpha=alpha, mu=mu, Delta=Delta, B_x=B_x, B_y=B_y, B_z=B_z)\n",
    "    ks = np.linspace(-0.9, 0.9, 201)\n",
    "    Es = wires.bands(lead_N.finalized(), params, ks)\n",
    "    kdims = [hv.Dimension(r\"$k$\", unit=r\"nm$^{-1}$\"), hv.Dimension(r\"$E$\", unit=r\"meV\")]\n",
    "    normal = hv.Path((ks, Es), kdims=kdims, label=\"normal\")\n",
    "    Es = wires.bands(lead_SC.finalized(), params, ks)\n",
    "    SC = hv.Path((ks, Es), kdims=kdims, label=\"SC\")\n",
    "    fermi_level = hv.HLine(0).opts(style=dict(linestyle=\"--\"))\n",
    "    return (normal * fermi_level).select(E=(-ylim, ylim)) + SC.select(E=(-ylim, ylim))\n",
    "\n",
    "\n",
    "kdims = [\n",
    "    hv.Dimension((\"ylim\"), range=(0.01, 50), default=5),\n",
    "    hv.Dimension((\"alpha\", r\"$\\alpha$\"), range=(0, 50), unit=\"meV·nm\", default=20),\n",
    "    hv.Dimension((\"mu\", r\"$\\mu$\"), range=(0.0, 30), unit=\"meV\", default=5),\n",
    "    hv.Dimension((\"Delta\", r\"$\\Delta$\"), range=(0, 1.0), unit=\"meV\", default=0.25),\n",
    "    hv.Dimension((\"B_x\", r\"$B_x$\"), range=(0, 10), unit=\"T\", default=1),\n",
    "    hv.Dimension((\"B_y\", r\"$B_y$\"), range=(0, 10), unit=\"T\", default=0),\n",
    "    hv.Dimension((\"B_z\", r\"$B_z$\"), range=(0, 0.5), unit=\"T\", default=0),\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1D band structure\n",
    "$$H(x) = \\left[ \\frac{p_x^2}{2m} - \\mu \\right]\\sigma_0\\otimes\\tau_z + \\alpha p_x \\sigma_y \\otimes \\tau_z + \\frac{1}{2}B_x\\sigma_x \\otimes \\tau_0  + \\Delta\\sigma_0\\otimes\\tau_x$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lead_SC = wires.make_lead_SC(10, dim=1)\n",
    "lead_N = wires.make_lead_normal(10, with_holes=False, dim=1)\n",
    "\n",
    "params = dict(V=0, g=50, **wires.constants)\n",
    "bands_1D = partial(bands, params=params, lead_SC=lead_SC, lead_N=lead_N)\n",
    "hv.DynamicMap(bands_1D, kdims=kdims)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D band structure\n",
    "\n",
    "$$H=\\left(\\frac{\\mathbf{p}^{2}}{2m}-\\mu\\right)\\sigma_{0}\\otimes\\tau_{z}+\\alpha\\left(p_{y}\\sigma_{x}\\otimes\\tau_{z}-p_{x}\\sigma_{y}\\otimes\\tau_{z}\\right)+\\frac{1}{2}B_{x}\\sigma_{x}\\otimes\\tau_{0}+\\frac{1}{2}B_{y}\\sigma_{y}\\otimes\\tau_{0}+\\frac{1}{2}B_{z}\\sigma_{z}\\otimes\\tau_{0}+\\Delta \\sigma_{0}\\otimes\\tau_{x}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "syst_pars = dict(a=10, r=30, dim=2)\n",
    "lead_SC = wires.make_lead_SC(**syst_pars)\n",
    "lead_N = wires.make_lead_normal(**syst_pars, with_holes=False)\n",
    "\n",
    "params = dict(a=syst_pars[\"a\"], V=0, g=50, **wires.constants)\n",
    "bands_2D = partial(bands, params=params, lead_SC=lead_SC, lead_N=lead_N)\n",
    "hv.DynamicMap(bands_2D, kdims=kdims)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D band structure\n",
    "\n",
    "\n",
    "$$H=\\left(\\frac{\\mathbf{p}^{2}}{2m}-\\mu\\right)\\sigma_{0}\\otimes\\tau_{z}+\\alpha\\left(p_{y}\\sigma_{x}\\otimes\\tau_{z}-p_{x}\\sigma_{y}\\otimes\\tau_{z}\\right)+\\frac{1}{2}B_{x}\\sigma_{x}\\otimes\\tau_{0}+\\frac{1}{2}B_{y}\\sigma_{y}\\otimes\\tau_{0}+\\frac{1}{2}B_{z}\\sigma_{z}\\otimes\\tau_{0}+\\Delta \\sigma_{0}\\otimes\\tau_{x}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "syst_pars = dict(a=10, r=30, dim=3)\n",
    "lead_SC = wires.make_lead_SC(**syst_pars)\n",
    "lead_N = wires.make_lead_normal(**syst_pars, with_holes=False)\n",
    "\n",
    "params = dict(a=syst_pars[\"a\"], V=0, g=50, **wires.constants)\n",
    "bands_3D = partial(bands, params=params, lead_SC=lead_SC, lead_N=lead_N)\n",
    "hv.DynamicMap(bands_3D, kdims=kdims)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Wave functions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wavefunction in the cross section of a 3D infinite lead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wavefunctions(lead, momentum, p):\n",
    "    h, t = lead.cell_hamiltonian(args=[p]), lead.inter_cell_hopping(args=[p])\n",
    "    h_k = lambda k: h + t * np.exp(1j * k) + t.T.conj() * np.exp(-1j * k)\n",
    "    vals, vecs = np.linalg.eigh(h_k(momentum))\n",
    "    indxs = np.argsort(abs(vals))\n",
    "    vecs = vecs[:, indxs]\n",
    "    vals = vals[indxs]\n",
    "    return vals, vecs\n",
    "\n",
    "\n",
    "syst_pars = dict(a=10, r=30, dim=3)\n",
    "lead_SC = wires.make_lead_SC(**syst_pars).finalized()\n",
    "\n",
    "params = dict(\n",
    "    a=syst_pars[\"a\"],\n",
    "    Delta=1,\n",
    "    B_x=1,  # in units of meV\n",
    "    mu=1,\n",
    "    V=0,\n",
    "    alpha=20,\n",
    "    g=1,\n",
    "    B_y=0,\n",
    "    B_z=0,\n",
    "    **wires.constants\n",
    ")\n",
    "params[\"mu_B\"] = 2\n",
    "\n",
    "\n",
    "wires.plot_wfs_in_cross_section(lead_SC, params, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wavefunction of finite system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "syst_pars = dict(a=10, r=30, dim=1, L=4000)\n",
    "syst = wires.make_wire(**syst_pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(  # should be topological\n",
    "    alpha=20, g=1, B_y=0, B_z=0, V=0, B_x=0.86, mu=0.42, Delta=0.384, **wires.constants\n",
    ")\n",
    "params[\"mu_B\"] = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = syst.hamiltonian_submatrix(params=params)\n",
    "Es, ψs = np.linalg.eigh(H)\n",
    "\n",
    "indices = np.abs(Es).argsort()\n",
    "\n",
    "rho = kwant.operator.Density(syst)\n",
    "wf = rho(ψs[:, indices[0]])\n",
    "\n",
    "hv.Curve(wf)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}