霍夫斯塔特蝴蝶（Hofstadter butterfly）的文献：Douglas R. Hofstadter, "Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields", Phys. Rev. B 14, 2239 (1976).
有限体系方格子的哈密顿量：方格子模型在实空间中的哈密顿量形式。
引入磁场：佩尔斯替换 Peierls substitution。
霍夫斯塔特蝴蝶的代码为：
"""
This code is supported by the website: https://www.guanjihuan.com
The newest version of this code is on the web page: https://www.guanjihuan.com/archives/8491
"""
import numpy as np
import cmath
import matplotlib.pyplot as plt
def main():
for n in np.arange(1, 11):
print('n=', n)
width = n
length = n
B_array = np.arange(0, 1, 0.001)
eigenvalue_all = np.zeros((B_array.shape[0], width*length))
i0 = 0
for B in B_array:
# print(B)
h = hamiltonian(width, length, B)
eigenvalue, eigenvector = np.linalg.eig(h)
eigenvalue_all[i0, :] = np.real(eigenvalue)
i0 += 1
plt.plot(B_array, eigenvalue_all, '.r', markersize=0.5)
plt.title('width=length='+str(n))
plt.xlabel('B*a^2/phi_0')
plt.ylabel('E')
plt.savefig('width=length='+str(n)+'.jpg', dpi=300)
plt.close('all') # 关闭所有plt，防止循环画图时占用内存
# plt.show()
def hamiltonian(width, length, B): # 方格子哈密顿量
h = np.zeros((width*length, width*length), dtype=complex)
# y方向的跃迁
for x in range(length):
for y in range(width-1):
h[x*width+y, x*width+y+1] = 1
h[x*width+y+1, x*width+y] = 1
# x方向的跃迁
for x in range(length-1):
for y in range(width):
h[x*width+y, (x+1)*width+y] = 1*cmath.exp(-2*np.pi*1j*B*y)
h[(x+1)*width+y, x*width+y] = 1*cmath.exp(2*np.pi*1j*B*y)
return h
if __name__ == "__main__":
main()
计算结果为：
参考资料：
[1] Landau levels, molecular orbitals, and the Hofstadter butterfly in finite systems
[2] 知乎：二维方格Hofstadter's butterfly
[3] 知乎：什么是凝聚态中的 Hofstadter butterfly（霍夫斯塔特蝴蝶）？
【说明：本站主要是个人的一些笔记和代码分享，内容可能会不定期修改。为了使全网显示的始终是最新版本，这里的文章未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】
关老师您好，对于两带体系，是否能像上述所写方式添加磁场呢？比如跃迁矩阵都是2×2泡利矩阵的形式Tx,Ty. 能否直接乘上e指数？
嗯，是可以的，如果和晶格位置无关，更多的维度当成内禀自由度就可以了，可直接乘。
关老师您好，请问如果把方格子模型Ny方向加周期性边界条件，也就是相当于把y方向的边卷起来，那么卷起来应该怎么加磁通才能保证磁通也是周期性的呢？
Ny取为磁周期的大小，这样才能保证连续性。需要说明的是：磁场越小，磁周期越大，Ny值也会越大。在计算中，Ny的值取决于最小磁场的大小，例如：磁场最小取为 0.001 就比磁场最小取为 0.01 的计算量大，因为 Ny 的值就有10倍的差距。
博主您好，您选取的规范为A_x = -B*y , t_x跃迁项会增加佩尔斯替代中的相位，但我发现您的程序中，在x方向的跃迁其实改变的是t_y的值,这里没看懂。这里的t_x和t_y指的是您“方格子模型在实空间中的哈密顿量形式”中的t_x和t_y。
是加在t_x上，代码片段如下，只是t_x的值是依赖于y的值。参考：磁场和磁势的选取。
关老师，那种没有边缘态能量的蝴蝶（就是最开始的那种蝴蝶翅膀里面是有空白的地方）该怎么画
原始文献好像是用连续模型进行计算的。这里是用紧束缚模型算的，可以在x和y方向简单地加上周期性边界条件，就可以得到没有边缘态的蝴蝶谱。
关老师您好，请问磁原胞的大小应该如何确定？
Peierls替换会导致在条带边界上不连续，只能通过扩胞来处理，这时候才可以应用周期性边界条件。元胞要取多大具体看Peierls替换中系数的周期。一般来说，磁场越小，需要的磁元胞越大。
关老师，您能否画一个2维电流分布? 谢谢
参考这篇：局域电流的计算（附Python代码）。
请问这样做不受gauge影响吗？好像这样做对应着连续场的微分，也就是 p=v+qA. 我们是在算v.
算能带和电流，计算结果跟gauge没关系吧
你好 这个模型有matalb程序嘛
没有的，可以根据语义进行改写。代码大概能看懂一点点吧。