{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Hofstadter Butterfly </center> ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate a fractal like-structure, called the Hofstadter Butterfly, which represents the energy\n",
    "levels of  an electron travelling through a periodic lattice under the influence of a\n",
    "magnetic field."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mathematical model related to the  Hamiltonian of an electron in a two dimensional lattice,\n",
    "subject to a perpendicular (uniform) magnetic field is the Almost Mathieu (AM) operator or Harper operator, which is \n",
    "a discrete one-dimensional operator that acts on the Hilbert space, $\\ell^2(\\mathbb{Z})$, of the infinite sequences. It is defined  by:\n",
    "\n",
    "$$(H_{\\Phi, K, \\theta}u)_n=u_{n+1}+u_{n-1}+K\\cos(n\\Phi +\\theta) u_n, \\quad\\Phi, K, \\theta\\in\\mathbb{R}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the magnetic flux  penetrating the lattice corresponds to a rational number $p/q$, i.e.  $\\Phi=2\\pi p/q$, with $p,q$  relative prime integers, the spectrum of the above operator  consists in $q$ bands  (closed intervals) separated by gaps  \n",
    "(*J. Avron, P. H. M. v. Mouche,  B. Simon*, *On the Measure of the Spectrum\n",
    "for the Almost Mathieu Operator, Commun Math Phys 132 (1990), 103-118*). \n",
    " \n",
    " For every  irrational $\\Phi$, and parameter $K>0$, the spectrum of the AM operator is a Cantor set \n",
    " (*A Avila, S Jitomirskaya, The Ten Martini Problem, Annals of math 170 (2009), 303-342*). \n",
    " \n",
    " \n",
    "For  a flux  $\\Phi=2\\pi n p/q$, corresponding to a rational number,  the potential $V_\\theta(n)=K\\cos(2\\pi n p/q+\\theta)$\n",
    "is periodic and the eigenvalue problem:\n",
    "\n",
    "$$(H_{\\Phi, K, \\theta}u)_n=E u_n$$\n",
    "\n",
    "reduces to a matrix eigenvalue problem associated to the following periodic Jacobi matrices, called Harper matrices:\n",
    "\n",
    "$$\n",
    "Ha(p, q, K, \\theta, s)=\\left(\\begin{array}{ccccccc}K\\cos(2\\pi 0 p/q+\\theta)&1 &0&\\ldots&0& 0& s\\\\\n",
    "                       1& K\\cos(2\\pi  p/q+\\theta)&1&\\ldots&0&0&0\\\\\n",
    "                        \\vdots&\\vdots&\\vdots&\\ldots&\\vdots&\\vdots&\\vdots\\\\\n",
    "                         0&0&0&\\ldots&1&K\\cos(2\\pi (q-2) p/q+\\theta)&1\\\\\n",
    "                         s&0&0&\\ldots&0&1&K\\cos(2\\pi (q-1) p/q+\\theta)\\end{array}\\right)$$\n",
    "                         \n",
    "with $s=\\pm 1$.\n",
    "\n",
    "More precisely, the spectrum of the  operator, $\\sigma(H_{2\\pi p/q, K, \\theta})$ is a union of intervals (bands) whose ends are the interlacing eigenvalues of the two Harper type matrices $Ha(p,q, K, \\theta, 1)$,  $Ha(p,q, K, \\theta, -1)$.\n",
    "\n",
    "The eigenvalues  $E_i$, respectively $E'_i$, $i=0, 1, \\ldots, q-1$, of the two matrices can be  ordered as follows:\n",
    "\n",
    "$$E_{2i} < E_{2i+1}\\leq E_{2i+2}, i\\geq 0$$\n",
    "respectively:\n",
    "$$E'_{2i}\\leq E'_{2i+1} < E'_{2i+2}, i\\geq 0$$\n",
    "and the two series are interlaced:\n",
    "\n",
    "$$E_0 < E'_0 \\leq  E'_1 < E_1\\leq E_2<E'_2\\leq \\cdots$$ \n",
    "\n",
    "\n",
    "The Hofstadter butterfly was defined and studied by the physicist Douglas Hofstadter  in 1976. It is a graphical representation   of all possible energies (eigenvalues) of the Harper matrices $H(p,q,s=1)$ corresponding to the  rational values $p/q$ in  [0,1). \n",
    "\n",
    "Hence to get the Hofstadter butterfly we have to plot  all points of coordinates, $(p/q, E_i)$, with $p/q\\in [0,1)$, $i=0, 1, \\ldots q-1$. For each $p/q$,   $E_i$  runs over the q eigenvalues of the Harper matrix,   $Ha(p, q, s)$.\n",
    "\n",
    "For  any $q<qmax$ we should compute the eigenvalues of all matrices $Ha(p, q, s=1)$, with $p\\in \\{1, 2, \\ldots q-1\\}$, such that\n",
    "$p, q$ are relative prime numbers. But since $\\cos$ is an odd $2\\pi$-periodic function,  we have that   \n",
    "\n",
    "$$\\cos(2\\pi n p/q)=\\cos(2\\pi n(q-p)/q),$$ and thus\n",
    " $$Ha(p, q, s)=Ha(q-p, s)$$.\n",
    " \n",
    " Hence  only the spectrum of the Harper matrices $Ha(p, q, s)$, with\n",
    "$p\\in\\{1, 2, \\ldots, q//2\\}$, if $q$ is odd, respectively $p\\in\\{1, 2, \\ldots, q/2-1\\}$, if $q$ is even, are calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python version: 3.6.4\n"
     ]
    }
   ],
   "source": [
    "import platform\n",
    "print(f'Python version: {platform.python_version()}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'3.1.1'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import plotly\n",
    "plotly.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objs as go\n",
    "import numpy as np\n",
    "from numpy import pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate the Hofstadter butterfly associated to the Harper operator corresponding to K=2, $\\theta=0$ and s=1 (these parameters were used by Hofstadter himself in the first numerical computation\n",
    "                                                                                                                  of data)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Gear(n, s=1): \n",
    "    # Generates a  Gear-type matrix, i.e. a periodic Jacobi matrix G=(0,..0; 1,...1; +-1), with 0 on the principal diagonal\n",
    "    #  1  in the positions G[i][i+1], G[i-1][i], and G[0][n-1], G[n-1][0]=s with s=1 or -1\n",
    "    \n",
    "    G=np.diag(np.ones(n - 1), -1) + np.diag(np.ones(n - 1), 1)\n",
    "    \n",
    "    G[0][n-1]=s\n",
    "    G[n-1][0]=s\n",
    "    return G\n",
    "\n",
    "\n",
    "def eigs_Harper(p, q, s, K):\n",
    "   \n",
    "    d=[K*np.cos(2*np.pi*m*p/q) for m in range(q)] #define the diagonal of the Harper matrix   Ha(p,q)    \n",
    "    Hd= np.diag(d)\n",
    "    G = Gear(q, s)\n",
    "    \n",
    "    return list(np.linalg.eigvalsh(Hd+G))#eigenvalues of the Harper matrix \n",
    "\n",
    "def gcd(a, b): # Greatest Common Divisor\n",
    "    if b == 0: return a\n",
    "    return gcd(b, a % b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_butterfly_points_even(qmax=101,  s=1, K=2):# for  qmax=101 value we define 1036 irreducible fractions p/q, \n",
    "                                             #and compute the eigvals for 1036/2=518 Harper matrices\n",
    " \n",
    "    phi=[]# the list of  of rational magnetic flux values, p/q\n",
    "    E=[]# the list of energies\n",
    "    text=[]# the list of hover strings\n",
    "    #take  all rational numbers p/q of even denominator, q<qmax\n",
    "    for q in range(4, qmax, 2):    \n",
    "        for p in range(1, q//2, 2):\n",
    "            if gcd(p, q) == 1:\n",
    "                phi.extend([p/q]*q+ [(q-p)/q]*q) #insert q copies of  p/q, respectively (q-p)/q, \n",
    "                                           #because the corresponding Harper matrix H(p,q), resp H(q-p, p), has q eigvals\n",
    "                eigs_pq=eigs_Harper(p, q, s, K) \n",
    "                E.extend(eigs_pq*2)           \n",
    "                p_text=[f\"(p, q) = {(p,q)}\"]*q+[f\"(p,q) = {(q-p, q)}\"]*q\n",
    "                text.extend([f\"{t}<br>E = {round(e, 3)}\" for t, e in zip(p_text, eigs_pq*2)])\n",
    "    return phi, E, text"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_butterfly_points_odd(qmax=70, s=1, K=2):\n",
    "\n",
    "    phi=[]\n",
    "    E=[]\n",
    "    text=[]\n",
    "    #take  all rational numbers p/q of odd denominator, q<qmax\n",
    "    for q in range(5, qmax, 1):    \n",
    "        for p in range(1, q//2+1, 1):\n",
    "            if gcd(p, q) == 1:\n",
    "                phi.extend([p/q]*q+ [(q-p)/q]*q) \n",
    "                eigs_pq=eigs_Harper(p, q, s, K) \n",
    "                E.extend(eigs_pq*2)           \n",
    "                p_text=[f\"(p, q) = {(p,q)}\"]*q+[f\"(p,q) = {(q-p, q)}\"]*q\n",
    "                text.extend([f\"{t}<br>E = {round(e, 3)}\" for t, e in zip(p_text, eigs_pq*2)])\n",
    "    return phi, E, text"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_butterfly_trace(phi, E, text, color='rgb(255,215, 0)', marker_size=1):\n",
    "    return dict(type='scatter',\n",
    "                x=phi,\n",
    "                y=E,\n",
    "                mode='markers',\n",
    "                text=text,  \n",
    "                marker=dict(color=color, size=marker_size),\n",
    "                hoverinfo='text')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the first sight the effective computation of the spectrum for many Harper matrices seems to be cumbersome. But with Anaconda Python packages it is very fast, because Anaconda packaged MKL-powered binary versions of some of the most popular numerical/scientific Python libraries into MKL Optimizations that improve performance. (MKL stands for Intel™ Math Kernel Library, a set of vectorized math routines that accelerate math functions).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us test `get_butterfly_points_even()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 598 ms\n"
     ]
    }
   ],
   "source": [
    "%time phi1, E1, text1=get_butterfly_points_even(qmax=101) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69836"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(E1)#points are plotted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generating the butterfly data,  with the function `get_butterfly_points_odd()` and then via `get_butterfly_points_even()`, corresponding both to the same s, we get two indistinguishable plots. \n",
    "\n",
    "The user can experiment plotting the Hofstadter butterfly associated to both data, succesively, or to their union. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "data=[get_butterfly_trace(phi1, E1, text1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis_style=dict(showline=True, \n",
    "                mirror=True, \n",
    "                zeroline=False, \n",
    "                showgrid=False, \n",
    "                ticklen=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "layout=dict(title='Hofstadter butterfly<br> K=2, s=1',\n",
    "            font=dict(family='Balto'),\n",
    "            width=600, height=675,\n",
    "            autosize=False,\n",
    "            showlegend=False,\n",
    "            xaxis=dict(axis_style, **dict( title='Phi (magnetic flux)', dtick=0.25)),\n",
    "            yaxis=dict(axis_style, **dict( title='E (Energy)')),\n",
    "            hovermode='closest',\n",
    "           plot_bgcolor='rgb(10,10,10)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "fw=go.FigureWidget(data=data, layout=layout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fw # running this cell  the FigureWidget is plotted in the next one"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative is to send the figure to Plotly cloud:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The draw time for this plot will be slow for all clients.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<iframe id=\"igraph\" scrolling=\"no\" style=\"border:none;\" seamless=\"seamless\" src=\"https://plot.ly/~empet/14979.embed\" height=\"675px\" width=\"600px\"></iframe>"
      ],
      "text/plain": [
       "<plotly.tools.PlotlyDisplay object>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import plotly.plotly as py\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "py.sign_in('empet', '')\n",
    "py.iplot(fw, filename='Hofstadter1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's naturally to ask yourself  how behaves the spectrum of the associated Harper-type matrix, $H(p,q, s=-1)$.\n",
    "Taking into account the interlacing property mentioned above we expect to get also a butterfly.\n",
    "\n",
    "Let us plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_m1, E_m1, text_m1=get_butterfly_points_even(qmax=101, s=-1)\n",
    "data1=[get_butterfly_trace(phi_m1, E_m1, text_m1, color='rgb(192,192,192)')]\n",
    "fw1=go.FigureWidget(data=data1, layout=layout)\n",
    "fw1.layout.title='Hofstadter butterfly<br> s = -1, K=2'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "fw1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The draw time for this plot will be slow for all clients.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<iframe id=\"igraph\" scrolling=\"no\" style=\"border:none;\" seamless=\"seamless\" src=\"https://plot.ly/~empet/14981.embed\" height=\"675px\" width=\"600px\"></iframe>"
      ],
      "text/plain": [
       "<plotly.tools.PlotlyDisplay object>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "py.iplot(fw1, filename='Hofstadterm1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two butterflies  look very similar. Even if we plot them in the same figure, we cannot distinguish the gold points from the silver ones,\n",
    "because the distance between two consecutive eigenvalues, one in H(p,q, s=1) and another in H(p,q, s=-1) is very small."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-2.698183697478866,\n",
       " -2.50205068512151,\n",
       " -2.4363187022725965,\n",
       " -0.6778990365402857,\n",
       " -0.2536753455832333,\n",
       " 0.25367534558323285,\n",
       " 0.6778990365402842,\n",
       " 2.4363187022725965,\n",
       " 2.5020506851215116,\n",
       " 2.698183697478866]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigs_Harper(3, 10, 1, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-2.685510165162846,\n",
       " -2.5485149479215528,\n",
       " -2.4002118113237896,\n",
       " -0.7089410427986523,\n",
       " -0.17173401423357063,\n",
       " 0.17173401423357224,\n",
       " 0.7089410427986518,\n",
       " 2.40021181132379,\n",
       " 2.5485149479215536,\n",
       " 2.6855101651628432]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigs_Harper(3, 10, -1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define data for plotting both above Hofstadter butterflies on the same figure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_global=[get_butterfly_trace(phi1, E1, text1, color='rgb(192,192,192)'), \n",
    "             get_butterfly_trace(phi_m1, E_m1, text_m1, color='rgb(255,215,0)')]\n",
    "fw_global=go.FigureWidget(data=data_global, layout=layout)\n",
    "fw_global.layout.title='Hofstadter Butterfly'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fw_global"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The most known and recommended algorithm on physics.stackexchange, math.stackexchange, for computing the points on the Hofstadter butterfly consists in reducing the Harper matrix to a tridiagonal form and then\n",
    "using a recursive formula to get its characteristic polynomial. This algorithm is presented in the book *W Kinzel/G Reents, Physics by Computer, Springer Verlag*. Moreover for its drawing as a raster image of resolution 500 x 500 or so, the authors suggest to map each pixel to a point $(a, E)$ in the rectangle $[0,1]\\times [-4,4]$, and then to evaluate two consecutive polynomials from the recursive formula  at this point. If the corresponding values  satisfy some condition  one decides  that the pixel mapped to $(a, E)$ is on  the Hofstadter butterfly. \n",
    "\n",
    "The Java, respectively the Python implementation of this algorithm can be found \n",
    "[here](https://web.archive.org/web/20010528223855/http://wptx15.physik.uni-wuerzburg.de/TP3/source_java/hofstadter.java), and \n",
    "[here](http://code.activestate.com/recipes/578670-hofstadter-butterfly-fractal/). \n",
    "\n",
    "This algorithm is computationally expensive. To ckeck it just run the above Python implementation.\n",
    "\n",
    "The method presented in this notebook is simple and straightforward, because  we generate a  SVG  image (i.e. a vector image), and no (approximative) mapping of a grid (raster image) to the real rectangle is needed. Moreover we get a performance benefit from Intel MKL provided by Anaconda."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}