I am trying to extend the Harper equation from two dimensional square lattice to monolayer graphene by using mathematica. For square lattice the code is given here "Poor rendering of fractals". The problem while dealing with graphene is that this equation has become matrix inside matrix. This is the eigen value equation that I am trying to plot,

$\begin{bmatrix} 0 & -t \\ -t & 0 \end{bmatrix}\begin{bmatrix}\psi_{m}\\\phi_{m}\end{bmatrix}=\begin{bmatrix}0 & 2*Cos[k-\frac{i\pi p}{q}(m-\frac{1}{6})]\\0&0\end{bmatrix}\begin{bmatrix}\psi_{m-1}\\\phi_{m-1}\end{bmatrix}+\begin{bmatrix}0 & 0\\2*Cos[k-\frac{i\pi p}{q}(m-\frac{5}{6})]&0\end{bmatrix}\begin{bmatrix}\psi_{m+1}\\\phi_{m+1}\end{bmatrix}$

One can take $t=1$ and $k=0$. However $\frac{p}{q}$ can take any rational fraction. Any help will be highly appreciated. This is required to see first the posted version of two dimensional square lattice case to understand this question,"Poor rendering of fractals"(Go for Hofstadter butterfly). Thanks

This is the code for two dimension square lattice copied from https://mathematica.stackexchange.com/questions/2392/poor-rendering-of-fractals.

        {{m_,m_} -> 2Cos[2Pi*m*p/q + nu], {i_,j_}/;Abs[i-j] == 1 -> 1},{q,q}]]


