# Hofstadter Butterfly for Graphene

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.

    ClearAll[matrix];
matrix[p_,q_,nu_:0]:=Module[
{sigma},
sigma=p/q;
N@SparseArray[
{{m_,m_} -> 2Cos[2Pi*m*p/q + nu], {i_,j_}/;Abs[i-j] == 1 -> 1},{q,q}]]

ClearAll[attachsigma]
attachsigma[sigma_,lst_]:={sigma,#}&/@lst

fracs = Table[p/q, {q, 2, 80}, {p, 2, q}] // Flatten //
DeleteDuplicates;
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
]


• 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. Nov 9, 2019 at 1:38
• What exactly are you trying to plot? Are $\psi_0,\psi_1,\phi_0,\phi_1$ given? Nov 9, 2019 at 9:13
• @anderstood This will be a solution of this matrix equation versus $\frac{p}{q}$. These $\psi$ and $\phi$ are not needed in secular equation. Nov 9, 2019 at 20:27
• 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$
– L.K.
Oct 23, 2021 at 15:18

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


• 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. Nov 27, 2023 at 17:10
• @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. Nov 27, 2023 at 23:08
• @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.
– Bart
Nov 28, 2023 at 1:10
• @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.
– Bart
Nov 28, 2023 at 1:13
• 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. Nov 28, 2023 at 5:03
tic
clear all
qmax = 15;
for i = 1:1:qmax
for j = 1:1:qmax
A(i,j) = i/j;
end
end
B(:) = unique(reshape(A',1,[]));
for i = 1:1:length(B)
if B(i)<=1
C1(i) = B(i);
end
end
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)));
end
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));
end
%             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);
end
end
end
end
figure(1);
hold on
plot(MF,(50/3)*EE,'.b','MarkerSize',4)
xlabel('p/q')
ylabel('E[meV]');
set(gca,'Fontsize',20);
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.

• Please provide Mathematica code. Nov 27, 2023 at 23:19