Hofstadter Butterfly for Graphene

Question

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
  ]

Answer 1

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

Answer 2

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.