fracs = Table[p/q, {q, 2, 80}, {p, 2, q}] // Flatten // 
pq = {Numerator@#, Denominator@#} & /@ fracs;
(ens = Eigenvalues[#] & /@ (matrix[#[[1]], #[[2]]] & /@ pq);) // Timing
pts = Flatten[#, 1] &@MapThread[attachsigma, {fracs, ens}];
plot = Graphics[
  {PointSize[0.001], Point[pts]},
  AspectRatio -> 1,
  ImageSize -> Full

  • $\begingroup$ For two dimensional square lattice the eigen value equation is $\psi_{m+1}+\psi_{m-1}+2Cos(2m\pi \frac{p}{q})\psi_{m}=0$. The link above plot this equation. $\endgroup$ Commented Nov 9, 2019 at 1:38
  • $\begingroup$ What exactly are you trying to plot? Are $\psi_0,\psi_1,\phi_0,\phi_1$ given? $\endgroup$
    – anderstood
    Commented Nov 9, 2019 at 9:13
  • $\begingroup$ @anderstood This will be a solution of this matrix equation versus $\frac{p}{q}$. These $\psi$ and $\phi$ are not needed in secular equation. $\endgroup$ Commented Nov 9, 2019 at 20:27
  • $\begingroup$ I think you can make the two could equations into one. Then it is possible to use the above code for any of the $\psi$ or $\phi$ $\endgroup$
    – L.K.
    Commented Oct 23, 2021 at 15:18

2 Answers 2


The most intuitive way to construct the Hofstadter butterfly for graphene is to write down two simultaneous Schrödinger equations for the basis sites, e.g.

$$E \Psi_{m,n}^A = -t e^{\mathrm{i}\theta_1}\Psi_{m,n}^B -te^{\mathrm{i}\theta_2}\Psi_{m-1,n}^B - t e^{\mathrm{i}\theta_3}\Psi_{m,n-1}^B,$$

$$E \Psi_{m,n}^B = -t e^{\mathrm{i}\theta_4}\Psi_{m+1,n}^A -te^{\mathrm{i}\theta_5}\Psi_{m,n+1}^A - t e^{\mathrm{i}\theta_6}\Psi_{m,n}^A,$$

where $\theta$ is the appropriate Peierls phase. You can then write these equations in matrix form, such that

$$\begin{pmatrix} H^{AA} & H^{AB} \\ H^{BA} & H^{BB} \end{pmatrix} \begin{pmatrix} \Psi^A \\ \Psi^B \end{pmatrix} = E \begin{pmatrix} \Psi^A \\ \Psi^B \end{pmatrix}.$$

For each block of this matrix, you can then construct the $q\times q$ Harper matrix, as before. Hence, you will ultimately be left with a $2q\times 2q$ matrix, which you can simply diagonalize using any programming language, e.g. Mathematica, Matlab, Python, etc.

The Hofstadter butterfly for nearest-neighbor hopping on the honeycomb lattice, as well as for any combination of hoppings on any lattice, can be easily computed using the open-source Python package HofstadterTools (https://hofstadter.tools). In the package, you can run the following command in the code directory to benchmark your results:

python butterfly.py -lat honeycomb -q 1999

enter image description here

  • 2
    $\begingroup$ Interesting, but out of scope for this site. I did not down-vote, by the way. But being that far out of scope should explain why that would happen. $\endgroup$ Commented Nov 27, 2023 at 17:10
  • $\begingroup$ @Daniel Lichtblau, the above Hofstadter butterfly is not correct. see R. Rammal, J. Phys. (Paris) 46, 1345 (1985). There should be no gap at n = 0.5. $\endgroup$ Commented Nov 27, 2023 at 23:08
  • $\begingroup$ @DanielLichtblau I apologize if my answer was off-topic for this site. Instead of focusing on the Mathematica code, I wanted to address the root cause of the problem. The OP said that they could solve the problem for the square lattice but were having problems with the honeycomb lattice because they got a "matrix inside matrix". My point was simply that if you formulate the problem correctly, it boils down to single matrix diagonalization, as in the square lattice case. $\endgroup$
    – Bart
    Commented Nov 28, 2023 at 1:10
  • $\begingroup$ @MuhammadImran Apologies, perhaps I should have clarified. The value of q that I used for the above plot is q=199, and I sweep through integer p. Since 199 is a prime number, there is no nphi=1/2 line in the above plot. $\endgroup$
    – Bart
    Commented Nov 28, 2023 at 1:13
  • $\begingroup$ I take no issue with the intent, which certainly appears to be in good faith. It’s a gray area here, since the guidance looks useful even if the tools used are otherwise outside the scope of the forum. Classic square-peg-round-hole problem. Or maybe the other way around. $\endgroup$ Commented Nov 28, 2023 at 5:03
clear all
qmax = 15;
for i = 1:1:qmax                                             
    for j = 1:1:qmax                                         
        A(i,j) = i/j;                                        
B(:) = unique(reshape(A',1,[]));                             
for i = 1:1:length(B)                                        
    if B(i)<=1
        C1(i) = B(i);
C2 = nonzeros(C1');                                          
[p,q1] = numden(sym(C2));                                                                           
Q(:) = double(q1(:));                                        
P(:) = double(p(:));                                         
aalpha(:) = C2(:); 
counter = 0;
a = 1;
k2 = 0:0.1:2*pi;
for i = 1:1:length(aalpha)
    k1 = 0:0.1/Q(i):2*pi/(Q(i));
    for ja = 1:1:length(k2)
        for l1 = 1:1:length(k1)
            h = zeros(Q(i),Q(i));
            Hf = zeros(2*Q(i),2*Q(i));
            for j = 1:2:2*Q(i)-1
%                 Hf(j,j+1) = (exp(1i*(k2(ja)*a))+exp(-1i*(2*pi*(aalpha(i))*(j+1)/2)));
                Hf(j,j+1) = (exp(1i*(4*pi/3*(aalpha(i))*(j+1)/2+k2(ja)*a))+exp(-1i*(2*pi/3*(aalpha(i))*(j+1)/2)));
            for j1 = 2:2:2*Q(i)-1
%                 Hf(j1,j1+1) = exp(-1i*(2*pi*(aalpha(i))*(j1+2)/2))*exp(-1i*(k1(l1)*a));
                  Hf(j1,j1+1) = exp(-1i*(2*pi/3*(aalpha(i))*(j1+2)/2));
%             Hf(1,2*Q(i)) = Hf(1,2*Q(i))+exp(1i*(k1(l1)*a));
            Hf(1,2*Q(i)) = exp(-1i*(2*pi/3*(aalpha(i))*(Q(i))-Q(i)*k1(l1)*a));
            L = eig(Hf+ctranspose(Hf));
            for ii = 1:1:2*Q(i)
                counter = counter+1;
                EE(counter) = L(ii);
                MF(counter) = aalpha(i);
hold on
box on
hold off
% save('Rammal.mat','EE','MF')
timeelapsed = toc

This is the right Hamiltonian for Hofstadter Butterfly for Monolayer Graphene. It reproduces the Rammal work as well. Its written in Matlab.

  • 9
    $\begingroup$ Please provide Mathematica code. $\endgroup$
    – bbgodfrey
    Commented Nov 27, 2023 at 23:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.