We generate a fractal like-structure, called the Hofstadter Butterfly, which represents the energy levels of an electron travelling through a periodic lattice under the influence of a magnetic field.

The mathematical model related to the Hamiltonian of an electron in a two dimensional lattice, subject to a perpendicular (uniform) magnetic field is the Almost Mathieu (AM) operator or Harper operator, which is a discrete one-dimensional operator that acts on the Hilbert space, $\ell^2(\mathbb{Z})$, of the infinite sequences. It is defined by:

$$(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}$$

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
(J. Avron, P. H. M. v. Mouche, B. Simon, On the Measure of the Spectrum for the Almost Mathieu Operator, Commun Math Phys 132 (1990), 103-118).

For every irrational $\Phi$, and parameter $K>0$, the spectrum of the AM operator is a Cantor set (A Avila, S Jitomirskaya, The Ten Martini Problem, Annals of math 170 (2009), 303-342).

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)$ is periodic and the eigenvalue problem:

$$(H_{\Phi, K, \theta}u)_n=E u_n$$

reduces to a matrix eigenvalue problem associated to the following periodic Jacobi matrices, called Harper matrices:

$$Ha(p, q, K, \theta, s)=\left(\begin{array}{ccccccc}K\cos(2\pi 0 p/q+\theta)&1 &0&\ldots&0& 0& s\\ 1& K\cos(2\pi p/q+\theta)&1&\ldots&0&0&0\\ \vdots&\vdots&\vdots&\ldots&\vdots&\vdots&\vdots\\ 0&0&0&\ldots&1&K\cos(2\pi (q-2) p/q+\theta)&1\\ s&0&0&\ldots&0&1&K\cos(2\pi (q-1) p/q+\theta)\end{array}\right)$$

with $s=\pm 1$.

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)$.

The eigenvalues $E_i$, respectively $E'_i$, $i=0, 1, \ldots, q-1$, of the two matrices can be ordered as follows:

$$E_{2i} < E_{2i+1}\leq E_{2i+2}, i\geq 0$$

respectively: $$E'_{2i}\leq E'_{2i+1} < E'_{2i+2}, i\geq 0$$ and the two series are interlaced:

$$E_0 < E'_0 \leq E'_1 < E_1\leq E_2<E'_2\leq \cdots$$

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).

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)$.

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 $p, q$ are relative prime numbers. But since $\cos$ is an odd $2\pi$-periodic function, we have that

$$\cos(2\pi n p/q)=\cos(2\pi n(q-p)/q),$$

and thus $$Ha(p, q, s)=Ha(q-p, s)$$.

Hence only the spectrum of the Harper matrices $Ha(p, q, s)$, with $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.

In [1]:
import platform
print(f'Python version: {platform.python_version()}')

Python version: 3.6.4

In [2]:
import plotly
plotly.__version__

Out[2]:
'3.1.1'
In [3]:
import plotly.graph_objs as go
import numpy as np
from numpy import pi


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 of data).

In [4]:
def Gear(n, s=1):
# Generates a  Gear-type matrix, i.e. a periodic Jacobi matrix G=(0,..0; 1,...1; +-1), with 0 on the principal diagonal
#  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

G=np.diag(np.ones(n - 1), -1) + np.diag(np.ones(n - 1), 1)

G[0][n-1]=s
G[n-1][0]=s
return G

def eigs_Harper(p, q, s, K):

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)
Hd= np.diag(d)
G = Gear(q, s)

return list(np.linalg.eigvalsh(Hd+G))#eigenvalues of the Harper matrix

def gcd(a, b): # Greatest Common Divisor
if b == 0: return a
return gcd(b, a % b)

In [5]:
def get_butterfly_points_even(qmax=101,  s=1, K=2):# for  qmax=101 value we define 1036 irreducible fractions p/q,
#and compute the eigvals for 1036/2=518 Harper matrices

phi=[]# the list of  of rational magnetic flux values, p/q
E=[]# the list of energies
text=[]# the list of hover strings
#take  all rational numbers p/q of even denominator, q<qmax
for q in range(4, qmax, 2):
for p in range(1, q//2, 2):
if gcd(p, q) == 1:
phi.extend([p/q]*q+ [(q-p)/q]*q) #insert q copies of  p/q, respectively (q-p)/q,
#because the corresponding Harper matrix H(p,q), resp H(q-p, p), has q eigvals
eigs_pq=eigs_Harper(p, q, s, K)
E.extend(eigs_pq*2)
p_text=[f"(p, q) = {(p,q)}"]*q+[f"(p,q) = {(q-p, q)}"]*q
text.extend([f"{t}<br>E = {round(e, 3)}" for t, e in zip(p_text, eigs_pq*2)])
return phi, E, text

In [6]:
def get_butterfly_points_odd(qmax=70, s=1, K=2):

phi=[]
E=[]
text=[]
#take  all rational numbers p/q of odd denominator, q<qmax
for q in range(5, qmax, 1):
for p in range(1, q//2+1, 1):
if gcd(p, q) == 1:
phi.extend([p/q]*q+ [(q-p)/q]*q)
eigs_pq=eigs_Harper(p, q, s, K)
E.extend(eigs_pq*2)
p_text=[f"(p, q) = {(p,q)}"]*q+[f"(p,q) = {(q-p, q)}"]*q
text.extend([f"{t}<br>E = {round(e, 3)}" for t, e in zip(p_text, eigs_pq*2)])
return phi, E, text

In [7]:
def get_butterfly_trace(phi, E, text, color='rgb(255,215, 0)', marker_size=1):
return dict(type='scatter',
x=phi,
y=E,
mode='markers',
text=text,
marker=dict(color=color, size=marker_size),
hoverinfo='text')


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).

Let us test get_butterfly_points_even():

In [8]:
%time phi1, E1, text1=get_butterfly_points_even(qmax=101)

Wall time: 598 ms

In [9]:
len(E1)#points are plotted

Out[9]:
69836

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.

The user can experiment plotting the Hofstadter butterfly associated to both data, succesively, or to their union.

In [10]:
data=[get_butterfly_trace(phi1, E1, text1)]

In [11]:
axis_style=dict(showline=True,
mirror=True,
zeroline=False,
showgrid=False,
ticklen=4)

In [12]:
layout=dict(title='Hofstadter butterfly<br> K=2, s=1',
font=dict(family='Balto'),
width=600, height=675,
autosize=False,
showlegend=False,
xaxis=dict(axis_style, **dict( title='Phi (magnetic flux)', dtick=0.25)),
yaxis=dict(axis_style, **dict( title='E (Energy)')),
hovermode='closest',
plot_bgcolor='rgb(10,10,10)')

In [13]:
fw=go.FigureWidget(data=data, layout=layout)

In [ ]:
fw # running this cell  the FigureWidget is plotted in the next one


An alternative is to send the figure to Plotly cloud:

In [15]:
import plotly.plotly as py

import warnings
warnings.filterwarnings("ignore")

py.sign_in('empet', '')

The draw time for this plot will be slow for all clients.

Out[15]:

It's naturally to ask yourself how behaves the spectrum of the associated Harper-type matrix, $H(p,q, s=-1)$. Taking into account the interlacing property mentioned above we expect to get also a butterfly.

Let us plot it:

In [21]:
phi_m1, E_m1, text_m1=get_butterfly_points_even(qmax=101, s=-1)
data1=[get_butterfly_trace(phi_m1, E_m1, text_m1, color='rgb(192,192,192)')]
fw1=go.FigureWidget(data=data1, layout=layout)
fw1.layout.title='Hofstadter butterfly<br> s = -1, K=2'

In [24]:
fw1

In [23]:
py.iplot(fw1, filename='Hofstadterm1')

The draw time for this plot will be slow for all clients.

Out[23]:

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, 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.

For example:

In [19]:
eigs_Harper(3, 10, 1, 2)

Out[19]:
[-2.698183697478866,
-2.50205068512151,
-2.4363187022725965,
-0.6778990365402857,
-0.2536753455832333,
0.25367534558323285,
0.6778990365402842,
2.4363187022725965,
2.5020506851215116,
2.698183697478866]
In [20]:
eigs_Harper(3, 10, -1, 2)

Out[20]:
[-2.685510165162846,
-2.5485149479215528,
-2.4002118113237896,
-0.7089410427986523,
-0.17173401423357063,
0.17173401423357224,
0.7089410427986518,
2.40021181132379,
2.5485149479215536,
2.6855101651628432]

Define data for plotting both above Hofstadter butterflies on the same figure.

In [ ]:
data_global=[get_butterfly_trace(phi1, E1, text1, color='rgb(192,192,192)'),
get_butterfly_trace(phi_m1, E_m1, text_m1, color='rgb(255,215,0)')]
fw_global=go.FigureWidget(data=data_global, layout=layout)

In [ ]:
fw_global


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 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.

The Java, respectively the Python implementation of this algorithm can be found here, and here.

This algorithm is computationally expensive. To ckeck it just run the above Python implementation.

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.

In [ ]: